From d3add008f7fad660aa0ef08011c9d757bb6fadec Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 13:20:44 +0100 Subject: [PATCH] P6 validated against GTCtoVCF too: zero unexplained genotype differences - 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 --- README.md | 9 +- docs/validation/P6.md | 35 +- docs/validation/P6_gsa36_vs_gtc2vcf.json | 1160 +++++++++++++++ docs/validation/P6_gsa36_vs_gtctovcf.json | 1641 +++++++++++++++++++++ scripts/validate_p6_vcf.py | 41 +- 5 files changed, 2868 insertions(+), 18 deletions(-) create mode 100644 docs/validation/P6_gsa36_vs_gtc2vcf.json create mode 100644 docs/validation/P6_gsa36_vs_gtctovcf.json diff --git a/README.md b/README.md index 89ceb05..0419982 100644 --- a/README.md +++ b/README.md @@ -16,8 +16,13 @@ Phase **P5 (QC)** is complete: control-probe metrics, thresholds/profiles, two-a KING-robust relatedness, Mendelian errors, locus summary and a self-contained HTML report (`simurg qc report`). On 2141 1000 Genomes samples all 717 pedigree parent–offspring pairs are recovered and trio Mendelian error rates are ~0.03 % (`docs/validation/P5.md`). +Phase **P6 (outputs)** is complete: `simurg gtc-to-vcf` (plus-strand harmonisation with audit, merged multi-probe +sites, per-sample VCF.gz + tabix or a cohort BCF + CSI) and `simurg export plink|bedgraph|finalreport|locus-summary`. +On the GSA demo, genotypes match gtc2vcf and GTCtoVCF with zero unexplained differences and PLINK files match plink2 +(`docs/validation/P6.md`). Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`, -`simurg qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`. +`simurg qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`, `simurg gtc-to-vcf`, +`simurg export plink|bedgraph|finalreport|locus-summary`. ## Quick start @@ -28,6 +33,8 @@ bash scripts/bootstrap_refs.sh # reference open-source repos → refs bash scripts/fetch_data.sh p0 # ~1.4 GB of public test data → data/ bash scripts/make_oracle_gtcs.sh # builds gtc2vcf, writes oracle GTCs bash scripts/make_oracle_qc.sh # locus-summary oracle (BeadArrayFiles, separate venv) +bash scripts/fetch_data.sh p6 # GRCh38 FASTA (~0.9 GB download, 3.1 GB unpacked) +bash scripts/make_oracle_vcf.sh # gtc2vcf oracle VCF for the P6 tests uv run pytest tests/oracle # real-data acceptance tests (add -m "not slow" to skip the long ones) uv run simurg inspect data/p0/gsa_v4_pgx/idats_Demo_36/idats_Demo_36/207895980019/207895980019_R03C01_Grn.idat ``` diff --git a/docs/validation/P6.md b/docs/validation/P6.md index 8bd6d10..33398e2 100644 --- a/docs/validation/P6.md +++ b/docs/validation/P6.md @@ -31,20 +31,37 @@ Both sides left-aligned with `bcftools norm -f`, compared per probe and sample a |---|---|---| | identical calls | 23 656 010 | | | identical no-calls | 102 077 | | -| our merged multi-allelic site differs from the single-probe call | 2 414 | e.g. an Infinium II [G/T] probe calls GG while another probe at the site shows G/C: C shares G's dye, so the probe cannot see it — our site genotype is consistent with every probe | -| our merged site missing, probe called | 3 945 | the probes at the site contradict each other (`W_MERGE_INCONSISTENT`) or leave several genotypes (`W_MERGE_AMBIGUOUS`) | -| our merged site called, probe no-call | 65 | other probes at the site resolve it | -| het in a haploid context → missing (OUT-018) | 14 | male X non-PAR and MT heterozygous calls | -| female Y → missing (OUT-018) | 11 | | -| indels gtc2vcf leaves unresolved (`REF=I/D`) that we resolve | 144 probes | checked against GTCtoVCF below | -| **unexplained differences** | **0** | | +| `grouping_differs` | 6 424 | the probe belongs to a SimurgArray multi-allelic site (kept merged even with `--unsquash-duplicates`, OUT-034); the site genotype uses every probe (e.g. an Infinium II [G/T] probe calls GG while the site is G/C: C shares G's dye) or is missing when the probes contradict each other | +| `ours_missing_het_in_haploid` | 16 | male X non-PAR and MT heterozygous calls → missing (OUT-018) | +| `ours_missing_haploid_Y` | 9 | Y of female samples → missing (OUT-018) | +| `representation_diff` | 5 184 (144 probes) | indels gtc2vcf leaves unresolved (`REF=I/D`); our representation agrees with GTCtoVCF (below) | +| **unexplained** | **0** | | + +Summary: `P6_gsa36_vs_gtc2vcf.json`. Locus sets: gtc2vcf writes 22 400 intensity-only loci that SimurgArray skips by design (OUT-019), and 49 indels it cannot resolve either (`W_SKIP_INDEL_REF` for us). Nothing we write is missing from gtc2vcf. -### Genotypes vs GTCtoVCF (Apache-2.0), 36 samples, merged sites +### Genotypes vs GTCtoVCF (Apache-2.0), 36 samples, merged sites (default squashing) -(filled in below) +GTCtoVCF 1.2.1 run under Python 3 with two local compatibility patches (K-022), alt-contig loci filtered. + +| Category | Count (probe × sample) | Explanation | +|---|---|---| +| identical calls | 23 656 575 | | +| identical no-calls | 109 760 | | +| `oracle_site_includes_probes_we_skip` | 2 917 | GTCtoVCF merges intensity-only CNV probes into SNV sites; their alleles make its call ambiguous. SimurgArray skips intensity-only loci (OUT-019) | +| `oracle_opposite_strand_indel_bug` | 177 | `[I/D]` and `[D/I]` duplicate indel probes all heterozygous: we write `0/1`, GTCtoVCF 1.2.1 `./.` — the reference bug named in OUT-020/EDGE-012 | +| `grouping_differs` | 155 | indel probes designed 1–3 bp apart normalise to the same variant: we merge them (OUT-030), GTCtoVCF groups by manifest position; the probes disagree, so the merged site is missing | +| `ours_missing_het_in_haploid` | 17 | OUT-018 | +| `ours_missing_haploid_Y` | 9 | OUT-018 | +| `oracle_diploid_candidates_in_haploid_context` | 2 | tri-allelic male-X site: we enumerate haploid candidates (OUT-031 step 3) and get one call; GTCtoVCF always enumerates diploid ones | +| representation differences | 0 | every indel (including the 144 gtc2vcf cannot resolve) is written identically after left-alignment | +| **unexplained** | **0** | | + +Probe sets: GTCtoVCF also writes the 22 400 intensity-only loci; it drops the 3 CYP2D6 indels whose `SourceSeq` +has the dash on the right (patched to skip instead of crash), which SimurgArray writes. Summary: +`P6_gsa36_vs_gtctovcf.json`. ### PLINK vs plink2 (GPL, subprocess oracle, built from source) diff --git a/docs/validation/P6_gsa36_vs_gtc2vcf.json b/docs/validation/P6_gsa36_vs_gtc2vcf.json new file mode 100644 index 0000000..4992dbf --- /dev/null +++ b/docs/validation/P6_gsa36_vs_gtc2vcf.json @@ -0,0 +1,1160 @@ +{ + "oracle": "gtc2vcf", + "samples": 36, + "probes_oracle": 682719, + "probes_ours": 660270, + "probes_in_merged_records": 5000, + "probes_only_in_oracle": 22449, + "only_in_oracle_examples": [ + "locus_0001", + "locus_0002", + "locus_0003", + "locus_0004", + "locus_0005", + "locus_0006", + "locus_0007", + "locus_0008", + "locus_0009", + "locus_0010", + "locus_0011", + "locus_0012", + "locus_0013", + "locus_0014", + "locus_0015", + "locus_0016", + "locus_0017", + "locus_0018", + "locus_0019", + "locus_0020" + ], + "categories": { + "grouping_differs": 6424, + "match_called": 23656010, + "match_missing": 102077, + "ours_missing_haploid_Y": 9, + "ours_missing_het_in_haploid": 16, + "representation_diff": 5184 + }, + "examples": { + "grouping_differs": [ + [ + "207895980019_R14C03", + "locus_0021", + [ + "1", + "19287570", + "G", + "T,C", + "0/2" + ], + "0/0", + [ + "1", + "19287570", + "G", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0022", + [ + "1", + "97515784", + "C", + "A,T", + "0/1" + ], + "./.", + [ + "1", + "97515784", + "C", + "A" + ] + ], + [ + "207895980019_R14C03", + "locus_0023", + [ + "1", + "97515784", + "C", + "A,T", + "0/1" + ], + "0/0", + [ + "1", + "97515784", + "C", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0024", + [ + "1", + "97515784", + "C", + "A,T", + "0/1" + ], + "./.", + [ + "1", + "97515784", + "C", + "A" + ] + ], + [ + "207895980019_R14C03", + "locus_0025", + [ + "1", + "108937356", + "G", + "C,T", + "./." + ], + "0/1", + [ + "1", + "108937356", + "G", + "C" + ] + ], + [ + "207895980019_R14C03", + "locus_0026", + [ + "1", + "108937356", + "G", + "C,T", + "./." + ], + "0/1", + [ + "1", + "108937356", + "G", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0027", + [ + "1", + "108937356", + "G", + "C,T", + "./." + ], + "0/1", + [ + "1", + "108937356", + "G", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0028", + [ + "1", + "186676356", + "A", + "C,G", + "./." + ], + "1/1", + [ + "1", + "186676356", + "A", + "G" + ] + ], + [ + "207895980019_R14C03", + "locus_0029", + [ + "1", + "186676356", + "A", + "C,G", + "./." + ], + "1/1", + [ + "1", + "186676356", + "A", + "G" + ] + ], + [ + "207895980019_R14C03", + "locus_0030", + [ + "1", + "186676356", + "A", + "C,G", + "./." + ], + "1/1", + [ + "1", + "186676356", + "A", + "G" + ] + ], + [ + "207895980019_R14C03", + "locus_0031", + [ + "1", + "192815708", + "T", + "A,C", + "0/2" + ], + "0/0", + [ + "1", + "192815708", + "T", + "A" + ] + ], + [ + "207895980019_R14C03", + "locus_0032", + [ + "10", + "33785635", + "G", + "A,T", + "./." + ], + "0/0", + [ + "10", + "33785635", + "G", + "A" + ] + ], + [ + "207895980019_R14C03", + "locus_0033", + [ + "10", + "33785635", + "G", + "A,T", + "./." + ], + "1/1", + [ + "10", + "33785635", + "G", + "A" + ] + ], + [ + "207895980019_R14C03", + "locus_0034", + [ + "10", + "33785635", + "G", + "A,T", + "./." + ], + "0/0", + [ + "10", + "33785635", + "G", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0035", + [ + "10", + "37501484", + "G", + "A,T", + "0/1" + ], + "0/0", + [ + "10", + "37501484", + "G", + "T" + ] + ], + [ + "207895980019_R14C03", + "locus_0036", + [ + "10", + "118836909", + "T", + "C,G", + "./." + ], + "1/1", + [ + "10", + "118836909", + "T", + "C" + ] + ], + [ + "207895980019_R14C03", + "locus_0037", + [ + "10", + "118836909", + "T", + "C,G", + "./." + ], + "0/0", + [ + "10", + "118836909", + "T", + "G" + ] + ], + [ + "207895980019_R14C03", + "locus_0038", + [ + "10", + "118836909", + "T", + "C,G", + "./." + ], + "1/1", + [ + "10", + "118836909", + "T", + "C" + ] + ], + [ + "207895980019_R14C03", + "locus_0039", + [ + "10", + "118836909", + "T", + "C,G", + "./." + ], + "1/1", + [ + "10", + "118836909", + "T", + "C" + ] + ], + [ + "207895980019_R14C03", + "locus_0040", + [ + "10", + "121586676", + "A", + "G,T", + "0/1" + ], + "0/0", + [ + "10", + "121586676", + "A", + "T" + ] + ] + ], + "representation_diff": [ + [ + "locus_0041", + [ + "1", + "20589096", + "TC", + "T" + ], + [ + "1", + "20589098", + "D", + "I" + ] + ], + [ + "locus_0042", + [ + "1", + "20589096", + "TC", + "T" + ], + [ + "1", + "20589098", + "D", + "I" + ] + ], + [ + "locus_0043", + [ + "1", + "20589096", + "TC", + "T" + ], + [ + "1", + "20589098", + "D", + "I" + ] + ], + [ + "locus_0044", + [ + "1", + "78303099", + "TG", + "T" + ], + [ + "1", + "78303101", + "D", + "I" + ] + ], + [ + "locus_0045", + [ + "1", + "78392093", + "AT", + "A" + ], + [ + "1", + "78392095", + "I", + "D" + ] + ], + [ + "locus_0046", + [ + "1", + "78411336", + "AC", + "A" + ], + [ + "1", + "78411338", + "D", + "I" + ] + ], + [ + "locus_0047", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0048", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0049", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0050", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0051", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0052", + [ + "1", + "97450065", + "TG", + "T" + ], + [ + "1", + "97450067", + "D", + "I" + ] + ], + [ + "locus_0053", + [ + "1", + "97593250", + "CT", + "C" + ], + [ + "1", + "97593252", + "I", + "D" + ] + ], + [ + "locus_0054", + [ + "1", + "97593250", + "CT", + "C" + ], + [ + "1", + "97593252", + "I", + "D" + ] + ], + [ + "locus_0055", + [ + "1", + "97593250", + "CT", + "C" + ], + [ + "1", + "97593252", + "I", + "D" + ] + ], + [ + "locus_0056", + [ + "1", + "97593250", + "CT", + "C" + ], + [ + "1", + "97593252", + "I", + "D" + ] + ], + [ + "locus_0057", + [ + "1", + "97593250", + "CT", + "C" + ], + [ + "1", + "97593252", + "I", + "D" + ] + ], + [ + "locus_0058", + [ + "1", + "147208902", + "CT", + "C" + ], + [ + "1", + "147208904", + "I", + "D" + ] + ], + [ + "locus_0059", + [ + "1", + "147208902", + "CT", + "C" + ], + [ + "1", + "147208904", + "I", + "D" + ] + ], + [ + "locus_0060", + [ + "1", + "147208902", + "CT", + "C" + ], + [ + "1", + "147208904", + "I", + "D" + ] + ] + ], + "ours_missing_haploid_Y": [ + [ + "207895980019_R11C02", + "locus_0061", + [ + "Y", + "19428214", + "G", + "C", + "." + ], + "0/0", + [ + "Y", + "19428214", + "G", + "C" + ] + ], + [ + "207895980035_R10C01", + "locus_0062", + [ + "Y", + "14081208", + "A", + "G", + "." + ], + "1/1", + [ + "Y", + "14081208", + "A", + "G" + ] + ], + [ + "207895980019_R12C01", + "locus_0063", + [ + "Y", + "12030916", + "C", + "T", + "." + ], + "0/0", + [ + "Y", + "12030916", + "C", + "T" + ] + ], + [ + "207895980019_R12C01", + "locus_0064", + [ + "Y", + "13331090", + "T", + "C", + "." + ], + "1/1", + [ + "Y", + "13331090", + "T", + "C" + ] + ], + [ + "207895980035_R04C03", + "locus_0065", + [ + "Y", + "9957467", + "A", + "T", + "." + ], + "0/0", + [ + "Y", + "9957467", + "A", + "T" + ] + ], + [ + "207895980035_R04C03", + "locus_0066", + [ + "Y", + "12816261", + "A", + "G", + "." + ], + "1/1", + [ + "Y", + "12816261", + "A", + "G" + ] + ], + [ + "207895980019_R03C01", + "locus_0067", + [ + "Y", + "19466559", + "T", + "G", + "." + ], + "1/1", + [ + "Y", + "19466559", + "T", + "G" + ] + ], + [ + "207895980019_R04C01", + "locus_0068", + [ + "Y", + "5680808", + "T", + "G", + "." + ], + "1/1", + [ + "Y", + "5680808", + "T", + "G" + ] + ], + [ + "207895980019_R03C02", + "locus_0069", + [ + "Y", + "12580293", + "T", + "C", + "." + ], + "0/0", + [ + "Y", + "12580293", + "T", + "C" + ] + ] + ], + "ours_missing_het_in_haploid": [ + [ + "207895980035_R13C03", + "locus_0070", + [ + "X", + "7467632", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "7467632", + "A", + "G" + ] + ], + [ + "207895980035_R13C03", + "locus_0071", + [ + "X", + "7478861", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "7478861", + "T", + "C" + ] + ], + [ + "207895980035_R13C03", + "locus_0072", + [ + "X", + "68391913", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "68391913", + "C", + "T" + ] + ], + [ + "207895980019_R06C01", + "locus_0073", + [ + "X", + "33006670", + "T", + "G", + "." + ], + "0/1", + [ + "X", + "33006670", + "T", + "G" + ] + ], + [ + "207895980035_R09C03", + "locus_0074", + [ + "MT", + "2690", + "G", + "A", + "." + ], + "0/1", + [ + "MT", + "2690", + "G", + "A" + ] + ], + [ + "207895980019_R15C02", + "locus_0075", + [ + "X", + "79535086", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "79535086", + "C", + "T" + ] + ], + [ + "207895980035_R14C01", + "locus_0076", + [ + "X", + "120442599", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "120442599", + "C", + "T" + ] + ], + [ + "207895980035_R03C03", + "locus_0077", + [ + "X", + "22216499", + "A", + "T", + "." + ], + "0/1", + [ + "X", + "22216499", + "A", + "T" + ] + ], + [ + "207895980019_R09C03", + "locus_0078", + [ + "Y", + "7044951", + "G", + "T", + "." + ], + "0/1", + [ + "Y", + "7044951", + "G", + "T" + ] + ], + [ + "207895980035_R12C03", + "locus_0079", + [ + "X", + "20830034", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "20830034", + "A", + "G" + ] + ], + [ + "207895980035_R11C03", + "locus_0080", + [ + "X", + "91616889", + "G", + "A", + "." + ], + "0/1", + [ + "X", + "91616889", + "G", + "A" + ] + ], + [ + "207895980035_R11C03", + "locus_0081", + [ + "X", + "154532408", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "154532408", + "T", + "C" + ] + ], + [ + "207895980042_R15C02", + "locus_0082", + [ + "Y", + "15734824", + "C", + "A", + "." + ], + "0/1", + [ + "Y", + "15734824", + "C", + "A" + ] + ], + [ + "207895980035_R05C03", + "locus_0083", + [ + "X", + "7462384", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "7462384", + "A", + "G" + ] + ], + [ + "207895980035_R05C03", + "locus_0084", + [ + "X", + "7475195", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "7475195", + "T", + "C" + ] + ], + [ + "207895980035_R05C03", + "locus_0072", + [ + "X", + "68391913", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "68391913", + "C", + "T" + ] + ] + ] + }, + "unexplained": 0, + "note": "Probe names from the vendor manifest are replaced by stable pseudonyms (locus_NNNN, shared by both oracle files); sample IDs, positions, alleles and all counts are unchanged." +} diff --git a/docs/validation/P6_gsa36_vs_gtctovcf.json b/docs/validation/P6_gsa36_vs_gtctovcf.json new file mode 100644 index 0000000..d8054bc --- /dev/null +++ b/docs/validation/P6_gsa36_vs_gtctovcf.json @@ -0,0 +1,1641 @@ +{ + "oracle": "GTCtoVCF", + "samples": 36, + "probes_oracle": 682667, + "probes_ours": 660270, + "probes_in_merged_records": 37785, + "probes_only_in_oracle": 22400, + "only_in_oracle_examples": [ + "locus_0001", + "locus_0002", + "locus_0003", + "locus_0004", + "locus_0005", + "locus_0006", + "locus_0007", + "locus_0008", + "locus_0009", + "locus_0010", + "locus_0011", + "locus_0012", + "locus_0013", + "locus_0014", + "locus_0015", + "locus_0016", + "locus_0017", + "locus_0018", + "locus_0019", + "locus_0020" + ], + "categories": { + "grouping_differs": 155, + "match_called": 23656575, + "match_missing": 109760, + "only_ours": 108, + "oracle_diploid_candidates_in_haploid_context": 2, + "oracle_opposite_strand_indel_bug": 177, + "oracle_site_includes_probes_we_skip": 2917, + "ours_missing_haploid_Y": 9, + "ours_missing_het_in_haploid": 17 + }, + "examples": { + "oracle_site_includes_probes_we_skip": [ + [ + "207895980019_R03C01", + "locus_0085", + [ + "1", + "201085496", + "G", + "A", + "0/0" + ], + "./.", + [ + "1", + "201085496", + "G", + "C,A" + ] + ], + [ + "207895980019_R03C01", + "locus_0086", + [ + "1", + "201085496", + "G", + "A", + "0/0" + ], + "./.", + [ + "1", + "201085496", + "G", + "C,A" + ] + ], + [ + "207895980019_R03C01", + "locus_0087", + [ + "10", + "94981305", + "C", + "G", + "0/0" + ], + "./.", + [ + "10", + "94981305", + "C", + "A,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0088", + [ + "10", + "94981305", + "C", + "G", + "0/0" + ], + "./.", + [ + "10", + "94981305", + "C", + "A,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0089", + [ + "10", + "94981305", + "C", + "G", + "0/0" + ], + "./.", + [ + "10", + "94981305", + "C", + "A,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0090", + [ + "10", + "95063364", + "C", + "G", + "0/0" + ], + "./.", + [ + "10", + "95063364", + "C", + "T,A,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0091", + [ + "12", + "21239885", + "T", + "C", + "0/0" + ], + "./.", + [ + "12", + "21239885", + "T", + "C,A" + ] + ], + [ + "207895980019_R03C01", + "locus_0092", + [ + "16", + "55822444", + "C", + "A", + "0/0" + ], + "./.", + [ + "16", + "55822444", + "C", + "A,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0093", + [ + "18", + "657095", + "C", + "T", + "1/1" + ], + "./.", + [ + "18", + "657095", + "C", + "T,A" + ] + ], + [ + "207895980019_R03C01", + "locus_0094", + [ + "18", + "664286", + "G", + "A", + "1/1" + ], + "./.", + [ + "18", + "664286", + "G", + "A,T" + ] + ], + [ + "207895980019_R03C01", + "locus_0095", + [ + "19", + "39244961", + "C", + "T", + "0/0" + ], + "./.", + [ + "19", + "39244961", + "C", + "G,T" + ] + ], + [ + "207895980019_R03C01", + "locus_0096", + [ + "19", + "39245004", + "C", + "G", + "1/1" + ], + "./.", + [ + "19", + "39245004", + "C", + "A,T,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0097", + [ + "19", + "39245004", + "C", + "G", + "1/1" + ], + "./.", + [ + "19", + "39245004", + "C", + "A,T,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0098", + [ + "19", + "39245004", + "C", + "G", + "1/1" + ], + "./.", + [ + "19", + "39245004", + "C", + "A,T,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0099", + [ + "19", + "39245004", + "C", + "G", + "1/1" + ], + "./.", + [ + "19", + "39245004", + "C", + "A,T,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0100", + [ + "19", + "40846021", + "G", + "A", + "0/0" + ], + "./.", + [ + "19", + "40846021", + "G", + "A,C" + ] + ], + [ + "207895980019_R03C01", + "locus_0101", + [ + "19", + "40849103", + "A", + "C", + "0/0" + ], + "./.", + [ + "19", + "40849103", + "A", + "T,C" + ] + ], + [ + "207895980019_R03C01", + "locus_0102", + [ + "2", + "85552428", + "C", + "T", + "0/0" + ], + "./.", + [ + "2", + "85552428", + "C", + "G,T" + ] + ], + [ + "207895980019_R03C01", + "locus_0103", + [ + "2", + "166073440", + "A", + "G", + "0/0" + ], + "./.", + [ + "2", + "166073440", + "A", + "T,G" + ] + ], + [ + "207895980019_R03C01", + "locus_0104", + [ + "22", + "42126481", + "G", + "T", + "0/0" + ], + "./.", + [ + "22", + "42126481", + "G", + "C,T" + ] + ] + ], + "grouping_differs": [ + [ + "207895980019_R03C01", + "locus_0105", + [ + "11", + "113904829", + "GAGA", + "G", + "./." + ], + "0/0", + [ + "11", + "113904829", + "GAGA", + "G" + ] + ], + [ + "207895980019_R03C01", + "locus_0106", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R03C01", + "locus_0107", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R03C02", + "locus_0105", + [ + "11", + "113904829", + "GAGA", + "G", + "./." + ], + "0/0", + [ + "11", + "113904829", + "GAGA", + "G" + ] + ], + [ + "207895980019_R03C02", + "locus_0106", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R03C02", + "locus_0107", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R04C01", + "locus_0105", + [ + "11", + "113904829", + "GAGA", + "G", + "./." + ], + "0/0", + [ + "11", + "113904829", + "GAGA", + "G" + ] + ], + [ + "207895980019_R04C01", + "locus_0106", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R04C01", + "locus_0107", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/1", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0105", + [ + "11", + "113904829", + "GAGA", + "G", + "./." + ], + "0/0", + [ + "11", + "113904829", + "GAGA", + "G" + ] + ], + [ + "207895980019_R05C01", + "locus_0106", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/0", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0107", + [ + "18", + "673443", + "TTTAAAG", + "T", + "./." + ], + "0/0", + [ + "18", + "673443", + "TTTAAAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0108", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T", + "0/1" + ], + "./.", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0109", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T", + "0/1" + ], + "./.", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0110", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T", + "0/1" + ], + "./.", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0111", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T", + "0/1" + ], + "./.", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T" + ] + ], + [ + "207895980019_R05C01", + "locus_0112", + [ + "6", + "160139865", + "CTGGTAAGT", + "C", + "./." + ], + "0/0", + [ + "6", + "160139865", + "CTGGTAAGT", + "C" + ] + ], + [ + "207895980019_R05C01", + "locus_0113", + [ + "6", + "160139865", + "CTGGTAAGT", + "C", + "./." + ], + "0/0", + [ + "6", + "160139865", + "CTGGTAAGT", + "C" + ] + ], + [ + "207895980019_R06C01", + "locus_0105", + [ + "11", + "113904829", + "GAGA", + "G", + "./." + ], + "0/1", + [ + "11", + "113904829", + "GAGA", + "G" + ] + ], + [ + "207895980019_R06C01", + "locus_0108", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T", + "0/1" + ], + "./.", + [ + "3", + "46373452", + "TACAGTCAGTATCAATTCTGGAAGAATTTCCAG", + "T" + ] + ] + ], + "oracle_opposite_strand_indel_bug": [ + [ + "207895980019_R03C01", + "locus_0114", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C01", + "locus_0115", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C01", + "locus_0116", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C01", + "locus_0117", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C01", + "locus_0118", + [ + "6", + "33087504", + "TC", + "T", + "0/1" + ], + "./.", + [ + "6", + "33087504", + "TC", + "T" + ] + ], + [ + "207895980019_R03C01", + "locus_0119", + [ + "6", + "33087504", + "TC", + "T", + "0/1" + ], + "./.", + [ + "6", + "33087504", + "TC", + "T" + ] + ], + [ + "207895980019_R03C01", + "locus_0120", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R03C01", + "locus_0121", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R03C01", + "locus_0122", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R03C02", + "locus_0114", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C02", + "locus_0115", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C02", + "locus_0116", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C02", + "locus_0117", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R03C02", + "locus_0118", + [ + "6", + "33087504", + "TC", + "T", + "0/1" + ], + "./.", + [ + "6", + "33087504", + "TC", + "T" + ] + ], + [ + "207895980019_R03C02", + "locus_0119", + [ + "6", + "33087504", + "TC", + "T", + "0/1" + ], + "./.", + [ + "6", + "33087504", + "TC", + "T" + ] + ], + [ + "207895980019_R03C02", + "locus_0120", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R03C02", + "locus_0121", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R03C02", + "locus_0122", + [ + "7", + "140836521", + "G", + "GT", + "0/1" + ], + "./.", + [ + "7", + "140836521", + "G", + "GT" + ] + ], + [ + "207895980019_R04C01", + "locus_0114", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ], + [ + "207895980019_R04C01", + "locus_0115", + [ + "12", + "21334835", + "C", + "CT", + "0/1" + ], + "./.", + [ + "12", + "21334835", + "C", + "CT" + ] + ] + ], + "ours_missing_haploid_Y": [ + [ + "207895980019_R03C01", + "locus_0067", + [ + "Y", + "19466559", + "T", + "G", + "." + ], + "1", + [ + "Y", + "19466559", + "T", + "G" + ] + ], + [ + "207895980019_R03C02", + "locus_0069", + [ + "Y", + "12580293", + "T", + "C", + "." + ], + "0", + [ + "Y", + "12580293", + "T", + "C" + ] + ], + [ + "207895980019_R04C01", + "locus_0068", + [ + "Y", + "5680808", + "T", + "G", + "." + ], + "1", + [ + "Y", + "5680808", + "T", + "G" + ] + ], + [ + "207895980019_R11C02", + "locus_0061", + [ + "Y", + "19428214", + "G", + "C", + "." + ], + "0", + [ + "Y", + "19428214", + "G", + "C" + ] + ], + [ + "207895980019_R12C01", + "locus_0063", + [ + "Y", + "12030916", + "C", + "T", + "." + ], + "0", + [ + "Y", + "12030916", + "C", + "T" + ] + ], + [ + "207895980019_R12C01", + "locus_0064", + [ + "Y", + "13331090", + "T", + "C", + "." + ], + "1", + [ + "Y", + "13331090", + "T", + "C" + ] + ], + [ + "207895980035_R04C03", + "locus_0065", + [ + "Y", + "9957467", + "A", + "T", + "." + ], + "0", + [ + "Y", + "9957467", + "A", + "T" + ] + ], + [ + "207895980035_R04C03", + "locus_0066", + [ + "Y", + "12816261", + "A", + "G", + "." + ], + "1", + [ + "Y", + "12816261", + "A", + "G" + ] + ], + [ + "207895980035_R10C01", + "locus_0062", + [ + "Y", + "14081208", + "A", + "G", + "." + ], + "1", + [ + "Y", + "14081208", + "A", + "G" + ] + ] + ], + "ours_missing_het_in_haploid": [ + [ + "207895980019_R06C01", + "locus_0073", + [ + "X", + "33006670", + "T", + "G", + "." + ], + "0/1", + [ + "X", + "33006670", + "T", + "G" + ] + ], + [ + "207895980019_R09C03", + "locus_0123", + [ + "Y", + "7044951", + "G", + "T", + "." + ], + "0/1", + [ + "Y", + "7044951", + "G", + "T" + ] + ], + [ + "207895980019_R09C03", + "locus_0078", + [ + "Y", + "7044951", + "G", + "T", + "." + ], + "0/1", + [ + "Y", + "7044951", + "G", + "T" + ] + ], + [ + "207895980019_R15C02", + "locus_0075", + [ + "X", + "79535086", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "79535086", + "C", + "T" + ] + ], + [ + "207895980035_R03C03", + "locus_0077", + [ + "X", + "22216499", + "A", + "T", + "." + ], + "0/1", + [ + "X", + "22216499", + "A", + "T" + ] + ], + [ + "207895980035_R05C03", + "locus_0083", + [ + "X", + "7462384", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "7462384", + "A", + "G" + ] + ], + [ + "207895980035_R05C03", + "locus_0084", + [ + "X", + "7475195", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "7475195", + "T", + "C" + ] + ], + [ + "207895980035_R05C03", + "locus_0072", + [ + "X", + "68391913", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "68391913", + "C", + "T" + ] + ], + [ + "207895980035_R09C03", + "locus_0074", + [ + "MT", + "2690", + "G", + "A", + "." + ], + "0/1", + [ + "MT", + "2690", + "G", + "A" + ] + ], + [ + "207895980035_R11C03", + "locus_0080", + [ + "X", + "91616889", + "G", + "A", + "." + ], + "0/1", + [ + "X", + "91616889", + "G", + "A" + ] + ], + [ + "207895980035_R12C03", + "locus_0079", + [ + "X", + "20830034", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "20830034", + "A", + "G" + ] + ], + [ + "207895980035_R13C03", + "locus_0070", + [ + "X", + "7467632", + "A", + "G", + "." + ], + "0/1", + [ + "X", + "7467632", + "A", + "G" + ] + ], + [ + "207895980035_R13C03", + "locus_0124", + [ + "X", + "7478861", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "7478861", + "T", + "C" + ] + ], + [ + "207895980035_R13C03", + "locus_0071", + [ + "X", + "7478861", + "T", + "C", + "." + ], + "0/1", + [ + "X", + "7478861", + "T", + "C" + ] + ], + [ + "207895980035_R13C03", + "locus_0072", + [ + "X", + "68391913", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "68391913", + "C", + "T" + ] + ], + [ + "207895980035_R14C01", + "locus_0076", + [ + "X", + "120442599", + "C", + "T", + "." + ], + "0/1", + [ + "X", + "120442599", + "C", + "T" + ] + ], + [ + "207895980042_R15C02", + "locus_0082", + [ + "Y", + "15734824", + "C", + "A", + "." + ], + "0/1", + [ + "Y", + "15734824", + "C", + "A" + ] + ] + ], + "oracle_diploid_candidates_in_haploid_context": [ + [ + "207895980035_R11C03", + "locus_0125", + [ + "X", + "22890860", + "A", + "G,T", + "1" + ], + ".", + [ + "X", + "22890860", + "A", + "G,T" + ] + ], + [ + "207895980035_R11C03", + "locus_0126", + [ + "X", + "22890860", + "A", + "G,T", + "1" + ], + ".", + [ + "X", + "22890860", + "A", + "G,T" + ] + ] + ] + }, + "unexplained": 0, + "note": "Probe names from the vendor manifest are replaced by stable pseudonyms (locus_NNNN, shared by both oracle files); sample IDs, positions, alleles and all counts are unchanged." +} diff --git a/scripts/validate_p6_vcf.py b/scripts/validate_p6_vcf.py index 9aa0f5f..3478c6a 100644 --- a/scripts/validate_p6_vcf.py +++ b/scripts/validate_p6_vcf.py @@ -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