Feature: tl.peptide_proximity() - #38
Merged
Merged
Conversation
Sequence proximity analysis for COPF proteoform clusters -- a reimplementation of CCprofiler's `evaluateProteoformLocation` (Bludau et al. 2021, branch `proteoformLocationMapping` @ 31a3043), which is the downstream characterisation COPF applies to the proteoform groups it detects: are a cluster's peptides closer together in the protein sequence than a random grouping of the same size? Peptide positions are resolved from a FASTA inside the call, reusing `pp.summarize_peptides_by_neighbourhood_union`'s locator so the two functions cannot drift apart on how a peptide is placed. Named `peptide_proximity`, not `proximity_analysis` and not the reference's `evaluateProteoformLocation`, for consistency with `pr.tl.differential_abundance` -- no `_analysis` suffix. Added to the Sphinx `tl.rst` COPF listing alongside the other COPF steps. ## What lands in `.var` Two families of three. Per cluster, repeated across every peptide of the `(protein_id, cluster_id)` group: - `peptide_proximity_pval` -- classical empirical p-value - `peptide_proximity_pseudo_pval` -- counts only strictly smaller permutations; the paper's "lowest possible p-value" criterion - `peptide_proximity_pval_adj` -- Benjamini-Hochberg over the tests Per protein, repeated across all its peptides, as `peptide_proximity_unified_*`: the minimum over the protein's clusters, which is the quantity `getProteoformStats` thresholds. The minimum skips clusters with no value, and a protein qualifies only with at least two distinct non-noise clusters -- below two there is no proteoform split to localise. Skipping rather than propagating is load-bearing: an untested single-peptide *noise* cluster is a placeholder the reference scores 1, and propagating it would suppress a real finding in a sibling cluster. `unified_pval_adj` is BH over `unified_pval` itself, one entry per protein, not the minimum of the per-cluster adjusted values. The per-test table and the parameters land in `.uns['peptide_proximity']`. ## Reference behaviours reproduced on purpose Documented in the module docstring: the noise cluster is tested like any other cluster; the normalisation mixes cluster size with protein-wide ranks; single-peptide clusters have no standard deviation; p-values are add-one empirical; unresolved peptide positions rank first, as `data.table`'s `setkey` puts them; and ineligible proteins are left `NaN`, so not every row carries a value (1,272 of 2,885 proteins are tested on the reference dataset). Permutations are re-seeded per peptide count, reproducing R's per-protein `set.seed(123)`; `random_state=None` means that seed, so the default run is the one comparable with the publication. Ties between the observed score and a permuted one decide the gap between the two p-values, and comparing float scores misjudges them, so `tie_arithmetic='exact'` compares an exact integer dispersion with the same ordering instead. A cluster of one peptide has no standard deviation, so it is left `NaN` rather than given the placeholder p-value of 1 the reference records. Nothing downstream changes: 1 could never be a protein's minimum, and the unified statistics skip a missing cluster. On the mouse tissue dataset the reference emits 283 such rows, every one of them a noise cluster, and the published counts are identical either way. `min_clusters_per_protein=2` is a guard the reference lacks, on by default and inert on the published analysis: it skips a protein whose single cluster holds every peptide, where the score is exactly 1 and the pseudo p-value hits its floor whatever the data say. Pass `1` for the reference's behaviour. ## Tests 39 tests pin the whole chain against CCprofiler's own output for the mouse tissue dataset: FASTA -> positions -> ranks -> per-cluster permutation p-values -> the counts the paper reports. - all 24,534 peptide positions equal to the reference's, including the one peptide it cannot locate and keeps - the same 1,272 tested proteins exactly, and every group tested here is one CCprofiler tested - the published split reproduced twice over, with identical protein sets: by aggregating the per-cluster columns, and read straight off the `unified_` columns - p-values within the sampling noise of two independent 1,000-draw estimates, which is all that is achievable without R's `sample()` stream Three reference assets added under `tests/data/mouse_tissue/`, with provenance and checksums in a new README there: the reference proximity p-values, the reference per-peptide proteoform assignment, and a subset FASTA whose derived positions are checked against the reference's for all 24,534 peptides rather than assumed equal. `tests/data/` is now excluded from the `trailing-whitespace`, `end-of-file-fixer` and `mixed-line-ending` hooks, as it already was from `black` and `pyupgrade`. It holds reference output where trailing whitespace is data: an untested proteoform group in the location table is a row with two empty trailing fields, and stripping them turns a 4-field row into a 2-field one. ## Two findings worth knowing `peptide_proximity_pval_adj` is 1.0 for every test on this dataset, and that is a property of the data rather than of the arithmetic: the BH curve `(n/j) * p_(j)` sits strictly above 1 for every `j < n` and touches 1 at `j = n`, i.e. `p_(j) > j/n` throughout, so the empirical distribution lies entirely below the uniform diagonal. Proximity is a sparse effect -- the publication flags 26 proteins among some three thousand clusters -- and no FDR level can declare an excess that thin. Raising `n_permutations` does not help. The docstring records why `pseudo_pval` is not a way round it either. The paper's Results text and Fig. 6C attach its two proximity numbers to the opposite criteria from the ones its own code computes: the classical criterion selects 7 proteins and the pseudo-only criterion 19, not the reverse. The direction is forced, since the pseudo p-value never exceeds the classical one. A run reporting 7 and 19 is correct; see `test_paper_sequence_proximity_result_vs_rcopf`. Full suite 488 passed. flake8 and pylint (errors only) clean. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 68ecf52adc
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Feature:
proteopy.tl.peptide_proximitytl;dr
Implements the proximity analysis from Bludau et al. 2021's COPF workflow: given a group of peptides that COPF has assigned to the same proteoform cluster, are those peptides closer together in the protein sequence than a random grouping of the same size would be? This is the second stage of the paper's proteoform-localization pipeline (COPF finds the groups;
peptide_proximityasks whether each group is spatially coherent), and ProteoPy previously had no equivalent — the mouse tissue benchmark could only be reproduced up to the paper's 63 significant proteins, not past them.The new function is validated against CCprofiler's own reference implementation rather than against ProteoPy's prior behavior, so the numbers in this PR are checked against the paper's actual pipeline, not against an assumption about what ProteoPy used to do.
Background: algorithm and assumptions
The analysis operates per
(protein_id, cluster_id)group produced by COPF. For each group, it computes a proximity score from the peptides' sequence positions and tests that score against a null distribution built by permuting peptide-to-cluster assignments within the protein. This produces an empirical p-value per group, which is then aggregated to a per-protein statistic — the quantity the paper's significance calls are actually based on.A few properties of the reference algorithm are preserved deliberately, because reproducing the paper means reproducing its exact conventions, not just its intent:
(hits + 1) / (n + 1)), not raw permutation fractions.data.table'ssetkeyin the reference R implementation.NaNrather than assigned a placeholder value.Two points where this implementation deliberately departs from the reference, both documented in the module docstring so they don't read as bugs later:
NaNhere instead of the reference's placeholder p-value of 1. This is inert for downstream results: a p-value of 1 could never be selected as a protein's minimum, and the per-protein aggregate simply skips a missing cluster.set.seed(123). Because R'ssample()stream isn't reproducible from Python, individual p-values differ by the sampling noise of two independent 1,000-draw estimates — but every discrete significance call matches exactly.Implementation
pr.tl.peptide_proximity()resolves peptide sequence positions from a FASTA file at call time, reusing the same locator thatpp.summarize_peptides_by_neighbourhood_unionalready uses. Sharing that locator (rather than reimplementing position lookup) is intentional: it's the only way to guarantee the two functions can't silently drift apart on how a peptide's position is defined.Output lands in six new
.varcolumns, two families of three:Per
(protein_id, cluster_id)group (repeated across the group's peptides):The FASTA subset itself is verified rather than trusted:
test_peptide_positions_from_fasta_vs_rcopfre-derives all 24,534 peptide positions from it and checks them against the reference's positions directly.tests/data/is now excluded from thetrailing-whitespace,end-of-file-fixer, andmixed-line-endingpre-commit hooks (it was already excluded fromblackandpyupgrade). This directory holds reference output, where trailing whitespace can be meaningful data — an untested proteoform group in the reference location table is represented as a row with two empty trailing fields, and stripping that whitespace would silently turn a 4-field row into a 2-field one.Full test results:
Two findings worth flagging for review
peptide_proximity_pval_adjis 1.0 for every test on this dataset, and this is a property of the data, not a bug in the arithmetic. BH's smallest adjusted value ismin_j (n/j)·p_(j); here the p-value curve stays strictly above thej/ndiagonal for everyj < nand only touches 1 atj = n, so the empirical distribution never crosses under the uniform diagonal. Proximity is a genuinely sparse effect in this dataset — the paper flags 26 proteins out of roughly three thousand clusters — and no FDR threshold can declare significance from an excess that thin. Increasingn_permutationsdoesn't change this. The module docstring explains whypseudo_pvalisn't a workaround either, and notes thatunified_pvalis a minimum over dependent tests, so it's a ranking statistic rather than an FDR-controlled quantity.The paper's Results text and Fig. 6C swap which number goes with which criterion. Running the reference workflow gives 7 significant proteins under the classical criterion and 19 under the pseudo-only criterion — not 19 and 7, as the text and figure state. This direction is forced: the pseudo p-value can never exceed the classical one, so 19-classical/7-pseudo is not an achievable outcome.
test_paper_sequence_proximity_result_vs_rcopfchecks for 7-classical/19-pseudo, and the module docstring flags this discrepancy explicitly so a future reader doesn't "fix" it into matching the paper's swapped numbers.