Skip to content

API reference

Top-level conveniences (import proteoform_regions as pfr) plus the module reference. Every public function is deterministic given its inputs; stochastic entry points take an explicit seed (default 42 = the paper's values).

Package

proteoform_regions

proteoform-regions.

Proteoform-aware region inference and detection-bias quantification from bottom-up proteomics peptide tables.

This package is the software implementation of the pipeline described in "Basic-peptide depletion and precursor masking in bottom-up proteomics" (Saadat et al., 2026). Modules:

  • schema -- canonical evidence schema, confidence tiers, pI classification
  • guard -- preview-input guards + sha256 manifest verification
  • physchem -- pI / MW / net-charge engine (multiple pK models) + benchmarks
  • digest -- in-silico tryptic digestion + digestion-matched backgrounds
  • regions -- single-linkage detected-region stitching (gap rule)
  • parsers -- search-engine peptide-table adapters + registry
  • study -- StudyConfig / cohort YAML
  • harmonize -- stage 1 orchestration (downloads + adapters + concat)
  • mapping -- UniProt accession mapping + protein properties
  • stats -- detection-bias models (permutation null, OR meta-analysis)
  • pipeline -- end-to-end orchestration (map/compute/stitch/bias + run)
  • cli -- command-line interface

FileSpec(name, url=None, sha256=None) dataclass

One raw result file belonging to a study.

name is the local filename under <data_dir>/<accession>/; url is the retrieval source (omitted for locally-curated files); sha256 pins the expected content when known.

RunReport(out_dir, stages, artifacts=dict(), counts=dict(), parameters=dict(), generated='', proteoform_regions_version='') dataclass

Run provenance: executed stages, artifacts, row counts, parameters.

StudyConfig(dataset_accession, adapter, study_id=None, study_year='2026', parser_family=None, sample_type='', disease_context='', assay_type='LC-MS/MS', acquisition_mode='', search_engine='', representative_file_arm='primary result file', representative_file_species='Homo sapiens', publication_doi=None, files=list(), build_usi=None, extract_protein_names=True, drop_zero_intensities=False, psm_count_zero_as_none=False) dataclass

Declarative description of one dataset's peptide-evidence extraction.

build_background(protein_sequences, max_missed=2, min_len=7, max_len=60)

Digestion-matched tryptic background over {accession: sequence} (S3).

Same digest protocol and same physicochemical columns as the detected peptides, so the RQ2 comparison is digestion-matched by construction.

Source code in src/proteoform_regions/pipeline.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def build_background(
    protein_sequences: dict[str, str],
    max_missed: int = 2,
    min_len: int = 7,
    max_len: int = 60,
) -> pd.DataFrame:
    """Digestion-matched tryptic background over {accession: sequence} (S3).

    Same digest protocol and same physicochemical columns as the detected
    peptides, so the RQ2 comparison is digestion-matched by construction.
    """
    digested = background_peptides(protein_sequences, max_missed, min_len, max_len)
    pi_cache: dict[str, tuple] = {}

    def props(seq: str) -> tuple:
        if seq not in pi_cache:
            pi_cache[seq] = (
                biopython_pi(seq),
                biopython_mw(seq),
                gravy(seq),
                peptide_features(seq),
            )
        return pi_cache[seq]

    rows = []
    for rec in digested:
        p_i, mw, gr, feats = props(rec["peptide_sequence"])
        rows.append(
            {
                "accession": rec["accession"],
                "peptide_sequence": rec["peptide_sequence"],
                "missed_cleavages": rec["missed_cleavages"],
                "peptide_pi": p_i,
                "peptide_mw_da": mw,
                "peptide_gravy": gr,
                "peptide_class_anl003": classify_pi(p_i),
                **feats,
            }
        )
    return pd.DataFrame(rows)

compute_features(evidence, mapping, proteins, calculator='bjellqvist')

Peptide feature table over exactly-mapped rows (S3 feature engine).

calculator selects the peptide pI engine: bjellqvist (Biopython, the paper's primary) or an HH table name (hh_rodwell_emboss_lineage / hh_lehninger) for the named sensitivity analysis. Net charge at pH 7.4 always uses the Rodwell/EMBOSS HH table, as in the source.

Source code in src/proteoform_regions/pipeline.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def compute_features(
    evidence: pd.DataFrame,
    mapping: pd.DataFrame,
    proteins: pd.DataFrame,
    calculator: str = "bjellqvist",
) -> pd.DataFrame:
    """Peptide feature table over exactly-mapped rows (S3 feature engine).

    ``calculator`` selects the peptide pI engine: ``bjellqvist`` (Biopython,
    the paper's primary) or an HH table name (``hh_rodwell_emboss_lineage`` /
    ``hh_lehninger``) for the named sensitivity analysis. Net charge at pH 7.4
    always uses the Rodwell/EMBOSS HH table, as in the source.
    """
    mapped = mapping[mapping["mapping_tier"].isin(_EXACT_TIERS)].copy()
    protein_lookup = proteins.set_index("accession")

    def _pi_map(seqs: list[str]) -> dict[str, float]:
        if calculator == "bjellqvist":
            return {s: biopython_pi(s) for s in seqs}
        if calculator in ("hh_rodwell_emboss_lineage", "hh_lehninger"):
            return dict(zip(seqs, hh_pi_batch(seqs, calculator), strict=True))
        raise ValueError(f"unknown calculator {calculator!r}")

    seqs_needed = pd.unique(mapped["peptide_sequence"]).tolist()
    pi_map = _pi_map([str(s) for s in seqs_needed])
    pi_cache: dict[str, tuple] = {}
    for seq in seqs_needed:
        pi_cache[seq] = (
            pi_map[seq],
            biopython_mw(seq),
            gravy(seq),
            peptide_features(seq),
        )

    feature_rows = []
    for record in mapped.itertuples(index=False):
        seq = record.peptide_sequence
        p_i, mw, gr, feats = pi_cache[seq]
        precursor = (
            protein_lookup.loc[record.mapped_accession]
            if record.mapped_accession in protein_lookup.index
            else None
        )
        row = {
            "peptide_row_id": record.peptide_row_id,
            "dataset_accession": record.dataset_accession,
            "peptide_sequence": seq,
            "mapped_accession": record.mapped_accession,
            "mapping_tier": record.mapping_tier,
            "peptide_pi": p_i,
            "peptide_mw_da": mw,
            "peptide_gravy": gr,
            "contaminant_flag": record.contaminant_flag,
            **feats,
            "peptide_class_anl003": classify_pi(p_i),
        }
        if precursor is not None:
            row.update(
                {
                    "precursor_pi": precursor["precursor_pi"],
                    "precursor_mw_da": precursor["precursor_mw_da"],
                    "precursor_length": precursor["sequence_length"],
                    "gene_name": precursor["gene_name"],
                    "protein_name": precursor["protein_name"],
                    "pi_divergence_peptide_minus_precursor": p_i - precursor["precursor_pi"],
                    "mw_ratio_peptide_over_precursor": (
                        mw / precursor["precursor_mw_da"] if precursor["precursor_mw_da"] else None
                    ),
                }
            )
        feature_rows.append(row)
    features = pd.DataFrame(feature_rows)

    # net charge at pH 7.4 (HH Rodwell lineage; batched as in the source)
    unique_seqs = [str(s) for s in features["peptide_sequence"].unique()]
    features["net_charge_ph74"] = features["peptide_sequence"].map(charge_at_ph(unique_seqs, 7.4))
    features["charge_density"] = features["net_charge_ph74"] / features["peptide_length"]
    return features

ensure_downloads(studies, data_dir, timeout=600.0)

Fetch missing raw files (skip-if-exists, one retry) and return a manifest.

Source code in src/proteoform_regions/harmonize.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def ensure_downloads(
    studies: list[StudyConfig], data_dir: str | Path, timeout: float = 600.0
) -> pd.DataFrame:
    """Fetch missing raw files (skip-if-exists, one retry) and return a manifest."""
    manifest_rows: list[dict] = []
    for study in studies:
        for spec in study.files:
            if spec.url is None:
                continue
            dest = Path(data_dir) / study.dataset_accession / spec.name
            dest.parent.mkdir(parents=True, exist_ok=True)
            if not dest.exists():
                _fetch(spec.url, dest, timeout)
            manifest_rows.append(
                {
                    "dataset_accession": study.dataset_accession,
                    "file_name": spec.name,
                    "source_url": spec.url,
                    "bytes_on_disk": dest.stat().st_size if dest.exists() else 0,
                    "sha256": sha256_of(dest) if dest.exists() else None,
                }
            )
    return pd.DataFrame(manifest_rows)

load_cohort(path)

Load a cohort file: a YAML list of study dicts, or a dict with 'studies'.

Source code in src/proteoform_regions/study.py 
122
123
124
125
126
127
128
129
def load_cohort(path: str | Path) -> list[StudyConfig]:
    """Load a cohort file: a YAML list of study dicts, or a dict with 'studies'."""
    import yaml

    payload = yaml.safe_load(Path(path).read_text(encoding="utf-8"))
    if isinstance(payload, dict):
        payload = payload["studies"]
    return [StudyConfig.from_dict(d) for d in payload]

map_to_uniprot(evidence, cache_dir='data/raw/uniprot_cache', snapshot_date='20260917', offline=False)

Full S2 stage: normalize -> resolve -> map. Returns (mapping, proteins).

proteins carries per-protein properties (precursor pI/MW via the same Biopython engine as peptide features); mapping carries per-row tiers, contaminant and degeneracy flags.

Source code in src/proteoform_regions/mapping.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
def map_to_uniprot(
    evidence: pd.DataFrame,
    cache_dir: str | Path = "data/raw/uniprot_cache",
    snapshot_date: str = "20260917",
    offline: bool = False,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Full S2 stage: normalize -> resolve -> map. Returns (mapping, proteins).

    ``proteins`` carries per-protein properties (precursor pI/MW via the same
    Biopython engine as peptide features); ``mapping`` carries per-row tiers,
    contaminant and degeneracy flags.
    """
    if "protein_accession_raw" not in evidence.columns:
        raise ValueError("evidence table lacks protein_accession_raw (not a canonical S1 table?)")
    all_bases: set[str] = set()
    for raw in evidence["protein_accession_raw"]:
        all_bases.update(token for token, _ in normalize_accessions(raw))
    entries = ensure_uniprot_entries(all_bases, cache_dir, snapshot_date, offline=offline)
    mapping = map_peptides(evidence, entries)
    proteins = protein_properties_table(entries)
    return mapping, proteins

normalize_accessions(raw_value)

Split a raw accession field into [(base_accession, isoform_suffix), ...].

Handles fasta sp|P02768|ALBU_HUMAN prefixes, ;/, separated groups, -2 isoform suffixes, and drops CON_/REV_ tokens. Verbatim from the publication's S2 (cell 4).

Source code in src/proteoform_regions/mapping.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def normalize_accessions(raw_value) -> list[tuple[str, str | None]]:
    """Split a raw accession field into [(base_accession, isoform_suffix), ...].

    Handles fasta ``sp|P02768|ALBU_HUMAN`` prefixes, ``;``/``,``
    separated groups, ``-2`` isoform suffixes, and drops ``CON_``/``REV_``
    tokens. Verbatim from the publication's S2 (cell 4).
    """
    if raw_value is None or (isinstance(raw_value, float) and np.isnan(raw_value)):
        return []
    tokens: list[tuple[str, str | None]] = []
    for token in str(raw_value).replace(",", ";").split(";"):
        token = token.strip()
        if not token:
            continue
        if "|" in token:
            parts = token.split("|")
            token = parts[1] if len(parts) >= 2 else parts[-1]
        if token.upper().startswith(("CON_", "REV_")):
            continue
        isoform = None
        match = ISOFORM_SUFFIX.search(token)
        if match and token[: match.start()].upper().startswith(("P", "Q", "O", "A")):
            isoform = token[match.start() :]
            token = token[: match.start()]
        tokens.append((token, isoform))
    return tokens

parse_study(config, data_dir)

Run a study's adapter over its files; returns (canonical rows, stats).

Source code in src/proteoform_regions/parsers/__init__.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def parse_study(config: StudyConfig, data_dir: str | Path) -> tuple[list[dict], dict]:
    """Run a study's adapter over its files; returns (canonical rows, stats)."""
    parser = get_parser(config.adapter)
    paths = config.file_paths(data_dir)
    result = parser(config, paths)
    rows, stats = result[0], result[1:]
    stats_out = {
        "dataset_accession": config.dataset_accession,
        "adapter": config.adapter,
        "n_files": len(paths),
        "n_evidence_rows": len(rows),
    }
    # adapters report heterogeneous leading stats (raw rows / SIRs / lines...)
    for i, value in enumerate(stats):
        stats_out[f"stat{i + 1}"] = value
    return rows, stats_out

quantify_bias(features, background, seed=st.SEED, b_perm=st.B_PERM, b_boot=st.B_BOOT, run_rq1=True)

Detection-bias quantification on feature + background tables (S8).

RQ2: per-study basic-vs-rest OR (Woolf, unconditional 0.5) against each study's digestion-matched parent-protein background; inverse-variance fixed effect; DerSimonian-Laird random effects; leave-one-study-out envelope. RQ1: per-protein divergence permutation null (seed-aware).

Set run_rq1=False to skip the permutation null (it dominates runtime on paper-scale tables: ~5,400 proteins x B=2,000 draws).

Source code in src/proteoform_regions/pipeline.py 
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
def quantify_bias(
    features: pd.DataFrame,
    background: pd.DataFrame,
    seed: int = st.SEED,
    b_perm: int = st.B_PERM,
    b_boot: int = st.B_BOOT,
    run_rq1: bool = True,
) -> dict:
    """Detection-bias quantification on feature + background tables (S8).

    RQ2: per-study basic-vs-rest OR (Woolf, unconditional 0.5) against each
    study's digestion-matched parent-protein background; inverse-variance
    fixed effect; DerSimonian-Laird random effects; leave-one-study-out
    envelope. RQ1: per-protein divergence permutation null (seed-aware).

    Set ``run_rq1=False`` to skip the permutation null (it dominates runtime
    on paper-scale tables: ~5,400 proteins x B=2,000 draws).
    """
    det = features.drop_duplicates(["dataset_accession", "peptide_sequence"])
    study_parents = (
        features.groupby("dataset_accession")["mapped_accession"].agg(lambda s: set(s.dropna())).to_dict()
    )

    study_rows = []
    for study in sorted(det["dataset_accession"].unique()):
        d_cls = det[det["dataset_accession"] == study]["peptide_class_anl003"].dropna()
        bg = background[background["accession"].isin(study_parents[study])].drop_duplicates(
            "peptide_sequence"
        )
        b_cls = bg["peptide_class_anl003"].dropna()
        r = st.study_or(d_cls.to_numpy(), b_cls.to_numpy(), "basic")
        r["dataset_accession"] = study
        study_rows.append(r)
    or_table = pd.DataFrame(study_rows)

    mu_fe, se_fe, p_fe, Q, I2 = st.fe_summary(or_table)
    dl = st.dl_random_effects(or_table)
    loo_lo, loo_hi = st.loo_envelope(or_table)

    result = {
        "k": int(len(or_table)),
        "or_table": or_table,
        "fe": {
            "log_or": mu_fe,
            "se": se_fe,
            "or": float(np.exp(mu_fe)),
            "ci_low": float(np.exp(mu_fe - 1.96 * se_fe)),
            "ci_high": float(np.exp(mu_fe + 1.96 * se_fe)),
            "p": p_fe,
            "Q": Q,
            "I2": I2,
        },
        "random_effects": dl,
        "loo_envelope": [loo_lo, loo_hi],
    }
    if run_rq1:
        result["rq1"] = st.rq1_permutation(features, background, seed=seed, b_perm=b_perm, b_boot=b_boot)
    return result

run(studies, out_dir, stages='all', data_dir='data/raw/peptide_tables', download=True, uniprot_cache='data/raw/uniprot_cache', snapshot_date='20260917', offline=False, calculator='bjellqvist', max_missed=2, gap=REGION_GAP_AA, terminal_fraction=TERMINAL_FRACTION, min_peptides=MIN_REGION_PEPTIDES, seed=st.SEED, b_perm=st.B_PERM, b_boot=st.B_BOOT, run_rq1=True)

Execute pipeline stages over a cohort, writing paper-vocabulary artifacts.

studies is a cohort YAML path or a list of StudyConfig. Artifacts land in out_dir under the paper's data/processed/ filenames; every run closes by writing a sha256 manifest (manifest.sha256) over its own outputs plus a machine-readable run_report.json.

Source code in src/proteoform_regions/pipeline.py 
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
def run(
    studies: str | Path | list[StudyConfig],
    out_dir: str | Path,
    stages: str | list[str] = "all",
    data_dir: str | Path = "data/raw/peptide_tables",
    download: bool = True,
    uniprot_cache: str | Path = "data/raw/uniprot_cache",
    snapshot_date: str = "20260917",
    offline: bool = False,
    calculator: str = "bjellqvist",
    max_missed: int = 2,
    gap: int = REGION_GAP_AA,
    terminal_fraction: float = TERMINAL_FRACTION,
    min_peptides: int = MIN_REGION_PEPTIDES,
    seed: int = st.SEED,
    b_perm: int = st.B_PERM,
    b_boot: int = st.B_BOOT,
    run_rq1: bool = True,
) -> RunReport:
    """Execute pipeline stages over a cohort, writing paper-vocabulary artifacts.

    ``studies`` is a cohort YAML path or a list of ``StudyConfig``. Artifacts
    land in ``out_dir`` under the paper's ``data/processed/`` filenames; every
    run closes by writing a sha256 manifest (``manifest.sha256``) over its own
    outputs plus a machine-readable ``run_report.json``.
    """
    import proteoform_regions

    if isinstance(studies, (str, Path)):
        studies = load_cohort(studies)
    selected = _selected_stages(stages)
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)

    report = RunReport(
        out_dir=out_dir,
        stages=selected,
        parameters={
            "calculator": calculator,
            "max_missed": max_missed,
            "gap": gap,
            "terminal_fraction": terminal_fraction,
            "min_peptides": min_peptides,
            "seed": seed,
            "b_perm": b_perm,
            "b_boot": b_boot,
            "snapshot_date": snapshot_date,
            "n_studies": len(studies),
        },
        generated=datetime.now().isoformat(timespec="minutes"),
        proteoform_regions_version=proteoform_regions.__version__,
    )
    evidence = mapping_df = features = background = regions = None
    proteins = None

    def _write(frame: pd.DataFrame, name: str, stage: str) -> Path:
        path = out_dir / name
        frame.to_csv(path, index=False)
        report.artifacts[stage] = name
        report.counts[stage] = int(len(frame))
        return path

    if "harmonize" in selected:
        evidence, _stats = harmonize(studies, data_dir, download=download)
        _write(evidence, "peptide_evidence_table.csv", "harmonize")
    elif "map" in selected or "compute" in selected or "stitch" in selected:
        evidence = pd.read_csv(out_dir / "peptide_evidence_table.csv", low_memory=False)

    if "map" in selected:
        mapping_df, proteins = _map_stage(evidence, uniprot_cache, snapshot_date, offline=offline)
        _write(mapping_df, "processed_peptide_protein_map.csv", "map")
        proteins_out = proteins.drop(columns=["sequence"], errors="ignore")
        _write(proteins_out, "processed_protein_properties.csv", "map_properties")
    elif "compute" in selected or "stitch" in selected:
        mapping_df = pd.read_csv(out_dir / "processed_peptide_protein_map.csv")

    if "compute" in selected:
        if proteins is None:
            proteins = pd.read_csv(out_dir / "processed_protein_properties.csv", low_memory=False)
        features = compute_features(evidence, mapping_df, proteins, calculator=calculator)
        _write(features, "processed_feature_table.csv", "compute")
    elif "stitch" in selected or "bias" in selected:
        features = pd.read_csv(out_dir / "processed_feature_table.csv", low_memory=False)

    if "background" in selected:
        if features is None:
            features = pd.read_csv(out_dir / "processed_feature_table.csv", low_memory=False)
        parents = set(features["mapped_accession"].dropna())
        sequences = _sequences_from_snapshot(parents, uniprot_cache, snapshot_date)
        background = build_background(sequences, max_missed=max_missed)
        _write(background, "processed_tryptic_background.csv", "background")
    elif "bias" in selected:
        background = pd.read_csv(out_dir / "processed_tryptic_background.csv")

    if "stitch" in selected:
        parents = set(mapping_df["mapped_accession"].dropna())
        precursor_sequences = _sequences_from_snapshot(parents, uniprot_cache, snapshot_date)
        regions = stitch_regions(
            features,
            mapping_df,
            precursor_sequences,
            gap=gap,
            terminal_fraction=terminal_fraction,
            min_peptides=min_peptides,
        )
        _write(regions, "processed_region_inference.csv", "stitch")

    if "bias" in selected:
        bias = quantify_bias(features, background, seed=seed, b_perm=b_perm, b_boot=b_boot, run_rq1=run_rq1)
        bias_payload = {
            "k": bias["k"],
            "fe": bias["fe"],
            "random_effects": bias["random_effects"],
            "loo_envelope": bias["loo_envelope"],
            "or_table": bias["or_table"].to_dict(orient="records"),
            "rq1": {k: v for k, v in bias.get("rq1", {}).items() if k != "per_protein"},
            "parameters": {"seed": seed, "b_perm": b_perm, "b_boot": b_boot},
        }
        (out_dir / "processed_inference_results.json").write_text(
            json.dumps(bias_payload, indent=2, default=str), encoding="utf-8"
        )
        report.artifacts["bias"] = "processed_inference_results.json"
        report.counts["bias"] = bias["k"]

    (out_dir / "run_report.json").write_text(json.dumps(report.to_dict(), indent=2), encoding="utf-8")
    outputs = [out_dir / name for name in report.artifacts.values()]
    if outputs:
        write_manifest(outputs, out_dir / "manifest.sha256")
    report.artifacts["run_report"] = "run_report.json"
    report.artifacts["manifest"] = "manifest.sha256"
    return report

stitch_regions(features, mapping, precursor_sequences=None, gap=REGION_GAP_AA, terminal_fraction=TERMINAL_FRACTION, min_peptides=MIN_REGION_PEPTIDES)

Detected-region inference (S4 walk): position spans -> stitched regions.

Start positions come from the mapping table (first exact-substring hit); rows without exact mapping never reach this stage. precursor_sequences (e.g. from the UniProt snapshot) enables region sequence / region pI/MW columns; without it they stay null (positions and tiers are unaffected).

Source code in src/proteoform_regions/pipeline.py 
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
def stitch_regions(
    features: pd.DataFrame,
    mapping: pd.DataFrame,
    precursor_sequences: dict[str, str] | None = None,
    gap: int = REGION_GAP_AA,
    terminal_fraction: float = TERMINAL_FRACTION,
    min_peptides: int = MIN_REGION_PEPTIDES,
) -> pd.DataFrame:
    """Detected-region inference (S4 walk): position spans -> stitched regions.

    Start positions come from the mapping table (first exact-substring hit);
    rows without exact mapping never reach this stage. ``precursor_sequences``
    (e.g. from the UniProt snapshot) enables region sequence / region pI/MW
    columns; without it they stay null (positions and tiers are unaffected).
    """
    precursor_sequences = precursor_sequences or {}
    mapped = mapping[mapping["mapping_tier"].isin(_EXACT_TIERS)].copy()
    if mapped.empty:
        return pd.DataFrame()
    mapped["start_position"] = mapped["positions"].astype(str).str.split(";").str[0].astype(float)
    mapped = mapped[mapped["start_position"].notna()]
    merged = mapped.merge(
        features[
            [
                "peptide_row_id",
                "peptide_pi",
                "peptide_length",
                "precursor_pi",
                "precursor_mw_da",
                "precursor_length",
                "gene_name",
            ]
        ],
        on="peptide_row_id",
        how="left",
        suffixes=("", "_feat"),
    )
    merged["end_position"] = merged["start_position"] + merged["peptide_length"] - 1

    region_rows = []
    for (study, accession), group in merged.groupby(["dataset_accession", "mapped_accession"]):
        precursor_length = group["precursor_length"].iloc[0]
        if not precursor_length or precursor_length != precursor_length:
            continue
        spans = list(
            zip(
                group["start_position"],
                group["end_position"],
                group["peptide_sequence"],
                group["peptide_pi"],
                group["peptide_row_id"],
                strict=True,
            )
        )
        for cluster in stitch_spans(spans, gap):
            region_rows.append(
                region_record(
                    study,
                    accession,
                    cluster,
                    precursor_length=int(precursor_length),
                    precursor_pi=group["precursor_pi"].iloc[0],
                    precursor_mw=group["precursor_mw_da"].iloc[0],
                    gene_name=group["gene_name"].iloc[0],
                    precursor_sequence=precursor_sequences.get(accession, ""),
                    gap=gap,
                    terminal_fraction=terminal_fraction,
                    min_peptides=min_peptides,
                )
            )
    return pd.DataFrame(region_rows)

Study configuration

proteoform_regions.study.StudyConfig(dataset_accession, adapter, study_id=None, study_year='2026', parser_family=None, sample_type='', disease_context='', assay_type='LC-MS/MS', acquisition_mode='', search_engine='', representative_file_arm='primary result file', representative_file_species='Homo sapiens', publication_doi=None, files=list(), build_usi=None, extract_protein_names=True, drop_zero_intensities=False, psm_count_zero_as_none=False) dataclass

Declarative description of one dataset's peptide-evidence extraction.

proteoform_regions.study.FileSpec(name, url=None, sha256=None) dataclass

One raw result file belonging to a study.

name is the local filename under <data_dir>/<accession>/; url is the retrieval source (omitted for locally-curated files); sha256 pins the expected content when known.

Harmonize (stage 1)

proteoform_regions.harmonize

Stage 1: harmonize heterogeneous peptide tables into the canonical evidence schema.

harmonize is a pure function over study configs + cached raw files: it downloads missing files (retry-once policy, failures documented never fabricated - the paper's S1 ensure_downloads semantics), runs each study's adapter, and concatenates canonical-schema rows. Write-out is the caller's job (see pipeline / cli).

ensure_downloads(studies, data_dir, timeout=600.0)

Fetch missing raw files (skip-if-exists, one retry) and return a manifest.

Source code in src/proteoform_regions/harmonize.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def ensure_downloads(
    studies: list[StudyConfig], data_dir: str | Path, timeout: float = 600.0
) -> pd.DataFrame:
    """Fetch missing raw files (skip-if-exists, one retry) and return a manifest."""
    manifest_rows: list[dict] = []
    for study in studies:
        for spec in study.files:
            if spec.url is None:
                continue
            dest = Path(data_dir) / study.dataset_accession / spec.name
            dest.parent.mkdir(parents=True, exist_ok=True)
            if not dest.exists():
                _fetch(spec.url, dest, timeout)
            manifest_rows.append(
                {
                    "dataset_accession": study.dataset_accession,
                    "file_name": spec.name,
                    "source_url": spec.url,
                    "bytes_on_disk": dest.stat().st_size if dest.exists() else 0,
                    "sha256": sha256_of(dest) if dest.exists() else None,
                }
            )
    return pd.DataFrame(manifest_rows)

harmonize(studies, data_dir, download=True, extraction_date=None)

Harmonize a cohort; returns (evidence DataFrame, run stats DataFrame).

The evidence frame carries exactly the 44 canonical columns; the stats frame records per-study row counts and adapter-reported counters.

Source code in src/proteoform_regions/harmonize.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def harmonize(
    studies: list[StudyConfig],
    data_dir: str | Path,
    download: bool = True,
    extraction_date: str | None = None,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Harmonize a cohort; returns (evidence DataFrame, run stats DataFrame).

    The evidence frame carries exactly the 44 canonical columns; the stats
    frame records per-study row counts and adapter-reported counters.
    """
    extraction_date = extraction_date or datetime.now().strftime("%Y-%m-%d")
    if download:
        ensure_downloads(studies, data_dir)
    all_rows: list[dict] = []
    stats_rows: list[dict] = []
    for study in studies:
        rows, stats = parse_study(study, data_dir)
        for row in rows:
            row["extraction_date"] = extraction_date
        all_rows.extend(rows)
        stats_rows.append(stats)
    evidence = pd.DataFrame(all_rows, columns=CANONICAL_COLUMNS)
    stats = pd.DataFrame(stats_rows)
    return evidence, stats

Mapping (stage 2)

proteoform_regions.mapping

Stage 2: UniProt accession mapping and protein properties.

Logic extracted from the source publication's S2 (uniprot_mapping.ipynb):

  • normalize_accessions sp|/tr| prefixes, semicolon groups, isoform suffixes, CON_/REV_ drops - verbatim
  • ensure_uniprot_entries chunked UniProtKB REST stream with per-chunk cache + dated snapshot, retry-once, 0.6 s politeness - verbatim semantics, snapshot path parameterized (no cwd-relative hack)
  • entry_properties / contaminant_flag protein metadata + the four contaminant rules (TRYPSIN/KERATIN/BSA/IG_CHAIN)
  • map_peptides exact-substring peptide->protein mapping with explicit failure tiers and NO silent drops (join-QC philosophy, review B9 part 1)

Contaminants are flagged, never dropped; degenerate peptides (one sequence mapping to >1 protein) are flagged per row for downstream ambiguity rules.

Verbatim-bug pin (divergence #32). entry_properties's reviewed flag reproduces the publication's entryType.endswith("Swiss-Prot") check, which is always False for real entries ("UniProtKB reviewed (Swiss-Prot)" ends with a paren). The frozen properties table therefore carries reviewed=False on all 5,777 rows; the package preserves this exactly for parity (pinned by test) - use entry_type directly if you need the review status.

contaminant_flag(accession, protein_name)

First matching contaminant rule over " ".

Source code in src/proteoform_regions/mapping.py 
193
194
195
196
197
198
199
def contaminant_flag(accession: str, protein_name) -> str | None:
    """First matching contaminant rule over "<accession> <protein_name>"."""
    text = f"{accession} {protein_name or ''}"
    for label, rule in CONTAMINANT_RULES:
        if rule.search(text):
            return label
    return None

ensure_uniprot_entries(base_accessions, cache_dir='data/raw/uniprot_cache', snapshot_date='20260917', chunk_size=100, politeness_seconds=0.6, offline=False)

Resolve base accessions against UniProtKB with cache + snapshot pinning.

Resolution order per chunk: dated snapshot -> per-chunk sha1-keyed cache -> live REST (skipped when offline). The merged snapshot is (re)written so later runs are network-free - the paper's reproducible-snapshot semantics.

Source code in src/proteoform_regions/mapping.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def ensure_uniprot_entries(
    base_accessions,
    cache_dir: str | Path = "data/raw/uniprot_cache",
    snapshot_date: str = "20260917",
    chunk_size: int = 100,
    politeness_seconds: float = 0.6,
    offline: bool = False,
) -> dict[str, dict]:
    """Resolve base accessions against UniProtKB with cache + snapshot pinning.

    Resolution order per chunk: dated snapshot -> per-chunk sha1-keyed cache ->
    live REST (skipped when ``offline``). The merged snapshot is (re)written so
    later runs are network-free - the paper's reproducible-snapshot semantics.
    """
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(parents=True, exist_ok=True)
    entries = load_snapshot(cache_dir, snapshot_date)
    ordered = sorted({a for a in base_accessions if a not in entries})
    chunks = [ordered[i : i + chunk_size] for i in range(0, len(ordered), chunk_size)]
    for chunk in chunks:
        key = hashlib.sha1(";".join(chunk).encode()).hexdigest()[:16]
        chunk_path = cache_dir / f"chunk_{key}.json"
        if chunk_path.exists():
            payload = json.loads(chunk_path.read_text(encoding="utf-8"))
        elif offline:
            continue
        else:
            payload = fetch_uniprot_chunk(chunk)
            if payload is None:
                continue
            chunk_path.write_text(json.dumps(payload), encoding="utf-8")
            time.sleep(politeness_seconds)
        for entry in payload:
            entries[entry["primaryAccession"]] = entry
    snapshot_path(cache_dir, snapshot_date).write_text(json.dumps(list(entries.values())), encoding="utf-8")
    return entries

entry_properties(entry)

Canonical sequence, gene/name, entryType, organism, version - verbatim S2.

Source code in src/proteoform_regions/mapping.py 
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
def entry_properties(entry: dict) -> dict:
    """Canonical sequence, gene/name, entryType, organism, version - verbatim S2."""
    sequence = entry.get("sequence", {}).get("value", "")
    description = entry.get("proteinDescription", {})
    name = (description.get("recommendedName") or {}).get("fullName", {}).get("value")
    if not name:
        names = description.get("submissionNames") or []
        name = names[0].get("fullName", {}).get("value") if names else None
    genes = [g.get("geneName", {}).get("value") for g in entry.get("genes", []) if g.get("geneName")]
    return {
        "sequence": sequence,
        "sequence_length": len(sequence),
        "gene_name": "; ".join(filter(None, genes)) or None,
        "protein_name": name,
        "entry_type": entry.get("entryType"),
        "organism": (entry.get("organism") or {}).get("scientificName"),
        "reviewed": entry.get("entryType", "").endswith("Swiss-Prot"),
        "sequence_version": entry.get("entryAudit", {}).get("sequenceVersion"),
    }

fetch_uniprot_chunk(accessions, timeout=120.0, retries=2)

GET one UniProtKB stream chunk (accession:A OR accession:B ...).

Retry-once; a hard failure returns None (documented, never fabricated).

Source code in src/proteoform_regions/mapping.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def fetch_uniprot_chunk(accessions: list[str], timeout: float = 120.0, retries: int = 2) -> list[dict] | None:
    """GET one UniProtKB stream chunk (``accession:A OR accession:B ...``).

    Retry-once; a hard failure returns ``None`` (documented, never fabricated).
    """
    query = " OR ".join(f"accession:{a}" for a in accessions)
    url = _UNIPROT_STREAM_URL + "?" + urllib.parse.urlencode({"query": query, "format": "json"})
    headers = {"Accept": "application/json", "User-Agent": _USER_AGENT}
    for attempt in range(1, retries + 1):
        try:
            with urllib.request.urlopen(
                urllib.request.Request(url, headers=headers), timeout=timeout
            ) as response:
                return json.loads(response.read().decode("utf-8")).get("results", [])
        except Exception:
            if attempt == retries:
                print(f"UNIPROT CHUNK FAILED after retry (documented): {len(accessions)} accessions")
                return None
            time.sleep(2.0)
    return None

load_snapshot(cache_dir, snapshot_date)

Load a dated snapshot as {primaryAccession: entry}; {} when absent.

Source code in src/proteoform_regions/mapping.py 
126
127
128
129
130
131
def load_snapshot(cache_dir: str | Path, snapshot_date: str) -> dict[str, dict]:
    """Load a dated snapshot as {primaryAccession: entry}; {} when absent."""
    path = snapshot_path(cache_dir, snapshot_date)
    if not path.exists():
        return {}
    return {e["primaryAccession"]: e for e in json.loads(path.read_text(encoding="utf-8"))}

map_peptides(evidence, entries)

Exact-substring peptide->protein mapping with explicit tiers (no silent drops).

Tiers: exact_unique / exact_multi_position on success; failures are no_accession_in_row / accession_unmapped_in_uniprot / not_found_in_group_sequences. Also flags isoform-suffixed rows, contaminant hits, and degenerate peptides (sequence mapping to >1 protein).

Source code in src/proteoform_regions/mapping.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
def map_peptides(evidence: pd.DataFrame, entries: dict[str, dict]) -> pd.DataFrame:
    """Exact-substring peptide->protein mapping with explicit tiers (no silent drops).

    Tiers: ``exact_unique`` / ``exact_multi_position`` on success; failures are
    ``no_accession_in_row`` / ``accession_unmapped_in_uniprot`` /
    ``not_found_in_group_sequences``. Also flags isoform-suffixed rows,
    contaminant hits, and degenerate peptides (sequence mapping to >1 protein).
    """
    prop_by_accession = {}
    for accession, entry in entries.items():
        props = entry_properties(entry)
        prop_by_accession[accession] = {
            "sequence": props["sequence"],
            "contaminant": contaminant_flag(accession, props["protein_name"]),
        }

    map_rows = []
    peptide_proteins: dict[str, set] = {}
    for row_index, record in enumerate(evidence.itertuples(index=False)):
        peptide = record.peptide_sequence
        tokens = normalize_accessions(record.protein_accession_raw)
        mapped_accession, positions, tier, flag = None, [], None, None
        resolved_any = False
        for token, _isoform in tokens:
            props = prop_by_accession.get(token)
            if props is None:
                continue
            resolved_any = True
            flag = flag or props["contaminant"]
            if peptide and peptide in props["sequence"]:
                found = [m.start() + 1 for m in re.finditer(f"(?={re.escape(peptide)})", props["sequence"])]
                if found:
                    mapped_accession, positions = token, found
                    tier = "exact_multi_position" if len(found) > 1 else "exact_unique"
                    break
        if tier is None:
            tier = (
                "no_accession_in_row"
                if not tokens
                else ("accession_unmapped_in_uniprot" if not resolved_any else "not_found_in_group_sequences")
            )
        if peptide and mapped_accession:
            peptide_proteins.setdefault(peptide, set()).add(mapped_accession)
        map_rows.append(
            {
                "peptide_row_id": row_index,
                "dataset_accession": record.dataset_accession,
                "peptide_sequence": peptide,
                "mapped_accession": mapped_accession,
                "positions": ";".join(str(p) for p in positions) if positions else None,
                "mapping_tier": tier,
                "row_had_isoform_suffix": any(i for _, i in tokens),
                "contaminant_flag": flag,
            }
        )
    mapping = pd.DataFrame(map_rows)
    mapping["degenerate_peptide"] = mapping["peptide_sequence"].isin(
        [p for p, accs in peptide_proteins.items() if len(accs) > 1]
    )
    return mapping

map_to_uniprot(evidence, cache_dir='data/raw/uniprot_cache', snapshot_date='20260917', offline=False)

Full S2 stage: normalize -> resolve -> map. Returns (mapping, proteins).

proteins carries per-protein properties (precursor pI/MW via the same Biopython engine as peptide features); mapping carries per-row tiers, contaminant and degeneracy flags.

Source code in src/proteoform_regions/mapping.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
def map_to_uniprot(
    evidence: pd.DataFrame,
    cache_dir: str | Path = "data/raw/uniprot_cache",
    snapshot_date: str = "20260917",
    offline: bool = False,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Full S2 stage: normalize -> resolve -> map. Returns (mapping, proteins).

    ``proteins`` carries per-protein properties (precursor pI/MW via the same
    Biopython engine as peptide features); ``mapping`` carries per-row tiers,
    contaminant and degeneracy flags.
    """
    if "protein_accession_raw" not in evidence.columns:
        raise ValueError("evidence table lacks protein_accession_raw (not a canonical S1 table?)")
    all_bases: set[str] = set()
    for raw in evidence["protein_accession_raw"]:
        all_bases.update(token for token, _ in normalize_accessions(raw))
    entries = ensure_uniprot_entries(all_bases, cache_dir, snapshot_date, offline=offline)
    mapping = map_peptides(evidence, entries)
    proteins = protein_properties_table(entries)
    return mapping, proteins

normalize_accessions(raw_value)

Split a raw accession field into [(base_accession, isoform_suffix), ...].

Handles fasta sp|P02768|ALBU_HUMAN prefixes, ;/, separated groups, -2 isoform suffixes, and drops CON_/REV_ tokens. Verbatim from the publication's S2 (cell 4).

Source code in src/proteoform_regions/mapping.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def normalize_accessions(raw_value) -> list[tuple[str, str | None]]:
    """Split a raw accession field into [(base_accession, isoform_suffix), ...].

    Handles fasta ``sp|P02768|ALBU_HUMAN`` prefixes, ``;``/``,``
    separated groups, ``-2`` isoform suffixes, and drops ``CON_``/``REV_``
    tokens. Verbatim from the publication's S2 (cell 4).
    """
    if raw_value is None or (isinstance(raw_value, float) and np.isnan(raw_value)):
        return []
    tokens: list[tuple[str, str | None]] = []
    for token in str(raw_value).replace(",", ";").split(";"):
        token = token.strip()
        if not token:
            continue
        if "|" in token:
            parts = token.split("|")
            token = parts[1] if len(parts) >= 2 else parts[-1]
        if token.upper().startswith(("CON_", "REV_")):
            continue
        isoform = None
        match = ISOFORM_SUFFIX.search(token)
        if match and token[: match.start()].upper().startswith(("P", "Q", "O", "A")):
            isoform = token[match.start() :]
            token = token[: match.start()]
        tokens.append((token, isoform))
    return tokens

protein_properties_table(entries)

Per-protein properties + contaminant flag + precursor pI/MW (Biopython).

Source code in src/proteoform_regions/mapping.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def protein_properties_table(entries: dict[str, dict]) -> pd.DataFrame:
    """Per-protein properties + contaminant flag + precursor pI/MW (Biopython)."""
    rows = []
    for accession in sorted(entries):
        props = entry_properties(entries[accession])
        rows.append(
            {
                "accession": accession,
                "contaminant_flag": contaminant_flag(accession, props["protein_name"]),
                **props,
            }
        )
    proteins = pd.DataFrame(rows)
    proteins["precursor_pi"] = [biopython_pi(seq) if seq else None for seq in proteins["sequence"]]
    proteins["precursor_mw_da"] = [biopython_mw(seq) if seq else None for seq in proteins["sequence"]]
    return proteins

Pipeline stages (3a/3b/4/8 + orchestration)

proteoform_regions.pipeline

End-to-end orchestration: harmonize -> map -> compute -> stitch -> bias.

Each stage function is a thin, faithful composition of the stage modules (the publication's S2-S4, S8 notebook flows), so run() reproduces the paper's pipeline as one call:

  • map_to_uniprot (stage 2) accession mapping + protein properties
  • compute_features (stage 3a) peptide physicochemical feature table
  • build_background (stage 3b) digestion-matched tryptic background
  • stitch_regions (stage 4) single-linkage detected-region inference
  • quantify_bias (stage 8) RQ1 permutation + RQ2 OR meta-analysis
  • run stage-selected execution + artifact writing + sha256 run manifest

Every artifact filename matches the paper's data/processed/ vocabulary; every stochastic knob (seed, B_perm, B_boot) is a parameter with the paper's defaults.

RunReport(out_dir, stages, artifacts=dict(), counts=dict(), parameters=dict(), generated='', proteoform_regions_version='') dataclass

Run provenance: executed stages, artifacts, row counts, parameters.

build_background(protein_sequences, max_missed=2, min_len=7, max_len=60)

Digestion-matched tryptic background over {accession: sequence} (S3).

Same digest protocol and same physicochemical columns as the detected peptides, so the RQ2 comparison is digestion-matched by construction.

Source code in src/proteoform_regions/pipeline.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def build_background(
    protein_sequences: dict[str, str],
    max_missed: int = 2,
    min_len: int = 7,
    max_len: int = 60,
) -> pd.DataFrame:
    """Digestion-matched tryptic background over {accession: sequence} (S3).

    Same digest protocol and same physicochemical columns as the detected
    peptides, so the RQ2 comparison is digestion-matched by construction.
    """
    digested = background_peptides(protein_sequences, max_missed, min_len, max_len)
    pi_cache: dict[str, tuple] = {}

    def props(seq: str) -> tuple:
        if seq not in pi_cache:
            pi_cache[seq] = (
                biopython_pi(seq),
                biopython_mw(seq),
                gravy(seq),
                peptide_features(seq),
            )
        return pi_cache[seq]

    rows = []
    for rec in digested:
        p_i, mw, gr, feats = props(rec["peptide_sequence"])
        rows.append(
            {
                "accession": rec["accession"],
                "peptide_sequence": rec["peptide_sequence"],
                "missed_cleavages": rec["missed_cleavages"],
                "peptide_pi": p_i,
                "peptide_mw_da": mw,
                "peptide_gravy": gr,
                "peptide_class_anl003": classify_pi(p_i),
                **feats,
            }
        )
    return pd.DataFrame(rows)

compute_features(evidence, mapping, proteins, calculator='bjellqvist')

Peptide feature table over exactly-mapped rows (S3 feature engine).

calculator selects the peptide pI engine: bjellqvist (Biopython, the paper's primary) or an HH table name (hh_rodwell_emboss_lineage / hh_lehninger) for the named sensitivity analysis. Net charge at pH 7.4 always uses the Rodwell/EMBOSS HH table, as in the source.

Source code in src/proteoform_regions/pipeline.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def compute_features(
    evidence: pd.DataFrame,
    mapping: pd.DataFrame,
    proteins: pd.DataFrame,
    calculator: str = "bjellqvist",
) -> pd.DataFrame:
    """Peptide feature table over exactly-mapped rows (S3 feature engine).

    ``calculator`` selects the peptide pI engine: ``bjellqvist`` (Biopython,
    the paper's primary) or an HH table name (``hh_rodwell_emboss_lineage`` /
    ``hh_lehninger``) for the named sensitivity analysis. Net charge at pH 7.4
    always uses the Rodwell/EMBOSS HH table, as in the source.
    """
    mapped = mapping[mapping["mapping_tier"].isin(_EXACT_TIERS)].copy()
    protein_lookup = proteins.set_index("accession")

    def _pi_map(seqs: list[str]) -> dict[str, float]:
        if calculator == "bjellqvist":
            return {s: biopython_pi(s) for s in seqs}
        if calculator in ("hh_rodwell_emboss_lineage", "hh_lehninger"):
            return dict(zip(seqs, hh_pi_batch(seqs, calculator), strict=True))
        raise ValueError(f"unknown calculator {calculator!r}")

    seqs_needed = pd.unique(mapped["peptide_sequence"]).tolist()
    pi_map = _pi_map([str(s) for s in seqs_needed])
    pi_cache: dict[str, tuple] = {}
    for seq in seqs_needed:
        pi_cache[seq] = (
            pi_map[seq],
            biopython_mw(seq),
            gravy(seq),
            peptide_features(seq),
        )

    feature_rows = []
    for record in mapped.itertuples(index=False):
        seq = record.peptide_sequence
        p_i, mw, gr, feats = pi_cache[seq]
        precursor = (
            protein_lookup.loc[record.mapped_accession]
            if record.mapped_accession in protein_lookup.index
            else None
        )
        row = {
            "peptide_row_id": record.peptide_row_id,
            "dataset_accession": record.dataset_accession,
            "peptide_sequence": seq,
            "mapped_accession": record.mapped_accession,
            "mapping_tier": record.mapping_tier,
            "peptide_pi": p_i,
            "peptide_mw_da": mw,
            "peptide_gravy": gr,
            "contaminant_flag": record.contaminant_flag,
            **feats,
            "peptide_class_anl003": classify_pi(p_i),
        }
        if precursor is not None:
            row.update(
                {
                    "precursor_pi": precursor["precursor_pi"],
                    "precursor_mw_da": precursor["precursor_mw_da"],
                    "precursor_length": precursor["sequence_length"],
                    "gene_name": precursor["gene_name"],
                    "protein_name": precursor["protein_name"],
                    "pi_divergence_peptide_minus_precursor": p_i - precursor["precursor_pi"],
                    "mw_ratio_peptide_over_precursor": (
                        mw / precursor["precursor_mw_da"] if precursor["precursor_mw_da"] else None
                    ),
                }
            )
        feature_rows.append(row)
    features = pd.DataFrame(feature_rows)

    # net charge at pH 7.4 (HH Rodwell lineage; batched as in the source)
    unique_seqs = [str(s) for s in features["peptide_sequence"].unique()]
    features["net_charge_ph74"] = features["peptide_sequence"].map(charge_at_ph(unique_seqs, 7.4))
    features["charge_density"] = features["net_charge_ph74"] / features["peptide_length"]
    return features

quantify_bias(features, background, seed=st.SEED, b_perm=st.B_PERM, b_boot=st.B_BOOT, run_rq1=True)

Detection-bias quantification on feature + background tables (S8).

RQ2: per-study basic-vs-rest OR (Woolf, unconditional 0.5) against each study's digestion-matched parent-protein background; inverse-variance fixed effect; DerSimonian-Laird random effects; leave-one-study-out envelope. RQ1: per-protein divergence permutation null (seed-aware).

Set run_rq1=False to skip the permutation null (it dominates runtime on paper-scale tables: ~5,400 proteins x B=2,000 draws).

Source code in src/proteoform_regions/pipeline.py 
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
def quantify_bias(
    features: pd.DataFrame,
    background: pd.DataFrame,
    seed: int = st.SEED,
    b_perm: int = st.B_PERM,
    b_boot: int = st.B_BOOT,
    run_rq1: bool = True,
) -> dict:
    """Detection-bias quantification on feature + background tables (S8).

    RQ2: per-study basic-vs-rest OR (Woolf, unconditional 0.5) against each
    study's digestion-matched parent-protein background; inverse-variance
    fixed effect; DerSimonian-Laird random effects; leave-one-study-out
    envelope. RQ1: per-protein divergence permutation null (seed-aware).

    Set ``run_rq1=False`` to skip the permutation null (it dominates runtime
    on paper-scale tables: ~5,400 proteins x B=2,000 draws).
    """
    det = features.drop_duplicates(["dataset_accession", "peptide_sequence"])
    study_parents = (
        features.groupby("dataset_accession")["mapped_accession"].agg(lambda s: set(s.dropna())).to_dict()
    )

    study_rows = []
    for study in sorted(det["dataset_accession"].unique()):
        d_cls = det[det["dataset_accession"] == study]["peptide_class_anl003"].dropna()
        bg = background[background["accession"].isin(study_parents[study])].drop_duplicates(
            "peptide_sequence"
        )
        b_cls = bg["peptide_class_anl003"].dropna()
        r = st.study_or(d_cls.to_numpy(), b_cls.to_numpy(), "basic")
        r["dataset_accession"] = study
        study_rows.append(r)
    or_table = pd.DataFrame(study_rows)

    mu_fe, se_fe, p_fe, Q, I2 = st.fe_summary(or_table)
    dl = st.dl_random_effects(or_table)
    loo_lo, loo_hi = st.loo_envelope(or_table)

    result = {
        "k": int(len(or_table)),
        "or_table": or_table,
        "fe": {
            "log_or": mu_fe,
            "se": se_fe,
            "or": float(np.exp(mu_fe)),
            "ci_low": float(np.exp(mu_fe - 1.96 * se_fe)),
            "ci_high": float(np.exp(mu_fe + 1.96 * se_fe)),
            "p": p_fe,
            "Q": Q,
            "I2": I2,
        },
        "random_effects": dl,
        "loo_envelope": [loo_lo, loo_hi],
    }
    if run_rq1:
        result["rq1"] = st.rq1_permutation(features, background, seed=seed, b_perm=b_perm, b_boot=b_boot)
    return result

run(studies, out_dir, stages='all', data_dir='data/raw/peptide_tables', download=True, uniprot_cache='data/raw/uniprot_cache', snapshot_date='20260917', offline=False, calculator='bjellqvist', max_missed=2, gap=REGION_GAP_AA, terminal_fraction=TERMINAL_FRACTION, min_peptides=MIN_REGION_PEPTIDES, seed=st.SEED, b_perm=st.B_PERM, b_boot=st.B_BOOT, run_rq1=True)

Execute pipeline stages over a cohort, writing paper-vocabulary artifacts.

studies is a cohort YAML path or a list of StudyConfig. Artifacts land in out_dir under the paper's data/processed/ filenames; every run closes by writing a sha256 manifest (manifest.sha256) over its own outputs plus a machine-readable run_report.json.

Source code in src/proteoform_regions/pipeline.py 
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
def run(
    studies: str | Path | list[StudyConfig],
    out_dir: str | Path,
    stages: str | list[str] = "all",
    data_dir: str | Path = "data/raw/peptide_tables",
    download: bool = True,
    uniprot_cache: str | Path = "data/raw/uniprot_cache",
    snapshot_date: str = "20260917",
    offline: bool = False,
    calculator: str = "bjellqvist",
    max_missed: int = 2,
    gap: int = REGION_GAP_AA,
    terminal_fraction: float = TERMINAL_FRACTION,
    min_peptides: int = MIN_REGION_PEPTIDES,
    seed: int = st.SEED,
    b_perm: int = st.B_PERM,
    b_boot: int = st.B_BOOT,
    run_rq1: bool = True,
) -> RunReport:
    """Execute pipeline stages over a cohort, writing paper-vocabulary artifacts.

    ``studies`` is a cohort YAML path or a list of ``StudyConfig``. Artifacts
    land in ``out_dir`` under the paper's ``data/processed/`` filenames; every
    run closes by writing a sha256 manifest (``manifest.sha256``) over its own
    outputs plus a machine-readable ``run_report.json``.
    """
    import proteoform_regions

    if isinstance(studies, (str, Path)):
        studies = load_cohort(studies)
    selected = _selected_stages(stages)
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)

    report = RunReport(
        out_dir=out_dir,
        stages=selected,
        parameters={
            "calculator": calculator,
            "max_missed": max_missed,
            "gap": gap,
            "terminal_fraction": terminal_fraction,
            "min_peptides": min_peptides,
            "seed": seed,
            "b_perm": b_perm,
            "b_boot": b_boot,
            "snapshot_date": snapshot_date,
            "n_studies": len(studies),
        },
        generated=datetime.now().isoformat(timespec="minutes"),
        proteoform_regions_version=proteoform_regions.__version__,
    )
    evidence = mapping_df = features = background = regions = None
    proteins = None

    def _write(frame: pd.DataFrame, name: str, stage: str) -> Path:
        path = out_dir / name
        frame.to_csv(path, index=False)
        report.artifacts[stage] = name
        report.counts[stage] = int(len(frame))
        return path

    if "harmonize" in selected:
        evidence, _stats = harmonize(studies, data_dir, download=download)
        _write(evidence, "peptide_evidence_table.csv", "harmonize")
    elif "map" in selected or "compute" in selected or "stitch" in selected:
        evidence = pd.read_csv(out_dir / "peptide_evidence_table.csv", low_memory=False)

    if "map" in selected:
        mapping_df, proteins = _map_stage(evidence, uniprot_cache, snapshot_date, offline=offline)
        _write(mapping_df, "processed_peptide_protein_map.csv", "map")
        proteins_out = proteins.drop(columns=["sequence"], errors="ignore")
        _write(proteins_out, "processed_protein_properties.csv", "map_properties")
    elif "compute" in selected or "stitch" in selected:
        mapping_df = pd.read_csv(out_dir / "processed_peptide_protein_map.csv")

    if "compute" in selected:
        if proteins is None:
            proteins = pd.read_csv(out_dir / "processed_protein_properties.csv", low_memory=False)
        features = compute_features(evidence, mapping_df, proteins, calculator=calculator)
        _write(features, "processed_feature_table.csv", "compute")
    elif "stitch" in selected or "bias" in selected:
        features = pd.read_csv(out_dir / "processed_feature_table.csv", low_memory=False)

    if "background" in selected:
        if features is None:
            features = pd.read_csv(out_dir / "processed_feature_table.csv", low_memory=False)
        parents = set(features["mapped_accession"].dropna())
        sequences = _sequences_from_snapshot(parents, uniprot_cache, snapshot_date)
        background = build_background(sequences, max_missed=max_missed)
        _write(background, "processed_tryptic_background.csv", "background")
    elif "bias" in selected:
        background = pd.read_csv(out_dir / "processed_tryptic_background.csv")

    if "stitch" in selected:
        parents = set(mapping_df["mapped_accession"].dropna())
        precursor_sequences = _sequences_from_snapshot(parents, uniprot_cache, snapshot_date)
        regions = stitch_regions(
            features,
            mapping_df,
            precursor_sequences,
            gap=gap,
            terminal_fraction=terminal_fraction,
            min_peptides=min_peptides,
        )
        _write(regions, "processed_region_inference.csv", "stitch")

    if "bias" in selected:
        bias = quantify_bias(features, background, seed=seed, b_perm=b_perm, b_boot=b_boot, run_rq1=run_rq1)
        bias_payload = {
            "k": bias["k"],
            "fe": bias["fe"],
            "random_effects": bias["random_effects"],
            "loo_envelope": bias["loo_envelope"],
            "or_table": bias["or_table"].to_dict(orient="records"),
            "rq1": {k: v for k, v in bias.get("rq1", {}).items() if k != "per_protein"},
            "parameters": {"seed": seed, "b_perm": b_perm, "b_boot": b_boot},
        }
        (out_dir / "processed_inference_results.json").write_text(
            json.dumps(bias_payload, indent=2, default=str), encoding="utf-8"
        )
        report.artifacts["bias"] = "processed_inference_results.json"
        report.counts["bias"] = bias["k"]

    (out_dir / "run_report.json").write_text(json.dumps(report.to_dict(), indent=2), encoding="utf-8")
    outputs = [out_dir / name for name in report.artifacts.values()]
    if outputs:
        write_manifest(outputs, out_dir / "manifest.sha256")
    report.artifacts["run_report"] = "run_report.json"
    report.artifacts["manifest"] = "manifest.sha256"
    return report

stitch_regions(features, mapping, precursor_sequences=None, gap=REGION_GAP_AA, terminal_fraction=TERMINAL_FRACTION, min_peptides=MIN_REGION_PEPTIDES)

Detected-region inference (S4 walk): position spans -> stitched regions.

Start positions come from the mapping table (first exact-substring hit); rows without exact mapping never reach this stage. precursor_sequences (e.g. from the UniProt snapshot) enables region sequence / region pI/MW columns; without it they stay null (positions and tiers are unaffected).

Source code in src/proteoform_regions/pipeline.py 
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
def stitch_regions(
    features: pd.DataFrame,
    mapping: pd.DataFrame,
    precursor_sequences: dict[str, str] | None = None,
    gap: int = REGION_GAP_AA,
    terminal_fraction: float = TERMINAL_FRACTION,
    min_peptides: int = MIN_REGION_PEPTIDES,
) -> pd.DataFrame:
    """Detected-region inference (S4 walk): position spans -> stitched regions.

    Start positions come from the mapping table (first exact-substring hit);
    rows without exact mapping never reach this stage. ``precursor_sequences``
    (e.g. from the UniProt snapshot) enables region sequence / region pI/MW
    columns; without it they stay null (positions and tiers are unaffected).
    """
    precursor_sequences = precursor_sequences or {}
    mapped = mapping[mapping["mapping_tier"].isin(_EXACT_TIERS)].copy()
    if mapped.empty:
        return pd.DataFrame()
    mapped["start_position"] = mapped["positions"].astype(str).str.split(";").str[0].astype(float)
    mapped = mapped[mapped["start_position"].notna()]
    merged = mapped.merge(
        features[
            [
                "peptide_row_id",
                "peptide_pi",
                "peptide_length",
                "precursor_pi",
                "precursor_mw_da",
                "precursor_length",
                "gene_name",
            ]
        ],
        on="peptide_row_id",
        how="left",
        suffixes=("", "_feat"),
    )
    merged["end_position"] = merged["start_position"] + merged["peptide_length"] - 1

    region_rows = []
    for (study, accession), group in merged.groupby(["dataset_accession", "mapped_accession"]):
        precursor_length = group["precursor_length"].iloc[0]
        if not precursor_length or precursor_length != precursor_length:
            continue
        spans = list(
            zip(
                group["start_position"],
                group["end_position"],
                group["peptide_sequence"],
                group["peptide_pi"],
                group["peptide_row_id"],
                strict=True,
            )
        )
        for cluster in stitch_spans(spans, gap):
            region_rows.append(
                region_record(
                    study,
                    accession,
                    cluster,
                    precursor_length=int(precursor_length),
                    precursor_pi=group["precursor_pi"].iloc[0],
                    precursor_mw=group["precursor_mw_da"].iloc[0],
                    gene_name=group["gene_name"].iloc[0],
                    precursor_sequence=precursor_sequences.get(accession, ""),
                    gap=gap,
                    terminal_fraction=terminal_fraction,
                    min_peptides=min_peptides,
                )
            )
    return pd.DataFrame(region_rows)

Statistics

proteoform_regions.stats

Detection-bias statistical models (Supplementary Methods Modules 4-6).

Pure, seed-aware functions extracted verbatim from the source publication's S5 (descriptive), S8 (inferential), and charge-state (J4) layers:

  • bh_q Benjamini-Hochberg q-values
  • or_2x2 / study_or 2x2 OR with Woolf SE (UNCONDITIONAL max(x, 0.5) zero-cell correction on every cell)
  • fe_summary inverse-variance fixed effect + Q + I2
  • dl_random_effects DerSimonian-Laird tau2 / RE summary
  • loo_envelope leave-one-study-out FE envelope
  • rq1_permutation digestion-matched per-protein permutation null (WITH replacement, B=2000, seed 42; two-stage global p via resampling null medians)
  • bootstrap_class_shares / enrichment_with_ci descriptive CIs
  • theoretical_charge_ph26 KRH+N-term protonation proxy (charge-state tests)

bh_q(pvals)

Benjamini-Hochberg q-values (step-up, monotone via reverse cummin).

Source code in src/proteoform_regions/stats.py 
43
44
45
46
47
48
49
50
51
52
53
def bh_q(pvals) -> np.ndarray:
    """Benjamini-Hochberg q-values (step-up, monotone via reverse cummin)."""
    pvals = np.asarray(pvals, dtype=float)
    order = np.argsort(pvals)
    ranked = pvals[order]
    m = len(ranked)
    q = ranked * m / (np.arange(m) + 1)
    q = np.minimum.accumulate(q[::-1])[::-1]
    out = np.empty_like(q)
    out[order] = np.minimum(q, 1.0)
    return out

bootstrap_class_shares(labels, n=2000, seed=SEED)

Bootstrap percentile CIs for class shares of a label array (S5).

Source code in src/proteoform_regions/stats.py 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def bootstrap_class_shares(labels, n: int = 2000, seed: int = SEED) -> tuple[dict, dict]:
    """Bootstrap percentile CIs for class shares of a label array (S5)."""
    rng = np.random.default_rng(seed)
    labels = np.asarray(labels)
    n_total = len(labels)
    masks = {c: labels == c for c in CLASSES}
    shares = {c: float(masks[c].mean()) if n_total else 0.0 for c in CLASSES}
    cis = {c: [0.0, 0.0] for c in CLASSES}
    if n_total == 0:
        return shares, cis
    draws = np.empty((n, len(CLASSES)))
    for i in range(n):
        sample = labels[rng.integers(0, n_total, n_total)]
        draws[i] = [np.mean(sample == c) for c in CLASSES]
    for j, c in enumerate(CLASSES):
        cis[c] = [float(np.percentile(draws[:, j], 2.5)), float(np.percentile(draws[:, j], 97.5))]
    return shares, cis

dl_random_effects(table)

DerSimonian-Laird random-effects summary (reported at k=10 per freeze spec).

Source code in src/proteoform_regions/stats.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
def dl_random_effects(table: pd.DataFrame) -> dict:
    """DerSimonian-Laird random-effects summary (reported at k=10 per freeze spec)."""
    w_inv = 1 / table["se"] ** 2
    mu_fe, _, _, Q, _ = fe_summary(table)
    k = len(table)
    tau2 = (
        max(0.0, (Q - (k - 1)) / (np.sum(w_inv) - np.sum(w_inv**2) / np.sum(w_inv))) if (k - 1) < Q else 0.0
    )
    w_re = 1 / (table["se"] ** 2 + tau2)
    mu_re = float(np.sum(w_re * table["log_or"]) / np.sum(w_re))
    se_re = float(np.sqrt(1 / np.sum(w_re)))
    p_re = float(2 * stats.norm.sf(abs(mu_re / se_re)))
    return dict(tau2=float(tau2), mu=float(mu_re), se=se_re, or_=float(np.exp(mu_re)), p=p_re)

enrichment_with_ci(detected_labels, background_labels, n=2000, seed=SEED)

Enrichment (detected share / background share) with multinomial bootstrap CIs (S5).

Source code in src/proteoform_regions/stats.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
def enrichment_with_ci(detected_labels, background_labels, n: int = 2000, seed: int = SEED) -> dict:
    """Enrichment (detected share / background share) with multinomial bootstrap CIs (S5)."""
    rng = np.random.default_rng(seed)
    d_arr = np.asarray(detected_labels)
    b_arr = np.asarray(background_labels)
    d_counts = np.array([(d_arr == c).sum() for c in CLASSES], dtype=float)
    b_counts = np.array([(b_arr == c).sum() for c in CLASSES], dtype=float)
    d_shares = d_counts / d_counts.sum() if d_counts.sum() else np.zeros(3)
    b_shares = b_counts / b_counts.sum() if b_counts.sum() else np.zeros(3)
    if d_counts.sum() == 0 or b_counts.sum() == 0:
        return {c: dict(enrichment=None, ci=[None, None]) for c in CLASSES}
    d_draws = rng.multinomial(int(d_counts.sum()), d_shares, size=n) / d_counts.sum()
    b_draws = rng.multinomial(int(b_counts.sum()), b_shares, size=n) / b_counts.sum()
    out = {}
    for j, c in enumerate(CLASSES):
        ratio = d_draws[:, j] / np.maximum(b_draws[:, j], 1e-12)
        out[c] = dict(
            enrichment=float(d_shares[j] / b_shares[j]) if b_shares[j] else None,
            ci=[float(np.percentile(ratio, 2.5)), float(np.percentile(ratio, 97.5))],
        )
    return out

fe_summary(table)

Inverse-variance fixed effect over a per-study OR table.

Requires columns log_or and se. Returns (mu, se, p, Q, I2).

Source code in src/proteoform_regions/stats.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def fe_summary(table: pd.DataFrame) -> tuple[float, float, float, float, float]:
    """Inverse-variance fixed effect over a per-study OR table.

    Requires columns ``log_or`` and ``se``. Returns (mu, se, p, Q, I2).
    """
    w = 1 / table["se"] ** 2
    mu = float(np.sum(w * table["log_or"]) / np.sum(w))
    se = float(np.sqrt(1 / np.sum(w)))
    z = mu / se
    p = float(2 * stats.norm.sf(abs(z)))
    Q = float(np.sum(w * (table["log_or"] - mu) ** 2))
    df = len(table) - 1
    I2 = max(0.0, (Q - df) / Q) if Q > 0 else 0.0
    return mu, se, p, Q, I2

loo_envelope(table, by='dataset_accession')

Leave-one-study-out FE summary ORs -> (min, max) envelope (OR scale).

Source code in src/proteoform_regions/stats.py 
123
124
125
126
127
128
129
130
def loo_envelope(table: pd.DataFrame, by: str = "dataset_accession") -> tuple[float, float]:
    """Leave-one-study-out FE summary ORs -> (min, max) envelope (OR scale)."""
    mus = []
    for key in table[by]:
        sub = table[table[by] != key]
        mu_l, _, _, _, _ = fe_summary(sub)
        mus.append(float(np.exp(mu_l)))
    return float(min(mus)), float(max(mus))

or_2x2(a, b, c, d)

(log OR, Woolf SE) for detected-vs-background x basic-vs-rest.

a=detect&basic, b=detect&rest, c=bg&basic, d=bg&rest. Zero-cell correction is UNCONDITIONAL: max(x, 0.5) on every cell always (always-on Haldane-Anscombe; simpler than the conditional variant).

Source code in src/proteoform_regions/stats.py 
56
57
58
59
60
61
62
63
64
65
66
def or_2x2(a: float, b: float, c: float, d: float) -> tuple[float, float]:
    """(log OR, Woolf SE) for detected-vs-background x basic-vs-rest.

    ``a``=detect&basic, ``b``=detect&rest, ``c``=bg&basic, ``d``=bg&rest.
    Zero-cell correction is UNCONDITIONAL: max(x, 0.5) on every cell always
    (always-on Haldane-Anscombe; simpler than the conditional variant).
    """
    a, b, c, d = max(a, 0.5), max(b, 0.5), max(c, 0.5), max(d, 0.5)
    lor = np.log((a / b) / (c / d))
    se = np.sqrt(1 / a + 1 / b + 1 / c + 1 / d)
    return lor, se

rq1_permutation(features, background, seed=SEED, b_perm=B_PERM, b_boot=B_BOOT, min_unique=2)

H1/RQ1: per-protein median peptide-pI divergence vs digestion-matched null.

features needs (dataset_accession, mapped_accession, peptide_sequence, peptide_pi, precursor_pi); background needs (accession, peptide_pi). Sampling is WITH replacement from each protein's tryptic background (multinomial bootstrap); the global null resamples the SET of per-protein null medians. Returns observed statistic, two-sided global p, 95% null CI, and the per-protein table with BH q.

Edge case (published): proteins with fewer background peptides than detected peptides are SKIPPED (documented absence, not corrected).

Source code in src/proteoform_regions/stats.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def rq1_permutation(
    features: pd.DataFrame,
    background: pd.DataFrame,
    seed: int = SEED,
    b_perm: int = B_PERM,
    b_boot: int = B_BOOT,
    min_unique: int = 2,
) -> dict:
    """H1/RQ1: per-protein median peptide-pI divergence vs digestion-matched null.

    ``features`` needs (dataset_accession, mapped_accession, peptide_sequence,
    peptide_pi, precursor_pi); ``background`` needs (accession, peptide_pi).
    Sampling is WITH replacement from each protein's tryptic background
    (multinomial bootstrap); the global null resamples the SET of per-protein
    null medians. Returns observed statistic, two-sided global p, 95% null CI,
    and the per-protein table with BH q.

    Edge case (published): proteins with fewer background peptides than
    detected peptides are SKIPPED (documented absence, not corrected).
    """
    rng = np.random.default_rng(seed)
    per_protein = (
        features.dropna(subset=["precursor_pi"])
        .groupby(["dataset_accession", "mapped_accession"])
        .agg(
            n_unique=("peptide_sequence", "nunique"),
            median_peptide_pi=("peptide_pi", "median"),
            precursor_pi=("precursor_pi", "first"),
        )
        .reset_index()
    )
    per_protein = per_protein[per_protein["n_unique"] >= min_unique]
    bg_by_accession = {acc: g["peptide_pi"].dropna().to_numpy() for acc, g in background.groupby("accession")}

    obs_divs, null_divs_global, per_protein_p = [], [], []
    for record in per_protein.itertuples(index=False):
        bg = bg_by_accession.get(record.mapped_accession)
        if bg is None or len(bg) < record.n_unique:
            continue
        obs = record.median_peptide_pi - record.precursor_pi
        obs_divs.append(obs)
        draws = bg[rng.integers(0, len(bg), size=(b_perm, record.n_unique))]
        null_medians = np.median(draws, axis=1)
        null_div = null_medians - record.precursor_pi
        null_divs_global.append(np.median(null_div))
        per_protein_p.append((1 + np.sum(np.abs(null_div) >= abs(obs))) / (b_perm + 1))

    obs_divs = np.array(obs_divs)
    T_obs = float(np.median(obs_divs))
    global_null = np.array(
        [
            np.median(rng.choice(null_divs_global, size=len(null_divs_global), replace=True))
            for _ in range(b_boot)
        ]
    )
    p_global = (
        1
        + np.sum(
            np.abs(global_null - np.median(null_divs_global)) >= abs(T_obs - np.median(null_divs_global))
        )
    ) / (b_boot + 1)
    per_protein_valid = per_protein.iloc[: len(per_protein_p)].copy()
    per_protein_valid["perm_p"] = per_protein_p
    per_protein_valid["perm_q"] = bh_q(per_protein_p)
    return dict(
        t_obs=T_obs,
        p_global=float(p_global),
        ci_low=float(np.percentile(global_null, 2.5)),
        ci_high=float(np.percentile(global_null, 97.5)),
        null_center=float(np.median(null_divs_global)),
        n_proteins=int(len(obs_divs)),
        n_significant=int((per_protein_valid["perm_q"] < 0.05).sum()),
        per_protein=per_protein_valid,
    )

study_or(det_classes, bg_classes, target='basic')

Per-study OR dictionary from class label arrays.

Source code in src/proteoform_regions/stats.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def study_or(det_classes, bg_classes, target: str = "basic") -> dict:
    """Per-study OR dictionary from class label arrays."""
    det_classes = pd.Series(det_classes).dropna()
    bg_classes = pd.Series(bg_classes).dropna()
    a = int((det_classes == target).sum())
    b = int((det_classes != target).sum())
    c = int((bg_classes == target).sum())
    d = int((bg_classes != target).sum())
    lor, se = or_2x2(a, b, c, d)
    return dict(
        a=a,
        b=b,
        c=c,
        d=d,
        log_or=lor,
        se=se,
        or_=float(np.exp(lor)),
        ci_low=float(np.exp(lor - 1.96 * se)),
        ci_high=float(np.exp(lor + 1.96 * se)),
        p=float(2 * stats.norm.sf(abs(lor / se))),
    )

theoretical_charge_ph26(seq)

Dominant-charge proxy at acidic LC pH (~2.6): count(K)+count(R)+count(H)+1 (N-term).

Upper-bound proxy; ESI charge-state distributions typically peak lower.

Source code in src/proteoform_regions/stats.py 
254
255
256
257
258
259
def theoretical_charge_ph26(seq: str) -> int:
    """Dominant-charge proxy at acidic LC pH (~2.6): count(K)+count(R)+count(H)+1 (N-term).

    Upper-bound proxy; ESI charge-state distributions typically peak lower.
    """
    return sum(1 for aa in str(seq) if aa in "KRH") + 1

Physicochemical engine

proteoform_regions.physchem

pI / MW / net-charge engine with multiple pK models.

Primary calculator: Biopython IsoelectricPoint (Bjellqvist pK set, with residue-specific terminal pKas). Sensitivity calculators: vectorized 60-iteration bisection over Henderson-Hasselbalch charge with the pyteomics Rodwell/EMBOSS and Lehninger pK tables.

Extraction note: the notebook version configured the Henderson-Hasselbalch model through module globals (hh_pi_configure); here the pK table is an explicit argument. Float operations are unchanged - equivalence against the original global-state implementation is pinned by tests/test_physchem.py (reference implementation comparison). Golden values from the publication's frozen ExPASy benchmark (16 sequences) live in tests/data/physchem_golden.csv.

biopython_mw(sequence)

Molecular weight: Biopython average-isotopic mass (X->A replacement before calculation).

Source code in src/proteoform_regions/physchem.py 
133
134
135
def biopython_mw(sequence: str) -> float:
    """Molecular weight: Biopython average-isotopic mass (X->A replacement before calculation)."""
    return ProteinAnalysis(str(sequence).replace("X", "A")).molecular_weight()

biopython_pi(sequence)

Primary pI calculator: Biopython IsoelectricPoint (Bjellqvist; X->A).

Source code in src/proteoform_regions/physchem.py 
128
129
130
def biopython_pi(sequence: str) -> float:
    """Primary pI calculator: Biopython IsoelectricPoint (Bjellqvist; X->A)."""
    return IsoelectricPoint(str(sequence).replace("X", "A")).pi()

charge_at_ph(sequences, ph=7.4, table='hh_rodwell_emboss_lineage')

Net charge per unique sequence at ph (paper default: pH 7.4, Rodwell lineage).

Source code in src/proteoform_regions/physchem.py 
115
116
117
118
119
120
121
122
123
124
125
def charge_at_ph(
    sequences: Sequence[str],
    ph: float = 7.4,
    table: str = "hh_rodwell_emboss_lineage",
) -> dict[str, float]:
    """Net charge per unique sequence at `ph` (paper default: pH 7.4, Rodwell lineage)."""
    uniq = list(dict.fromkeys(str(s) for s in sequences))
    pks = PK_TABLES[table] if isinstance(table, str) else table
    counts = count_matrix(uniq)
    charges = hh_charge(ph, counts, pks)
    return dict(zip(uniq, (float(c) for c in charges), strict=True))

classify_pi(value)

Pre-specified pI classification: acidic < 7.0 <= near_neutral < 7.5 <= basic.

Returns None for None/NaN input (the S3 variant of the paper's helper, which had the most defensive NaN handling of the four notebook copies).

Source code in src/proteoform_regions/schema.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def classify_pi(value) -> str | None:
    """Pre-specified pI classification: acidic < 7.0 <= near_neutral < 7.5 <= basic.

    Returns None for None/NaN input (the S3 variant of the paper's helper, which
    had the most defensive NaN handling of the four notebook copies).
    """
    if value is None:
        return None
    if isinstance(value, float) and math.isnan(value):
        return None
    if value < ACIDIC_BOUND:
        return "acidic"
    if value < NEUTRAL_BOUND:
        return "near_neutral"
    return "basic"

count_matrix(sequences)

Per-ionizable-residue count vectors for a list of sequences (verbatim).

Source code in src/proteoform_regions/physchem.py 
66
67
68
69
def count_matrix(sequences: Sequence[str]) -> dict[str, np.ndarray]:
    """Per-ionizable-residue count vectors for a list of sequences (verbatim)."""
    seqs = [str(s).upper() for s in sequences]
    return {res: np.array([seq.count(res) for seq in seqs], dtype=float) for res in IONIZABLE}

fetch_expasy(sequence, timeout=60.0, retries=2)

Theoretical pI/MW from live ExPASy ProtParam (network; benchmark use only).

Returns (pi, mw) floats or (None, None) on failure. Marked network-reliant: tests that use it are skipped in offline CI.

Source code in src/proteoform_regions/physchem.py 
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def fetch_expasy(sequence: str, timeout: float = 60.0, retries: int = 2):
    """Theoretical pI/MW from live ExPASy ProtParam (network; benchmark use only).

    Returns (pi, mw) floats or (None, None) on failure. Marked network-reliant:
    tests that use it are skipped in offline CI.
    """
    import time

    import requests

    url = "https://web.expasy.org/cgi-bin/protparam/protparam"
    for attempt in range(1, retries + 1):
        try:
            response = requests.get(
                url,
                params={"sequence": str(sequence)},
                timeout=timeout,
                headers={"User-Agent": "proteoform-regions benchmark (validation)"},
            )
            response.raise_for_status()
            import re

            text = re.sub(r"<[^>]+>", " ", response.text)
            p_i = re.search(r"Theoretical pI[/:]?\s*([\d.]+)", text)
            mw = re.search(r"Molecular weight[/:]?\s*([\d.]+)", text)
            return (
                float(p_i.group(1)) if p_i else None,
                float(mw.group(1)) if mw else None,
            )
        except Exception:
            if attempt == retries:
                return None, None
            time.sleep(2.0)

gravy(seq)

Mean Kyte-Doolittle hydropathicity over recognized residues (None if empty).

Source code in src/proteoform_regions/physchem.py 
143
144
145
146
def gravy(seq: str) -> float | None:
    """Mean Kyte-Doolittle hydropathicity over recognized residues (None if empty)."""
    values = [_GRAVY_TABLE[aa] for aa in str(seq) if aa in _GRAVY_TABLE]
    return float(np.mean(values)) if values else None

hh_charge(ph, counts, table='hh_rodwell_emboss_lineage')

Net charge at ph from residue-count vectors under a pK table.

table may be a PK_TABLES key or a {'n_term','c_term','side'} dict (the notebook's global-state configuration becomes an argument).

Source code in src/proteoform_regions/physchem.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def hh_charge(ph: float, counts: Mapping[str, np.ndarray], table: str | dict = "hh_rodwell_emboss_lineage"):
    """Net charge at `ph` from residue-count vectors under a pK table.

    `table` may be a PK_TABLES key or a {'n_term','c_term','side'} dict
    (the notebook's global-state configuration becomes an argument).
    """
    pks = PK_TABLES[table] if isinstance(table, str) else table
    p_k_n, p_k_c, side = pks["n_term"], pks["c_term"], pks["side"]
    charge = 1.0 / (1.0 + 10 ** (ph - p_k_n)) + -1.0 / (1.0 + 10 ** (p_k_c - ph))
    for res, p_k in side.items():
        n = counts[res]
        if p_k is None or not np.any(n):
            continue
        if res in ("D", "E", "C", "Y"):
            charge += -n / (1.0 + 10 ** (p_k - ph))
        else:
            charge += n / (1.0 + 10 ** (ph - p_k))
    return charge

hh_pi_batch(sequences, table='hh_rodwell_emboss_lineage', iterations=60)

Vectorized bisection pI for all sequences with the given pK table.

60 iterations of bisection on [0, 14]; convergence implicit (~14/2**60). Sign convention: pI is where charge crosses from positive (low pH) to non-positive (high pH).

Source code in src/proteoform_regions/physchem.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def hh_pi_batch(
    sequences: Sequence[str],
    table: str = "hh_rodwell_emboss_lineage",
    iterations: int = 60,
) -> np.ndarray:
    """Vectorized bisection pI for all sequences with the given pK table.

    60 iterations of bisection on [0, 14]; convergence implicit (~14/2**60).
    Sign convention: pI is where charge crosses from positive (low pH) to
    non-positive (high pH).
    """
    pks = PK_TABLES[table] if isinstance(table, str) else table
    counts = count_matrix(sequences)
    lo = np.full(len(sequences), 0.0)
    hi = np.full(len(sequences), 14.0)
    for _ in range(iterations):
        mid = (lo + hi) / 2.0
        positive = hh_charge(mid, counts, pks) > 0
        lo = np.where(positive, mid, lo)
        hi = np.where(positive, hi, mid)
    return (lo + hi) / 2.0

peptide_features(seq)

Composition features (verbatim): length, acidic/basic fractions, Cys/Met, tryptic-terminus flag, missed-cleavage count from internal KR|not-P sites.

Source code in src/proteoform_regions/physchem.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def peptide_features(seq: str) -> dict:
    """Composition features (verbatim): length, acidic/basic fractions, Cys/Met,
    tryptic-terminus flag, missed-cleavage count from internal KR|not-P sites."""
    seq = str(seq)
    n = len(seq)
    acidic = sum(seq.count(r) for r in "DE")
    basic = sum(seq.count(r) for r in "KRH")
    return {
        "peptide_length": n,
        "acidic_frac": acidic / n,
        "basic_frac": basic / n,
        "cys_count": seq.count("C"),
        "met_count": seq.count("M"),
        "tryptic_terminus": seq[-1] in "KR",
        "missed_cleavages": sum(1 for i, aa in enumerate(seq[:-1]) if aa in "KR" and seq[i + 1] != "P"),
    }

pi_for_calculator(sequences, calculator='bjellqvist')

pI per sequence under a named calculator: 'bjellqvist' (Biopython primary) or any PK_TABLES key ('hh_rodwell_emboss_lineage', 'hh_lehninger').

Source code in src/proteoform_regions/physchem.py 
205
206
207
208
209
210
211
212
def pi_for_calculator(sequences: Iterable[str], calculator: str = "bjellqvist") -> dict[str, float]:
    """pI per sequence under a named calculator: 'bjellqvist' (Biopython primary)
    or any PK_TABLES key ('hh_rodwell_emboss_lineage', 'hh_lehninger')."""
    uniq = list(dict.fromkeys(str(s) for s in sequences))
    if calculator == "bjellqvist":
        return {s: biopython_pi(s) for s in uniq}
    values = hh_pi_batch(uniq, table=calculator)
    return dict(zip(uniq, (float(v) for v in values), strict=True))

Regions

proteoform_regions.regions

Single-linkage detected-region stitching from peptide evidence.

Algorithm (Supplementary Methods Module 3 of the source publication):

  • Sort peptide spans by start position (tuple sort: start, end, then payload).
  • Walk spans in order; join a span to the current cluster when next.start - last_member.end <= gap (paper default gap = 25).
  • Region bounds are min(member starts) .. max(member ends); the region sequence is the contiguous precursor slice between the bounds (unobserved intermediate residues included).
  • A region with fewer than 2 distinct member peptide sequences stays peptide_only; >= 2 is region_supported; a numerically consistent MW constraint window upgrades to fragment_constrained (ANL-020).

Exact-semantics note (pinned by tests). The gap check uses the END OF THE LAST-APPENDED MEMBER, not the running maximum end of the cluster. Because spans are sorted by start, a later-starting nested span (smaller end) becomes the reference end for the next comparison, which can SPLIT a span that textbook running-max single linkage would join. This is the published algorithm's exact behavior ("D1 last-member-end"); it is intentionally preserved, not normalized. The join condition next.start - last.end <= 25 corresponds to at most 24 unobserved residues between the last member and the next span.

confidence_label(n_unique_sequences, region_mw=None, constraint_window=None, min_peptides=MIN_REGION_PEPTIDES)

Tiered region confidence with ANL-020 MW-window gating.

Returns (label, mw_constraint_applied):

  • peptide_only : fewer than min_peptides distinct sequences, OR an available numeric window that the region MW violates (an out-of-window region cannot exceed the evidence of its peptides).
  • region_supported : >= min_peptides sequences and no numeric window.
  • fragment_constrained : >= min_peptides sequences and MW inside the window.

The paper's pilots exposed no numeric windows, so every pilot region is peptide_only/region_supported with mw_constraint_applied=False - the unreachable path is exercised only by tests.

Source code in src/proteoform_regions/regions.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def confidence_label(
    n_unique_sequences: int,
    region_mw: float | None = None,
    constraint_window: tuple[float, float] | None = None,
    min_peptides: int = MIN_REGION_PEPTIDES,
) -> tuple[str, bool]:
    """Tiered region confidence with ANL-020 MW-window gating.

    Returns (label, mw_constraint_applied):

    - ``peptide_only``          : fewer than `min_peptides` distinct sequences, OR an
                                  available numeric window that the region MW violates
                                  (an out-of-window region cannot exceed the evidence
                                  of its peptides).
    - ``region_supported``      : >= `min_peptides` sequences and no numeric window.
    - ``fragment_constrained``  : >= `min_peptides` sequences and MW inside the window.

    The paper's pilots exposed no numeric windows, so every pilot region is
    peptide_only/region_supported with ``mw_constraint_applied=False`` - the
    unreachable path is exercised only by tests.
    """
    if n_unique_sequences < min_peptides:
        return "peptide_only", False
    if constraint_window is not None and region_mw is not None:
        low, high = constraint_window
        if not (low <= region_mw <= high):
            return "peptide_only", True
        return "fragment_constrained", True
    return "region_supported", False

region_bounds(cluster)

(min member start, max member end) -> 1-based inclusive region bounds.

Source code in src/proteoform_regions/regions.py 
73
74
75
def region_bounds(cluster: Sequence[tuple]) -> tuple[int, int]:
    """(min member start, max member end) -> 1-based inclusive region bounds."""
    return int(min(s[0] for s in cluster)), int(max(s[1] for s in cluster))

region_physchem(sequence)

(pI, MW) of a region sequence (X->A; None inputs -> (None, None)).

Source code in src/proteoform_regions/regions.py 
129
130
131
132
133
134
def region_physchem(sequence: str | None) -> tuple[float | None, float | None]:
    """(pI, MW) of a region sequence (X->A; None inputs -> (None, None))."""
    if not sequence:
        return None, None
    cleaned = sequence.replace("X", "A")
    return biopython_pi(cleaned), biopython_mw(cleaned)

region_record(study, accession, cluster, precursor_length, precursor_pi=None, precursor_mw=None, gene_name=None, precursor_sequence='', gap=REGION_GAP_AA, terminal_fraction=TERMINAL_FRACTION, min_peptides=MIN_REGION_PEPTIDES, constraint_window=None, extra_flags=None)

Build one canonical region row (S4 column set) from a stitched cluster.

cluster spans are (start, end, sequence, peptide_pi, row_id) tuples, as produced by the paper's S4 walk. extra_flags carries RQ3 UniProt-feature flags when a feature index is supplied upstream.

Source code in src/proteoform_regions/regions.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def region_record(
    study: str,
    accession: str,
    cluster: Sequence[tuple],
    precursor_length: int,
    precursor_pi: float | None = None,
    precursor_mw: float | None = None,
    gene_name: str | None = None,
    precursor_sequence: str = "",
    gap: int = REGION_GAP_AA,
    terminal_fraction: float = TERMINAL_FRACTION,
    min_peptides: int = MIN_REGION_PEPTIDES,
    constraint_window: tuple[float, float] | None = None,
    extra_flags: str | None = None,
) -> dict:
    """Build one canonical region row (S4 column set) from a stitched cluster.

    `cluster` spans are (start, end, sequence, peptide_pi, row_id) tuples, as
    produced by the paper's S4 walk. `extra_flags` carries RQ3 UniProt-feature
    flags when a feature index is supplied upstream.
    """
    region_start, region_end = region_bounds(cluster)
    n_peptides = len({s[2] for s in cluster})
    span_length = region_end - region_start + 1
    seq = region_sequence(precursor_sequence, region_start, region_end)
    region_pi, region_mw = region_physchem(seq)
    n_term, c_term = terminal_flags(region_start, region_end, precursor_length, terminal_fraction)
    label, constraint_applied = confidence_label(n_peptides, region_mw, constraint_window, min_peptides)
    return {
        "dataset_accession": study,
        "mapped_accession": accession,
        "gene_name": gene_name,
        "region_start": region_start,
        "region_end": region_end,
        "region_span_aa": span_length,
        "region_n_peptides": n_peptides,
        "region_unique_sequences": n_peptides,
        "member_peptide_rows": ";".join(str(int(s[4])) for s in cluster),
        "member_sequences": ";".join(sorted({s[2] for s in cluster})),
        "precursor_length": int(precursor_length),
        "coverage_fraction": span_length / precursor_length,
        "n_terminal_inclusive": n_term,
        "c_terminal_inclusive": c_term,
        "region_pi": region_pi,
        "region_mw_da": region_mw,
        "precursor_pi": precursor_pi,
        "precursor_mw_da": precursor_mw,
        "region_pi_divergence": (region_pi - precursor_pi) if region_pi is not None else None,
        "region_pi_class_anl003": classify_pi(region_pi),
        "precursor_pi_class_anl003": classify_pi(precursor_pi),
        "class_flip_region_vs_precursor": (
            classify_pi(region_pi) != classify_pi(precursor_pi)
            if region_pi is not None and precursor_pi is not None
            else None
        ),
        "confidence_label": label,
        "mw_constraint_applied": constraint_applied,
        "rq3_class_flags": extra_flags,
        "entity_level": "detected_region",
    }

region_sequence(precursor_sequence, region_start, region_end)

Contiguous precursor slice covering the region (1-based inclusive bounds).

Source code in src/proteoform_regions/regions.py 
122
123
124
125
126
def region_sequence(precursor_sequence: str, region_start: int, region_end: int) -> str | None:
    """Contiguous precursor slice covering the region (1-based inclusive bounds)."""
    if not precursor_sequence:
        return None
    return precursor_sequence[region_start - 1 : region_end]

stitch_spans(spans, gap=REGION_GAP_AA)

Cluster sorted-by-start spans into regions via the paper's exact loop.

Each span is a tuple whose first two elements are (start, end); remaining elements are carried through untouched (the paper carries (sequence, peptide_pi, row_id)). Tuples are sorted exactly as the notebook does (lexicographic), so input order does not matter - determinism is total.

Source code in src/proteoform_regions/regions.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def stitch_spans(spans: Sequence[tuple], gap: int = REGION_GAP_AA) -> list[list[tuple]]:
    """Cluster sorted-by-start spans into regions via the paper's exact loop.

    Each span is a tuple whose first two elements are (start, end); remaining
    elements are carried through untouched (the paper carries
    (sequence, peptide_pi, row_id)). Tuples are sorted exactly as the notebook
    does (lexicographic), so input order does not matter - determinism is total.
    """
    spans = sorted(spans)
    clusters: list[list[tuple]] = [[spans[0]]]
    current = clusters[0]
    for span in spans[1:]:
        if span[0] - current[-1][1] <= gap:
            current.append(span)
        else:
            current = [span]
            clusters.append(current)
    return clusters

terminal_flags(region_start, region_end, precursor_length, terminal_fraction=TERMINAL_FRACTION)

(n_terminal_inclusive, c_terminal_inclusive) per the 10% rule.

Source code in src/proteoform_regions/regions.py 
78
79
80
81
82
83
84
85
86
87
88
def terminal_flags(
    region_start: int,
    region_end: int,
    precursor_length: int,
    terminal_fraction: float = TERMINAL_FRACTION,
) -> tuple[bool, bool]:
    """(n_terminal_inclusive, c_terminal_inclusive) per the 10% rule."""
    return (
        region_start <= terminal_fraction * precursor_length,
        region_end >= (1 - terminal_fraction) * precursor_length,
    )