An Innovative Data Analysis Strategy For Accurate NGS Detection of Tumor mtDNA Mutations

Background: Next generation sequencing (NGS) technology has been commonly applied to detect mtDNA mutations, which are reported to be strongly associated with cancers. However, several key challenges still exist in the bioinformatic analysis of mtDNA sequencing data, which greatly affect the detection accuracy of mtDNA mutations. Methods: Here, we systematically evaluated several key analysis procedures, including trimming, mapping and ltering, in mtDNA mutation detection of FFPE tissues, fresh tissues and plasma samples from cancer patients. Furthermore, the innovative bioinformatics pipeline integrating a newly-developed ltering algorithm was established. Results: We found that trimming procedure was essential for improving mtDNA mapping performance in plasma but not tissue samples. Mapping with rCRS-hg19 reference was strongly suggested for mtDNA mutation detection in plasma samples due to the extreme abundance of NUMTs. In addition, our results showed that the setting of 3 mismatches was most appropriate for mtDNA mutation calling. More importantly, we revealed the presence of a negative logarithmic relationship between mtDNA site sequencing depth and minimum detectable mutation frequency and thus build up an innovative and ecient ltering strategy to increase the accuracy and sensitivity of mutation detection. Finally, we veried that higher sequencing depth was required for PCR-based than capture-based enrichment strategy. Conclusions: Collectively, we established an innovative data analysis strategy, which is of great signicance for improving the accuracy of mtDNA mutation detection for different types of tumor samples.

importantly, we revealed the presence of a negative logarithmic relationship between mtDNA site sequencing depth and minimum detectable mutation frequency and thus build up an innovative and e cient ltering strategy to increase the accuracy and sensitivity of mutation detection. Finally, we veri ed that higher sequencing depth was required for PCR-based than capture-based enrichment strategy.
Conclusions: Collectively, we established an innovative data analysis strategy, which is of great signi cance for improving the accuracy of mtDNA mutation detection for different types of tumor samples.

Background
Human mitochondria possess their own genome, which is a double-stranded, maternally inherited, circular DNA of 16,569 base pairs (bp), with up to 10 3 -10 4 copies in each cell (1). Mitochondrial DNA (mtDNA) exists with a plenty of sequence variants, which are commonly observed either as germline or somatic mutations (2). Numerous studies have demonstrated that heteroplasmy, termed as the presence of a mixed mtDNA mutation genotypes within a cell, is strongly associated with lots of disease phenotypes, especially when the percentage of mutations exceeds a critical threshold (3). Thus, it is of great importance to accurately detect mtDNA mutations for better understanding mitochondrial biology and cancers.
Traditional techniques, such as Sanger sequencing and high-resolution melt analysis, have been used for mtDNA mutation detection, which still face the disadvantages of low throughput and sensitivity.
Fortunately, the advent of next-generation sequencing (NGS) technologies provides an opportunity for high throughput and low-cost detection of mitochondrial genome-wide mutations (4). Two major methods have been used for the NGS-based detection of mtDNA mutations (5). One is direct data extraction from whole genome sequencing (WGS) or whole exome sequencing (WES) for mutation analysis, which is low-throughput and not cost-effective (6). Another is to rst enrich mtDNA from total genomic DNA and then detect by NGS, mainly including PCR-based and capture-based enrichment strategies (7). They can not only achieve the detection of low minor allele frequency (MAF), but also permit customized sequencing strategies, allowing segment-speci ed or whole-genome sequencing of mtDNA. Combined with NGS technology, such strategies have been proved to be powerful and e cient when depicting spectrum of mtDNA somatic mutation in many types of cancers (3).
Although NGS has been extensively applied to effectively detect mtDNA mutations, there are still several big challenges existing in data analysis of mtDNA sequencing (8). Among them, most critical consideration is to decrease the number of false-positive and false-negative mutations, which are greatly affected by mtDNA sequencing depth, especially for those with low heteroplasmy level (9). To date, the quantitative relationship between sequencing depth and the detection accuracy of mtDNA mutations remains to be determined. Therefore, most of previous studies arbitrarily selected 5%, 2% or 1% as the minimal heteroplasmy level to lter mtDNA mutations (10)(11)(12). In addition, owing to the presence of homologous sequences within the nuclear genome, referred as "nuclear sequences of mitochondrial origin" (NUMTs), suitable mtDNA mapping strategy should be carefully assessed to reduce the effect of NUMTs in different sample types (13). Considering the great decrease of NGS cost in recent years, high sequencing depth (usually more than thousands) is becoming practical for tumor mtDNA mutation detection. It is urgently needed to develop an innovative data analysis strategy for more accurate NGS detection of mtDNA mutations, especially for those with low heteroplasmy levels in cancer cells.
In the present study, mainly based on capture-based NGS data, we systematically evaluated several key analysis procedures, including trimming, mapping and ltering, in mtDNA mutation detection of FFPE tissues, fresh tissues and plasma samples from cancer patients. Finally, an innovative bioinformatics pipeline integrating a newly-developed ltering algorithm was established. Furthermore, the application of our analysis strategy in different sample types has also been evaluated, providing a exible selection of key analysis procedures.

