kemal/SimurgArray
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

P6 validated against GTCtoVCF too: zero unexplained genotype differences

Browse Source 
- GTCtoVCF (Apache-2.0, 36 GSA samples, merged sites): 23.66 M identical calls; every difference named
  (intensity-only probes merged by the oracle, its opposite-strand indel bug EDGE-012, normalisation-based
  grouping, OUT-018 ploidy); indel representations identical, including 144 gtc2vcf cannot resolve
- validate_p6_vcf.py: rule-based categories (--bpm for [I/D]/[D/I] detection)
- docs: validation P6 complete, README status
This commit is contained in:
kemal committed 2026-09-28 13:20:44 +01:00
1 parent e65d2426bb
commit d3add008f7
5 files changed
+2868 -18
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

scripts/validate_p6_vcf.py
+33 -8
View File
@@ -49,14 +49,27 @@ def main() -> int:
ap.add_argument("--bcftools", required=True)
ap.add_argument("--out", required=True, type=Path)
ap.add_argument("--name", default="oracle")
ap.add_argument("--bpm", type=Path, help="manifest, to recognise opposite-strand indel designs (EDGE-012)")
a = ap.parse_args()
samples = subprocess.run([a.bcftools, "query", "-l", str(a.oracle)], capture_output=True, text=True,
check=True).stdout.split() # fmt: skip
oracle: dict[str, tuple[str, str, str, str, list[str]]] = {}
oracle_group: dict[str, frozenset[str]] = {}
for f in query(a.bcftools, a.fasta, a.oracle):
for pid in f[0].replace(";", ",").split(","):
ids = f[0].replace(";", ",").split(",")
for pid in ids:
oracle[pid] = (f[1], f[2], f[3], f[4], f[5:])
oracle_group[pid] = frozenset(ids)
snp: dict[str, str] = {}
if a.bpm:
sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
from simurg.formats.bpm import read_bpm
m = read_bpm(a.bpm)
snp = dict(zip(m.names, m.snps, strict=True))
first = a.ours / f"{samples[0]}.snv.vcf.gz"
ours_all = {pid for f in query(a.bcftools, a.fasta, first) for pid in f[0].replace(";", ",").split(",")}
cats: Counter[str] = Counter()
examples: dict[str, list[Any]] = defaultdict(list)
ours_ids: set[str] = set()
@@ -71,6 +84,8 @@ def main() -> int:
merged = len(ids) > 1
ours_al = alleles_of(f[5], f[3], f[4])
haploid = "/" not in f[5] and f[5] != "."
group = frozenset(ids)
mixed_indel = {snp.get(x) for x in ids} >= {"[I/D]", "[D/I]"}
for pid in ids:
if si == 0:
ours_ids.add(pid)
@@ -88,17 +103,27 @@ def main() -> int:
if ours_al == theirs:
cats["match_called" if ours_al else "match_missing"] += 1
continue
if ours_al is None:
if merged:
key = "ours_missing_merged_site"
og = oracle_group[pid]
if og - ours_all:
key = "oracle_site_includes_probes_we_skip" # e.g. intensity-only CNV probes (OUT-019)
elif og != group and (merged or len(og) > 1):
key = "grouping_differs" # same normalised variant, different manifest positions (OUT-030)
elif ours_al is not None and theirs is None and mixed_indel:
key = "oracle_opposite_strand_indel_bug" # EDGE-012: [I/D] + [D/I] both het → 0/1
elif ours_al is None:
if theirs is not None and theirs[0] != theirs[1] and f[5] == ".":
key = "ours_missing_het_in_haploid" # OUT-018
elif f[1] in ("Y", "chrY") and f[5] == ".":
key = "ours_missing_haploid_Y"
elif theirs is not None and theirs[0] != theirs[1] and (f[5] == "." or haploid):
key = "ours_missing_het_in_haploid"
key = "ours_missing_haploid_Y" # OUT-018: Y of female samples
elif merged:
key = "ours_missing_merged_site"
else:
key = "ours_missing_UNEXPLAINED"
elif theirs is None:
key = "oracle_missing_ours_called_merged" if merged else "oracle_missing_ours_called"
if merged and haploid:
key = "oracle_diploid_candidates_in_haploid_context" # OUT-031 enumerates with the ploidy
else:
key = "oracle_missing_ours_called_merged" if merged else "oracle_missing_ours_called"
else:
key = "diff_merged_site" if merged else "diff_UNEXPLAINED"
cats[key] += 1