P4 validated on 2141 vendor GTCs; legacy AutoCall 1.6.2 preset
- vendor comparison (1000 Genomes Omni2.5, AutoCall 1.6.2.2): every computed field identical in 2138/2141 samples; the 3 others trace back to the single-locus raw-intensity rule (K-009); sex rule 1 matches all - AutoCall 1.6.2.x leaves GC10/GC50/#calls/#no-calls at 0 in GTC v3: new `--preset autocall-1.6.2` (1.1.2 normalization, sex rule 1, GTC v3, legacy statistics); the sample summary keeps the real statistics - NaN sign canonicalised to the x86 pattern (K-014); report P4.md Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
1 parent
eb19ac7f91
commit
8b1e5d8dab
9 files changed
+792
-29
No files matched your search
@@ -39,6 +39,27 @@ class GenotypeParams:
|
||||
gtc_version: int = 5
|
||||
autocall_version: str = f"{PRODUCT_NAME} {__version__}"
|
||||
autocall_date: str | None = None # None: run start time (invariant format, OUT-001)
|
||||
legacy_sample_stats: bool = False # AutoCall 1.6.2.x left GC10, GC50, #calls and #no-calls at 0 (K-015)
|
||||
|
||||
|
||||
PRESETS: dict[str, dict[str, object]] = {
|
||||
# default: current vendor software (normalization 1.2.0 / GenTrain 3, sex rule 2, GTC v5)
|
||||
"current": {},
|
||||
# AutoCall 1.6.2.x as used for the public 1000 Genomes Omni2.5 GTCs (validated on 2141 samples, P4)
|
||||
"autocall-1.6.2": {
|
||||
"norm_version": "1.1.2",
|
||||
"sex": SexParams(version=1),
|
||||
"gtc_version": 3,
|
||||
"legacy_sample_stats": True,
|
||||
},
|
||||
}
|
||||
|
||||
|
||||
def params_for_preset(name: str, **overrides: object) -> GenotypeParams:
|
||||
if name not in PRESETS:
|
||||
raise ValueError(f"unknown preset {name!r}; choose from {sorted(PRESETS)}")
|
||||
values = {**PRESETS[name], **{k: v for k, v in overrides.items() if v is not None}}
|
||||
return GenotypeParams(**values) # type: ignore[arg-type]
|
||||
|
||||
|
||||
@dataclass
|
||||
@@ -193,4 +214,23 @@ def genotype_sample(
|
||||
g.sample_data["num_calls"] = calls.num_calls
|
||||
g.sample_data["num_no_calls"] = m.num_loci - calls.num_intensity_only - calls.num_calls
|
||||
g.sample_data["num_intensity_only"] = calls.num_intensity_only
|
||||
if params.legacy_sample_stats:
|
||||
g.gc10 = 0.0
|
||||
g.sample_data["p50gc"] = 0.0
|
||||
g.sample_data["num_calls"] = 0
|
||||
g.sample_data["num_no_calls"] = 0
|
||||
return g
|
||||
|
||||
|
||||
def summary_stats(g: Gtc, intensity_only: npt.NDArray[np.bool_], cutoff: float) -> dict[str, float | int]:
|
||||
"""Sample statistics recomputed from calls and scores (independent of legacy GTC field conventions)."""
|
||||
gs = np.sort(g.genotype_scores[(g.genotype_scores > F32(cutoff)) & ~intensity_only])
|
||||
n_calls = int((g.genotypes > 0).sum())
|
||||
n_io = int(intensity_only.sum())
|
||||
return {
|
||||
"gc10": float(percentile(gs, 10)),
|
||||
"gc50": float(percentile(gs, 50)),
|
||||
"n_calls": n_calls,
|
||||
"n_nocalls": g.num_snps - n_io - n_calls,
|
||||
"n_intensity_only": n_io,
|
||||
}
|
||||
Reference in new issue
Block a user