P6: gtc-to-vcf (VCF/BCF + indexes), harmonisation, exports (plink, bedgraph, finalreport)
- formats/fasta: plain/BGZF FASTA, .fai/.gzi built identical to htslib, index check, build detection, contig aliases, mmap genome cache (IN-060..064) - formats/bgzf, formats/bcf: pure-Python BGZF, tabix .tbi, BCF 2.2 and CSI (OUT-010, OUT-038, OUT-040) - harmonize/sites: plus-strand mapping with audit rows, indel resolution from SourceSeq (port of GTCtoVCF, Apache-2.0) and left normalisation (OUT-014..021, OUT-030, OUT-042) - harmonize/combine: OUT-031 genotype combination (ploidy-aware, Infinium I/II signal classes, OUT-037) - pipeline/vcf: `simurg gtc-to-vcf` (per-sample VCF.gz+tbi, --merge cohort BCF+csi, audit parquet) - pipeline/export: `simurg export bedgraph|finalreport|plink|locus-summary` (OUT-043, OUT-050..052, OUT-071) - validation: 36 GSA samples vs gtc2vcf: 23.66 M identical calls, 0 unexplained differences; PLINK identical to plink2 on 637,397 variants; bcftools validates VCF/BCF and our indexes answer region queries correctly - fetch_data.sh p6 (Ensembl GRCh38 with checksum), make_oracle_vcf.sh, validate_p6_vcf.py; K-022/K-023
This commit is contained in:
kemal committed
1 parent
1b6f834583
commit
e65d2426bb
19 files changed
+3184
-9
No files matched your search
@@ -0,0 +1,121 @@
|
||||
"""Genotypes of VCF sites from probe calls (spec OUT-018, OUT-031, OUT-032, OUT-037; plan §8.3).
|
||||
|
||||
Single-probe sites map the A/B call straight to allele indices. Sites with several probes use the exact candidate
|
||||
algorithm of OUT-031: enumerate every genotype over {REF} ∪ ALT with the site's ploidy, predict each probe's call
|
||||
from the alleles' signal classes (Infinium II: dye group of the extended base; Infinium I and indels: exact allele,
|
||||
no signal otherwise) and keep the candidates consistent with every called probe.
|
||||
[OSS: Illumina GTCtoVCF ``GenotypeFormat.RecordCombiner`` (Apache-2.0) — the consistency idea; the class model and
|
||||
the ploidy-aware enumeration are ours.]
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from itertools import combinations_with_replacement
|
||||
|
||||
from simurg.harmonize.sites import GREEN, RED, ProbeSite, Site
|
||||
|
||||
NC, AA, AB, BB = 0, 1, 2, 3
|
||||
CALL_TEXT = {NC: "NC", AA: "AA", AB: "AB", BB: "BB"}
|
||||
|
||||
|
||||
def signal_class(p: ProbeSite, allele: str) -> str | None:
|
||||
"""OUT-031 step 1."""
|
||||
a, b = p.alleles
|
||||
if p.kind == "snv" and p.infinium == 2 and len(allele) == 1:
|
||||
grp = RED if allele in RED else GREEN
|
||||
if a in grp:
|
||||
return "A"
|
||||
if b in grp:
|
||||
return "B"
|
||||
return None
|
||||
if allele == a:
|
||||
return "A"
|
||||
if allele == b:
|
||||
return "B"
|
||||
return None
|
||||
|
||||
|
||||
def predicted_call(p: ProbeSite, genotype: tuple[str, ...]) -> int:
|
||||
classes = sorted(c for c in (signal_class(p, z) for z in genotype) if c is not None)
|
||||
if not classes:
|
||||
return NC
|
||||
u = set(classes)
|
||||
if u == {"A"}:
|
||||
return AA
|
||||
if u == {"B"}:
|
||||
return BB
|
||||
return AB
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Outcome:
|
||||
alleles: tuple[int, ...] | None # allele indices (sorted) or None for missing
|
||||
reason: str # "" | "het_in_haploid" | W_MERGE_* | W_INF1_NC_CONFLICT
|
||||
|
||||
|
||||
def single(site: Site, call: int, ploidy: int) -> Outcome:
|
||||
p = site.probes[0]
|
||||
if call == NC:
|
||||
return Outcome(None, "")
|
||||
ia, ib = site.allele_index(p.alleles[0]), site.allele_index(p.alleles[1])
|
||||
if ploidy == 1:
|
||||
if call == AB:
|
||||
return Outcome(None, "het_in_haploid")
|
||||
return Outcome((ia if call == AA else ib,), "")
|
||||
pair = {AA: (ia, ia), AB: (ia, ib), BB: (ib, ib)}[call]
|
||||
return Outcome(tuple(sorted(pair)), "")
|
||||
|
||||
|
||||
def merged(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome:
|
||||
"""OUT-031 steps 3–6."""
|
||||
if all(c == NC for c in calls):
|
||||
return Outcome(None, "W_MERGE_ALL_NC")
|
||||
alleles = [site.ref, *site.alts]
|
||||
cands = list(combinations_with_replacement(range(len(alleles)), ploidy))
|
||||
keep = [
|
||||
g
|
||||
for g in cands
|
||||
if all(
|
||||
c == NC or predicted_call(p, tuple(alleles[i] for i in g)) == c
|
||||
for p, c in zip(site.probes, calls, strict=True)
|
||||
)
|
||||
]
|
||||
reason = ""
|
||||
if inf1_nc_informative:
|
||||
for p, c in zip(site.probes, calls, strict=True):
|
||||
if c == NC and p.infinium == 1:
|
||||
sub = [g for g in keep if all(signal_class(p, alleles[i]) is None for i in g)]
|
||||
if sub:
|
||||
keep = sub
|
||||
else:
|
||||
reason = "W_INF1_NC_CONFLICT"
|
||||
if len(keep) == 1:
|
||||
return Outcome(keep[0], reason)
|
||||
return Outcome(None, "W_MERGE_INCONSISTENT" if not keep else "W_MERGE_AMBIGUOUS")
|
||||
|
||||
|
||||
def genotype(site: Site, calls: tuple[int, ...], ploidy: int, inf1_nc_informative: bool = False) -> Outcome:
|
||||
if len(site.probes) == 1:
|
||||
return single(site, calls[0], ploidy)
|
||||
key = (calls, ploidy, inf1_nc_informative)
|
||||
hit = site.memo.get(key)
|
||||
if hit is None:
|
||||
hit = site.memo[key] = merged(site, calls, ploidy, inf1_nc_informative)
|
||||
assert isinstance(hit, Outcome)
|
||||
return hit
|
||||
|
||||
|
||||
def gt_text(o: Outcome, ploidy: int) -> str:
|
||||
if o.alleles is None:
|
||||
return "." if ploidy == 1 else "./."
|
||||
return "/".join(str(i) for i in o.alleles)
|
||||
|
||||
|
||||
def failure_message(site: Site, calls: tuple[int, ...], reason: str) -> str:
|
||||
what = {"W_MERGE_INCONSISTENT": "inconsistent", "W_MERGE_AMBIGUOUS": "ambiguous", "W_MERGE_ALL_NC": "all no-calls"}
|
||||
probes = ", ".join(
|
||||
f"{p.name} ({'InfiniumI' if p.infinium == 1 else 'InfiniumII'}): {CALL_TEXT[c]}"
|
||||
for p, c in zip(site.probes, calls, strict=True)
|
||||
)
|
||||
return f"Failed to combine genotypes ({what.get(reason, reason)}) - {probes}"
|
||||
@@ -0,0 +1,307 @@
|
||||
"""Manifest loci → plus-strand VCF sites, with an audit row per locus (spec IN-046, OUT-014 … OUT-021, OUT-030,
|
||||
OUT-042; plan §8.1–§8.2).
|
||||
|
||||
[DOC: Illumina "TOP/BOT strand and A/B allele" technical note — strand conventions]
|
||||
[OSS: Illumina GTCtoVCF ``BPMRecord.py`` (Apache-2.0) — indel source-sequence handling: design-strand orientation,
|
||||
left shift, IUPAC-aware flank matching and the deletion-versus-insertion context score; ported below.]
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import re
|
||||
from collections.abc import Callable, Sequence
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
COMPLEMENT = dict(zip("ABCDGHKMRTVYNID", "TVGHCDMKYABRNID", strict=True))
|
||||
DEGENERACY = {"A": "A", "C": "C", "G": "G", "T": "T", "R": "AG", "Y": "CT", "S": "GC", "W": "AT", "K": "GT",
|
||||
"M": "AC", "B": "CGT", "D": "AGT", "H": "ACT", "V": "ACG", "N": "ACGT"} # fmt: skip
|
||||
RED, GREEN = set("AT"), set("CG")
|
||||
SEQ_RE = re.compile(r"^([^\[]*)\[([^/\]]*)/([^\]]*)\](.*)$")
|
||||
|
||||
Fetch = Callable[[str, int, int], str] # (contig, start0, end0) -> upper-case bases
|
||||
|
||||
|
||||
def revcomp(s: str) -> str:
|
||||
return "".join(COMPLEMENT.get(c, "N") for c in reversed(s))
|
||||
|
||||
|
||||
def complement_base(b: str) -> str:
|
||||
return COMPLEMENT.get(b, "N")
|
||||
|
||||
|
||||
@dataclass
|
||||
class ProbeSite:
|
||||
"""One manifest locus mapped to the plus strand."""
|
||||
|
||||
index: int
|
||||
name: str
|
||||
contig: str
|
||||
pos: int # VCF POS (1-based)
|
||||
ref: str
|
||||
alleles: tuple[str, str] # VCF allele strings of the probe's A and B alleles
|
||||
kind: str # "snv" | "indel"
|
||||
infinium: int # 1 or 2
|
||||
manifest_chrom: str
|
||||
par: bool = False
|
||||
confidence: str = "high"
|
||||
flipped: bool = False # design strand is minus
|
||||
reason: str = "ok"
|
||||
|
||||
|
||||
@dataclass
|
||||
class AuditRow:
|
||||
locus: str
|
||||
original_alleles: str
|
||||
strand_in: str
|
||||
transformed_alleles: str
|
||||
strand_out: str
|
||||
ref: str
|
||||
alt: str
|
||||
flipped: bool
|
||||
swapped: bool
|
||||
reason: str
|
||||
build: str
|
||||
confidence: str
|
||||
|
||||
|
||||
@dataclass
|
||||
class Site:
|
||||
"""A VCF record: one or more probes that share CHROM/POS/REF (and the indel span)."""
|
||||
|
||||
contig: str
|
||||
pos: int
|
||||
ref: str
|
||||
alts: list[str]
|
||||
probes: list[ProbeSite] = field(default_factory=list)
|
||||
kind: str = "snv"
|
||||
memo: dict[object, object] = field(default_factory=dict, repr=False, compare=False) # merged-genotype cache
|
||||
|
||||
@property
|
||||
def id(self) -> str:
|
||||
return (";" if self.kind == "indel" else ",").join(p.name for p in self.probes)
|
||||
|
||||
def allele_index(self, allele: str) -> int:
|
||||
return 0 if allele == self.ref else 1 + self.alts.index(allele)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ indel helpers
|
||||
def split_source(seq: str) -> tuple[str, str, str, str] | None:
|
||||
"""``left[a/b]right`` → (left, a, b, right); None when malformed."""
|
||||
m = SEQ_RE.match(seq.strip().upper())
|
||||
return (m.group(1), m.group(2), m.group(3), m.group(4)) if m else None
|
||||
|
||||
|
||||
def left_shift(five: str, indel: str, three: str) -> tuple[str, str]:
|
||||
"""Shift an indel fully 5′ within its context (GTCtoVCF ``determine_left_shift``)."""
|
||||
while indel and five.endswith(indel):
|
||||
five, three = five[: -len(indel)], indel + three
|
||||
while indel and five and len(indel) * five[-1] == indel:
|
||||
three, five = five[-1] + three, five[:-1]
|
||||
return five, three
|
||||
|
||||
|
||||
def _match(a: str, b: str) -> bool:
|
||||
return a in DEGENERACY.get(b, "")
|
||||
|
||||
|
||||
def max_suffix_match(genome: str, design: str) -> int:
|
||||
n = 0
|
||||
for x, y in zip(genome[::-1], design[::-1], strict=False):
|
||||
if not _match(x, y):
|
||||
break
|
||||
n += 1
|
||||
return n
|
||||
|
||||
|
||||
def max_prefix_match(genome: str, design: str) -> int:
|
||||
n = 0
|
||||
for x, y in zip(genome, design, strict=False):
|
||||
if not _match(x, y):
|
||||
break
|
||||
n += 1
|
||||
return n
|
||||
|
||||
|
||||
def classify_indel(fetch: Fetch, contig: str, pos: int, five: str, indel: str, three: str) -> tuple[bool | None, float]:
|
||||
"""Is the plus-strand ``indel`` sequence present in the reference at ``pos`` (a reference deletion)?
|
||||
Returns (is_deletion or None when undecidable, best context score). Context scoring as GTCtoVCF."""
|
||||
start = pos - 1
|
||||
L = len(indel)
|
||||
seq_match = fetch(contig, start, start + L) == indel
|
||||
d5, d3 = left_shift(fetch(contig, start - len(five), start), indel,
|
||||
fetch(contig, start + L, start + L + len(three))) # fmt: skip
|
||||
i5, i3 = left_shift(fetch(contig, start - len(five) + 1, start + 1), indel,
|
||||
fetch(contig, start + 1, start + len(three) + 1)) # fmt: skip
|
||||
dm = (max_suffix_match(d5, five), max_prefix_match(d3, three))
|
||||
dmax = min(len(d5), len(five)) + min(len(d3), len(three)) + L
|
||||
dscore = (sum(dm) + L if seq_match else 0) / dmax if dmax else 0.0
|
||||
im = (max_suffix_match(i5, five), max_prefix_match(i3, three))
|
||||
imax = min(len(i5), len(five)) + min(len(i3), len(three))
|
||||
iscore = sum(im) / imax if imax else 0.0
|
||||
is_del = seq_match and dscore > iscore and min(dm) >= 1
|
||||
is_ins = iscore > dscore and min(im) >= 1
|
||||
if is_del == is_ins:
|
||||
return None, max(dscore, iscore)
|
||||
return is_del, dscore if is_del else iscore
|
||||
|
||||
|
||||
def normalize_indel(fetch: Fetch, contig: str, pos: int, indel: str, is_deletion: bool) -> tuple[int, str, str]:
|
||||
"""Left-normalised, anchored VCF representation. ``pos`` is the manifest position: the first deleted base for
|
||||
a deletion, the base 5′ of the insertion for an insertion (IN-046)."""
|
||||
if is_deletion:
|
||||
p0 = pos - 1 # 0-based first deleted base
|
||||
while p0 > 0 and fetch(contig, p0 - 1, p0) == indel[-1]:
|
||||
indel = indel[-1] + indel[:-1]
|
||||
p0 -= 1
|
||||
anchor = fetch(contig, p0 - 1, p0)
|
||||
return p0, anchor + indel, anchor # VCF POS = p0 (1-based anchor)
|
||||
a0 = pos - 1 # 0-based anchor base
|
||||
while a0 > 0 and fetch(contig, a0, a0 + 1) == indel[-1]:
|
||||
indel = indel[-1] + indel[:-1]
|
||||
a0 -= 1
|
||||
anchor = fetch(contig, a0, a0 + 1)
|
||||
return a0 + 1, anchor, anchor + indel
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ strand
|
||||
def strand_from_flanks(fetch: Fetch, contig: str, pos: int, left: str, right: str) -> tuple[str, float]:
|
||||
"""'+' when the design flanks match the reference as written, '-' when reverse-complemented; with a score."""
|
||||
left, right = left[-30:], right[:30]
|
||||
g_left = fetch(contig, pos - 1 - len(left), pos - 1)
|
||||
g_right = fetch(contig, pos, pos + len(right))
|
||||
|
||||
def score(lf: str, rt: str) -> float:
|
||||
n = len(lf) + len(rt)
|
||||
if not n:
|
||||
return 0.0
|
||||
hits = sum(_match(a, b) for a, b in zip(g_left[::-1], lf[::-1], strict=False))
|
||||
hits += sum(_match(a, b) for a, b in zip(g_right, rt, strict=False))
|
||||
return hits / n
|
||||
|
||||
fwd, rev = score(left, right), score(revcomp(right), revcomp(left))
|
||||
return ("+", fwd) if fwd >= rev else ("-", rev)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ mapping
|
||||
@dataclass
|
||||
class ManifestView:
|
||||
"""The manifest columns the mapper needs (BPM and/or CSV)."""
|
||||
|
||||
names: Sequence[str]
|
||||
snps: Sequence[str]
|
||||
chroms: Sequence[str]
|
||||
positions: Sequence[int]
|
||||
ilmn_strands: Sequence[str]
|
||||
source_strands: Sequence[str]
|
||||
ref_strands: Sequence[str]
|
||||
source_seqs: Sequence[str]
|
||||
assay_types: Sequence[int] # 0 Infinium II, 1/2 Infinium I
|
||||
intensity_only: Sequence[bool]
|
||||
|
||||
|
||||
def map_locus(
|
||||
i: int, mv: ManifestView, fetch: Fetch, contig_of: Callable[[str], str | None], build: str = ""
|
||||
) -> tuple[ProbeSite | None, AuditRow]:
|
||||
name, snp, chrom = mv.names[i], mv.snps[i], mv.chroms[i]
|
||||
pos = int(mv.positions[i])
|
||||
ilmn, ref_strand = mv.ilmn_strands[i] or "", (mv.ref_strands[i] or "").strip()
|
||||
audit = AuditRow(name, snp, f"{ilmn}/{ref_strand or '?'}", "", "+", "", "", False, False, "ok", build, "high")
|
||||
|
||||
def skip(code: str) -> tuple[None, AuditRow]:
|
||||
audit.reason, audit.strand_out = code, ""
|
||||
return None, audit
|
||||
|
||||
if chrom in ("0", "") or pos <= 0:
|
||||
return skip("W_SKIP_UNMAPPED")
|
||||
if mv.intensity_only[i]:
|
||||
return skip("W_SKIP_INTENSITY_ONLY")
|
||||
contig = contig_of(chrom)
|
||||
if contig is None:
|
||||
return skip("W_SKIP_CONTIG")
|
||||
m = re.match(r"^\[([^/\]]+)/([^\]]+)\]$", snp.strip())
|
||||
if not m:
|
||||
return skip("W_SKIP_UNMAPPED")
|
||||
a, b = m.group(1).upper(), m.group(2).upper()
|
||||
infinium = 2 if int(mv.assay_types[i]) == 0 else 1
|
||||
par = chrom == "XY"
|
||||
src = mv.source_seqs[i] if i < len(mv.source_seqs) else ""
|
||||
parts = split_source(src) if src else None
|
||||
if a in ("I", "D") or b in ("I", "D"):
|
||||
# ---------------------------------------------------------------- indel
|
||||
if parts is None or "-" not in (parts[1], parts[2]):
|
||||
return skip("W_SKIP_INDEL_CONTEXT")
|
||||
five, x, y, three = parts
|
||||
indel_seq = y if x == "-" else x
|
||||
s_strand = (mv.source_strands[i] or "")[:1].upper()
|
||||
d_strand = ilmn[:1].upper()
|
||||
if not ref_strand or s_strand in ("U", "") or d_strand in ("U", ""):
|
||||
return skip("W_SKIP_INDEL_CONTEXT")
|
||||
on_design = s_strand == d_strand
|
||||
if on_design != (ref_strand == "+"):
|
||||
five, indel_seq, three = revcomp(three), revcomp(indel_seq), revcomp(five)
|
||||
five, three = left_shift(five, indel_seq, three)
|
||||
is_del, score = classify_indel(fetch, contig, pos, five, indel_seq, three)
|
||||
if is_del is None:
|
||||
return skip("W_SKIP_INDEL_REF")
|
||||
vpos, vref, valt = normalize_indel(fetch, contig, pos, indel_seq, is_del)
|
||||
if not vref or "N" in vref:
|
||||
return skip("W_SKIP_INDEL_REF")
|
||||
# I (insertion allele) carries the sequence; D lacks it
|
||||
allele_of = {"I": vref if is_del else valt, "D": valt if is_del else vref}
|
||||
alleles = (allele_of[a], allele_of[b])
|
||||
audit.reason = "ok" if score >= 1.0 else "W_INDEL_PARTIAL_MATCH"
|
||||
audit.transformed_alleles = f"{alleles[0]}/{alleles[1]}"
|
||||
audit.ref, audit.alt = vref, valt
|
||||
audit.swapped = alleles[0] != vref
|
||||
return ProbeSite(i, name, contig, vpos, vref, alleles, "indel", infinium, chrom, par,
|
||||
"high", ref_strand == "-", audit.reason), audit # fmt: skip
|
||||
# -------------------------------------------------------------------- SNV
|
||||
confidence = "high"
|
||||
strand = ref_strand
|
||||
if strand not in ("+", "-"):
|
||||
if parts is None:
|
||||
return skip("W_SKIP_NO_STRAND") # no RefStrand and no sequence to align
|
||||
five, x, y, three = parts
|
||||
src_plus, score = strand_from_flanks(fetch, contig, pos, five, three)
|
||||
same = (mv.source_strands[i] or "")[:1].upper() == ilmn[:1].upper()
|
||||
strand = src_plus if same else ("-" if src_plus == "+" else "+")
|
||||
if score < 0.9:
|
||||
confidence = "low"
|
||||
if {a, b} in ({"A", "T"}, {"C", "G"}):
|
||||
confidence = "low" # plan §8.2: ambiguous without RefStrand
|
||||
pa, pb = (a, b) if strand == "+" else (complement_base(a), complement_base(b))
|
||||
ref = fetch(contig, pos - 1, pos)
|
||||
if not ref or ref not in "ACGT":
|
||||
return skip("W_SKIP_UNMAPPED")
|
||||
audit.transformed_alleles = f"{pa}/{pb}"
|
||||
audit.ref = ref
|
||||
audit.alt = ",".join(x for x in (pa, pb) if x != ref)
|
||||
audit.flipped = strand == "-"
|
||||
audit.swapped = pa != ref
|
||||
audit.confidence = confidence
|
||||
if ref not in (pa, pb):
|
||||
audit.reason = "W_REF_NOT_QUERIED"
|
||||
return ProbeSite(i, name, contig, pos, ref, (pa, pb), "snv", infinium, chrom, par, confidence,
|
||||
strand == "-", audit.reason), audit # fmt: skip
|
||||
|
||||
|
||||
def group_sites(probes: Sequence[ProbeSite], contig_order: dict[str, int], *, unsquash: bool = False) -> list[Site]:
|
||||
"""OUT-030/OUT-033/OUT-034: probes sharing CHROM/POS/REF (and, for indels, the ALT) form one site; ALT alleles in
|
||||
order of first appearance over BPM order; records sorted by contig (.fai order), POS, REF, ALT, ID."""
|
||||
groups: dict[tuple[str, int, str, str, str], list[ProbeSite]] = {}
|
||||
for p in sorted(probes, key=lambda p: p.index):
|
||||
extra = ",".join(sorted(x for x in p.alleles if x != p.ref)) if p.kind == "indel" else ""
|
||||
groups.setdefault((p.contig, p.pos, p.ref, p.kind, extra), []).append(p)
|
||||
sites: list[Site] = []
|
||||
for (contig, pos, ref, kind, _), ps in groups.items():
|
||||
alts: list[str] = []
|
||||
for p in ps:
|
||||
for x in p.alleles:
|
||||
if x != ref and x not in alts:
|
||||
alts.append(x)
|
||||
if unsquash and len(ps) > 1 and len(alts) <= 1:
|
||||
sites += [Site(contig, pos, ref, [x for x in p.alleles if x != ref], [p], kind) for p in ps]
|
||||
else:
|
||||
sites.append(Site(contig, pos, ref, alts, ps, kind))
|
||||
sites.sort(key=lambda s: (contig_order.get(s.contig, 1 << 30), s.pos, s.ref, ",".join(s.alts), s.id))
|
||||
return sites
|
||||
Reference in new issue
Block a user