T-cell atlas: nomination and reuse of reliability weights¶

Start with a response gene and obtain candidate perturbations for follow-up. This offline tutorial uses authentic archived sufficient statistics from the primary CD4+ T-cell atlas. The reliability statistics retain information from the full original gene axis, while the effect matrix stores the queried gene coordinates. All candidates in this archived task are retained.

Install the accompanying wheel or run python -m pip install . from the package directory. No atlas download is required for this tutorial. Its computations recompute weights and nominations from those statistics. The fixed scientific parameters are five control directions and attenuation 0.9, as in manuscript v25.2. See scientific design for how the reliability calculation and gene readout are separated.

In [1]:
import numpy as np
import pandas as pd
from pathlib import Path
import pertrank as pr
from pertrank.datasets import load_example

data, provenance = load_example("tcell")
print(provenance["dataset"])
print("Candidates:", len(data["target_ids"]))
print("Weight genes:", int(data["weight_gene_count"]))
print("Stored effect columns:", len(data["gene_ids"]))
queries = data["selected_query_ids"].tolist()
print("Example queries:", queries)
Primary CD4+ T-cell atlas, donor D3, resting state
Candidates: 52
Weight genes: 13390
Stored effect columns: 293
Example queries: ['ENSG00000008517', 'ENSG00000010278', 'ENSG00000010610']

Generate candidate lists¶

The example queries are the first three archived IDs in lexicographic order. They were chosen without consulting rankings or evaluation outcomes. The effect_gene_subset declaration preserves the distinction between the full reliability geometry and the displayed gene effects.

In [2]:
result = pr.from_statistics(
    data, queries, top_k=10, effect_gene_subset=True,
    weight_gene_count=int(data["weight_gene_count"]),
)
for method in ["Euclidean", "PertRank"]:
    np.testing.assert_allclose(
        result.weights[method].to_numpy(), data[("archived_weight_" if "archived_weight_" + method in data else "frozen_formula_weight_") + method],
        rtol=0, atol=1e-12,
    )
print("Both weight vectors match the saved implementation within 1e-12.")
query = queries[0]
shortlist = result.nominations.query(
    "method == 'PertRank' and query_id == @query and direction == 'up'"
)
display(shortlist[["rank", "target_id", "effect", "weight", "signed_score"]])
Both weight vectors match the saved implementation within 1e-12.
rank target_id effect weight signed_score
20 1 ENSG00000196843 0.310750 0.142306 0.044222
21 2 ENSG00000185811 0.181581 0.236137 0.042878
22 3 ENSG00000112182 0.355635 0.084425 0.030025
23 4 ENSG00000126456 0.280062 0.104852 0.029365
24 5 ENSG00000181894 0.165483 0.168286 0.027848
25 6 ENSG00000143379 0.156661 0.176335 0.027625
26 7 ENSG00000176472 0.299976 0.082506 0.024750
27 8 ENSG00000077150 0.105026 0.224170 0.023544
28 9 ENSG00000125347 0.104918 0.221143 0.023202
29 10 ENSG00000162676 0.099980 0.195021 0.019498

Inspect how reliability changes the shortlist¶

The scalar changes each candidate's competitive strength while the readout uses its original gene-specific effect. Compare PertRank with Unweighted to inspect the complete weighting workflow, and with Euclidean to inspect the added control geometry. Candidate movement describes this archived list; accuracy requires a separately defined reference.

In [3]:
raw = result.nominations.query(
    "method == 'Unweighted' and query_id == @query and direction == 'up'"
)
changes = pd.DataFrame({
    "added_by_PertRank": pd.Series(sorted(set(shortlist.target_id) - set(raw.target_id))),
    "removed_from_raw": pd.Series(sorted(set(raw.target_id) - set(shortlist.target_id))),
})
display(changes)
display(result.weights[["n_input_cells", "Euclidean", "PertRank"]].head())
added_by_PertRank removed_from_raw
0 ENSG00000077150 ENSG00000010818
1 ENSG00000125347 ENSG00000127528
2 ENSG00000162676 ENSG00000130844
n_input_cells Euclidean PertRank
target_id
ENSG00000105127 20.0 0.038611 0.055316
ENSG00000196843 20.0 0.323435 0.142306
ENSG00000166454 20.0 0.104323 0.090079
ENSG00000087152 20.0 0.126278 0.104914
ENSG00000112182 20.0 0.138561 0.084425

Reuse the weights on saved FR-Perturb estimates¶

The external estimates below come from a completed archived FR-Perturb fit. Labels align rows even when an upstream estimator returns them in a different order. The source weights and full candidate pool remain fixed, so the three returned arms compare weighting rules on the same FR effects. This reuses the reliability calculation across estimators; it does not transform the gene coordinate or train a new FR model.

fr_result.effects contains the FR readout. fr_result.statistics preserves the original source-cell statistics that generated the weights, keeping effect provenance and reliability provenance distinct.

In [4]:
external = pd.DataFrame(
    data["external_fr_effects"], index=data["target_ids"], columns=data["gene_ids"]
).iloc[::-1]
fr_result = result.rank_effects(external, queries=queries, top_k=10)
np.testing.assert_array_equal(fr_result.weights["PertRank"], result.weights["PertRank"])
display(fr_result.nominations.query(
    "method == 'PertRank' and query_id == @query and direction == 'up'"
)[["rank", "target_id", "effect", "weight", "signed_score"]])
rank target_id effect weight signed_score
20 1 ENSG00000077150 0.109191 0.224170 0.024477
21 2 ENSG00000185811 0.095389 0.236137 0.022525
22 3 ENSG00000125347 0.070550 0.221143 0.015602
23 4 ENSG00000126456 0.117504 0.104852 0.012321
24 5 ENSG00000176472 0.119042 0.082506 0.009822
25 6 ENSG00000010818 0.097811 0.088280 0.008635
26 7 ENSG00000112182 0.094404 0.084425 0.007970
27 8 ENSG00000118513 0.032298 0.232325 0.007504
28 9 ENSG00000143379 0.041232 0.176335 0.007271
29 10 ENSG00000196843 0.047757 0.142306 0.006796

Export for review¶

Up and down are signed expression directions. Boundary-tie inclusion is the chance of inclusion under uniform tie breaking, not a biological confidence probability. The complete score table retains eligible candidates with zero weight as well as candidates outside the shortlist.

In [5]:
out = Path("outputs/tcell")
out.mkdir(parents=True, exist_ok=True)
result.nominations.to_csv(out / "nominations.csv", index=False)
result.weights.to_csv(out / "weights.csv", index=False)
result.all_candidate_scores.to_csv(out / "all_candidate_scores.csv", index=False)
print("Exported", len(result.nominations), "nomination rows to", out)
Exported 180 nomination rows to outputs/tcell

Use your own AnnData¶

pr.rank(adata, queries, input_scale="counts", full_gene_axis=True, strata_keys=["batch"], top_k=25) is the cell-input entry point. Its input must contain the full measured raw gene axis, selected perturbation cells and matched controls in one declared stratum. Supply the metadata columns shown in the README. Already log-normalized data use input_scale="log" with an explicitly declared analysis axis.

The executable 00_anndata_quickstart.py demonstrates the complete cell-input route on a synthetic software fixture. It is separate from these authentic archived atlas statistics. Consult input contracts before using a new screen, and the full manuscript archive for scientific evaluation. The provenance JSON packaged with each fixture records its source files, their hashes, units and selection rules.