GLnexus: joint variant calling for large cohort sequencing

Key results

1. ShownRun on identical gVCFs, GLnexus and GATK made nearly identical genotype calls at corresponding sites, but GLnexus kept more alleles, especially rare indels; its indel sites with no GATK counterpart were nearly all seen in ten or fewer copies.

   Chromosome 2 exome targets, 50K DiscovEHR exomes: genotypes disagreed in 0.00073% of 12.8 billion matrix cells, or 0.12% of the 80 million cells where either pVCF showed variation; over 33% more rare indels; of 3,228 biallelic indel sites absent from GATK's pVCF, 98.7% had ten or fewer copies and 45.2% were singletons (median genotype quality, GQ, 74; the rest median 57). · 3. VALIDATION STUDY WITH 50,000 EXOMES

2. ShownTrio checks indicate the added sensitivity did not cost genotype quality: GLnexus calls showed fewer apparent Mendelian violations than GATK, most clearly for indels, and a more unified site representation.

   877 child–mother–father trios, violations used as an error proxy: GATK flagged a violation in 8.2% of 44,107 indel-variant genotype triplets vs 6.8% for GLnexus. GLnexus had more multiallelic sites (5.9% vs 4.5%) and fewer overlapping sites (1.5% vs 2.8%); monoallelic fallback sites were 0.35% of sites. · 3. VALIDATION STUDY WITH 50,000 EXOMES

3. ShownGLnexus, with its ramped allele-frequency prior, kept rare-variant sensitivity stable as the cohort grew, whereas GATK increasingly dropped very rare sites of below-average call quality when the same samples were joint-called inside larger cohorts.

   Nested random subsets: of GATK's 77,138 sites at 5K, 309 (0.40%) were absent at 10K, 465 (0.60%) at 25K and 492 (0.64%) at 50K (median GQ 71 for the 492). GLnexus lost fewer than 10 of its 5K sites, attributed to low-complexity-region filtering. · Joint calling scaling from N=1K to N=50K; Empirical frequency prior

4. ShownGLnexus consumed several-fold fewer core-hours than GATK 3.8, and the advantage grew with cohort size even under conservative timing.

   Nested 1K–50K cohorts: 3.8-fold fewer core-hours at 1K rising to 6.6-fold fewer at 50K. GLnexus timings began from the original gVCFs and included a final Spark stage writing one pVCF file; the authors expect incremental, sharded operation to cost relatively less. · Runtime comparison

5. ShownThe distributed, cohort-sharded deployment has already produced production pVCFs for hundreds of thousands of exomes and tens of thousands of whole genomes.

   243,953 exomes: 33 TB compressed gVCF to 6.8 TB pVCF (2.7 trillion entries) in about 36 h on 1,600 threads, roughly 14 thread-minutes per exome. 22,609 genomes (36.6x median coverage): 6.7 trillion-entry matrix from 481 TB of CRAM in about 35 h, under 7 core-hours per genome, for the combined CRAM-to-pVCF workflow. · 5. APPLICATIONS