Sample collection
A total of 50 FFPE (Formalin-Fixed and Para n-Embedded) tissue samples, 50 fresh tissue samples and 50 plasma samples were collected from liver cancer patients in Xijing Hospital, Fourth Military Medical University (FMMU) in Xi'an, China. This study was approved by the Ethical Committees of FMMU and written consent was obtained from each patient.

DNA extraction
Genomic DNA was extracted from fresh tissue, FFPE tissue and plasma using ENZA Tissue DNA Kit (Omega), QIAamp DNA FFPE Kit (Qiagen) and QIAamp Circulating Nucleic Acid Kit (Qiagen) according to manufacturer's protocols, respectively. DNA quality and concentration were assessed using 2100 Bioanalyzer (Agilent) and Qubit (Invitrogen).

Library construction and mtDNA enrichment
Whole genome sequencing (WGS) library for Illumina platform was constructed as previously described (14). In brief, 1μg genomic DNA from FFPE and fresh tissues were randomly sonicated by focusedultrasonicator (Scientz98, Ningbo, China) to obtain fragments mainly distributed between 300 and 500 bp in length. DNA fragments were end-repaired, ligated with sequencing adapters and slightly PCR-ampli ed (9 cycles). For plasma samples, 20 ng genomic DNA were used to construct the sequencing library using NEB ultra v2 kit (NEB). Then, WGS libraries were mixed with home-made biotinylated mtDNA capture probes for hybridization. Furthermore, to examine the effect of different enrichment strategy on mtDNA mutation detection, we also constructed PCR-based mtDNA enrichment library using QIAseq Targeted DNA Human Mitochondrial Panel (Qiagen) for 6 plasma samples following the manufacture's protocol.

Sequencing platforms and next generation sequencing
The capture-and PCR-based mtDNA libraries were sequenced on Illumina XTen platform using paired-end runs (2 x 150 cycles). To further evaluate the suitability of the optimized mtDNA bioinformatics pipeline in different sequencing platforms, 10 capture-based mtDNA libraries from fresh tissue were also sequenced on MGISEQ-2000 platform (BGI) using paired-end runs (2 x 100 cycles). A summarization of the mtDNA sequencing data used in this study was shown in table S1.

Bioinformatics pipeline for mtDNA mutation calling
We systematically evaluated the analysis pipeline for mtDNA deep-sequencing data. Brie y, raw mtDNA sequencing data rst encounter two options, trimming or without trimming for quality control. The mtDNA reads were then mapped to either rCRS (revised Cambridge Reference Sequence) or combined rCRS and human genome reference (hg19) using BWA software. After sorting and removing duplicated reads with Picard, the Genome Analysis Toolkit 4 (GATK4) was used for local realignment. Finally, we applied a series of ltering conditions (removing false positive mutations) to detect mtDNA mutations and analyze heteroplasmy levels.

