MitoTracer facilitates the identification of informative mitochondrial mutations for precise lineage reconstruction

Abstract

Mitochondrial (MT) mutations serve as natural genetic markers for inferring clonal relationships using single cell sequencing data. However, the fundamental challenge of MT mutation-based lineage tracing is automated identification of informative MT mutations. Here, we introduced an open-source computational algorithm called "MitoTracer", which accurately identified clonally informative MT mutations and inferred evolutionary lineage from scRNA-seq or scATAC-seq samples. We benchmarked MitoTracer using the ground-truth experimental lineage sequencing data and demonstrated its superior performance over the existing methods measured by high sensitivity and specificity. MitoTracer is compatible with multiple single cell sequencing platforms. Its application to a cancer evolution dataset revealed the genes related to primary BRAF-inhibitor resistance from scRNA-seq data of BRAF-mutated cancer cells. Overall, our work provided a valuable tool for capturing real informative MT mutations and tracing the lineages among cells.

Fig 1. Overview of MitoTracer algorithm.

(A) The whole analysis process of identifying informative MT mutations. MERCI-mtSNP calls MT mutations from single cell RNA or DNA sequencing data. The VAF matrix is generated for MitoTracer. (B) Informative MT mutation selection. MitoTracer firstly removes mutations caused by sequencing errors, and filters variants with extremely low/high cell- or sample-level frequency. Dirichlet process Gaussian mixture model is conducted on each MT mutation to find out informative MT mutation. We define the informative MT mutation as the absolute difference between the mean in the top two Gaussian distributions larger than the cutoff. MitoTracer uses the VAF matrix of these informative mutations to calculate the similarity matrix based on several distance methods, including the mitochondrial distance defined by ourselves, Euclidean, and correlation.

Fig 2. Performance comparison.

(A-D) Performance comparison on a gold-standard dataset with 15 clones among four methods, including (A) VAF_cutoff, (B) MitoTracer, (C) MQuad, and (D) SClineager. Clone information is labeled by the “Clone” annotation bar on the right side of the heatmap. (E-F) ROC of the manual selection method and the above four methods under between-clone and within-clone. (G-H) PRC of the manual selection method and the above four methods under between-clone and within-clone.

Fig 3. Validation of the MitoTracer algorithm.

(A) Unsupervised hierarchical clustering of the VAF matrix of 24 informative MT mutations showed a clear clustering of BT142 and K562 cells for the MAESTER dataset. (B) Unsupervised hierarchical clustering of the VAF matrix of 47 informative MT mutations showed a clear clustering of GM11906 and TF1 cells for mtscATAC-seq dataset. (C) Unsupervised hierarchical clustering of the VAF matrix of 23 informative MT mutations showed a clear clustering of C9, D6, and G10 hematopoietic cells for the scATAC-seq dataset. (D) Unsupervised hierarchical clustering of the VAF matrix of 31 informative MT mutations showed a clear clustering of cells by their patient origin for the SMART-seq2 dataset. (E) Unsupervised hierarchical clustering of the VAF matrix of 82 informative MT mutations showed a clear clustering of cells by their patient origin for the 10X scRNA-seq dataset.

Fig 4. MitoTracer identified essential primary BRAF inhibitor-resistant genes.

(A)The reconstructed lineage was visualized by heatmap from MitoTracer. All these cells were clustered into two major clusters labeled “Cluster_1” and “Cluster_2”. Cells were labeled according to their original BRAF inhibitor resistance status, “Resistant” and “Parental”. We also labeled the primary resistant cells predicted by MitoTracer with “Resistant Cells”. (B) Graphic description of positive selection for primary resistant-related MT mutation, MT_16389_G-A. (C) Gene ontology biological process enrichment results of 674 differentially expressed genes. (D) GSEA enrichment results of 674 differentially expressed genes. (E) The overall expression level across resistant and MT_16389_G-A status. We defined the group “MT0” which presented the cells with MT_16389_G-A mutation, and MT1 indicated the cells without MT_16389_G-A mutation. (F-G) the expression level of PSMB3, SNRPD2, and UBL5 across resistant and MT_16389_G-A status. The definition of the group was the same as (E). ***P < 0.0001 and **** P < 0.00001.