`.peptide_proximity`

proteopy.tl.peptide_proximity(adata, annotator, *, protein_col='protein_id', peptide_col='peptide_id', cluster_col='cluster_id', noise=1000000.0, n_permutations=1000, min_peptides_per_protein=4, min_peptides_per_proteoform=2, min_clusters_per_protein=2, mod_regex='\\\\(UniMod:[0-9]+\\\\)', alphabet='ACDEFGHIKLMNPQRSTVWYBJOUXZ', on_unknown_protein='raise', on_unlocated_peptide='raise', tie_arithmetic='exact', random_state=None, key_added='peptide_proximity', inplace=True, verbose=False)[source]

Test whether a proteoform’s peptides cluster in the sequence.

Reimplements CCprofiler’s evaluateProteoformLocation [1], the downstream characterisation COPF applies to the proteoform groups it detects. Peptide positions are resolved from annotator, peptides are ranked by position within their protein, and each proteoform cluster’s rank dispersion is compared against n_permutations random groupings of the same size. See the module docstring for the formula and for the six reference behaviours reproduced on purpose.

Six columns are written to .var, in two families of three. The first is a property of the (protein_id, cluster_id) group and is therefore repeated across every peptide of that group:

`{key_added}_pval`	classical empirical p-value
`{key_added}_pval_adj`	the same, Benjamini-Hochberg adjusted over the tests performed
`{key_added}_pseudo_pval`	pseudo p-value, counting only strictly smaller permutations

The second is a property of the protein: the minimum over its clusters, repeated across all its peptides, and named {key_added}_unified_pval, {key_added}_unified_pval_adj and {key_added}_unified_pseudo_pval. This is the quantity CCprofiler’s getProteoformStats thresholds, so it reproduces the published split directly.

The minimum skips clusters with no value, and what decides whether a protein gets one at all is its structure: it needs at least two distinct non-noise clusters, and is NaN otherwise. Below two there is no proteoform split to localise, so a minimum would say nothing; at two or more, every non-noise cluster carries a value anyway, because COPF never leaves a real cluster with fewer than two peptides. Skipping rather than propagating is what keeps an untested noise cluster – one peptide, no standard deviation, a placeholder the reference scores 1 – from suppressing a real finding in a sibling cluster.

The two families are corrected separately, because they answer different questions. pval_adj corrects over the tests, one entry per cluster. unified_pval_adj is Benjamini-Hochberg applied to unified_pval itself, one entry per protein – not the minimum of the per-cluster adjusted values, which would drag the cluster-level n into a protein-level statement. In both families n counts only the entries that carry a p-value, so a family is never inflated by members that were not tested.

⚠️ unified_pval is a minimum over several dependent tests, so it is not a calibrated p-value: its null is stochastically smaller than uniform, and on the reference dataset its median is 0.30 rather than 0.50. Benjamini-Hochberg over it is therefore anti-conservative. Read unified_pval_adj as a ranking statistic, and if a calibrated per-protein p-value is needed, correct within the protein first – Šidák, 1 - (1 - p) ** k over its k clusters – rather than treating the minimum as though it were one test.

P10852 is the case that fixes the rule. Its evidence is a two-peptide cluster whose pseudo p-value reaches the floor, while its noise cluster holds a single peptide no test can score. It is one of the 19 the publication counts, and propagating the missing cluster instead of skipping it would have dropped it.

🔴 On the reference dataset every adjusted p-value comes out at 1.0, and that is a property of these data rather than of the arithmetic. Benjamini-Hochberg’s smallest adjusted value is min_j (n / j) * p_(j) over the ascending p-values – the j = 1 term is the largest member of that family, not a bound on it, so a floor on the raw p-value saturates nothing by itself. Here the curve sits strictly above 1 for every j < n and touches 1 at j = n; equivalently p_(j) > j / n throughout, so the empirical p-value distribution lies entirely below the uniform diagonal and no FDR level could declare an excess of small p-values.

Two reasons, and only the second is decisive. The classical p-value is conservative by construction, because it counts permutations scoring exactly the observed value as evidence against, and rank dispersion is discrete enough that ties are abundant – a mean of 89 tied permutations per 1,000, and 160 for the two-peptide clusters that are 41 % of the family. But the calibrated mid-p variant removes that bias exactly, its median landing on 0.500, and BH still returns nothing below 0.95. So the real reason is that proximity is a sparse effect: the publication flags 26 proteins, which is 26 of some three thousand clusters – far too few to move a procedure controlling the false discovery rate over the whole family. Raising n_permutations does not help, and neither does restricting the family: over just the 63 significant proteins’ clusters the smallest adjusted value is 0.28.