Trimming procedure
The FASTQ preprocessor fastp (version 0.20.0) (15) was used for trimming of mtDNA sequencing data with three parameters: rst, all sequencing adaptors were initially removed; second, a sliding window (4 bp in length) approach was used to scan reads from front (5') to tail (3'). If average base quality in the window was below Q30, these bases and downstream part were dropped; third, the reads with length below 50 bp were discarded to avoid ambiguous mapping of short reads.

Mapping strategies and mismatch selection
Two mapping strategies were compared: the rst strategy was to only map sequencing reads to rCRS reference (rCRS-alone); the second strategy was to map sequencing reads to both rCRS and hg19 reference (combined rCRS-hg19), but keep only the reads uniquely mapped to rCRS. In addition, sequencing reads with mapping quality below Q20 were removed from subsequent analysis. Different mismatch lter ranging from 1-4 were evaluated in mtDNA mutation calling. We used two error-evaluating parameters (assumed false positive and assumed false negative mutants) to identify the optimum mismatch number selection. The assumed false positive (AFP) mutations were de ned as those only detected in one mismatch group. The assumed false negative (AFN) mutations were de ned as those absent in this group but presented in at least two other mismatch groups.
Identi cation of minimum detectable mutation frequency based on site sequencing depth To identify the minimum detectable mutation threshold, genomic DNA from different sample types were used for independent library construction and sequencing. These repeated sequencing datasets enabled us to experimentally analyze the consistency of mtDNA mutation calling. Here, consistency between two repeated experiments was de ned as: Where A is the collection of mtDNA mutations in experiment 1, B is the collection of mtDNA mutations in experiment 2. Setting the consistency level at 95%, a dynamic minimum detectable mutation threshold was identi ed for sites with different sequencing depth. Moreover, the relationship between site sequencing depth and the minimum detectable mutation threshold were explored by logistic regression analysis using R software.

Development of novel ltering criteria for mtDNA mutation calling
For each site, we rst counted the respective reads number of the major and minor allele and calculated site-speci c minor allele frequency (MAF). Then, mtDNA mutations were initially called using the default settings as we previously described (14,16,17): 1) at least 3 reads on each strand have the mutation site; 2) minimum MAF cutoff 1%; 3) remove heterogeneity sites in rCRS repeat regions (66-71, 303-311, 514-523, 12418-12425, 16184-16193).
In addition to the default settings, we further explored the impact of three extra ltering strategies on mtDNA mutation calling: lter 1 removes C > A/G > T mutations with low MAF (1%) and strong sequence context bias (at CpCpN>CpApN, most frequently CpCpG>CpApG), which is known to arise from arti cial guanine oxidation during sequencing library preparation (18,19); lter 2 removes mtDNA mutations if the mutant-rate and mutant base-quality does not pass binominal test (P > 0.0001) (20,21); lter 3 removes mtDNA mutations if the MAF is smaller than the site-speci c threshold determined by site sequencing depth.

Statistical analysis
Graphpad Prism 7.0 (GraphPad Software, USA) was used for statistical analysis. Mann-Whitney U test was used to compare the difference between two groups. All P values were two-tailed and reported using a signi cance level of 0.05.

