This blog is part of our Paper Spotlight series, which features peer-reviewed research publications involving work done in Terra and highlights how the analysis methods were applied.
Mitochondrial DNA variation across 56,434 individuals in gnomAD
By Kristen M. Laricchia, Nicole J. Lake, […] and Sarah E. Calvo, 2022
Genome Research 32(3):569-582 (2022) https://doi.org/10.1101/gr.276013.121
Abstract: Genomic databases of allele frequency are extremely helpful for evaluating clinical variants of unknown significance; however, until now, databases such as the Genome Aggregation Database (gnomAD) have focused on nuclear DNA and have ignored the mitochondrial genome (mtDNA). Here, we present a pipeline to call mtDNA variants that addresses three technical challenges: (1) detecting homoplasmic and heteroplasmic variants, present, respectively, in all or a fraction of mtDNA molecules; (2) circular mtDNA genome; and (3) misalignment of nuclear sequences of mitochondrial origin (NUMTs). We observed that mtDNA copy number per cell varied across gnomAD cohorts and influenced the fraction of NUMT-derived false-positive variant calls, which can account for the majority of putative heteroplasmies. To avoid false positives, we excluded contaminated samples, cell lines, and samples prone to NUMT misalignment due to few mtDNA copies. Furthermore, we report variants with heteroplasmy ≥10%. We applied this pipeline to 56,434 whole-genome sequences in the gnomAD v3.1 database that includes individuals of European (58%), African (25%), Latino (10%), and Asian (5%) ancestry. Our gnomAD v3.1 release contains population frequencies for 10,850 unique mtDNA variants at more than half of all mtDNA bases. Importantly, we report frequencies within each nuclear ancestral population and mitochondrial haplogroup. Homoplasmic variants account for most variant calls (98%) and unique variants (85%). We observed that 1/250 individuals carry a pathogenic mtDNA variant with heteroplasmy above 10%. These mtDNA population allele frequencies are freely accessible and will aid in diagnostic interpretation and research studies.
What part of the work was done in Terra?
Excerpts from the paper’s Methods section:
Mitochondrial variant calling pipeline in single samples
WGS data were aligned to reference genome GRCh38, which includes Chr M […] using BWA-MEM version 0.7.15-r1140 […]. For each sample CRAM, Terra MitochondrialPipeline version 25 was run.
Briefly, GATK version 4.1.2.0 tools were used to estimate the median nuclear genome coverage, to exclude duplicates, to pull reads from Chr M, and to call variants. […] Mutect2 variants were then filtered; multi-allelic sites were split into different variants; and HaploGrep/HaploCheck was run to assign haplogroup and estimate mtDNA contamination. The min_vaf_threshold was set to 0.01 and calls below 0.01 VAF were later set to homoplasmic reference. For each input sample, a VCF with mtDNA variants was produced. […]
How did they do it?
To call mitochondrial variants from whole genome samples, the authors developed a new pipeline based on the GATK Best Practices, using the Mutect2 somatic variant caller with certain adaptations designed to address the idiosyncrasies of mitochondrial genomes. The methodology is described in full in the paper; there is also step-by-step documentation of each tool involved on the GATK website.
The authors implemented their new pipeline in the Workflow Description Language (WDL) and shared it in the Broad Methods Repository under mitochondria/MitochondriaPipeline (requires sign-in). In the paper, they used version 25 of the pipeline (note that in the repository, versions are called “snapshots”).
Terra also supports importing workflows from Dockstore, a free and open-source platform for sharing reusable and scalable analytical tools and workflows.
The authors ran the workflows at scale using Terra’s workflow execution service. To try your hand at running a workflow in Terra, check out this Quickstart Tutorial Workspace.
You can get started with the MitochondriaPipeline in this Terra workspace preconfigured to run on example data.
Appendix: Data and code availability
- Variants and population frequencies generated by this study are available in gnomAD v3.1. A user-friendly website provides variant annotations, including distributions across heteroplasmy levels, populations, haplogroups, and age.
- Data are available for download in multiple formats, including VCF, Hail Table, and simple tab-delimited files.
- The Mutect2 pipeline is available through the GATK repository in Github and through the Terra Methods Repository (version 25 was used in this paper).
- The authors also provide a simplified version of the pipeline (MitochondriaPipelineSlim.wdl) that omits the realignment step to reduce complexity, input parameters, and cost, without substantial differences to the output VCF.
Other materials related to the paper:
- The Hail scripts used for combining the VCFs, filtering samples and variants, adding annotations, and performing analyses can be found in the gnomad project Github repository.
- Scripts are also available for download in the Supplemental Code archive provided alongside the paper on the Genome Research website.