CCprofiler computes no adjusted p-value at all – evaluateProteoformLocation accepts an adj.method argument and never uses it – and the published counts are raw, so compare against the publication on {key_added}_pval.

⚠️ Do not reach for {key_added}_pseudo_pval as a way round this. It is the opposite extreme – it discards the ties rather than counting them – and is anti-conservative to the point of being unusable under FDR control: 383 tests reach its floor, 284 of them two-peptide clusters whose classical p-value has a median of 0.23. Applying BH to it would “reject” 390 tests on the strength of that discreteness alone. The publication uses it as a secondary “lowest possible p-value” criterion, never as a p-value to correct.

Peptides of a protein that was not eligible, or that carry no cluster, get NaN – as does every peptide when on_unknown_protein or on_unlocated_peptide excluded it. Do not expect a value on every row: with the reference defaults a protein needs at least four peptides and a median cluster size of two, so on the mouse tissue dataset 1,272 of 2,885 proteins are tested.

.X is neither read nor written, so layer, zero_to_na and fill_na have no meaning here and are absent.

Parameters:

adata (AnnData) – Peptide-level data. Only .var is read.
annotator (str | pathlib.Path | dict) – Path to a FASTA file, or a pre-parsed {accession: sequence} mapping, supplying the protein sequences that peptide positions are resolved against. This is the load-bearing sequence input: the ranks, and hence every p-value, come from it.
protein_col (str, optional) – Columns in .var holding the protein and peptide identifiers.
peptide_col (str, optional) – Columns in .var holding the protein and peptide identifiers.
cluster_col (str, optional) – Column in .var holding the peptide cluster, as written by peptide_clusters_from_dendograms().
noise (float, optional) – Cluster label marking peptides that joined no multi-peptide cluster. The noise cluster is tested like any other – it is one of the protein’s proteoform groups in the reference – but is flagged in the per-test table.
n_permutations (int, optional) – Random groupings drawn per protein. The empirical p-value cannot go below 1 / (n_permutations + 1).
min_peptides_per_protein (int, optional) – Proteins with fewer peptides are not tested. The reference’s 4, which is deliberately stricter than the 2 the main COPF pipeline filters on.
min_peptides_per_proteoform (int, optional) – Proteins whose median cluster size is below this are not tested.
min_clusters_per_protein (int, optional) – Proteins with fewer clusters are not tested. A single cluster holding every peptide of a protein has rank vector 1..n, and so does every permutation of it, so its score is exactly 1 and its pseudo p-value hits the floor whatever the data say. Pass 1 for the reference’s behaviour; on the published analysis both select the identical 1,272 proteins.
mod_regex (str, optional) – Pattern for annotations to disregard when locating a peptide in its protein sequence. Defaults to CCprofiler’s \(UniMod:[0-9]+\).
alphabet (Iterable[str], optional) – Characters allowed to remain after mod_regex is applied, which is what makes mod_regex self-checking.
on_unknown_protein ({'raise', 'keep', 'skip'}, optional) – What to do about a protein absent from annotator. 'keep' gives its peptides unresolved positions, which ranks them first and is the CCprofiler behaviour; 'skip' excludes them from the analysis. 'raise' is the default because silently proceeding is the reference’s blind spot, and reproducing it should be an explicit choice.
on_unlocated_peptide ({'raise', 'keep', 'skip'}, optional) – The same, for a peptide whose protein is present but whose sequence is not found in it.
tie_arithmetic ({'exact', 'reference'}, optional) – How a permutation scoring equal to the observation is detected – the comparison that separates the two p-values. 'exact' compares an exact integer dispersion with the same ordering as the score; 'reference' compares float scores in R’s arithmetic order, which misjudges ties.
random_state (int | None, optional) – Seed for the permutations. None uses CCprofiler’s hard-coded 123, so the default run is the one comparable with the publication.
key_added (str, optional) – Prefix of the three .var columns, and the .uns key holding the per-test table.
inplace (bool, optional) – If False, return a modified copy instead of writing to adata.
verbose (bool, optional) – Print what was read, how many tests ran, and where results were stored.

Returns:

None when inplace is True, otherwise the modified copy.

Return type:

AnnData | None

Raises:

KeyError – If protein_col, peptide_col or cluster_col is absent from .var.
ValueError – On an out-of-range or unknown argument, or – under the default policies – if a protein or peptide sequence cannot be resolved from annotator.

.peptide_proximity

`.peptide_proximity`