Results
Effect of trimming procedure on mapped mtDNA sequencing data In NGS data analysis, trimming procedure is commonly used to remove the adaptor sequences and bases with low quality from sequencing reads. Therefore, we rst evaluated the effect of trimming procedure on mtDNA mapping of sequencing reads in three types of samples, including FFPE tissues, fresh tissues and plasma. As shown in Fig.1A and B, the trimming procedure had no signi cant in uence on both the mtDNA mapping rate of sequencing reads and average sequencing depth in FFPE and fresh tissue samples. In contrast, the mtDNA mapping rate and average sequencing depth were signi cantly increased in plasma samples after trimming ( Fig. 1A and B). Furthermore, our data indicated the very consistent size distribution of mtDNA insert fragments in FFPE and fresh tissue samples (Fig. 1C). However, the mapped mtDNA sequencing reads were signi cantly increased in plasma samples after trimming (Fig. 1C). These results indicate that trimming procedure is essential for improving the mtDNA mappability in plasma but not tissue samples.
Evaluation of mtDNA mapping strategy in mtDNA mutation calling Next, we evaluated the application of two commonly used mapping strategies in three different sample types. As shown in Fig. 2A, when compared with the rst mapping strategy (rCRS-alone), the second mapping strategy (combined rCRS-hg19) clearly removed the sequencing reads mapped to both hg19 and rCRS, leading to 11 more remarkable peaks. Among them, the most signi cant peak occurred around mtDNA site 8500, with 45% -80% removed reads in three different sample types. Compared to tissue samples, plasma DNA exhibited a clearly different distribution pattern of removed reads. Then, the mtDNA mutations detected by the two mapping strategies were depicted in Fig. 2B, with mutation site and heteroplasmy level (MAF≥1%). The venn diagram showed that all mtDNA mutations detected by the second mapping strategy were included in those detected by the rst one, while 1, 5 and 1256 of mtDNA mutations only detected by mapping to rCRS were observed in FFPE tissue, fresh tissue and plasma samples, respectively (Fig. 2C). The site and heteroplasmy level of mtDNA mutations only detected by mapping to rCRS were shown in Fig. 2D. To explore the potential reason explaining extremely high number of mtDNA mutations only detected by mapping to rCRS in plasma samples, we further analyzed the mtDNA copy number and percentage of mtDNA sequencing reads mapped to both hg19 and rCRS in three sample types. Our results showed that the mtDNA copy number in plasma was signi cantly lower than tissues (Fig. S1), and the percentage of mtDNA sequencing reads mapped to both hg19 and rCRS was signi cantly higher in plasma (17.18%) compared to FFPE and fresh tissues (3.46% and 4.03%, respectively) (Fig. S2). To investigate the NUMTs-derived possibility of those mutations, we determined the percentage of mutant reads mapped to rCRS only or both hg19 and rCRS. As shown in Fig. 2E, 93.33%, 92.86% and 91.91% of the mutant reads were mapped to both hg19 and rCRS in three sample types, respectively, while 6.67%, 7.14% and 9.09% were only mapped to rCRS. Further analysis showed that 82.57%, 91.50% and 92.72% of the mutant reads had higher alignment scores when mapped to hg19 than rCRS in three sample types, respectively (Fig. 2F). These data suggest that these mutations may be introduced by NUMTs. Collectively, when mtDNA mutations are detected in tissue samples, both the two mapping strategies can be used, although rare mutations may be introduced by NUMTs. In comparison, the second mapping strategy is strongly suggested for mtDNA mutation detection in plasma samples due to the extreme abundance of NUMTs.
Effect of mismatch number selection on mtDNA mutation calling Then, we investigated the effect of mismatch number selection in the second mapping strategy on mtDNA mutation calling by using the repeated sequencing data in three sample types. As shown in Fig.  3A, the total repeated mtDNA mutation numbers in two repeat experiments of 10 samples gradually increased when the mismatch number changed from 1 to 4 in three sample types. We then found that the average number of repeatable mutations in the group with 3 or 4 mismatches was signi cantly higher than that in the group with 1 mismatch in all three sample types. The average number of repeatable mutations varied from 35.5-42.6, 32.3-36.6, and 37.7-49.2 under different mismatches in three sample types, respectively. Furthermore, the assumed false positive (AFP) mutations were de ned as those only detected in one mismatch group. The assumed false negative (AFN) mutations were de ned as those not detected in this group but detected in at least two other mismatch groups. We found no signi cant difference of AFP mutation number among the groups with 1, 2, or 3 mismatches in three sample types, while it was signi cantly increased in the group with 4 mismatches when compared to the groups with 1, 2, or 3 mismatches (Fig. 3B). In contrast, the AFN mutation number in groups with 1 or 2 mismatches were signi cantly higher than that in the groups with 3 or 4 mismatches among three sample types (Fig.  3C). Collectively, these results demonstrate the impacts to some extent of mismatch number selection on mtDNA mutation detection, and strongly suggests the setting of 3 mismatches for mtDNA mutation calling.
Relationship between the site sequencing depth and consistency of mtDNA mutations at different heteroplasmy levels.
Next, we comprehensively investigated the effect of mtDNA site sequencing depth on the detection accuracy of mutations at low heteroplasmy levels by using the repeated sequencing data in three sample types. As shown in Fig. 4A, a very consistent trend was observed in all three sample types, indicating a faster rise of the consistency when the site sequencing depth was relatively low. To achieve the consistency of 90%, the site sequencing depth of greater than 2700X, 2300X, and 3200X was needed in FFPE tissues, fresh tissues and plasma, respectively, when the heteroplasmy level was set above 0.5%. In comparison, the consistency of mtDNA mutations was notably decreased in all three sample types when the heteroplasmy level was set at 0.1%, suggesting that the system used in this study may not be suitable for mtDNA mutation detection at the heteroplasmy level of 0.1%. Furthermore, the negative logarithmic relationship was established between mtDNA site sequencing depth and minimum detectable mutation frequency when the consistency was set at 90% (Fig. 4B). When the minimum detectable mutation frequency was 1% and 0.5%, the site sequencing depth was required to be higher than 1000X and 4000X, 1200X and 3600X, 1700X and 4700X in FFPE tissue, fresh tissue and plasma, respectively.
Commonly, the average sequencing depth was calculated and provided based on mtDNA sequencing data. Considering the accuracy of mtDNA mutation calling was actually affected by site sequencing depth, we thus determined the relationship between them. As shown in Fig. 4C, with the increasing average sequencing depth, the percentage of mtDNA sites with site depth greater than average sequencing depth were gradually increased. When the average sequencing depth was greater than 8000X, there were 80% of sites with site depth higher than the average.

