Identifying Top Mitochondrial Variants from mgatk
Using previous authors’ own confirmed thresholds, with the field-standard VMR ranking and mitochondrial read depth thresholding.
Rachael White
September 17, 2026
Overview
This analysis investigates candidate mitochondrial variants across our 3-sample celltype-annotated object, using thresholds confirmed directly from two sources: Mullin et al.’s own executed code, and the Lareau et al. Nature Protocols paper (the mtscATAC-seq/mgatk method’s own origin, https://doi.org/10.1038/s41596-022-00795-3).
The works above both used 10X as a minimum cutoff for mitochondrial read depth. This cutoff was used accordingly when generating the data analyzed here.
Thresholds used, and why:
Minimum cutoff for mitochondrial read depth: 10X – confirmed quality threshold used by Mullin & Lareau
n_cells_conf_detected >= 5– confirmed active in Mullin et al.’s codestrand_correlation >= 0.65– confirmed active in Mullin et al.’s codevmris NOT used as a hard filter (Mullin’s own code considered it, but it was commented out in the code and not mentioned numerically in the paper.) It IS used below for ranking, since the Lareau protocol paper specifically describes strand correlation + VMR together as the standard way to identify informative variants.
A threshold deliberately NOT used:
Some clonal-tracing literature adds
mean heteroplasmy < 0.5, to exclude near-homoplasmic
variants that don’t help distinguish cell lineages. Mullin et al. didn’t
do clonal tracing, so they didn’t discriminate between high and low
heteroplasmy values. We actually need that low-heteroplasmy information
here to help distinguish interesting variants.
Setup
annotated_object_path <- "outputs/checkpoint_annotated.rds"
# Any ONE sample's mgatk final/ directory -- only needed to pull the
# reference allele calls, which are genome-based, not sample-specific.
# Doesn't matter which sample you pick.
mgatk_final_dir_for_refallele <- "mgatk_out/melas_1_rpe_choroid_rep1_mgatk/final"
# Confirmed thresholds (see Overview above)
MIN_CELLS_DETECTED <- 5
MIN_STRAND_CORR <- 0.65mine <- readRDS(annotated_object_path)
message("Loaded object: ", ncol(mine), " cells")1. Run IdentifyVariants() and inspect the raw output
DefaultAssay(mine) <- "mito"
variable_sites <- IdentifyVariants(mine, assay = "mito", refallele = mito_ref_data$refallele)#sanity checks
#message(nrow(variable_sites), " total candidate variant positions/alleles found")
#head(variable_sites)
#Confirm the columns we expect are actually #there- this should include at minimum #`variant`, `position`,#
#`n_cells_conf_detected`, `strand_correlation`, #`vmr`, and `mean`.
#colnames(variable_sites)3. Apply the confirmed QC thresholds
high_conf <- variable_sites %>%
filter(n_cells_conf_detected >= MIN_CELLS_DETECTED,
strand_correlation >= MIN_STRAND_CORR)
message(nrow(high_conf), " variants pass QC thresholds ",
"(out of ", nrow(variable_sites), " total candidates)")4. Sanity check: does m.3243A>G actually appear in this list?
check_3243 <- high_conf %>% filter(position == 3243)
if (nrow(check_3243) == 0) {
warning("m.3243A>G did NOT pass the QC thresholds above.")
} else {
message("m.3243A>G present in the high-confidence variant list.")
}## m.3243A>G present in the high-confidence variant list.print(check_3243)## position nucleotide variant vmr mean variance
## 3243A>G 3243 A>G 3243A>G 0.01042155 0.9518696 0.00992
## n_cells_conf_detected n_cells_over_5 n_cells_over_10 n_cells_over_20
## 3243A>G 347 541 540 538
## strand_correlation mean_coverage
## 3243A>G 0.9415496 9.7607895. Rank and view the top variants
Ranked by VMR, matching the Lareau protocol paper’s own described approach for identifying the most informative variants.
We’ll treat this as the standard way to prioritize which high-confidence variants are most worth a closer look. We also explore a few alternative methods later in this report.
top_variants <- high_conf %>%
arrange(desc(vmr)) %>%
select(variant, position, n_cells_conf_detected, strand_correlation, vmr, mean)
head(top_variants, 15)## variant position n_cells_conf_detected strand_correlation
## 72T>C 72T>C 72 6 0.9971011
## 3243A>G 3243A>G 3243 347 0.9415496
## 9448A>G 9448A>G 9448 312 0.9134332
## 310T>C 310T>C 310 177 0.7662445
## 16189T>C 16189T>C 16189 236 0.8092604
## 263A>G 263A>G 263 302 0.8949466
## 8860A>G 8860A>G 8860 327 0.9147965
## 16519T>C 16519T>C 16519 360 0.9366273
## 750A>G 750A>G 750 372 0.9496424
## 11893A>G 11893A>G 11893 349 0.9234525
## 15326A>G 15326A>G 15326 340 0.9353755
## 4769A>G 4769A>G 4769 325 0.8841710
## 3010G>A 3010G>A 3010 374 0.9236634
## 1438A>G 1438A>G 1438 373 0.9203289
## 12519T>C 12519T>C 12519 321 0.9184452
## vmr mean
## 72T>C 1.131044e-01 0.0251883
## 3243A>G 1.042155e-02 0.9518696
## 9448A>G 4.954301e-03 0.9535470
## 310T>C 3.263383e-03 0.7442912
## 16189T>C 7.605553e-05 0.9944897
## 263A>G 5.746567e-05 0.9974504
## 8860A>G 5.126319e-05 0.9960120
## 16519T>C 3.612037e-05 0.9979421
## 750A>G 3.357966e-05 0.9982984
## 11893A>G 3.262667e-05 0.9983553
## 15326A>G 1.953978e-05 0.9977827
## 4769A>G 1.179316e-05 0.9982204
## 3010G>A 9.946766e-06 0.9981163
## 1438A>G 9.791110e-06 0.9992124
## 12519T>C 3.258128e-06 0.9989945Visualize strand correlation vs. VMR
This is the standard visualization described directly in the Lareau protocol paper for distinguishing high- from low-quality variants.
How to read this: real, informative variants typically show both high strand correlation and high VMR, clustering toward the upper right. Noise tends to cluster with low strand correlation, scattered VMR. Note that the gold star (m.3243A>G) falls in the high-quality range, as we expect.
6. Save variants
write.csv(top_variants, "outputs/top_mitochondrial_variants.csv", row.names = FALSE)
message("Full ranked list saved to outputs/top_mitochondrial_variants.csv")## Full ranked list saved to outputs/top_mitochondrial_variants.csv7. Plot top 5 variants (ranked by VMR) and analyze diversion trends of those variants between celltypes
Let’s investigate the top highest-confident variants in more detail.
top5_variants <- head(top_variants$variant, 5)
top5_variants## [1] "72T>C" "3243A>G" "9448A>G" "310T>C" "16189T>C"
Same plot excluding the known variant (m.3243A>G) and the really low, essentially homoplasmic one across celltypes (72T>C):
UMAPs (can be re-generated for any one variant of interest):
8. A key result:
16189T>C looks interesting!
It’s high in RPE and low in Endo, just like the variant they studied
in the paper. We could realistically say we identified this one
additional variant as an RPE /Endo separator in addition to
mt.3243A>G.**
Stats comparison to mt.3243A>G:
## variant n_cells_conf_detected strand_correlation vmr
## 16189T>C 16189T>C 236 0.8092604 7.605553e-05
## 3243A>G 3243A>G 347 0.9415496 1.042155e-02
## mean
## 16189T>C 0.9944897
## 3243A>G 0.9518696Reminder: This variant was identified as a top-5 when variants are ranked by VMR. We could also assess their rankings by the sister QC metric, strand correlation. We’ll do that next.
9. Now plot top 5 variants ranked by Strand Correlation and analyze representation of those variants across celltypes
# Re-rank by strand correlation instead of VMR
top_variants_by_strand <- high_conf %>%
arrange(desc(strand_correlation)) %>%
select(variant, position, n_cells_conf_detected, strand_correlation, vmr, mean)
head(top_variants_by_strand, 6)## variant position n_cells_conf_detected strand_correlation
## 72T>C 72T>C 72 6 0.9971011
## 750A>G 750A>G 750 372 0.9496424
## 3243A>G 3243A>G 3243 347 0.9415496
## 16519T>C 16519T>C 16519 360 0.9366273
## 15326A>G 15326A>G 15326 340 0.9353755
## 3010G>A 3010G>A 3010 374 0.9236634
## vmr mean
## 72T>C 1.131044e-01 0.0251883
## 750A>G 3.357966e-05 0.9982984
## 3243A>G 1.042155e-02 0.9518696
## 16519T>C 3.612037e-05 0.9979421
## 15326A>G 1.953978e-05 0.9977827
## 3010G>A 9.946766e-06 0.9981163top6_variants <- head(top_variants_by_strand$variant, 6)
top6_variants## [1] "72T>C" "750A>G" "3243A>G" "16519T>C" "15326A>G" "3010G>A"Mean heteroplasmy plots by celltype, same pattern, now for 6 variants since we want to get the top 5 AFTER 3243A>G
Same plot excluding the known variant (m.3243A>G) and the really low, essentially homoplasmic one across celltypes (72T>C):
Slightly more detailed view instead of just averages per celltype
Faceted violin/jitter:
Basically, we can see significant overlap in our top variants,
whether we use VMR or strand concordance as
our index for variant quality.
Irrespective of quality, it might be interesting to find out which variants differ most drastically between celltypes in general. Next, we calculate that as a ranking metric, like VMR or strand correlation. Then we can cross-check if any “high-variance” candidates we find overlap with our other top candidates.
10. Rank variants by pairwise difference between any 2 celltypes (to find the ones that are most distinct in their variance)
Here we calculate the MAXIMUM pairwise difference in mean heteroplasmy between any 2 celltypes represented, for each variant:
# 1. Compute allele frequencies for the ENTIRE high-confidence list, not just top N
all_high_conf_variants <- high_conf$variant
mine <- AlleleFreq(mine, variants = all_high_conf_variants, assay = "mito")
DefaultAssay(mine) <- "alleles"
# confirm all variants landed
SeuratObject::Assays(mine)## [1] "RNA" "ATAC" "mito" "integrated" "alleles"length(rownames(mine[["alleles"]]))## [1] 15# 2. Reshape to long format: one row per cell per variant
allele_freq_data_all <- FetchData(mine, vars = c(all_high_conf_variants, "cell_type"))
allele_freq_long_all <- allele_freq_data_all %>%
pivot_longer(cols = all_of(all_high_conf_variants), names_to = "variant", values_to = "allele_freq")
# 3. Per-cell-type mean for each variant
celltype_means <- allele_freq_long_all %>%
group_by(variant, cell_type) %>%
summarise(mean_freq = mean(allele_freq, na.rm = TRUE), .groups = "drop")Top 10 variants by their variance between celltypes:
# 4. The actual metric: max pairwise difference = max(means) - min(means), per variant
variant_max_diff <- celltype_means %>%
group_by(variant) %>%
summarise(
max_diff = max(mean_freq) - min(mean_freq),
celltype_max = cell_type[which.max(mean_freq)],
celltype_min = cell_type[which.min(mean_freq)],
.groups = "drop"
) %>%
arrange(desc(max_diff))
#show all
variant_max_diff## # A tibble: 15 × 4
## variant max_diff celltype_max celltype_min
## <chr> <dbl> <fct> <fct>
## 1 3243A>G 0.537 Fibroblast Endothelial
## 2 9448A>G 0.419 RPE Endothelial
## 3 16189T>C 0.370 Pericyte/SMC Endothelial
## 4 8860A>G 0.344 Pericyte/SMC Endothelial
## 5 12519T>C 0.338 Pericyte/SMC Endothelial
## 6 15326A>G 0.336 Pericyte/SMC Endothelial
## 7 16519T>C 0.327 Pericyte/SMC Endothelial
## 8 750A>G 0.314 Pericyte/SMC Endothelial
## 9 3010G>A 0.313 Pericyte/SMC Endothelial
## 10 11893A>G 0.299 RPE Melanocyte
## 11 4769A>G 0.269 Pericyte/SMC Melanocyte
## 12 1438A>G 0.263 Pericyte/SMC Melanocyte
## 13 310T>C 0.262 Fibroblast Endothelial
## 14 263A>G 0.256 Pericyte/SMC Endothelial
## 15 72T>C 0.0442 RPE EndothelialThere’s our 16189T>C cropping up again! It ranks #3 for variance as well as high on our QC metrics!
Let’s scope out its biological significance.
11. Biological deep-dive of a key identified variant: 16189T>C
We cross-reference variants of interest against known mtDNA disease-variant databases (e.g., MITOMAP) to see if there’s anything known about them that we can connect to our own research interests.
From the MITOMAP entry: this specific variant has 198 associated citations. It’s a known population polymorphism studied almost entirely in association-study contexts. Based on inspection of the citation list, it’s been repeatedly discussed in the context of metabolic and cardiovascular conditions across many ethnic populations.
It could be interesting to say that we also found celltype-specific separation for this more common variant, as well as the more niche one obviously linked to a pathogenic state. Maybe it even fits with Mullin’s overall hypothesis about the heteroplasmy of variants being non-random between celltypes.