Conversation
Add the corrected Pairwise Homoplasy Score (cPHS) as a TreeData-first functional tool. Rewrites the starter implementation to match package conventions: standard argument names, resolution via shared utils helpers, reuse of fraction_mutated for the mutation-rate estimate, and a clean per-tree core. - Returns a dict of scores keyed by tree name when multiple trees are present and tree_key is None; a float otherwise. - Raises TreeMetricError with a pointer to ancestral_characters (and a note on reconstruction sensitivity) when internal-node character states are missing. - Docstring warns that cPHS is strongly influenced by ancestral reconstruction accuracy and should not be used to compare reconstruction algorithms. - Adds pytest coverage and docs/API entry. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## 3.0.0 #289 +/- ##
==========================================
+ Coverage 85.61% 85.72% +0.10%
==========================================
Files 57 57
Lines 5937 6017 +80
==========================================
+ Hits 5083 5158 +75
- Misses 854 859 +5
🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull request overview
Adds a new tree metric, cassiopeia.tl.calculate_cPHS, implementing the corrected Pairwise Homoplasy Score (cPHS) as a TreeData-first tool with package-consistent parameter resolution, multi-tree handling, and accompanying documentation/tests.
Changes:
- Implement
calculate_cPHSintools/tree_metrics.py, including collision-probability estimation from priors and ultrametric/ancestral-state validation. - Export the new tool from
cassiopeia.toolsand document it in the API docs. - Add pytest coverage for scalar vs dict return, error cases, and collision-probability estimation behavior.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
| tests/tools_tests/tree_metrics_test.py | Adds tests covering cPHS behavior (returns, errors, priors/default q). |
| src/cassiopeia/tools/tree_metrics.py | Implements _collision_probability, _calculate_cphs, and calculate_cPHS. |
| src/cassiopeia/tools/init.py | Exposes calculate_cPHS via imports and __all__. |
| docs/api/tools.rst | Adds tl.calculate_cPHS to the documented Metrics API list. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| alpha = np.exp(-mutation_rate * lca_heights) | ||
| beta = 1 - np.exp(-mutation_rate * (1 - lca_heights)) | ||
| prob = alpha * beta**2 * collision_probability | ||
| prob[np.isclose(lca_heights, 1)] = 1 | ||
| pvalues = 1 - scipy.stats.binom.cdf(phs - 1, k, prob) |
There was a problem hiding this comment.
Short version: at lca_height == 1 the pair carries no information, and
prob = 1 is what removes it from the score. prob = 0 would do the
opposite.
The reason is what prob feeds into on the next line:
pvalues = 1 - scipy.stats.binom.cdf(phs - 1, k, prob)
prob is the null probability of seeing a homoplasy in one character.
Setting it to 0 tells the binomial that a homoplasy there is impossible.
So if the pair has even one homoplasy, the p-value is exactly 0, which
the following line converts to eps (~1e-16). The cPHS is the minimum
adjusted p-value over all pairs, so that single degenerate pair would
become the score for the whole tree.
prob = 1 goes the other way: the p-value is 1, the largest possible, so
the pair can never be the minimum and drops out of the calculation. That
is the intent. At lca_height == 1 the LCA sits at leaf level, meaning
zero elapsed time between it and the leaves, so there was no opportunity
for independent mutations to arise. The pair tells us nothing and should
not influence the result.
In practice it never triggers. IIDExponentialMLE enforces a minimum
branch length, so internal nodes cannot reach leaf depth. Across 74
KPTracer trees the deepest internal node sits at 0.99 on a unit-scaled
tree, a gap of 1e-2, six orders of magnitude larger than the np.isclose
tolerance. The guard only matters if branch lengths come from elsewhere
and permit zero-length terminal edges.
| # Benjamini-Hochberg style adjustment; the cPHS is the minimum adjusted p-value. | ||
| pvalues_sorted = np.sort(pvalues) | ||
| adjusted_pvalues = pvalues_sorted * len(pvalues) / np.arange(1, len(pvalues) + 1) | ||
| return float(np.min(adjusted_pvalues)) |
There was a problem hiding this comment.
The minimum cannot exceed 1, so the cap is not needed.
The adjusted values are p_sorted[i] * N / i. The last entry (i = N) is
multiplied by N/N = 1, so it is just p_max unchanged, which is at most 1
since it is a p-value. There is therefore always at least one entry <= 1,
and the function returns the minimum, so the result is always <= 1.
Individual entries can exceed 1: p = [0.9, 0.95, 1.0] gives adjusted
values [2.7, 1.425, 1.0]. But the 2.7 is never selected, since 1.0 is
smaller. Checked numerically over 200k random configurations, the maximum
returned value was 0.9999977, and the all-p-equal-1 case returns exactly
1.0.
Adding np.minimum(..., 1.0) would be harmless but never has an effect, so
I would leave it out to avoid implying a failure mode that cannot occur.
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
_collision_probability estimates q = sum_s p_s^2, the chance that two independent mutations at a character produce the same state. This value scales the null probability of a homoplasy in cPHS, so any error in q propagates to every pair's p-value. The previouse implementation squares the prior values as stored, which is correct only when each character's priors already sum to 1. That does not hold for real biological datasets hold like the KPTracer priors files (Yang et al. 2022, Zenodo 5847462), which store unnormalized allele-frequency weights. In the real data each character's weights sum to roughly 3.4 rather than 1. Measured on tumor 3435_NT_T3 (30 characters): correct q (renormalized per character, then averaged): 0.0966 q returned by the previous code: 1.2778 By definition q<=1, so the old calculated that satisfies q>1 is outside the valid range for a probability. The null probability is computed as alpha * beta**2 * q, so an inflated q inflates the per-pair binomial tail probability and therefore every p-value. As cPHS is the minimum BH-adjusted p-value across pairs, the score is driven toward 1. The failure mode is systematic and directional: reconstructions that should be flagged as containing impossible homoplasies would instead pass. On a representative pair (k=30, phs=3, lambda=2.0, tau=0.3) the inflated q raised the p-value by roughly 60x, from 1.6e-2 to 9.6e-1. The bug does not surface on simulated data, where priors are generated already normalized and identical across characters, which is why the existing tests do not catch it. It affects only real datasets with unnormalized priors, i.e. precisely the data the metric is intended for. This commit adds a _normalized_collision helper that renormalizes a single character's priors to sum to 1 before computing sum_s p_s^2, and routes both branches of _collision_probability (the per-character dict-of-dicts case and the flat state-to-probability case) through it. The helper raises TreeMetricError if a character's weights sum to zero or less, which would otherwise produce a silent division by zero. For priors that already sum to 1 the renormalization divides by 1 and is, so results on simulated data and on any correctly normalized input are unchanged. Only previously incorrect results change. Verified on both unnormalized (weights summing to ~3.4) and normalized inputs: both now yield the same q, matching the reference implementation used to produce the published KPTracer cPHS values.
for more information, see https://pre-commit.ci
The existing tests for _collision_probability only use priors that already sum to 1, which is why the missing normalization was not caught before. This adds a test covering the unnormalized case. The test checks that the same distribution yields the same q whether it is given as probabilities or as raw weights, that q stays within [0, 1], that per-character distributions are averaged correctly when they differ from one another, and that priors summing to zero raise TreeMetricError rather than dividing by zero.
for more information, see https://pre-commit.ci
Summary
Adds
cassiopeia.tl.calculate_cPHS, the corrected Pairwise Homoplasy Score (cPHS; Zilber et al., 2026), as a TreeData-first functional tool. The starter implementation intools/tree_metrics.pyis rewritten to match the package's functional conventions.Details
tree_key,characters_key,time_key), inputs resolved via the sharedutilshelpers (_get_parameter,_get_characters,_get_digraph,_get_root,get_leaves), and mutation-rate estimation reusesparameter_estimators.fraction_mutated.dict[str, float]keyed by tree name whentree_key is Noneandtdataholds more than one tree; afloatotherwise (an explicittree_keyalways yields a scalar).TreeMetricErrorpointing tocassiopeia.tl.ancestral_characterswhen internal-node character states are missing, and when the tree is not ultrametric undertime_key.q) is estimated from priors (per-character dict-of-dicts or a flat state→prob dict), falling back to a uniform1/mwith a warning when no priors are available.tl.calculate_cPHStotools/__init__(__all__) anddocs/api/tools.rst.Testing
tests/tools_tests/tree_metrics_test.py(scalar/dict return, explicit-key scalar, ancestral-states and ultrametric errors, priors-based and default collision probability).61 passed.🤖 Generated with Claude Code