Assessment of different mtDNA mutation ltering strategies
To reduce false positive mutations, several ltering strategies has commonly been applied (18). First, we assessed the effect of two ltering strategies previously reported on mtDNA mutation calling. One ltering strategy ( lter 1) is to remove C > A/G > T mutations with low MAF (1%) and strong sequence context bias (at CpCpN>CpApN; most frequently at CpCpG>CpApG), which is known to arise from arti cial guanine oxidation during sequencing library preparation steps (18,19). As shown in Fig. 5A, no difference of mtDNA mutation numbers was observed between two groups ltered with and without lter 1 in FFPE tissue, fresh tissue and plasma. Furthermore, second ltering strategy ( lter 2) is to remove mtDNA mutations where the quality of mutant bases does not t well with binominal distribution (P > 0.0001) (20,21). Those may be false positive mutations introduced by sequencing errors. Very similarly with lter 1, lter 2 exhibited no signi cant effect on mtDNA mutation numbers in all three sample types (Fig. 5B). Considering the great in uence of site sequencing depth on mtDNA mutation detection, we built up a novel ltering strategy ( lter 3) based on negative logarithmic functions presented in Fig. 4B to remove mtDNA mutations with heteroplasmy level lower than minimum detectable mutation frequency. As shown in Fig 5C, the number of mtDNA mutations was signi cantly decreased after ltering in all three sample types. Moreover, the repeated mtDNA sequencing data were used to analyze whether those ltered mutations were repeatable or not. We found that 91.84%, 89.4%, and 93.3% of these ltered mutations were unrepeatable in FFPE tissue, fresh tissue and plasma, respectively. In comparison, only 9.4%, 5.7% and 15% of the unaffected mtDNA mutations was unrepeatable (Fig. S3). These ndings indicate that both lter 1 and lter 2 are not necessary in our analysis pipeline, whereas lter 3 greatly contribute to improve the detection accuracy.
Comparison between PCR-based and capture-based enrichment strategies or two sequencing platforms To reduce sequencing cost, both capture-based and PCR-based strategies are commonly used for mtDNA enrichment from total genomic DNA. Therefore, we systematically compared the performance of two enrichment strategies in NGS-based mtDNA mutation detection of plasma samples. We found that capture-based mtDNA sequencing exhibited a more uniform distribution of the coverage across the whole mitochondrial genome and the coe cient of variation (CV) signi cantly decreased when compared to PCR-based approach (Fig. 6A). Then, six plasma DNA samples were sequenced three times, twice by capture-based approach and one time by PCR-based approach. The number of mtDNA mutation sites with more than 1% heteroplasmy level was depicted in Fig. 6B. Moreover, we found that the consistency of the mtDNA mutations between two capture-based sequencing data was signi cantly higher than those between each capture-based sequencing data and PCR-based sequencing data (Fig. 6C). With the increasing of mtDNA site sequencing depth, the mtDNA mutation consistency between capture-based and PCR-based sequencing data was gradually elevated (Fig. 6D). To get more accurate mtDNA mutations, our results suggest higher sequencing depth for PCR-based sequencing approach when compared with capture-based sequencing.
Additionally, to test the robust application of our innovative bioinformatics pipeline, we evaluated the consistence of mtDNA mutation pro ling between Illumina and BGI sequencing platforms, which were most widely used for next-generation sequencing. The mtDNA sequencing data of 10 fresh tissue samples from two platforms were analyzed. As shown in Fig. 6E, all samples exhibited the good consistency of mtDNA mutation numbers. We further compared the heteroplasmy level of mtDNA mutations. Our results revealed a remarkable correlation between two sequencing platforms (r = 0.9974, P < 0.01), indicating the widespread applicability of our innovative data analysis pipeline.

