P4 (WIP): GenCall-equivalent calling, BAF/LRR, sex, GTC v5, simurg genotype
- genotype/clusters.py: EGT aligned to manifest order (float32 mean R) - genotype/gencall.py: flair, six-case likelihood, priors/posteriors, fuzzy sub-scores and score map, vectorised with the reference's float/double evaluation order and libm transcendentals - genotype/baf_lrr.py, sex.py (rules v1/v2), sample.py (IDAT -> Gtc) - pipeline/inputs.py shared by ingest and genotype - pipeline/genotype.py + CLI: GTC per sample, sample_summary.tsv/json (OUT-060..063), run manifest - GTCs from IDATs are byte-identical to the oracle for 36 GSA samples x 2 normalization versions (NaN sign canonicalised to the x86 pattern) Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
kemalandClaude Opus 5.5 committed
1 parent
0cd7eb8aa0
commit
febdd121ef
12 files changed
+1270
-53
No files matched your search
@@ -0,0 +1,77 @@
|
||||
"""B-allele frequency, log R ratio and LogR deviation — plan §7.6 (Peiffer et al. 2006).
|
||||
|
||||
Exact mode follows gtc2vcf ``get_baf_lrr``/``calculate_baf_lrr`` (MIT): float32 interpolation between the cluster
|
||||
θ/R means, ``log2`` in double of a float32 ratio, NaN where θ is undefined, BAF clamped to [0, 1] with NaN kept.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.genotype.clusters import ClusterArrays
|
||||
from simurg.genotype.gencall import _libm
|
||||
|
||||
F32 = np.float32
|
||||
F32Arr = npt.NDArray[np.float32]
|
||||
NEG_NAN = np.frombuffer(np.uint32(0xFFC00000).tobytes(), dtype=np.float32)[0] # C's -NAN
|
||||
|
||||
|
||||
def baf_lrr(
|
||||
theta: F32Arr, r: F32Arr, raw_zero: npt.NDArray[np.bool_], c: ClusterArrays, intensity_only: npt.NDArray[np.bool_]
|
||||
) -> tuple[F32Arr, F32Arr]:
|
||||
t = np.asarray(theta, np.float32)
|
||||
rr = np.asarray(r, np.float32)
|
||||
at, bt, ct = c.theta_mean[:, 0], c.theta_mean[:, 1], c.theta_mean[:, 2]
|
||||
ar, br, cr = c.r_mean[:, 0], c.r_mean[:, 1], c.r_mean[:, 2]
|
||||
with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
|
||||
eq = t == bt
|
||||
lo = t < bt
|
||||
hi = t > bt
|
||||
r_lo = (ar + ((t - at) * (ar - br)) / (at - bt)).astype(np.float32)
|
||||
b_lo = (((t - at) / (bt - at)) * F32(0.5)).astype(np.float32)
|
||||
r_hi = (br + ((t - bt) * (cr - br)) / (ct - bt)).astype(np.float32)
|
||||
b_hi = (F32(0.5) + ((t - bt) / (ct - bt)) * F32(0.5)).astype(np.float32)
|
||||
r_ref = np.where(eq, br, np.where(lo, r_lo, np.where(hi, r_hi, np.nan))).astype(np.float32)
|
||||
baf = np.where(eq, F32(0.5), np.where(lo, b_lo, np.where(hi, b_hi, NEG_NAN))).astype(np.float32)
|
||||
denom = np.where(intensity_only, c.mean_r, r_ref).astype(np.float32)
|
||||
ratio = (rr / denom).astype(np.float32).astype(np.float64)
|
||||
lrr = _libm(math.log2, ratio).astype(np.float32)
|
||||
lrr = np.where(rr != F32(0.0), lrr, -np.finfo(np.float32).max).astype(np.float32)
|
||||
undefined = ~(eq | lo | hi)
|
||||
lrr[undefined] = NEG_NAN
|
||||
baf = np.where(baf < 0.0, F32(0.0), np.where(baf > 1.0, F32(1.0), baf)).astype(np.float32)
|
||||
skip = raw_zero | ~c.present
|
||||
baf[skip] = NEG_NAN
|
||||
lrr[skip] = NEG_NAN
|
||||
# canonical NaN bit pattern (0xFFC00000), as written by the reference on every platform
|
||||
baf[np.isnan(baf)] = NEG_NAN
|
||||
lrr[np.isnan(lrr)] = NEG_NAN
|
||||
return baf, lrr
|
||||
|
||||
|
||||
def chrom_core(chrom: str) -> str:
|
||||
return chrom[3:] if chrom[:3].upper() == "CHR" else chrom
|
||||
|
||||
|
||||
def autosomal_mask(chroms: list[str]) -> npt.NDArray[np.bool_]:
|
||||
"""LogR-dev locus set: first character (after an optional 'chr') not X, Y or M (so '0' and alt contigs count)."""
|
||||
return np.array([(chrom_core(ch)[:1].upper() not in ("X", "Y", "M")) for ch in chroms], dtype=np.bool_)
|
||||
|
||||
|
||||
def logr_dev(
|
||||
lrr: F32Arr, autosomal: npt.NDArray[np.bool_], intensity_only: npt.NDArray[np.bool_], present: npt.NDArray[np.bool_]
|
||||
) -> np.float32:
|
||||
"""Population SD of LRR over autosomal, non-intensity-only, finite values (sequential double sums)."""
|
||||
sel = autosomal & ~intensity_only & present & np.isfinite(lrr)
|
||||
v = lrr[sel].astype(np.float64)
|
||||
count = len(v)
|
||||
if count == 0:
|
||||
return F32(np.nan)
|
||||
s = float(np.cumsum(v)[-1])
|
||||
s2 = float(np.cumsum(v * v)[-1])
|
||||
mean = s / count
|
||||
var = s2 / count - mean * mean
|
||||
return F32(math.sqrt(var)) if var >= 0 else F32(np.nan)
|
||||
@@ -0,0 +1,75 @@
|
||||
"""Cluster-file records aligned to manifest (GTC) locus order — input of genotype calling (plan §7.5–§7.6)."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import SimurgError
|
||||
from simurg.formats.bpm import Manifest
|
||||
from simurg.formats.egt import ClusterFile
|
||||
|
||||
F32Arr = npt.NDArray[np.float32]
|
||||
|
||||
|
||||
@dataclass
|
||||
class ClusterArrays:
|
||||
"""Per-locus cluster parameters; (n, 3) arrays are ordered (AA, AB, BB)."""
|
||||
|
||||
present: npt.NDArray[np.bool_] # locus has a record in the cluster file
|
||||
n: npt.NDArray[np.int32]
|
||||
theta_mean: F32Arr
|
||||
theta_dev: F32Arr
|
||||
r_mean: F32Arr
|
||||
r_dev: F32Arr # floored at 0.1 as the caller does (plan §7.5 step 2)
|
||||
intensity_threshold: F32Arr # NaN where the cluster file has none
|
||||
gentrain: F32Arr
|
||||
mean_r: F32Arr # count-weighted mean R in float32 (expected R of intensity-only loci)
|
||||
|
||||
@property
|
||||
def num_loci(self) -> int:
|
||||
return len(self.present)
|
||||
|
||||
|
||||
def align_clusters(egt: ClusterFile, manifest: Manifest, *, allow_missing: bool = False) -> ClusterArrays:
|
||||
"""Join cluster records to the manifest by locus name (IN-043)."""
|
||||
idx = egt.name_index()
|
||||
rows = np.array([idx.get(nm, -1) for nm in manifest.names], dtype=np.int64)
|
||||
present = rows >= 0
|
||||
if not present.all() and not allow_missing:
|
||||
missing = [manifest.names[i] for i in np.flatnonzero(~present)[:10]]
|
||||
raise SimurgError(
|
||||
"E_EGT_MANIFEST_MISMATCH",
|
||||
f"{int((~present).sum())} manifest loci are missing from cluster file {egt.path}, "
|
||||
f"e.g. {', '.join(missing)}",
|
||||
hint="use the cluster file built for this manifest",
|
||||
)
|
||||
safe = np.where(present, rows, 0)
|
||||
|
||||
def take(a: npt.NDArray, fill: float = np.nan) -> npt.NDArray: # type: ignore[type-arg]
|
||||
out = a[safe].copy()
|
||||
out[~present] = fill
|
||||
return out
|
||||
|
||||
n = take(egt.n, 0).astype(np.int32)
|
||||
r_mean = take(egt.r_mean).astype(np.float32)
|
||||
r_dev = take(egt.r_dev).astype(np.float32)
|
||||
r_dev = np.where(r_dev < np.float32(0.1), np.float32(0.1), r_dev).astype(np.float32)
|
||||
# float32 weighted mean exactly as the reference: ((Naa*Raa + Nab*Rab) + Nbb*Rbb) / (Naa+Nab+Nbb)
|
||||
nf = n.astype(np.float32)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
num = (nf[:, 0] * r_mean[:, 0] + nf[:, 1] * r_mean[:, 1]) + nf[:, 2] * r_mean[:, 2]
|
||||
mean_r = (num / (n.sum(axis=1)).astype(np.float32)).astype(np.float32)
|
||||
return ClusterArrays(
|
||||
present=present,
|
||||
n=n,
|
||||
theta_mean=take(egt.theta_mean).astype(np.float32),
|
||||
theta_dev=take(egt.theta_dev).astype(np.float32),
|
||||
r_mean=r_mean,
|
||||
r_dev=r_dev,
|
||||
intensity_threshold=take(egt.intensity_threshold).astype(np.float32),
|
||||
gentrain=take(egt.gentrain_score).astype(np.float32),
|
||||
mean_r=mean_r,
|
||||
)
|
||||
@@ -0,0 +1,270 @@
|
||||
"""GenCall-equivalent genotype calling from a cluster file — plan §7.5.
|
||||
|
||||
Vectorised port of gtc2vcf ``gen_std_flair``/``modilik``/``compute_score_call_prelim``/``gencall_score_map``/
|
||||
``make_calls`` (MIT, © Giulio Genovese), which reproduce the vendor's GenTrain 6.x calling (concepts: Kermani
|
||||
patents, expired; Illumina GenCall technical note). Every expression keeps the reference's float/double types and
|
||||
evaluation order; scalar transcendental functions come from libm through ``math`` where NumPy's SIMD kernels could
|
||||
round differently.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.genotype.clusters import ClusterArrays
|
||||
|
||||
F32 = np.float32
|
||||
F32Arr = npt.NDArray[np.float32]
|
||||
F64Arr = npt.NDArray[np.float64]
|
||||
|
||||
M_2_SQRTPI = 1.12837916709551257390
|
||||
M_SQRT1_2 = 0.707106781186547524401
|
||||
PEAK = 0.5 * M_2_SQRTPI * M_SQRT1_2 # constant sub-expression of the normal pdf (evaluated left to right)
|
||||
TINY_D = float(np.finfo(np.float64).tiny) * float(np.finfo(np.float64).eps) # DBL_MIN * DBL_EPSILON
|
||||
DEFAULT_CUTOFF = 0.15
|
||||
|
||||
|
||||
def _c_call(fn, v: float) -> float: # type: ignore[no-untyped-def]
|
||||
"""Call a ``math`` function with C semantics: domain/range errors give NaN/±inf instead of raising."""
|
||||
if v != v:
|
||||
return v
|
||||
try:
|
||||
return float(fn(v))
|
||||
except ValueError:
|
||||
if fn is math.log10 and v == 0.0:
|
||||
return -math.inf
|
||||
return math.nan
|
||||
except OverflowError:
|
||||
return math.inf
|
||||
|
||||
|
||||
def _libm(fn, x: F64Arr) -> F64Arr: # type: ignore[no-untyped-def]
|
||||
"""Element-wise libm call (bit-identical to C)."""
|
||||
return np.fromiter((_c_call(fn, v) for v in x.tolist()), dtype=np.float64, count=len(x))
|
||||
|
||||
|
||||
def _pow(e: float): # type: ignore[no-untyped-def]
|
||||
def f(v: float) -> float:
|
||||
return math.pow(v, e)
|
||||
|
||||
return f
|
||||
|
||||
|
||||
def smf(x: F32Arr, a: F32Arr | np.float32, b: F32Arr | np.float32) -> F32Arr:
|
||||
"""MATLAB ``smf`` in float32, including the reference's ``a >= b`` rule."""
|
||||
x = np.asarray(x, np.float32)
|
||||
a = np.broadcast_to(np.asarray(a, np.float32), x.shape)
|
||||
b = np.broadcast_to(np.asarray(b, np.float32), x.shape)
|
||||
with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
|
||||
mid = ((a + b) / F32(2.0)).astype(np.float32)
|
||||
u = ((x - a) / (b - a)).astype(np.float32)
|
||||
v = ((x - b) / (b - a)).astype(np.float32)
|
||||
rise = (F32(2.0) * (u * u)).astype(np.float32)
|
||||
top = (F32(1.0) - F32(2.0) * (v * v)).astype(np.float32)
|
||||
out = np.where(
|
||||
x <= a, F32(0), np.where((a < x) & (x <= mid), rise, np.where((mid < x) & (x <= b), top, F32(1)))
|
||||
).astype(np.float32)
|
||||
degenerate = a >= b
|
||||
out = np.where(degenerate, np.where(x >= mid, F32(1), F32(0)), out).astype(np.float32)
|
||||
return out
|
||||
|
||||
|
||||
def zmf(x: F32Arr, a: F32Arr | np.float32, b: F32Arr | np.float32) -> F32Arr:
|
||||
"""MATLAB ``zmf`` in float32 (callers guarantee a < b)."""
|
||||
x = np.asarray(x, np.float32)
|
||||
a = np.broadcast_to(np.asarray(a, np.float32), x.shape)
|
||||
b = np.broadcast_to(np.asarray(b, np.float32), x.shape)
|
||||
with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
|
||||
mid = ((a + b) / F32(2.0)).astype(np.float32)
|
||||
u = ((x - a) / (b - a)).astype(np.float32)
|
||||
v = ((x - b) / (b - a)).astype(np.float32)
|
||||
fall = (F32(1.0) - F32(2.0) * (u * u)).astype(np.float32)
|
||||
tail = (F32(2.0) * (v * v)).astype(np.float32)
|
||||
return np.where(
|
||||
x <= a, F32(1), np.where((a < x) & (x <= mid), fall, np.where((mid < x) & (x <= b), tail, F32(0)))
|
||||
).astype(np.float32)
|
||||
|
||||
|
||||
@dataclass
|
||||
class _Pdf:
|
||||
"""Normal pdf pieces for one cluster and one dimension: the full density and the flat peak."""
|
||||
|
||||
x: F32Arr
|
||||
mu: F32Arr
|
||||
full: F64Arr
|
||||
peak: F64Arr
|
||||
bad: npt.NDArray[np.bool_]
|
||||
|
||||
def get(self, side: str) -> F64Arr:
|
||||
if side == "left":
|
||||
out = np.where(self.x < self.mu, self.peak, self.full)
|
||||
elif side == "right":
|
||||
out = np.where(self.x > self.mu, self.peak, self.full)
|
||||
else:
|
||||
out = self.full
|
||||
return np.where(self.bad, np.nan, out)
|
||||
|
||||
|
||||
def _pdf(x: F32Arr, mu: F32Arr, sigma: F32Arr) -> _Pdf:
|
||||
"""``exp(-0.5*((x-mu)/sigma)^2) * 0.5 * M_2_SQRTPI * M_SQRT1_2 / sigma`` with float32 z, double arithmetic."""
|
||||
with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
|
||||
z = ((x - mu) / sigma).astype(np.float32).astype(np.float64)
|
||||
e = _libm(math.exp, -0.5 * (z * z))
|
||||
sd = sigma.astype(np.float64)
|
||||
full = e * 0.5 * M_2_SQRTPI * M_SQRT1_2 / sd
|
||||
peak = PEAK / sd
|
||||
return _Pdf(x, mu, full, peak, sigma <= F32(0))
|
||||
|
||||
|
||||
@dataclass
|
||||
class Calls:
|
||||
genotypes: npt.NDArray[np.uint8] # 0 NC, 1 AA, 2 AB, 3 BB (after the cutoff)
|
||||
scores: F32Arr # GenCall scores (0 where not computable)
|
||||
num_calls: int
|
||||
num_intensity_only: int
|
||||
|
||||
|
||||
def call_genotypes(
|
||||
theta: F32Arr,
|
||||
r: F32Arr,
|
||||
raw_zero: npt.NDArray[np.bool_],
|
||||
c: ClusterArrays,
|
||||
intensity_only: npt.NDArray[np.bool_],
|
||||
cutoff: float = DEFAULT_CUTOFF,
|
||||
) -> Calls:
|
||||
"""Bayesian calls and GenCall scores for one sample (plan §7.5)."""
|
||||
n = len(theta)
|
||||
t = np.asarray(theta, np.float32)
|
||||
rr = np.asarray(r, np.float32)
|
||||
tm, td, rm, rd = c.theta_mean, c.theta_dev, c.r_mean, c.r_dev
|
||||
with np.errstate(invalid="ignore", divide="ignore", over="ignore"):
|
||||
# --- gen_std_flair ------------------------------------------------------------------------
|
||||
z1 = (F32(0.5) * (tm[:, 1] - tm[:, 0]) / (td[:, 1] + td[:, 0])).astype(np.float32)
|
||||
z2 = (F32(0.5) * (tm[:, 2] - tm[:, 1]) / (td[:, 2] + td[:, 1])).astype(np.float32)
|
||||
mz = np.fmin(z1, z2)
|
||||
third = F32(1.0) / F32(3.0)
|
||||
alpha = np.fmax(third * mz, F32(1.0))
|
||||
beta = np.fmin(third * mz, F32(1.0))
|
||||
eta = (alpha * beta).astype(np.float32)
|
||||
ftd = (td * eta[:, None]).astype(np.float32)
|
||||
ftd = np.where(ftd < F32(0.02), F32(0.02), ftd).astype(np.float32)
|
||||
med = np.fmax(np.fmin(td[:, 0], td[:, 1]), np.fmin(np.fmax(td[:, 0], td[:, 1]), td[:, 2]))
|
||||
min_r = np.where(F32(0.2) < med, med, F32(0.2)).astype(np.float32)
|
||||
min_rk = (rm / F32(7.0)).astype(np.float32)
|
||||
min_rk = np.where(min_rk < min_r[:, None], min_r[:, None], min_rk).astype(np.float32)
|
||||
frd = np.where(rd < min_rk, min_rk, rd).astype(np.float32)
|
||||
|
||||
# --- modilik: six ordered cases, later ones overwrite -----------------------------------------
|
||||
d = [
|
||||
100.0 * np.abs(t - tm[:, k]).astype(np.float64) + np.abs(rr - rm[:, k]).astype(np.float64) for k in range(3)
|
||||
]
|
||||
daa, dab, dbb = d
|
||||
ok_t = ~np.isnan(t)
|
||||
|
||||
pt = [_pdf(t, tm[:, k], ftd[:, k]) for k in range(3)]
|
||||
pr = [_pdf(rr, rm[:, k], frd[:, k]) for k in range(3)]
|
||||
|
||||
def pdf2(k: int, st: str, sr: str) -> F64Arr:
|
||||
return pt[k].get(st) * pr[k].get(sr)
|
||||
|
||||
L = np.zeros((n, 3), dtype=np.float64)
|
||||
covered = np.zeros(n, dtype=np.bool_)
|
||||
cases = [
|
||||
((daa <= dbb) & (dab <= dbb) & (rm[:, 0] <= rm[:, 1]), ((0, "left", "left"), (1, "right", "right"))),
|
||||
((daa <= dab) & (dbb <= dab) & (rm[:, 0] <= rm[:, 2]), ((0, "left", "left"), (2, "right", "right"))),
|
||||
((dab <= daa) & (dbb <= daa) & (rm[:, 1] <= rm[:, 2]), ((1, "left", "left"), (2, "right", "right"))),
|
||||
((daa <= dbb) & (dab <= dbb) & (rm[:, 0] > rm[:, 1]), ((0, "left", "right"), (1, "right", "left"))),
|
||||
((dab <= daa) & (dbb <= daa) & (rm[:, 1] > rm[:, 2]), ((1, "left", "right"), (2, "right", "left"))),
|
||||
((daa <= dab) & (dbb <= dab) & (rm[:, 0] > rm[:, 2]), ((0, "left", "right"), (2, "right", "left"))),
|
||||
]
|
||||
for cond, assigns in cases:
|
||||
m = cond & ok_t
|
||||
if m.any():
|
||||
for k, st, sr in assigns:
|
||||
L[m, k] = pdf2(k, st, sr)[m]
|
||||
covered |= m
|
||||
nc = ~covered
|
||||
if nc.any():
|
||||
for k in range(3):
|
||||
L[nc, k] = pdf2(k, "full", "full")[nc]
|
||||
L = np.where(np.isnan(L), 0.0, L)
|
||||
L[~ok_t] = np.nan
|
||||
|
||||
# --- priors and posteriors -----------------------------------------------------------------
|
||||
ntot = c.n.sum(axis=1)
|
||||
P = (
|
||||
(c.n.astype(np.float32) / ntot[:, None].astype(np.float32)).astype(np.float32).astype(np.float64) + 1.0
|
||||
) / 4.0
|
||||
LP = L * P
|
||||
ev = (LP[:, 0] + LP[:, 1]) + LP[:, 2]
|
||||
ap = LP / ev[:, None]
|
||||
aa, ab, bb = ap[:, 0], ap[:, 1], ap[:, 2]
|
||||
iap = np.where(
|
||||
(aa >= ab) & (aa >= bb), 1, np.where((ab >= aa) & (ab >= bb), 2, np.where((bb >= aa) & (bb >= ab), 3, 0))
|
||||
).astype(np.uint8)
|
||||
mx = np.where(aa > ab, aa, ab)
|
||||
scnd = np.where(aa > ab, ab, aa)
|
||||
b_gt_mx = bb > mx
|
||||
scnd = np.where(b_gt_mx, mx, np.where(bb > scnd, bb, scnd))
|
||||
mx = np.where(b_gt_mx, bb, mx)
|
||||
ap_lod = _libm(math.log10, mx / (TINY_D + scnd))
|
||||
score_ap = smf(ap_lod.astype(np.float32), F32(0.0), F32(2.0)).astype(np.float32)
|
||||
|
||||
# --- sub-scores --------------------------------------------------------------------------
|
||||
score_r1 = smf(rr, F32(0.0), F32(0.1))
|
||||
nclus = (c.n > 0).sum(axis=1)
|
||||
misclust = np.where(
|
||||
(nclus == 1) & (c.n[:, 1] == 0), F32(0.7), np.where(nclus != 3, F32(0.95), F32(1.0))
|
||||
).astype(np.float32)
|
||||
g = np.clip(iap.astype(np.int64) - 1, 0, 2)
|
||||
gm_t, gd_t = tm[np.arange(n), g], ftd[np.arange(n), g]
|
||||
gm_r, gd_r = rm[np.arange(n), g], frd[np.arange(n), g]
|
||||
st_aa = zmf(t, gm_t + F32(2.0) * gd_t, gm_t + F32(6.0) * gd_t)
|
||||
st_ab = zmf(np.abs((t - gm_t) / gd_t).astype(np.float32), F32(2.0), F32(6.0))
|
||||
st_bb = smf(t, gm_t - F32(6.0) * gd_t, gm_t - F32(2.0) * gd_t)
|
||||
score_t = np.where(iap == 1, st_aa, np.where(iap == 2, st_ab, np.where(iap == 3, st_bb, F32(1.0))))
|
||||
r2 = smf(rr, F32(0.0), (gm_r / F32(10.0)).astype(np.float32))
|
||||
r3 = smf(rr, gm_r - F32(6.0) * gd_r, gm_r - F32(2.0) * gd_r)
|
||||
r4 = zmf(rr, gm_r + F32(6.0) * gd_r, gm_r + F32(12.0) * gd_r)
|
||||
called = iap > 0
|
||||
score_r = np.where(called, ((score_r1 * r2) * r3) * r4, F32(0.0)).astype(np.float32)
|
||||
score_t = score_t.astype(np.float32)
|
||||
score_n = np.ones(n, dtype=np.float32)
|
||||
nan_gt = np.isnan(c.gentrain)
|
||||
score_t[nan_gt] = np.nan
|
||||
score_r[nan_gt] = np.nan
|
||||
score_t[~ok_t] = np.nan
|
||||
score_r[~ok_t] = np.nan
|
||||
score_n[~ok_t] = np.nan
|
||||
prelim = (((((score_ap * c.gentrain) * score_t) * score_r) * score_n) * misclust).astype(np.float32)
|
||||
|
||||
# intensity threshold: below it the call is NC with a NaN score
|
||||
below = rr < c.intensity_threshold
|
||||
prelim[below] = np.nan
|
||||
iap[below] = 0
|
||||
|
||||
# --- gencall_score_map -------------------------------------------------------------------
|
||||
x = prelim.astype(np.float64)
|
||||
q = (x - 1.08) / 0.504
|
||||
gbell = 1.0 / (1.0 + _libm(_pow(1.71), q * q))
|
||||
score = (_libm(_pow(0.35), x) * gbell).astype(np.float32)
|
||||
score = np.where(np.isnan(score), F32(0.0), score).astype(np.float32)
|
||||
|
||||
# loci without intensity or without a cluster record
|
||||
no_call = raw_zero | ~c.present
|
||||
score[no_call] = F32(0.0)
|
||||
iap[no_call] = 0
|
||||
geno = np.where(score < F32(cutoff), 0, iap).astype(np.uint8)
|
||||
geno[intensity_only] = 0
|
||||
return Calls(
|
||||
genotypes=geno, scores=score, num_calls=int((geno > 0).sum()), num_intensity_only=int(intensity_only.sum())
|
||||
)
|
||||
|
||||
|
||||
def call_rate(calls: Calls, num_snps: int) -> np.float32:
|
||||
tiny = F32(np.finfo(np.float32).tiny) * F32(np.finfo(np.float32).eps)
|
||||
return F32(F32(calls.num_calls) / ((F32(num_snps) - F32(calls.num_intensity_only)) + tiny))
|
||||
@@ -0,0 +1,196 @@
|
||||
"""One sample from IDAT pair to a complete GTC record — plan §7.1–§7.8, spec OUT-001/OUT-002.
|
||||
|
||||
Order of operations follows the reference AutoCall pipeline (gtc2vcf ``gtc_init``, MIT): raw X/Y → normalization →
|
||||
calls → sex → BAF/LRR/LogR dev → GC10/GC50 → intensity percentiles → controls → sample metadata.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import datetime as _dt
|
||||
import re
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg import PRODUCT_NAME, __version__
|
||||
from simurg.formats.bpm import Manifest
|
||||
from simurg.formats.gtc import XFORM_DTYPE, Gtc, ScannerData
|
||||
from simurg.formats.idat import Idat
|
||||
from simurg.genotype.baf_lrr import autosomal_mask, baf_lrr, logr_dev
|
||||
from simurg.genotype.clusters import ClusterArrays
|
||||
from simurg.genotype.gencall import DEFAULT_CUTOFF, call_genotypes, call_rate
|
||||
from simurg.genotype.sex import SexParams, chromosome_classes, estimate_sex
|
||||
from simurg.signal.normalization.apply import normalize, polar
|
||||
from simurg.signal.normalization.fit import XForm, fit_transforms
|
||||
from simurg.signal.normalization.numerics import percentile
|
||||
from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy
|
||||
|
||||
F32 = np.float32
|
||||
_COMPLEMENT = bytes.maketrans(b"ACGTDI", b"TGCADI")
|
||||
|
||||
|
||||
@dataclass
|
||||
class GenotypeParams:
|
||||
cutoff: float = DEFAULT_CUTOFF
|
||||
norm_version: str = "1.2.0"
|
||||
sex: SexParams = field(default_factory=SexParams)
|
||||
gtc_version: int = 5
|
||||
autocall_version: str = f"{PRODUCT_NAME} {__version__}"
|
||||
autocall_date: str | None = None # None: run start time (invariant format, OUT-001)
|
||||
|
||||
|
||||
@dataclass
|
||||
class ManifestContext:
|
||||
"""Per-manifest data reused for every sample (computed once)."""
|
||||
|
||||
manifest: Manifest
|
||||
clusters: ClusterArrays
|
||||
norm_bins: npt.NDArray[np.uint8]
|
||||
lookup: npt.NDArray[np.int32]
|
||||
allele_a: npt.NDArray[np.uint8]
|
||||
allele_b: npt.NDArray[np.uint8]
|
||||
autosomal: npt.NDArray[np.bool_]
|
||||
sex_classes: tuple[npt.NDArray[np.bool_], ...]
|
||||
cache: LocusMapCache
|
||||
cluster_file_name: str = ""
|
||||
|
||||
|
||||
def top_alleles(manifest: Manifest) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint8]]:
|
||||
"""TOP-strand A/B base per locus (gtc2vcf ``get_base_call``): SNP column on TOP designs, complemented otherwise.
|
||||
|
||||
Unknown bases complement to NUL, as in the reference's lookup table.
|
||||
"""
|
||||
a = np.zeros(manifest.num_loci, dtype=np.uint8)
|
||||
b = np.zeros(manifest.num_loci, dtype=np.uint8)
|
||||
for i, (snp, strand) in enumerate(zip(manifest.snps, manifest.ilmn_strands, strict=True)):
|
||||
raw = snp.encode("latin-1", errors="replace")
|
||||
ca = raw[1:2] or b"\x00"
|
||||
cb = raw[3:4] or b"\x00"
|
||||
if strand[:1].upper() != "T":
|
||||
ca = ca.translate(_COMPLEMENT) if ca in b"ACGTDI" else b"\x00"
|
||||
cb = cb.translate(_COMPLEMENT) if cb in b"ACGTDI" else b"\x00"
|
||||
a[i], b[i] = ca[0], cb[0]
|
||||
return a, b
|
||||
|
||||
|
||||
def manifest_context(manifest: Manifest, clusters: ClusterArrays, cluster_file_name: str = "") -> ManifestContext:
|
||||
a, b = top_alleles(manifest)
|
||||
return ManifestContext(
|
||||
manifest=manifest,
|
||||
clusters=clusters,
|
||||
norm_bins=manifest.norm_bin_ids(),
|
||||
lookup=manifest.normalization_lookups(),
|
||||
allele_a=a,
|
||||
allele_b=b,
|
||||
autosomal=autosomal_mask(manifest.chroms),
|
||||
sex_classes=chromosome_classes(manifest.chroms),
|
||||
cache=LocusMapCache(manifest),
|
||||
cluster_file_name=cluster_file_name,
|
||||
)
|
||||
|
||||
|
||||
def base_calls(
|
||||
genotypes: npt.NDArray[np.uint8], a: npt.NDArray[np.uint8], b: npt.NDArray[np.uint8]
|
||||
) -> npt.NDArray[np.uint8]:
|
||||
dash = np.uint8(ord("-"))
|
||||
first = np.where(genotypes == 1, a, np.where(genotypes == 2, a, np.where(genotypes == 3, b, dash)))
|
||||
second = np.where(genotypes == 1, a, np.where(genotypes == 2, b, np.where(genotypes == 3, b, dash)))
|
||||
return np.stack([first, second], axis=1).astype(np.uint8)
|
||||
|
||||
|
||||
def invariant_datetime(t: _dt.datetime) -> str:
|
||||
"""``M/d/yyyy h:mm:ss tt`` with ASCII digits (OUT-001), independent of the process locale."""
|
||||
hour = t.hour % 12 or 12
|
||||
return f"{t.month}/{t.day}/{t.year} {hour}:{t.minute:02d}:{t.second:02d} {'AM' if t.hour < 12 else 'PM'}"
|
||||
|
||||
|
||||
def _param(pars: str, key: str) -> str:
|
||||
m = re.search(rf"(?:^|\|){re.escape(key)}=([^|]*)", pars)
|
||||
return m.group(1) if m else ""
|
||||
|
||||
|
||||
def xforms_to_records(xforms: list[XForm]) -> npt.NDArray[np.void]:
|
||||
rec = np.zeros(len(xforms), dtype=XFORM_DTYPE)
|
||||
for i, x in enumerate(xforms):
|
||||
rec[i] = (
|
||||
x.version,
|
||||
x.offset_x,
|
||||
x.offset_y,
|
||||
x.scale_x,
|
||||
x.scale_y,
|
||||
x.shear,
|
||||
x.theta,
|
||||
(x.cvx, x.cvy, x.nn12, x.rr12, x.taa, x.tbb),
|
||||
)
|
||||
return rec
|
||||
|
||||
|
||||
def genotype_sample(
|
||||
grn: Idat, red: Idat, ctx: ManifestContext, params: GenotypeParams, run_started: _dt.datetime | None = None
|
||||
) -> Gtc:
|
||||
m, c = ctx.manifest, ctx.clusters
|
||||
raw = assemble_raw_xy(grn, red, m, cache=ctx.cache)
|
||||
xforms = fit_transforms(raw.raw_x, raw.raw_y, ctx.norm_bins, params.norm_version)
|
||||
nx, ny = normalize(raw.raw_x, raw.raw_y, ctx.lookup, xforms)
|
||||
theta, r = polar(nx, ny)
|
||||
raw_zero = (raw.raw_x == 0) & (raw.raw_y == 0)
|
||||
theta_c = np.where(raw_zero, np.float32(np.nan), theta).astype(np.float32)
|
||||
r_c = np.where(raw_zero, np.float32(np.nan), r).astype(np.float32)
|
||||
|
||||
calls = call_genotypes(theta_c, r_c, raw_zero, c, m.intensity_only, params.cutoff)
|
||||
sex = estimate_sex(calls.genotypes, r, c.gentrain, c.present, ctx.sex_classes, params.sex)
|
||||
baf, lrr = baf_lrr(theta_c, r_c, raw_zero, c, m.intensity_only)
|
||||
ldev = logr_dev(lrr, ctx.autosomal, m.intensity_only, c.present)
|
||||
|
||||
gs = np.sort(calls.scores[(calls.scores > F32(params.cutoff)) & ~m.intensity_only])
|
||||
pct_x = [int(percentile(np.sort(raw.raw_x.astype(np.float32)), p)) for p in (5, 50, 95)]
|
||||
pct_y = [int(percentile(np.sort(raw.raw_y.astype(np.float32)), p)) for p in (5, 50, 95)]
|
||||
|
||||
scan = red.scan_entry()
|
||||
pars = scan.block_pars if scan else ""
|
||||
started = run_started or _dt.datetime.now()
|
||||
g = Gtc(
|
||||
version=params.gtc_version,
|
||||
num_snps=m.num_loci,
|
||||
ploidy=2,
|
||||
ploidy_type=1,
|
||||
sample_name=red.sample_id or Path(grn.path).name.split("_Grn.idat")[0].split("_grn.idat")[0],
|
||||
sample_plate=red.plate,
|
||||
sample_well=red.well,
|
||||
cluster_file=ctx.cluster_file_name,
|
||||
snp_manifest=Path(m.path).name,
|
||||
imaging_date=scan.run_time if scan else "",
|
||||
autocall_date=params.autocall_date if params.autocall_date is not None else invariant_datetime(started),
|
||||
autocall_version=params.autocall_version,
|
||||
normalization_transforms=xforms_to_records(xforms),
|
||||
controls_x=raw.controls_x,
|
||||
controls_y=raw.controls_y,
|
||||
raw_x=raw.raw_x,
|
||||
raw_y=raw.raw_y,
|
||||
genotypes=calls.genotypes,
|
||||
base_calls=base_calls(calls.genotypes, ctx.allele_a, ctx.allele_b),
|
||||
genotype_scores=calls.scores,
|
||||
scanner_data=ScannerData(
|
||||
name=_param(pars, "sherlockID"),
|
||||
pmt_green=int(_param(pars, "PMTGainCY3") or 0).to_bytes(4, "little", signed=True),
|
||||
pmt_red=int(_param(pars, "PMTGainCY5") or 0).to_bytes(4, "little", signed=True),
|
||||
version=scan.code_version if scan else "",
|
||||
user=_param(pars, "Username"),
|
||||
),
|
||||
call_rate=float(call_rate(calls, m.num_loci)),
|
||||
gender=sex.encode(),
|
||||
logr_dev=float(ldev),
|
||||
gc10=float(percentile(gs, 10)),
|
||||
b_allele_freqs=baf,
|
||||
logr_ratios=lrr,
|
||||
percentiles_x=np.array(pct_x, dtype=np.uint16),
|
||||
percentiles_y=np.array(pct_y, dtype=np.uint16),
|
||||
slide_identifier=red.barcode,
|
||||
)
|
||||
g.sample_data["p50gc"] = percentile(gs, 50)
|
||||
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
|
||||
return g
|
||||
@@ -0,0 +1,88 @@
|
||||
"""Sex estimate from X heterozygosity and X/Y intensities — plan §7.8.
|
||||
|
||||
Port of gtc2vcf ``estimate_gender`` (MIT), which reproduces the vendor's two generations:
|
||||
|
||||
* version 1 (AutoConvert 1.x): male unless the X heterozygosity rate among called X loci exceeds 0.1;
|
||||
* version 2 (AutoConvert 2.0, IAAP, current CLI): median normalized R of Y loci > 0.3 → M; otherwise median R of
|
||||
X loci > 0.9 → F (or F when there are too few X loci), else U. Falls back to version 1 when there are too few Y or
|
||||
autosomal loci.
|
||||
|
||||
Only loci whose cluster record has a non-zero GenTrain score count.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.genotype.baf_lrr import chrom_core
|
||||
from simurg.signal.normalization.numerics import median_sorted
|
||||
|
||||
F32 = np.float32
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class SexParams:
|
||||
version: int = 2
|
||||
min_loci: int = 100
|
||||
min_x_loci: int = 20
|
||||
min_y_loci: int = 20
|
||||
call_rate_threshold: float = 0.0
|
||||
y_threshold: float = 0.3
|
||||
x_threshold: float = 0.9
|
||||
x_het_rate_threshold: float = 0.1
|
||||
|
||||
|
||||
def chromosome_classes(chroms: list[str]) -> tuple[npt.NDArray[np.bool_], ...]:
|
||||
core = [chrom_core(c) for c in chroms]
|
||||
is_x = np.array([c == "X" for c in core])
|
||||
is_y = np.array([c == "Y" for c in core])
|
||||
is_auto = np.array([c not in ("X", "Y", "XY", "MT") for c in core])
|
||||
return is_x, is_y, is_auto
|
||||
|
||||
|
||||
def estimate_sex(
|
||||
genotypes: npt.NDArray[np.uint8],
|
||||
r_all: npt.NDArray[np.float32],
|
||||
gentrain: npt.NDArray[np.float32],
|
||||
present: npt.NDArray[np.bool_],
|
||||
classes: tuple[npt.NDArray[np.bool_], ...],
|
||||
p: SexParams | None = None,
|
||||
) -> str:
|
||||
"""``r_all``: normalized R of every locus, computed even where raw intensities are zero (as the reference)."""
|
||||
p = p or SexParams()
|
||||
is_x, is_y, is_auto = classes
|
||||
|
||||
def thr(v: float) -> float: # thresholds are C floats promoted to double in the comparisons
|
||||
return float(F32(v))
|
||||
|
||||
use = present & (gentrain != F32(0.0)) # NaN != 0 is True, as in C
|
||||
x = use & is_x
|
||||
y = use & is_y
|
||||
auto = use & is_auto & ~is_x & ~is_y
|
||||
x_count, y_count, auto_count = int(x.sum()), int(y.sum()), int(auto.sum())
|
||||
called = genotypes > 0
|
||||
x_non_missing = int((x & called).sum())
|
||||
x_hets = int((x & (genotypes == 2)).sum())
|
||||
auto_non_missing = int((auto & called).sum())
|
||||
|
||||
if p.version == 1 or y_count < p.min_y_loci or auto_count < p.min_loci:
|
||||
if x_non_missing > p.min_x_loci:
|
||||
rate = float(F32(F32(x_hets) / F32(x_non_missing)))
|
||||
return "M" if not rate > thr(p.x_het_rate_threshold) else "F"
|
||||
return "U"
|
||||
if auto_count > 0 and auto_non_missing / auto_count > thr(p.call_rate_threshold):
|
||||
ry = r_all[y].astype(np.float32)
|
||||
ry = np.where(np.isfinite(ry), ry, F32(0.0)).astype(np.float32)
|
||||
y_med = median_sorted(np.sort(ry))
|
||||
if float(y_med) > thr(p.y_threshold):
|
||||
return "M"
|
||||
if x_count < p.min_x_loci:
|
||||
return "F"
|
||||
rx = r_all[x].astype(np.float32)
|
||||
rx = np.where(np.isfinite(rx), rx, F32(0.0)).astype(np.float32)
|
||||
x_med = median_sorted(np.sort(rx))
|
||||
return "F" if float(x_med) > thr(p.x_threshold) else "U"
|
||||
return "U"
|
||||
Reference in new issue
Block a user