PGx CN: scan-date batch warning (PCN-006), CN-consistent BAF evidence (PCN-010, --bafs-in-cn),
noise-outlier exclusion + robustness test (PCN-017/EDGE-026), platform from the CSV manifest (CLI-032) Also: star-allele calls are checked for one manifest/build per run (PGX-011), duplications at the CN cap are named xN (PGX-022), and tandem names attach supporting variants to their full copy.
This commit is contained in:
kemal committed
1 parent
ed5e92b859
commit
1cb7b59559
8 files changed
+216
-25
No files matched your search
@@ -935,6 +935,7 @@ def pgx_cn_train_cmd(
|
||||
mask: Annotated[Path | None, typer.Option("--mask", help="Probe mask TSV (IN-052).")] = None,
|
||||
regions: Annotated[Path | None, typer.Option("--regions", help="Region BED (default: built-in v1).")] = None,
|
||||
summary: Annotated[Path | None, typer.Option("--summary", help="sample_summary.tsv (assay control).")] = None,
|
||||
csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (Assay Format → platform).")] = None,
|
||||
platform: Annotated[str, typer.Option("--platform", help="auto|LCG|EX")] = "auto",
|
||||
min_samples: Annotated[int, typer.Option("--min-samples", min=1)] = 96,
|
||||
min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 90,
|
||||
@@ -954,7 +955,8 @@ def pgx_cn_train_cmd(
|
||||
raise typer.BadParameter("--platform auto|LCG|EX")
|
||||
rep = _reporter("pgx copy-number train", out, log_json, debug, quiet)
|
||||
o = PgxTrainOptions(gtcs=gtcs, bpm=bpm, fasta=fasta, out=out, mask=mask, regions=regions, summary=summary,
|
||||
platform=platform, min_samples=min_samples, min_passing=min_passing, max_samples=max_samples,
|
||||
platform=platform, csv=csv, min_samples=min_samples, min_passing=min_passing,
|
||||
max_samples=max_samples,
|
||||
lrr_sd_max=train_lrr_sd_max, call_rate_min=train_call_rate_min,
|
||||
assay_control_min=train_assay_control_min, force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_train(o, rep))
|
||||
@@ -971,6 +973,9 @@ def pgx_cn_call_cmd(
|
||||
recommended_passing: Annotated[int, typer.Option("--recommended-passing", min=1)] = 96,
|
||||
lrr_sd_max: Annotated[float, typer.Option("--lrr-sd-max")] = 0.20,
|
||||
cluster_mismatch: Annotated[str, typer.Option("--cluster-mismatch", help="error|warn")] = "error",
|
||||
bafs_in_cn: Annotated[
|
||||
bool, typer.Option("--bafs-in-cn", help="Experimental: add CN-consistent BAF evidence to the CN posterior.")
|
||||
] = False,
|
||||
no_bgzip: Annotated[bool, typer.Option("--no-bgzip")] = False,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
@@ -987,7 +992,8 @@ def pgx_cn_call_cmd(
|
||||
rep = _reporter("pgx copy-number call", out, log_json, debug, quiet)
|
||||
o = PgxCallOptions(gtcs=gtcs, model=model, fasta=fasta, out=out, min_batch=min_batch, min_passing=min_passing,
|
||||
recommended_passing=recommended_passing, lrr_sd_max=lrr_sd_max,
|
||||
cluster_mismatch=cluster_mismatch, bgzip=not no_bgzip, force=force) # fmt: skip
|
||||
cluster_mismatch=cluster_mismatch, bafs_in_cn=bafs_in_cn, bgzip=not no_bgzip,
|
||||
force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_call(o, rep))
|
||||
|
||||
|
||||
|
||||
@@ -176,6 +176,7 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("E_PGX_VCF_UNSUITABLE", "error", "SNV VCF unsuitable for PGx (squashing/locus filters/foreign)", "PGX-011"),
|
||||
_m("E_PGX_GUIDELINE", "error", "Guideline tables not present in the PGx DB", "PGX-025"),
|
||||
_m("W_PGX_NO_PROBE_MAP", "warning", "The PGx DB has no probe map for this manifest", "PGX-035"),
|
||||
_m("E_PGX_PLATFORM_UNKNOWN", "error", "PGx CN platform (LCG/EX) unknown: give --csv or --platform", "CLI-032"),
|
||||
_m("E_PGX_SOURCE", "error", "PGx source tables missing; run scripts/pgx/fetch_sources.sh", "PGX-001"),
|
||||
# --- methylation control QC (MET-001 .. MET-018) -------------------------------------------------
|
||||
_m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"),
|
||||
|
||||
+30
-2
@@ -195,9 +195,37 @@ def posterior(x: float, fit: RegionFit) -> F64:
|
||||
return out
|
||||
|
||||
|
||||
def call_cn(x: float, fit: RegionFit) -> tuple[int, int, F64]:
|
||||
"""PCN-009 steps 5–6: (CN, QUAL, posterior)."""
|
||||
def baf_loglik(baf: F64, pfb: F64 | None = None, sd: float = 0.05) -> F64:
|
||||
"""PCN-010: per-CN-state log-likelihood of a region's SNP BAFs: for CN c the BAF is k/c with binomial weights
|
||||
C(c, k) p^k (1 − p)^(c − k) from the probe's population B-allele frequency p (0.5 when unknown); CN 0 → uniform
|
||||
[LIT: Wang et al. 2007 (PennCNV BAF emissions); OURS: region-level use]."""
|
||||
ok = np.isfinite(baf)
|
||||
b = baf[ok]
|
||||
p = np.clip(np.nan_to_num(pfb[ok] if pfb is not None else np.full(len(b), 0.5), nan=0.5), 0.01, 0.99)
|
||||
out = np.zeros(len(CNS))
|
||||
if len(b) == 0:
|
||||
return out
|
||||
for c in CNS:
|
||||
if c == 0:
|
||||
out[c] = 0.0 # uniform on [0, 1]
|
||||
continue
|
||||
k = np.arange(c + 1)
|
||||
means = k / c
|
||||
w = np.array([math.comb(c, int(j)) for j in k])[None, :] * p[:, None] ** k[None, :] * (1 - p[:, None]) ** (
|
||||
c - k[None, :]) # fmt: skip
|
||||
dens = np.exp(-0.5 * ((b[:, None] - means[None, :]) / sd) ** 2) / (sd * math.sqrt(2 * math.pi))
|
||||
out[c] = float(np.log(0.98 * (w * dens).sum(axis=1) + 0.02).sum())
|
||||
return out
|
||||
|
||||
|
||||
def call_cn(x: float, fit: RegionFit, baf_ll: F64 | None = None) -> tuple[int, int, F64]:
|
||||
"""PCN-009 steps 5–6: (CN, QUAL, posterior); ``baf_ll`` (PCN-010, experimental ``--bafs-in-cn``) is added to the
|
||||
log posterior."""
|
||||
p = posterior(x, fit)
|
||||
if baf_ll is not None:
|
||||
lp = np.log(np.maximum(p, 1e-300)) + baf_ll
|
||||
p = np.exp(lp - lp.max())
|
||||
p = p / p.sum()
|
||||
cn = int(np.argmax(p))
|
||||
pmax = float(p[cn])
|
||||
qual = 60 if pmax >= 1.0 else int(min(60, max(0, round(-10 * math.log10(max(1e-300, 1 - pmax))))))
|
||||
|
||||
@@ -31,7 +31,7 @@ from simurg import __version__
|
||||
from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.pgx.cn import CN_MAX, Region, RegionFit, call_cn, fit_mixture, parse_regions
|
||||
from simurg.pgx.cn import CN_MAX, Region, RegionFit, baf_loglik, call_cn, fit_mixture, parse_regions
|
||||
|
||||
MODEL_SCHEMA = 1
|
||||
REGIONS_FILE = "pgx_cn_regions.GRCh38.v1.bed"
|
||||
@@ -60,6 +60,7 @@ class PgxCnModel:
|
||||
probe_region: list[str] # region name, "" for scale probes
|
||||
probe_ref: F32 # training reference (median sample-centred LRR)
|
||||
fits: dict[str, RegionFit] = field(default_factory=dict)
|
||||
probe_pfb: F32 | None = None # training B-allele frequency per probe (NaN: no called genotypes), PCN-010
|
||||
|
||||
def region_sel(self, name: str) -> npt.NDArray[np.int64]:
|
||||
return np.flatnonzero(np.asarray(self.probe_region) == name)
|
||||
@@ -111,8 +112,10 @@ def region_values(d: npt.NDArray[np.float64], model: PgxCnModel) -> dict[str, tu
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ model file
|
||||
def write_model(path: Path, m: PgxCnModel) -> None:
|
||||
table = pa.table({"index": m.probe_idx, "region": m.probe_region, "ref": m.probe_ref,
|
||||
"name": m.header.pop("_names")}) # fmt: skip
|
||||
cols = {"index": m.probe_idx, "region": m.probe_region, "ref": m.probe_ref, "name": m.header.pop("_names")}
|
||||
if m.probe_pfb is not None:
|
||||
cols["pfb"] = m.probe_pfb
|
||||
table = pa.table(cols)
|
||||
buf = io.BytesIO()
|
||||
pq.write_table(table, buf, compression="zstd")
|
||||
mixture = {name: {"delta": f.delta, "gamma": f.gamma, "sigma": f.sigma, "weights": f.weights, "n": f.n,
|
||||
@@ -151,7 +154,9 @@ def read_model(path: Path) -> PgxCnModel:
|
||||
}
|
||||
return PgxCnModel(header, parse_regions(regions_text), regions_text,
|
||||
table.column("index").to_numpy().astype(np.int64), table.column("region").to_pylist(),
|
||||
table.column("ref").to_numpy().astype(np.float32), fits) # fmt: skip
|
||||
table.column("ref").to_numpy().astype(np.float32), fits,
|
||||
table.column("pfb").to_numpy().astype(np.float32)
|
||||
if "pfb" in table.column_names else None) # fmt: skip
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ train
|
||||
@@ -165,6 +170,7 @@ class PgxTrainOptions:
|
||||
regions: Path | None = None
|
||||
summary: Path | None = None
|
||||
platform: str = "auto"
|
||||
csv: Path | None = None # CSV manifest: source of the Assay Format heading (CLI-032)
|
||||
min_samples: int = 96
|
||||
min_passing: int = 90
|
||||
max_samples: int = 3000
|
||||
@@ -219,6 +225,17 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int:
|
||||
from simurg.pipeline.qc import write_tsv
|
||||
|
||||
man = read_bpm(o.bpm)
|
||||
csv_fmt = o.platform
|
||||
if o.platform == "auto":
|
||||
from simurg.formats.csv_manifest import assay_format, read_csv_manifest
|
||||
|
||||
csv_fmt = assay_format(read_csv_manifest(o.csv)) if o.csv else ""
|
||||
if csv_fmt not in ("LCG", "EX"):
|
||||
raise SimurgError(
|
||||
"E_PGX_PLATFORM_UNKNOWN",
|
||||
"cannot tell the platform (LCG/EX)" + (f": Assay Format {csv_fmt!r}" if csv_fmt else ": no --csv"),
|
||||
hint="give --csv or --platform LCG|EX",
|
||||
)
|
||||
out_model = o.out / f"{o.bpm.stem}.simurg-pgxcn"
|
||||
check_outputs_absent([out_model, o.out / "train_exclusions.tsv"], force=o.force)
|
||||
genome = Genome(o.fasta)
|
||||
@@ -249,6 +266,8 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int:
|
||||
rep.info(f"training on a seeded random subset of {o.max_samples} of {len(gtcs)} samples (PCN-005)")
|
||||
summ = _summary(o.summary or (o.gtcs / "sample_summary.tsv"))
|
||||
rows = np.full((len(keys), len(probe_idx)), np.nan, dtype=np.float32)
|
||||
alt_count = np.zeros(len(probe_idx))
|
||||
n_called = np.zeros(len(probe_idx))
|
||||
exclusions: list[dict[str, str]] = []
|
||||
passing = np.zeros(len(keys), dtype=np.bool_)
|
||||
cluster_files = set()
|
||||
@@ -259,6 +278,9 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int:
|
||||
f"{man.num_loci}") # fmt: skip
|
||||
cluster_files.add(g.cluster_file)
|
||||
rows[i] = _centred(g.logr_ratios, sc_idx)[probe_idx]
|
||||
gt = np.asarray(g.genotypes)[probe_idx]
|
||||
n_called += gt > 0
|
||||
alt_count += np.where(gt == 2, 1, np.where(gt == 3, 2, 0))
|
||||
reasons = []
|
||||
if not g.logr_dev <= o.lrr_sd_max:
|
||||
reasons.append(f"logr_dev {g.logr_dev:.3f} > {o.lrr_sd_max}")
|
||||
@@ -280,22 +302,32 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int:
|
||||
d = rows - ref
|
||||
sc_sel = np.flatnonzero(np.asarray(probe_region) == "")
|
||||
sc = np.nanmedian(d[:, sc_sel], axis=1)
|
||||
mad = float(np.nanmedian(np.abs(sc[passing] - np.nanmedian(sc[passing])))) * 1.4826 or 1e-9
|
||||
z = np.abs(sc - np.nanmedian(sc[passing])) / mad
|
||||
use = passing & (z <= 5) # PCN-009: outlier samples out of the reference
|
||||
|
||||
def robust_z(v: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
|
||||
med = float(np.nanmedian(v[passing]))
|
||||
mad = float(np.nanmedian(np.abs(v[passing] - med))) * 1.4826 or 1e-9
|
||||
out: npt.NDArray[np.float64] = np.abs(v - med) / mad
|
||||
return out
|
||||
|
||||
# PCN-009 / EDGE-026: outlier samples out of the reference — genome-wide level or genome-wide noise
|
||||
noise = 1.4826 * np.nanmedian(np.abs(d[:, sc_sel] - sc[:, None]), axis=1)
|
||||
z, zn = robust_z(sc), robust_z(noise)
|
||||
use = passing & (z <= 5) & (zn <= 5)
|
||||
for i in np.flatnonzero(passing & ~use).tolist():
|
||||
exclusions.append({"sample_id": keys[i], "reason": f"genome-wide NR outlier (|z| = {z[i]:.1f} > 5)"})
|
||||
exclusions.append({"sample_id": keys[i], "reason": f"genome-wide NR outlier (level |z| = {z[i]:.1f}, "
|
||||
f"noise |z| = {zn[i]:.1f}; limit 5)"}) # fmt: skip
|
||||
if (~use & passing).any():
|
||||
ref = np.nanmedian(rows[use], axis=0)
|
||||
d = rows - ref
|
||||
model = PgxCnModel({}, regions, reg_text, probe_idx, probe_region, ref.astype(np.float32))
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
pfb = np.where(n_called > 0, alt_count / (2 * n_called), np.nan).astype(np.float32)
|
||||
model = PgxCnModel({}, regions, reg_text, probe_idx, probe_region, ref.astype(np.float32), probe_pfb=pfb)
|
||||
values = {}
|
||||
for r in regions:
|
||||
if len(rp[r.name]) < MIN_REGION_PROBES:
|
||||
continue
|
||||
values[r.name] = np.array([region_values(d[i], model)[r.name][0] for i in np.flatnonzero(use)])
|
||||
model.fits = fit_mixture(values)
|
||||
csv_fmt = "unknown" if o.platform == "auto" else o.platform
|
||||
model.header = {
|
||||
"schema": "simurg.pgx_cn_model", "schema_version": MODEL_SCHEMA, "software_version": __version__,
|
||||
"manifest_name": man.manifest_name, "manifest_revision": getattr(man, "version", None),
|
||||
@@ -331,6 +363,7 @@ class PgxCallOptions:
|
||||
recommended_passing: int = 96
|
||||
lrr_sd_max: float = 0.20
|
||||
cluster_mismatch: str = "error"
|
||||
bafs_in_cn: bool = False # PCN-010 experimental
|
||||
bgzip: bool = True
|
||||
force: bool = False
|
||||
argv: list[str] = field(default_factory=list)
|
||||
@@ -359,8 +392,14 @@ def _vcf_header(sample: str, contigs: list[tuple[str, int]], meta: dict[str, str
|
||||
return "\n".join(h) + "\n"
|
||||
|
||||
|
||||
def call_sample_regions(d: npt.NDArray[np.float64], model: PgxCnModel, passing: bool) -> list[dict[str, Any]]:
|
||||
"""PCN-002/009: one result per region."""
|
||||
def call_sample_regions(
|
||||
d: npt.NDArray[np.float64],
|
||||
model: PgxCnModel,
|
||||
passing: bool,
|
||||
baf: npt.NDArray[np.float64] | None = None,
|
||||
bafs_in_cn: bool = False,
|
||||
) -> list[dict[str, Any]]:
|
||||
"""PCN-002/009/010: one result per region; ``baf`` = the sample's BAF at ``model.probe_idx``."""
|
||||
vals = region_values(d, model)
|
||||
out = []
|
||||
for r in model.regions:
|
||||
@@ -368,23 +407,44 @@ def call_sample_regions(d: npt.NDArray[np.float64], model: PgxCnModel, passing:
|
||||
n_all = len(model.region_sel(r.name))
|
||||
fit = model.fits.get(r.name)
|
||||
filters: list[str] = []
|
||||
bll = None
|
||||
n_snp = 0
|
||||
if baf is not None:
|
||||
sel = model.region_sel(r.name)
|
||||
rb = baf[sel]
|
||||
pf = model.probe_pfb[sel] if model.probe_pfb is not None else None
|
||||
if pf is not None:
|
||||
rb = np.where(np.isfinite(pf) & (pf > 0.01) & (pf < 0.99), rb, np.nan) # informative SNPs only
|
||||
n_snp = int(np.isfinite(rb).sum())
|
||||
bll = baf_loglik(rb, pf)
|
||||
if fit is None or n_all < MIN_REGION_PROBES or nf < MIN_FINITE_FRAC * n_all or not math.isfinite(x):
|
||||
cn, qual, nr = None, 0, math.nan
|
||||
filters.append("NoCall")
|
||||
else:
|
||||
cn, qual, _ = call_cn(x, fit)
|
||||
cn, qual, _ = call_cn(x, fit, bll if bafs_in_cn and n_snp else None)
|
||||
# NR ≈ CN / 2: the continuous CN that the region's mixture means imply for x [OURS]
|
||||
nr = max(0.0, 3.0 * 2.0 ** ((x - fit.delta) / fit.gamma) - 1.0) / 2.0 if cn else 0.0
|
||||
if qual < 7:
|
||||
filters.append("Q7")
|
||||
if not passing:
|
||||
filters.append("SampleQuality")
|
||||
out.append(
|
||||
{"region": r, "cn": cn, "qual": qual, "nr": nr, "x": x, "probes": nf, "filters": filters or ["PASS"]}
|
||||
)
|
||||
out.append({"region": r, "cn": cn, "qual": qual, "nr": nr, "x": x, "probes": nf, "filters": filters or ["PASS"],
|
||||
"n_snp": n_snp,
|
||||
"baf_ll": None if bll is None or not n_snp else [round(float(v), 3) for v in bll]}) # fmt: skip
|
||||
return out
|
||||
|
||||
|
||||
def _parse_date(text: str) -> Any:
|
||||
from datetime import datetime as _dt
|
||||
|
||||
for fmt in ("%m/%d/%Y %I:%M:%S %p", "%m/%d/%Y %H:%M:%S", "%Y-%m-%dT%H:%M:%S", "%Y-%m-%d %H:%M:%S"):
|
||||
try:
|
||||
return _dt.strptime(text.strip(), fmt)
|
||||
except ValueError:
|
||||
continue
|
||||
return None
|
||||
|
||||
|
||||
def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
"""``simurg pgx copy-number call`` (CLI-031)."""
|
||||
from simurg.core.atomic import atomic_write_text, check_outputs_absent
|
||||
@@ -406,13 +466,20 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
mcl = model.header.get("cluster_file_name")
|
||||
# pass 1: batch QC from GTC headers (PCN-004)
|
||||
heads = {}
|
||||
scan_dates: list[Any] = []
|
||||
for k, p in gtcs.items():
|
||||
g = read_gtc(p)
|
||||
heads[k] = (g.logr_dev, g.sample_name or k, g.cluster_file, len(g.logr_ratios))
|
||||
scan_dates.append(_parse_date(g.imaging_date))
|
||||
n_pass = sum(1 for v in heads.values() if v[0] < o.lrr_sd_max)
|
||||
if len(heads) < o.min_batch or n_pass < o.min_passing:
|
||||
raise SimurgError("E_PGX_BATCH_TOO_SMALL", f"{n_pass} sample(s) passed QC. Requires at least "
|
||||
f"{max(o.min_passing, 1)} samples to proceed.") # fmt: skip
|
||||
dates = [x for x in scan_dates if x is not None]
|
||||
if dates and (max(dates) - min(dates)).days > 14: # PCN-006 (sample-sheet lot/batch columns: no sheet input here)
|
||||
span = (max(dates) - min(dates)).days
|
||||
rep.warn("W_PGX_BATCH_COMPOSITION", f"scan dates span {span} days ({min(dates):%Y-%m-%d} … "
|
||||
f"{max(dates):%Y-%m-%d}); batch statistics may mix processing batches") # fmt: skip
|
||||
if n_pass < o.recommended_passing:
|
||||
rep.warn("W_PGX_BATCH_BELOW_RECOMMENDED", f"{n_pass} passing sample(s); {o.recommended_passing} recommended")
|
||||
for k, (_, _, cf, nl) in heads.items():
|
||||
@@ -445,6 +512,7 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
}
|
||||
centre_idx = model.probe_idx[model.region_sel("")] # the scale probes (autosomal), as in training
|
||||
o.out.mkdir(parents=True, exist_ok=True)
|
||||
baf_rows: list[dict[str, Any]] = []
|
||||
srows, per_region = [], {r.name: {"n_gain": 0, "n_loss": 0, "n_neutral": 0, "n_nocall": 0, "n_filtered": 0}
|
||||
for r in model.regions} # fmt: skip
|
||||
failed = 0
|
||||
@@ -453,13 +521,25 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
g = read_gtc(p)
|
||||
d = _centred(g.logr_ratios, centre_idx)[model.probe_idx] - model.probe_ref
|
||||
passing = heads[k][0] < o.lrr_sd_max
|
||||
res = call_sample_regions(d, model, passing)
|
||||
bafs = None
|
||||
if len(
|
||||
g.b_allele_freqs
|
||||
): # PCN-010: SNP probes only (called genotypes; intensity-only CNV probes never are)
|
||||
bafs = np.asarray(g.b_allele_freqs, dtype=np.float64)[model.probe_idx]
|
||||
bafs[np.asarray(g.genotypes)[model.probe_idx] == 0] = np.nan
|
||||
res = call_sample_regions(d, model, passing, bafs, o.bafs_in_cn)
|
||||
except SimurgError:
|
||||
raise
|
||||
except Exception as e: # PCN: one sample never stops the batch
|
||||
failed += 1
|
||||
rep.error("E_PGX_CN_FAILED", f"{k}: {e}", sample_id=k)
|
||||
continue
|
||||
for x in res:
|
||||
if x["baf_ll"] is not None:
|
||||
bl = x["baf_ll"]
|
||||
baf_rows.append({"sample_id": k, "region": x["region"].name, "n_snp": x["n_snp"],
|
||||
**{f"ll_cn{c}": bl[c] for c in range(len(bl))},
|
||||
"best_cn_by_baf": int(np.argmax(bl[1:]) + 1)}) # fmt: skip
|
||||
called = [(x["region"], x["cn"]) for x in res if x["cn"] is not None]
|
||||
tot = sum(r.end - r.start0 for r, _ in called)
|
||||
ploidy = sum((r.end - r.start0) * min(cn, CN_MAX) for r, cn in called) / tot if tot else math.nan
|
||||
@@ -502,6 +582,8 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
scols = ["sample_id", "sample_name", "sample_folder", "logr_dev", "passed_batch_qc", "used_in_reference"]
|
||||
scols += [f"cn_{r.name}" for r in model.regions] + ["flags"]
|
||||
write_tsv(o.out / "pgx_cn_sample_summary.tsv", scols, srows)
|
||||
write_tsv(o.out / "pgx_cn_baf_evidence.tsv", ["sample_id", "region", "n_snp"] + [f"ll_cn{c}" for c in range(6)]
|
||||
+ ["best_cn_by_baf"], baf_rows) # PCN-010: for the star-allele caller and review # fmt: skip
|
||||
run = {"software": f"SimurgArray {__version__}", **meta, "samples": len(gtcs), "failed": failed}
|
||||
atomic_write_text(o.out / "pgx_cn_run.json", json.dumps(run, indent=1))
|
||||
rep.info(f"PGx CN: {len(srows)}/{len(gtcs)} sample(s) called")
|
||||
|
||||
@@ -11,6 +11,7 @@ import contextlib
|
||||
import gzip
|
||||
import json
|
||||
import math
|
||||
import re
|
||||
import sqlite3
|
||||
from collections import defaultdict
|
||||
from dataclasses import dataclass, field
|
||||
@@ -205,7 +206,9 @@ def solution_detail(m: GeneModel, obs: Observations, recs: list[dict[str, Any] |
|
||||
hap_of = {m.haplotypes[i]: i for i in s.copies}
|
||||
reason = {NO_CALL: "no_call", MASKED: "masked", NOT_COVERED: "not_on_array"}
|
||||
for name in names:
|
||||
core = name.split("x")[0] if "x" in name and name.split("x")[-1].isdigit() else name
|
||||
unit = name.split("+")[-1] # "*68+*4" → the full copy "*4"; "*2x2"/"*1xN" → "*2"/"*1"
|
||||
mx = re.fullmatch(r"(.+?)x(\d+|N)", unit)
|
||||
core = mx.group(1) if mx else unit
|
||||
h = hap_of.get(core)
|
||||
sup, unres = [], []
|
||||
if h is not None:
|
||||
@@ -300,6 +303,7 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int:
|
||||
unres_rows: list[dict[str, Any]] = []
|
||||
not_cov_rows: dict[tuple[str, str], dict[str, Any]] = {}
|
||||
calls_json = []
|
||||
first_ident: tuple[str, str] | None = None
|
||||
for sid, path in sorted(snv.items()):
|
||||
recs, meta = read_vcf_sites(path, wanted)
|
||||
opts = meta.get("simurg_options", "")
|
||||
@@ -313,6 +317,12 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int:
|
||||
or kv.get("include_loci", "none").lower() != "none"
|
||||
): # PGX-011
|
||||
raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: made with unsquashed duplicates or a locus filter")
|
||||
ident = (meta.get("simurg_manifest", "").split(";")[0], meta.get("simurg_genome_build", ""))
|
||||
if first_ident is None:
|
||||
first_ident = ident
|
||||
elif ident != first_ident and not o.allow_foreign_vcf: # PGX-011 (4)
|
||||
raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: manifest/build {ident} differs from {first_ident}",
|
||||
hint="call one manifest and build per run") # fmt: skip
|
||||
cn = {} if o.no_cn or sid not in cnv else read_cn_vcf(cnv[sid])
|
||||
masked = masks.get(meta.get("simurg_manifest", "").split(";")[0], set())
|
||||
sample_name = sid
|
||||
@@ -375,6 +385,12 @@ def call_one_gene(
|
||||
feasible_ratio=o.feasible_ratio, structure=structure) # fmt: skip
|
||||
if not sols:
|
||||
return None, [], [], [], []
|
||||
if gcn is not None and gcn >= 5: # PGX-022: CN 5 means "5 or more" → the duplication count is unknown (xN)
|
||||
import re as _re
|
||||
|
||||
for s_ in sols:
|
||||
s_.haplotypes = tuple(_re.sub(r"x\d+$", "xN", h) for h in s_.haplotypes)
|
||||
s_.members = [tuple(_re.sub(r"x\d+$", "xN", h) for h in mm) for mm in s_.members]
|
||||
top = sols[0]
|
||||
details = [solution_detail(m, obs, vrec, s_, top, o.low_gs, MAX_EDIT_DISTANCE) for s_ in sols]
|
||||
flags = []
|
||||
|
||||
Reference in new issue
Block a user