"""Reanalyse published splice-event peptides against GENCODE v50. Python 3.12; openpyxl 3.1.5; pyahocorasick 2.2.0. Inputs are downloaded only with --fetch; large reference files stay in --input-dir. """ import argparse import collections import csv import gzip import hashlib import importlib.metadata import json import pathlib import re import sys import urllib.request import ahocorasick import openpyxl SOURCES = { "supplement.xlsx": "https://media.springernature.com/original/springer-static/esm/art%3A10.1038%2Fs41587-023-01714-x/MediaObjects/41587_2023_1714_MOESM3_ESM.xlsx", "gencode.v50.pc_translations.fa.gz": "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_50/gencode.v50.pc_translations.fa.gz", "gencode.v50.annotation.gtf.gz": "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_50/gencode.v50.annotation.gtf.gz", } def write_csv(path, rows, fields): with gzip.open(path, "wt", encoding="utf-8", newline="") as out: writer = csv.DictWriter(out, fieldnames=fields) writer.writeheader() writer.writerows(rows) def fasta(path): header, parts = None, [] with gzip.open(path, "rt") as handle: for line in handle: if line.startswith(">"): if header is not None: yield header, "".join(parts) header, parts = line[1:].strip().split("|"), [] else: parts.append(line.strip()) if header is not None: yield header, "".join(parts) def nonnested_pair(peptides): peptides = sorted({p for p in peptides if len(p) >= 9}) return any(a not in b and b not in a for i, a in enumerate(peptides) for b in peptides[i + 1:]) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--input-dir", type=pathlib.Path, required=True) parser.add_argument("--output-dir", type=pathlib.Path, required=True) parser.add_argument("--fetch", action="store_true") args = parser.parse_args() args.input_dir.mkdir(parents=True, exist_ok=True) args.output_dir.mkdir(parents=True, exist_ok=True) manifest = {} expected_path = pathlib.Path(__file__).with_name("sources.json") expected = json.loads(expected_path.read_text()) if expected_path.exists() else None for name, url in SOURCES.items(): path = args.input_dir / name if args.fetch and not path.exists(): urllib.request.urlretrieve(url, path) digest = hashlib.file_digest(path.open("rb"), "sha256").hexdigest() manifest[name] = {"url": url, "sha256": digest, "bytes": path.stat().st_size} if expected: assert digest == expected["inputs"][name]["sha256"], f"Input changed: {name}" workbook = openpyxl.load_workbook(args.input_dir / "supplement.xlsx", read_only=True, data_only=True) iterator = workbook["Table S3"].values headers = next(iterator) assert headers[14] == "Proteomics|aspn|Peptides|Path0" assert headers[37] == "Proteomics|trypsin|MSMScount|Path1" observations, peptide_context, events = [], collections.defaultdict(set), {} max_spectrum_count = collections.Counter() rows = 0 listed, zero_only, malformed_peptides = set(), set(), set() exclusions = [] cells = collections.Counter() for excel_row, row in enumerate(iterator, 2): rows += 1 cells[row[0]] += 1 # A local event, pooled over cell lines; do not rely on row-local Id. event_key = tuple(str(row[i]) for i in (2, 3, 4, 5, 7, 8, 10)) event = events.setdefault(event_key, {"paths": set(), "rna_paths": set(), "type": row[9]}) for path in (0, 1): if float(row[12 + path] or 0) > 0: event["rna_paths"].add(path) for column in range(14, 38, 2): if row[column] is None: assert row[column + 1] is None continue peptides = str(row[column]).split(";") counts = str(row[column + 1]).split(";") if len(peptides) != len(counts): malformed_peptides.update(peptides) exclusions.append({"excel_row": excel_row, "cell_line": row[0], "column": headers[column], "peptides": str(row[column]), "counts": str(row[column + 1]), "reason": "Peptide and spectrum-count lists have different lengths"}) continue _, enzyme, _, path_name = headers[column].split("|") for peptide, count in zip(peptides, counts): assert re.fullmatch("[ACDEFGHIKLMNPQRSTVWYU]+", peptide), peptide count = float(count) assert count >= 0 and count.is_integer() listed.add(peptide) if count == 0: zero_only.add(peptide) continue event["paths"].add(int(path_name[-1])) max_spectrum_count[peptide] = max(max_spectrum_count[peptide], int(count)) peptide_context[peptide].add((str(row[5]), str(row[6]), row[0], enzyme)) observations.append({"excel_row": excel_row, "cell_line": row[0], "gene_id": row[5], "gene_symbol": row[6], "event_type": row[9], "path": path_name, "enzyme": enzyme, "peptide": peptide, "spectrum_count": int(count)}) workbook.close() print(f"Read {rows} event/cell rows, {len(peptide_context)} observed peptide sequences", flush=True) # Restrict targets to full-length chromosome-only protein_coding models. target_transcripts = set() target_gene_ids = set() all_pc_transcripts = 0 attrs_re = re.compile(r'(\w+) "([^"]*)"') with gzip.open(args.input_dir / "gencode.v50.annotation.gtf.gz", "rt") as handle: for line in handle: if line.startswith("#"): continue columns = line.rstrip().split("\t") if columns[2] != "transcript": continue attrs = dict(attrs_re.findall(columns[8])) if attrs.get("transcript_type") == "protein_coding": all_pc_transcripts += 1 if '"cds_start_NF"' not in columns[8] and '"cds_end_NF"' not in columns[8]: target_transcripts.add(attrs["transcript_id"]) target_gene_ids.add(attrs["gene_id"]) # Keep all FASTA regions/biotypes as background competitors. Group identical # complete sequences; separate transcripts encoding the same sequence are not # separate sequence identities. I/L grouping reflects MS indistinguishability. sequences, sequence_index, records, target_ids = [], {}, collections.defaultdict(list), set() fasta_records = 0 matched_transcripts = set() for header, sequence in fasta(args.input_dir / "gencode.v50.pc_translations.fa.gz"): fasta_records += 1 sid = sequence_index.get(sequence) if sid is None: sid = len(sequences) sequence_index[sequence] = sid sequences.append(sequence) protein_id, transcript_id, gene_id = header[:3] records[sid].append((protein_id, transcript_id, gene_id)) if transcript_id in target_transcripts: target_ids.add(sid) matched_transcripts.add(transcript_id) il_sequences = {} il_id = {} for sid, sequence in enumerate(sequences): normalized = sequence.replace("I", "L") il_id[sid] = il_sequences.setdefault(normalized, len(il_sequences)) peptides = sorted(peptide_context) automaton = ahocorasick.Automaton() normalized_peptides = collections.defaultdict(list) for peptide in peptides: normalized_peptides[peptide.replace("I", "L")].append(peptide) for peptide in normalized_peptides: automaton.add_word(peptide, peptide) automaton.make_automaton() exact_hits = {p: set() for p in peptides} il_hits = {p: set() for p in peptides} for sid, sequence in enumerate(sequences): for end, normalized in automaton.iter(sequence.replace("I", "L")): for peptide in normalized_peptides[normalized]: il_hits[peptide].add(sid) if sequence[end + 1 - len(peptide):end + 1] == peptide: exact_hits[peptide].add(sid) print(f"Mapped against {len(sequences)} distinct reference sequences", flush=True) peptide_rows, target_any, target_unique = [], set(), collections.defaultdict(set) category_counts = collections.Counter() exact_category_counts = collections.Counter() for peptide in peptides: hits = il_hits[peptide] category = "no_match" if not hits else "one_sequence" if len(hits) == 1 else "multiple_sequences" category_counts[category] += 1 exact_category_counts["no_match" if not exact_hits[peptide] else "one_sequence" if len(exact_hits[peptide]) == 1 else "multiple_sequences"] += 1 target_any.update(hits & target_ids) if len(hits) == 1: sid = next(iter(hits)) if sid in target_ids: target_unique[sid].add(peptide.replace("I", "L")) contexts = peptide_context[peptide] peptide_rows.append({"peptide": peptide, "length": len(peptide), "il_category": category, "exact_sequence_matches": len(exact_hits[peptide]), "il_sequence_matches": len(hits), "il_equivalence_classes": len({il_id[sid] for sid in hits}), "target_sequence_matches": len(hits & target_ids), "max_spectrum_count_in_source_cell": max_spectrum_count[peptide], "reported_genes": ";".join(sorted({x[0] for x in contexts})), "reported_symbols": ";".join(sorted({x[1] for x in contexts})), "cell_lines": ";".join(sorted({x[2] for x in contexts})), "enzymes": ";".join(sorted({x[3] for x in contexts})), "matched_protein_ids": ";".join(sorted({r[0] for sid in hits for r in records[sid]}))}) sequence_rows = [] for sid in sorted(target_any): unique = target_unique.get(sid, set()) sequence_rows.append({"sequence_sha256": hashlib.sha256(sequences[sid].encode()).hexdigest(), "length": len(sequences[sid]), "protein_ids": ";".join(sorted({r[0] for r in records[sid]})), "gene_ids": ";".join(sorted({r[2] for r in records[sid]})), "any_matching_splice_peptide": int(sid in target_any), "unique_peptides": ";".join(sorted(unique)), "unique_peptide_count": len(unique), "two_nonnested_unique_9mers": int(nonnested_pair(unique))}) summary = { "table_s3_rows": rows, "cell_line_rows": dict(cells), "listed_peptides_including_zero_counts": len(listed), "peptides_only_ever_with_zero_count": len(zero_only - peptide_context.keys()), "malformed_cells_excluded": len(exclusions), "malformed_cell_peptides_without_independent_positive_count": len(malformed_peptides - peptide_context.keys()), "observed_distinct_peptides": len(peptides), "positive_peptide_observations": len(observations), "observed_genes": len({x[0] for contexts in peptide_context.values() for x in contexts}), "pooled_local_events": len(events), "pooled_events_by_observed_paths": dict(collections.Counter(len(e["paths"]) for e in events.values())), "fasta_records": fasta_records, "background_distinct_sequences": len(sequences), "background_il_equivalence_classes": len(il_sequences), "chromosome_protein_coding_transcripts": all_pc_transcripts, "full_length_chromosome_protein_coding_transcripts": len(target_transcripts), "target_transcripts_with_fasta_sequence": len(matched_transcripts), "target_distinct_sequences": len(target_ids), "exact_peptide_categories": dict(exact_category_counts), "il_peptide_categories": dict(category_counts), "target_sequences_with_any_match": len(target_any), "target_sequences_without_match": len(target_ids - target_any), "target_sequences_with_unique_peptide": len(target_unique), "target_sequences_with_two_nonnested_unique_9mers": sum(nonnested_pair(x) for x in target_unique.values()), "sensitivity_9mers_with_at_least_two_spectra_in_one_source_cell": dict(collections.Counter( row["il_category"] for row in peptide_rows if row["length"] >= 9 and row["max_spectrum_count_in_source_cell"] >= 2)), "sensitivity_full_length_targets_only": dict(collections.Counter( "no_match" if row["target_sequence_matches"] == 0 else "one_sequence" if row["target_sequence_matches"] == 1 else "multiple_sequences" for row in peptide_rows)), "APP_example_peptide_categories": dict(collections.Counter( row["il_category"] for row in peptide_rows if "APP" in row["reported_symbols"].split(";"))), } assert sum(category_counts.values()) == len(peptides) assert set(target_unique) <= target_any <= target_ids assert len(listed) == len(peptides) + len(zero_only - peptide_context.keys()) write_csv(args.output_dir / "observations.csv.gz", observations, list(observations[0])) write_csv(args.output_dir / "excluded-cells.csv.gz", exclusions, list(exclusions[0])) write_csv(args.output_dir / "peptides.csv.gz", peptide_rows, list(peptide_rows[0])) write_csv(args.output_dir / "sequences.csv.gz", sequence_rows, list(sequence_rows[0])) (args.output_dir / "summary.json").write_text(json.dumps(summary, indent=2) + "\n") (args.output_dir / "sources.json").write_text(json.dumps({ "retrieved": "2026-09-19", "inputs": manifest, "python": sys.version.split()[0], "packages": {p: importlib.metadata.version(p) for p in ("openpyxl", "pyahocorasick")}, }, indent=2) + "\n") print(json.dumps(summary, indent=2), flush=True) if __name__ == "__main__": main()