Discussion
Here, based on systematical evaluation of key analysis procedures, we established an innovative data analysis strategy for improving the accuracy of NGS-based mtDNA mutation detection, which is mainly integrated with trimming procedure, mapping strategy with rCRS and hg19 as reference genomes and a new-developed ltering approach. Therefore, we for the rst time reported several key ndings about the application of data analysis pipeline in three different sample types. First, we demonstrated that trimming procedure is essential for improving mtDNA mapping performance in plasma but not tumor tissue samples. Second, our systematic analysis of mtDNA reference genomes clearly showed the great impact of NUMTs in plasma samples and provided the optimal choice for reference genomes in different sample types. Third, we demonstrated that the setting of 3 mismatches is most suitable for mtDNA mutation calling. More importantly, we found the negative logarithmic relationship between mtDNA site sequencing depth and minimum detectable mutation frequency and thus innovatively built up an e cient ltering strategy to increase accuracy of mutation detection. All these efforts greatly contribute to the establishment of an innovative and versatile bioinformatics pipeline for the accuracy of mtDNA mutation detection, laying a foundation for translational application of mitochondrial mutations in clinical practices.
Quality control and preprocessing of sequencing data are critical to obtain highly accurate mtDNA mutations in downstream data analysis, especially important for detecting low-MAF mutations. The trimming procedure, which refers to eliminate adaptors and poor-quality bases of the sequencing reads, is an initial step of the analysis that is heterogeneously applied in previous mtDNA mutation analyses (22).
Whether the trimming procedure is indispensable to mtDNA mutation analysis remains to be con rmed, especially for different sample types. In the present study, we demonstrated that trimming procedure is essential for improving mtDNA mapping performance in plasma but not tumor tissue samples.
Since the true origin of homologous reads are di cult to discriminate, NUMTs and mtDNA cross-mapping occurs, leading to the detection of false positive (NUMTs reads aligning to chrM) or false negative (mtDNA reads aligning to NUMTs loci) mtDNA mutations. To address this concerning source of error, Li et al. have previously created a database of NUMTs containing mismatches from the mitochondrial genome, which may appear as false heteroplasmies if aligned to chrM (23). Complementing this approach, two main mapping strategies have been commonly used at present. The rst is mapping to the human reference mtDNA sequence rCRS, and the second is mapping to rCRS and hg19 (human genome 19) that contributes to remove NUMTs during analysis. Therefore, the accuracy of mtDNA mutation calling can be largely affected by mapping strategies with different selection of reference genomes and mismatch number, owing to existence of sequence similarity between NUMTs and mtDNA (a known source of confounding in mtDNA NGS studies) (24). However, there are few studies focusing on the applicability of the two mapping strategies under different circumstances. In the present study, we performed a systematic performance comparison of mtDNA mutation calling pipelines with different mapping strategies (rCRS-alone or combined rCRS-hg19) or mismatch settings (changed from 1 to 4) in three different sample types. Our data indicate that both the two mapping strategies can be effectively used in both FFPE and fresh tissue samples, although only mapping to rCRS may introduce a very small amount of NUMTs. However, the mapping strategy with both hg19 and rCRS as reference genomes is strongly suggested for mtDNA mutation detection in plasma samples due to the extreme abundance of NUMTs. Furthermore, our study strongly suggests the setting of 3 mismatches for mtDNA mutation calling in all three sample types.
The detectable level of mtDNA heteroplasmy is heavily dependent on the depth of NGS coverage. Furthermore, with the sequencing depth becoming deeper, the accuracy of mtDNA mutation detection is increasing. Recent studies have detected the somatic mtDNA mutations at a very low heteroplasmic level (MAF > 1%), with the NGS depth for the mitochondrial genome 9959 × (18). However, as we lower the detectable threshold of heteroplasmy, it becomes increasingly di cult to distinguish between ultra-low MAF mutations and sequencing errors. Thus, addressing the technical question of correctly standardizing the sequencing depth in order to obtain con dent and reproducible detections of low MAF mutations is of great importance. However, which sequencing depth can provide su cient information to detect lowfrequency mtDNA mutations remains to be investigated up to now. Since there are no consensus criteria for determining mtDNA heteroplasmy threshold, it is very common to arbitrarily select a constant heteroplasmy level (mainly 2% or 5%) in previous NGS analyses for mtDNA mutations (12,17). Here, we innovatively revealed the relationship of site sequencing depth and the minimum mutation threshold via establishing logarithmic functions in different sample types, which contributes to the discrimination of false-positive and false-negative mutations with ultra-low heteroplasmy level (up to 0.5%). By providing a standardized option of proper mutation frequency threshold, it can largely increase the sensitivity and accuracy of mtDNA mutation detection during NGS data analysis.
One unique aspect of our study is the integrative analysis of mtDNA mutation calling based on the assessment of key procedures. We found that the ltering strategy based on site sequencing depth greatly contributes to improve the detection consistency of repeated data. Previous studies on mtDNA mutation calling have also applied several ltering conditions, such as at least 3 reads per mutation allele and strand bias correction (3), to reduce false positive mtDNA mutations. Besides, several studies have also focused on developing lter-based methods to remove oxidation-mediated mutations during DNA shearing or artifacts arising from Illumina sequence errors (19,21). However, those ltering strategies may not be necessary for the mutation analysis in ultra-deep sequencing (with sequencing depth is deep enough). For example, the minimum number of mutant reads must be greater than 3 as long as the site depth is greater than 600X when detecting mutations at the threshold greater than 1%. Moreover, we also observed that the ltering strategies considering DNA oxidation and binominal distribution exhibited no signi cant impact in our analysis system, which may at least partially be due to the strict removing of low-quality bases during the trimming procedure. In addition, the negative logarithmic function established in our study is based on capture-based sequencing data, whose applicability to all sequencing data may indeterminate, but can be applied for providing a reference for the range of sequencing depth and heteroplasmy level in NGS data.
Both capture-based and PCR-based methods are commonly used for mtDNA enrichment before sequencing (16). Therefore, selecting the right enriching approach and sequencing platform is of great importance to accurately identify mtDNA mutations, especially those with ultra-low heteroplasmy level. A previous study has compared the detection performance of different enrichment strategies in fresh tissue samples and illustrated that DNA quality was a great challenge for PCR-based approach, which would not work with highly fragmented DNA samples such as FFPE tissues and plasma (25). In the present study, we explored different enriching approaches and sequencing platforms in three different sample types. Our results suggest higher sequencing depth for PCR-based sequencing approach when compared with capture-based sequencing, and the innovative data analysis pipeline is applicable for both Illumina and BGI sequencing platforms in the detection of mtDNA mutations.

