SimurgArray/docs/validation/P6.md

Phase P6 validation — outputs (VCF/BCF, PLINK, bedGraph, final report)

Date: 2026-09-28. All work on the Raspberry Pi (no cloud cost). Data: the 36 GSA v4 PGx oracle GTCs (P4, byte-identical to gtc2vcf's), the GSA BPM + CSV manifest, and the Ensembl GRCh38 primary assembly (release 110, fetch_data.sh p6).

Delivered

Item	Module
FASTA access: plain/BGZF, .fai/.gzi read and built (identical to htslib's), index check, build detection, contig aliases, mmap genome cache	formats/fasta.py
BGZF writer, tabix .tbi, BCF 2.2 encoder, CSI index (pure Python)	formats/bgzf.py, formats/bcf.py
Plus-strand mapping with audit row per locus; indel resolution from SourceSeq (context score), left normalisation	harmonize/sites.py
Site grouping and the OUT-031 genotype-combination algorithm (ploidy-aware, Infinium I/II signal classes)	harmonize/combine.py
simurg gtc-to-vcf (per-sample VCF.gz + .tbi, --merge cohort BCF + .csi, audit Parquet)	pipeline/vcf.py
`simurg export bedgraph	finalreport

Acceptance (plan §18 P6)

bcftools validation and indexes

Every VCF passes bcftools view; bcftools index --nrecords on our .tbi equals the record count (635 796); region queries through our index (1:1-2M, X:150-156M, MT, Y, 22:30.0-30.1M) return exactly the records a full scan finds. The cohort BCF reads back identically to the per-sample VCFs (635 796 records × 3 samples, 0 differences in GT/GS/BAF/LRR) and its CSI index answers region queries.

Genotypes vs gtc2vcf (MIT), 36 samples, one record per probe (--unsquash-duplicates)

Both sides left-aligned with bcftools norm -f, compared per probe and sample as allele sets.

Category	Count (probe × sample)	Explanation
identical calls	23 656 010
identical no-calls	102 077
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 (default squashing)

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)

simurg export plink and plink2 --bcf <bcftools norm -m- of our BCF> --make-bed --split-par hg38 agree on all 637 397 variants × 3 samples (bed codes identical; the 1 601 extra variants are split multi-allelic sites).

BAF orientation, bedGraph

• EDGE-017: at every called homozygous single-probe site where the probe's A allele is ALT, VCF BAF = 1 − GTC BAF (float32), otherwise unchanged (> 500 sites of each kind checked per sample).
• EDGE-023: raw LRR bedGraph values equal the GTC values after %.6g formatting for every single-probe position.

Optional outputs (completed 2026-09-28)

Option	Check	Result
--info-fields ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS (--egt for the cluster scores)	INFO vs gtc2vcf ASSAY_TYPE, GenTrain_Score, Cluster_Sep, 3 samples	identical
--extra-format X,Y,NORMX,NORMY,THETA,R	vs gtc2vcf FORMAT tags, 1 965 810 probe × sample values	X, Y, NORMX, NORMY identical (within the 6-digit print of bcftools query); THETA/R identical except 148 values where both normalized intensities are 0: missing for SimurgArray (plan §7.4, as the caller sees them), 0 for gtc2vcf
both, in --merge BCF	BCF vs per-sample VCFs, 1 970 064 records × all fields	0 differences
--merge-snv-indel (OUT-021, EDGE-013)	1 GSA sample; bcftools norm -c e -f	677 SNV+indel positions merged (635 796 → 635 119 records), 0 REF mismatches; unit test of the combined genotypes
--auxiliary-loci (OUT-035, EDGE-015)	unit tests for every rule (a)–(e) and a minus-strand design; 3 real probes (plus and minus designs) redefined as MNVs	records replaced, genotypes preserved, 0 REF mismatches
`export plink --format pgen	both` (plan §8.4)	plink2 --pfile on our .pgen → .bed

The example rows in P6_gsa36_vs_gtc2vcf.json and P6_gsa36_vs_gtctovcf.json show probe names as stable pseudonyms (locus_NNNN); sample IDs, positions, alleles and all counts are the original values.