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
This commit is contained in:
1 parent
e6cf273ab4
commit
d11226ca31
9 files changed
+792
-29
No files matched your search
+19
-5
@@ -284,9 +284,18 @@ def genotype_cmd(
|
||||
sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None,
|
||||
csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None,
|
||||
cutoff: Annotated[float, typer.Option("--cutoff", help="GenCall score below which a call is a no-call.")] = 0.15,
|
||||
norm_version: Annotated[str, typer.Option("--norm-version", help="1.2.0 (current) or 1.1.2")] = "1.2.0",
|
||||
sex_version: Annotated[int, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule (1 or 2).")] = 2,
|
||||
gtc_version: Annotated[int, typer.Option("--gtc-version", help="GTC format version (3 or 5).")] = 5,
|
||||
preset: Annotated[
|
||||
str, typer.Option("--preset", help="current (default) or autocall-1.6.2 (legacy vendor behaviour)")
|
||||
] = "current",
|
||||
norm_version: Annotated[
|
||||
str | None, typer.Option("--norm-version", help="1.2.0 or 1.1.2 (default: from --preset)")
|
||||
] = None,
|
||||
sex_version: Annotated[
|
||||
int | None, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule 1 or 2 (default: preset)")
|
||||
] = None,
|
||||
gtc_version: Annotated[
|
||||
int | None, typer.Option("--gtc-version", help="GTC format 3 or 5 (default: preset)")
|
||||
] = None,
|
||||
autocall_date: Annotated[
|
||||
str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).")
|
||||
] = None,
|
||||
@@ -303,9 +312,13 @@ def genotype_cmd(
|
||||
|
||||
if layout not in ("nested", "flat"):
|
||||
raise typer.BadParameter("--layout must be nested or flat")
|
||||
if gtc_version not in (3, 5):
|
||||
from simurg.genotype.sample import PRESETS
|
||||
|
||||
if preset not in PRESETS:
|
||||
raise typer.BadParameter(f"--preset must be one of {sorted(PRESETS)}")
|
||||
if gtc_version not in (None, 3, 5):
|
||||
raise typer.BadParameter("--gtc-version must be 3 or 5")
|
||||
if norm_version not in ("1.2.0", "1.1.2"):
|
||||
if norm_version not in (None, "1.2.0", "1.1.2"):
|
||||
raise typer.BadParameter("--norm-version must be 1.2.0 or 1.1.2")
|
||||
rep = Reporter(
|
||||
command="genotype",
|
||||
@@ -327,6 +340,7 @@ def genotype_cmd(
|
||||
layout=layout,
|
||||
on_duplicate_idat=on_duplicate_idat,
|
||||
cutoff=cutoff,
|
||||
preset=preset,
|
||||
norm_version=norm_version,
|
||||
sex_version=sex_version,
|
||||
gtc_version=gtc_version,
|
||||
|
||||
@@ -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,
|
||||
}
|
||||
@@ -27,7 +27,14 @@ from simurg.formats.gtc import write_gtc
|
||||
from simurg.formats.idat import check_pair, read_idat
|
||||
from simurg.genotype.baf_lrr import chrom_core
|
||||
from simurg.genotype.clusters import align_clusters
|
||||
from simurg.genotype.sample import GenotypeParams, ManifestContext, genotype_sample, manifest_context
|
||||
from simurg.genotype.sample import (
|
||||
GenotypeParams,
|
||||
ManifestContext,
|
||||
genotype_sample,
|
||||
manifest_context,
|
||||
params_for_preset,
|
||||
summary_stats,
|
||||
)
|
||||
from simurg.genotype.sex import SexParams
|
||||
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
|
||||
from simurg.signal.normalization.apply import normalize, polar
|
||||
@@ -77,9 +84,10 @@ class GenotypeOptions:
|
||||
layout: str = "nested"
|
||||
on_duplicate_idat: str = "error"
|
||||
cutoff: float = 0.15
|
||||
norm_version: str = "1.2.0"
|
||||
sex_version: int = 2
|
||||
gtc_version: int = 5
|
||||
preset: str = "current"
|
||||
norm_version: str | None = None # None: from the preset
|
||||
sex_version: int | None = None
|
||||
gtc_version: int | None = None
|
||||
autocall_date: str | None = None
|
||||
assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
|
||||
argv: list[str] = field(default_factory=list)
|
||||
@@ -133,11 +141,7 @@ def _process(job: Job) -> dict[str, Any]:
|
||||
call_rate=float(g.call_rate),
|
||||
logr_dev=float(g.logr_dev),
|
||||
sex_est=g.gender_char,
|
||||
gc10=float(g.gc10),
|
||||
gc50=float(g.gc50),
|
||||
n_calls=g.num_calls,
|
||||
n_nocalls=g.num_no_calls,
|
||||
n_intensity_only=g.num_intensity_only,
|
||||
**summary_stats(g, ctx.manifest.intensity_only, params.cutoff),
|
||||
p05_x=int(g.percentiles_x[0]),
|
||||
p50_x=int(g.percentiles_x[1]),
|
||||
p95_x=int(g.percentiles_x[2]),
|
||||
@@ -261,13 +265,17 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
auto_1_22 = np.array([chrom_core(c) in {str(i) for i in range(1, 23)} for c in manifest.chroms])
|
||||
name_idx = manifest.name_index()
|
||||
probes = {p: name_idx[p] for p in opts.assay_control_probes if p in name_idx}
|
||||
params = GenotypeParams(
|
||||
params = params_for_preset(
|
||||
opts.preset,
|
||||
cutoff=opts.cutoff,
|
||||
norm_version=opts.norm_version,
|
||||
sex=SexParams(version=opts.sex_version),
|
||||
sex=SexParams(version=opts.sex_version) if opts.sex_version else None,
|
||||
gtc_version=opts.gtc_version,
|
||||
autocall_date=opts.autocall_date,
|
||||
)
|
||||
rm.set_option("effective_norm_version", params.norm_version, f"preset:{opts.preset}")
|
||||
rm.set_option("effective_sex_version", params.sex.version, f"preset:{opts.preset}")
|
||||
rm.set_option("effective_gtc_version", params.gtc_version, f"preset:{opts.preset}")
|
||||
started = _dt.datetime.now()
|
||||
|
||||
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(jobs) or 1))
|
||||
@@ -335,7 +343,7 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
"warnings": rep.warning_counts(),
|
||||
}
|
||||
rm.add_output(summary_tsv, SUMMARY_SCHEMA)
|
||||
rm.schema_versions["*.gtc"] = f"gtc/{opts.gtc_version}"
|
||||
rm.schema_versions["*.gtc"] = f"gtc/{params.gtc_version}"
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / "run_manifest.json")
|
||||
rep.info(f"done: {ok}/{len(out_rows)} samples genotyped; summary {summary_tsv}")
|
||||
|
||||
Reference in new issue
Block a user