Conclusions
Taken together, based on systematic evaluation of key analysis procedures, we established an innovative data analysis pipeline for different tumor sample types, which is of great signi cance for improving the accuracy of mtDNA mutation detection. These efforts lay a foundation for broader biomedical applicability for accurate investigation of mitochondrial genome in cancer cells. Effect of trimming procedure on mapped mtDNA sequencing data A. Comparison of mtDNA mapping rate between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. B. Comparison of mtDNA average sequencing depth between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. C. Distribution of mtDNA insert size between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. ns, no signi cance; *P < 0.05.

Figure 1
Effect of trimming procedure on mapped mtDNA sequencing data A. Comparison of mtDNA mapping rate between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. B. Comparison of mtDNA average sequencing depth between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. C. Distribution of mtDNA insert size between the trimming group and without trimming group in FFPE tissue, fresh tissue, and plasma samples. ns, no signi cance; *P < 0.05.       Assessment of different mtDNA mutation ltering strategies A. Comparison of mtDNA mutation calling with or without lter 1 in FFPE tissue, fresh tissue and plasma samples. Filter 1 is to remove C > A/G > T mutations with low MAF (1%) which is known to arise from arti cial guanine oxidation during sequencing library preparation. B. Comparison of mtDNA mutation calling with or without lter 2 in FFPE tissue, fresh tissue and plasma samples. Filter 2 is to remove mtDNA mutations where mutation-rate and mutation base-quality does not pass binominal test (P > 0.0001). C. Comparison of mtDNA mutation calling with or without lter 3 in FFPE tissue, fresh tissue and plasma samples. Filter 3 is the novel lter that removes mtDNA mutations if the MAF is smaller than the site-speci c threshold determined by site sequencing depth. D. Percentage of the repeatable and unrepeatable mtDNA mutations that were removed by lter 3 in FFPE tissue, fresh tissue and plasma samples. ns, no signi cance; **P < 0.01.

Figure 5
Assessment of different mtDNA mutation ltering strategies A. Comparison of mtDNA mutation calling with or without lter 1 in FFPE tissue, fresh tissue and plasma samples. Filter 1 is to remove C > A/G > T mutations with low MAF (1%) which is known to arise from arti cial guanine oxidation during sequencing library preparation. B. Comparison of mtDNA mutation calling with or without lter 2 in FFPE tissue, fresh tissue and plasma samples. Filter 2 is to remove mtDNA mutations where mutation-rate and mutation base-quality does not pass binominal test (P > 0.0001). C. Comparison of mtDNA mutation calling with or without lter 3 in FFPE tissue, fresh tissue and plasma samples. Filter 3 is the novel lter that removes mtDNA mutations if the MAF is smaller than the site-speci c threshold determined by site sequencing depth. D. Percentage of the repeatable and unrepeatable mtDNA mutations that were removed by lter 3 in FFPE tissue, fresh tissue and plasma samples. ns, no signi cance; **P < 0.01.