Mitochondrial DNA quality control

Abstract

The present disclosure provides a method for identifying unreliable biological samples that may be mislabeled or contaminated by determining the heteroplasmy and homoplasmy, and therefore the mtDNA haplogroup, of the mitochondrial DNA (mtDNA) present in the biological sample.

Description

[Technical field]

[0001] The present disclosure is directed, in part, to methods for identifying unreliable biological samples that may be mislabeled or contaminated. [Background technology]

[0002] For over a decade, next-generation sequencing (NGS) has become an important component of biological and biomedical research because it allows for the sequencing of large batches of DNA or RNA samples. NGS has a wide range of applications, such as whole-genome and whole-exome sequencing for large-scale cohort genetic surveys, bulk RNA-seq for disease gene expression signature identification in clinical evaluation, tissue biopsy sequencing in tumor research / diagnosis, and the recently emerged single-cell sequencing studies, providing answers and solutions to many different challenges and questions. However, in studies involving large samples, sample identity complications are a common and almost unavoidable challenge. Estimates of sample identification error rates can range from 0.2% to 6% in practice (Non-Patent Document 1, Non-Patent Document 2, Non-Patent Document 3, and Non-Patent Document 4). Errors can occur to various degrees: 1) complete substitution between samples, and / or 2) contamination of one sample with one or more other samples. Errors can occur at various steps during sample processing, such as mislabeling of samples during sample collection, spillage of material during pipetting, index swapping in pooled libraries when performing sequencing, and many other unexpected situations. Sample swapping / contamination will subsequently reduce the quality and accuracy of downstream analyses. For example, sample swapping in whole transcriptome analysis may lead to false discoveries or loss of ability to detect differentially expressed genes. In cancer research, the identification of somatic mutations is routinely used given that many of them were present at very low frequencies (<5%), and therefore even low levels of contamination (1%-5%) may result in false positive mutation calls. For those reasons, accurate detection of sample swapping and contamination is a critical quality control step in large-scale NGS studies.

[0003] Mitochondria are essential organelles in most eukaryotic cells. Human mitochondrial DNA (mtDNA) is a 16.5 kb circular DNA molecule located within mitochondria that encodes gene products essential for mitochondrial function. There are hundreds to thousands of mtDNA copies in a single cell. mtDNA is maternally inherited, with negligible recombination. Because mtDNA is uniparentally inherited, with negligible recombination at the population level, mutations acquired over time have led to the subdivision of human populations into several distinct mtDNA haplogroups. On average, two random individuals will have 30-40 nucleotide differences in their mitochondrial genomes (Non-Patent Document 5, Non-Patent Document 6, and Non-Patent Document 7). Due to their multicopy nature, mtDNA mutations are often present in only a small fraction of a cell's mtDNA (a condition called heteroplasmy). The percentage of mtDNA carrying a mutation is called the heteroplasmy frequency. In contrast, if a mutation is found in all mtDNA molecules, the mutation is called homoplasmy. Previous studies have demonstrated that in generally healthy populations, most individuals carry less than five heteroplasmies (frequency >1-2%) in the mitochondrial genome (Non-Patent Documents 8 and 9). For a batch of samples, all samples collected from the same individual should belong to the same haplogroup. [Prior art documents] [Non-patent literature]

[0004] [Non-Patent Document 1] Pfeifer et al.,Amer.J.Clin.Pathol.,2013,139,93-100 [Non-Patent Document 2] Costello et al.,BMC Genomics,2018,19,332 [Non-Patent Document 3] Lerner et al.,Cancer Res.,2015,75,Abstract P5-02-08 [Non-Patent Document 4] Sehn et al.,Amer.J.Clin.Pathol.,2015.144,667-674 [Non-Patent Document 5] Gunnarsdottir et al.,Nature Commun.,2011,2,228 [Non-Patent Document 6] Slatkin et al.,Genetics,1991,129,555-562 [Non-Patent Document 7] Ye et al.,Proc.Nat'l Acad.Sci.USA,2014,111,E4548-E4550 [Non-Patent Document 8] Zhang et al.,BMC Genomics,2017,18,890 [Non-Patent Document 9] Ye et al.,Proc.Nat'l Acad.Sci.USA,2014,111,10654-10659 Summary of the Invention [Means for solving the problem]

[0005] The present disclosure provides a method for identifying unreliable biological samples, the method including: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) sequencing reads for each biological sample; b) identifying heteroplasmy and homoplasmy in the mtDNA sequencing reads from the previous step for each of the biological samples; and c) assigning a primary mtDNA haplogroup to each biological sample, wherein any biological sample having an assigned primary mtDNA haplogroup that differs from the primary mtDNA haplogroup assigned to a majority of biological samples from the same individual is an unreliable biological sample that is a mislabeled biological sample.

[0006] The present disclosure also provides a method for identifying unreliable biological samples, the method including: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) raw sequencing reads for each biological sample; b) identifying heteroplasmy and homoplasmy in the mtDNA sequencing reads from the previous step for each of the biological samples; c) assigning a primary mtDNA haplogroup to each biological sample; and d) determining a total heteroplasmy count for each biological sample, wherein if the heteroplasmy count of the biological sample is high, the biological sample is assigned a secondary mtDNA haplogroup based on minor alleles at the heteroplasmy sites, and determining that a biological sample having an assigned secondary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroup is a contaminated unreliable sample.

[0007] The disclosure also provides a method for identifying an unreliable biological sample, the method comprising: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) raw sequencing reads for each biological sample; b) processing the mtDNA raw sequencing reads for quality control and adapter sequence removal to generate quality controlled mtDNA sequencing reads; c) mapping the quality controlled mtDNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; and d) determining whether or not the candidate mtDNA sequencing read is a candidate mtDNA sequencing read if i) the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome, and ii) the candidate mtDNA sequencing read is a candidate mtDNA sequencing read. The method includes: remapping candidate mtDNA sequencing reads to a human reference genome and retaining the candidate mtDNA sequencing reads if the alignment mismatch count of the sequencing reads is less than 5; e) performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and removing duplicates; f) identifying heteroplasmy and homoplasmy in the retained candidate mtDNA sequencing reads for each of the biological samples; and g) assigning a primary mtDNA haplogroup to each biological sample, wherein any biological sample having an assigned primary mtDNA haplogroup that is different from the primary mtDNA haplogroup assigned to the majority of biological samples from the same individual is an unreliable biological sample that is a mislabeled biological sample.

[0008] The present disclosure also provides a method for identifying an unreliable biological sample, the method comprising: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) raw sequencing reads for each biological sample; b) processing the mtDNA raw sequencing reads for quality control and adapter sequence removal to generate quality controlled mtDNA sequencing reads; c) mapping the quality controlled mtDNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; and d) determining whether or not the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome, and ii) determining whether or not the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome. e) performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and removing duplicates; f) identifying heteroplasmy and homoplasmy in the retained candidate mtDNA sequencing reads for each of the biological samples; g) assigning a primary mtDNA haplogroup to each biological sample; and h) determining a total heteroplasmy count for each biological sample, wherein if the heteroplasmy count of a biological sample is high, the biological sample is assigned a secondary mtDNA haplogroup and determining that a biological sample having an assigned secondary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroup is a contaminated unreliable sample.

[0009] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the U.S. Patent and Trademark Office upon request and payment of the necessary fee. [Brief description of the drawings]

[0010] [Figure 1] A representative schematic showing suitable steps for performing the quality control analysis described herein is shown. In box 1, mtDNA homoplasmy and heteroplasmy are identified from fastq files. An optional downsampling step can be applied for samples with high mtDNA coverage. After QC of the reads, mtDNA reads are selected by a two-stage mapping strategy. mtDNA variants are identified from the mtDNA mapping results, and primary and secondary mtDNA haplogroups are assigned to each sample based on the variant information. In box 2, sample swapping / mislabeling can be detected by comparing the haplogroup assignments of samples from a given individual. In box 3, sample contamination can be detected by abnormally high mtDNA heteroplasmy counts and unmatched primary and secondary haplogroups. [Figure 2A] The performance of the method described herein on virtual contaminated samples is shown. Virtual contaminated samples were created by mixing two samples from the 1000 Genomes Project in different ratios. The X-axis shows the theoretical contamination level, and the Y-axis shows the heteroplasmy frequency determined from each virtual contaminated sample. Each colored dot represents one heteroplasmy, the black dot represents the average heteroplasmy frequency in the samples, and the error bars represent the standard error of the frequency. The average frequency is significantly correlated with the theoretical contamination level (Pearson correlation=0.996781, P-value=6.212e-09). [Figure 2B]The performance of the method described herein on virtual contaminated samples is shown. Virtual contaminated samples were created by mixing two samples from the 1000 Genomes Project in different ratios. The X-axis shows the theoretical contamination level, and the Y-axis shows the heteroplasmy frequency determined from each virtual contaminated sample. Each colored dot represents one heteroplasmy, the black dot represents the average heteroplasmy frequency in the samples, and the error bars represent the standard error of the frequency. The average frequency is significantly correlated with the theoretical contamination level (Pearson correlation=0.996781, P-value=6.212e-09). [Diagram 3] 1 shows the results of detection of contamination in a hypothetical contaminated sample. [Figure 4-1] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 4-2] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 4-3] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 4-4] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 4-5] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 4-6] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Diagram 4-7] Results from detection of sample swapping and contamination in RNA-seq data 1 are shown. [Figure 5-1] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-2] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-3] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-4]Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-5] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-6] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-7] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. [Figure 5-8] Results from detection of sample swapping and contamination in RNA-seq data 2 are shown. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0011] Herein, we present a method to leverage mtDNA sequence information to detect potential sample mislabeling and contamination in NGS data. mtDNA polymorphisms and mutations can be used to estimate the identity of a particular biological sample and can serve as indicators of sample mislabeling. In addition, if a biological sample is contaminated by DNA / RNA from another biological sample, aberrant mtDNA mutation patterns will be revealed, which can help identify and further quantify the contaminant. Compared to approaches based on nuclear DNA mutations, the method described herein allows for higher sensitivity even in low-coverage sequencing data.

[0012] The methods described herein can take as input any NGS data with sufficient mtDNA reads, identify mtDNA variants (heteroplasmy and homoplasmy) from the data, and use the variant information to assign haplogroups to each sample to detect potential sample shuffling or mislabeling. By assessing the heteroplasmy information of samples, the methods described herein can further detect contamination between individuals.

[0013] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting. A human can have wild-type (or reference) mtDNA molecules and mutant mtDNA molecules. If a human does not have mutant mtDNA molecules, then such a human is considered to be homoplasmic wild-type (or homoplasmic reference). If a human does not have wild-type mtDNA molecules (i.e., has only mutant mtDNA), then such a human is considered to be homoplasmic mutant. Thus, homoplasmy is a measure of having all or none of the copies of mutant mtDNA.

[0014] If a human has a mixture of wild-type and mutant mtDNA molecules, the human is considered to have heteroplasmy. The proportion of mutant copies is referred to herein as "heteroplasmy frequency." For example, assuming that a human has eight copies of mtDNA molecules and a single copy of the eight mtDNA molecules with a particular mutation in gene A, the heteroplasmy frequency of such a human is considered to be 12.5% ​​(i.e., 1 / 8). Heteroplasmy can be determined for each mutation in the mtDNA genome of a particular individual. Thus, an individual with two mtDNA mutations (compared to wild-type mtDNA) can have two heteroplasmies. Each heteroplasmy is associated with its own heteroplasmy frequency.

[0015] The present disclosure provides a method for identifying an unreliable biological sample. The method includes performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) sequencing reads for each biological sample. The method also includes identifying the presence of one or more heteroplasmies and homoplasmies in the mtDNA sequencing reads for each of the biological samples. The method also includes assigning a primary mtDNA haplogroup to each biological sample. A biological sample having an assigned primary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroups of the majority of biological samples from the same individual is an unreliable biological sample. Such an unreliable biological sample may, for example, be mislabeled or swapped with another biological sample.

[0016] The nucleic acid sequencing assay is any nucleic acid sequencing protocol. In some embodiments, the sequencing assay comprises next generation sequencing (NGS). In some embodiments, the NGS comprises whole genome sequencing. In some embodiments, the NGS comprises whole exome sequencing. In some embodiments, the NGS comprises RNA sequencing. In some embodiments, the NGS comprises bisulfite sequencing.

[0017] A nucleic acid sequencing assay is performed on each biological sample of the multiple biological samples obtained from a single individual. In some embodiments, the multiple biological samples can be as few as two to several thousand samples. In some embodiments, the multiple biological samples can be as few as two to several hundred samples. In some embodiments, the multiple biological samples are obtained from one or more clinical studies. In some embodiments, the multiple biological samples of a single individual may be mixed or batched with multiple biological samples from another individual. mtDNA sequencing reads are obtained for each biological sample.

[0018] The presence of one or more heteroplasmies and homoplasmies in the mtDNA sequencing reads for each of the biological samples is determined. Thus, for each mutation identified in the mtDNA sequencing reads, a heteroplasmy and homoplasmy analysis is performed. The sum of all heteroplasmies is represented by the total heteroplasmy count for a particular biological sample. The mtDNA sequencing read information for each mtDNA mutation site is compiled to provide a summary of the sequencing information of the mapped reads at each single site. In some embodiments, the compiling can be performed, for example, using the samtools mpileup function (Li et al., Bioinformatics, 2009, 25, 2078-2079). The mtDNA sequencing read information for each mtDNA mutation site is filtered by sequence quality, for example, to remove sequencing bases with low sequencing quality to reduce sequencing errors. In some embodiments, a sequence quality score (Q) is determined. This is a property that is logarithmically related to the sequencing error probability (Q=-10 * log 10 (P), where P is the sequencing error probability. In some embodiments, the sequence quality Q is ≧20. When Q is 20, the sequencing error probability is 1%.

[0019] In some embodiments, heteroplasmy is identified by determining sequence coverage, the presence of minor alleles, and minor allele frequency. Sequencing coverage represents the number of reads that align to a known mtDNA reference base. In some embodiments, the sequencing coverage is ≧50. Sequencing coverage is generated by mpileup, excluding bases with Q<20. In some embodiments, the minor allele frequency is ≧1% for DNA sequencing data and ≧5% for RNA sequencing data. In some embodiments, the minor allele is observed at least twice from each DNA strand, or the minor allele is observed at least three times for RNA. For example, the following mtDNA sequencing reads may be obtained (the first sequence is the reference sequence):

[0020] [ka]

[0021] [ka]

[0022] The heteroplasmy frequency for this candidate mtDNA heteroplasmy site is 25% (2 / 8). In this particular analysis, the sequencing quality is >20 and the sequencing coverage is >50. The minor allele is observed in both strands of DNA. Therefore, this particular mutation site (i.e., the candidate mtDNA heteroplasmy site) is mtDNA heteroplasmic.

[0023] In some embodiments, homoplasmy is identified by determining the sequencing coverage and the presence of one or more alleles. Homoplasmy exists if i) the sequence coverage is ≧10, and ii) only one allele is observed at a particular nucleic acid mutation site and is different from the corresponding reference allele, or multiple alleles are observed at a particular nucleic acid mutation site and the major allele is different from the corresponding reference allele, and the particular nucleic acid site is not heteroplasmic and does not meet the heteroplasmy identification criteria.

[0024] In some embodiments, assigning a primary mtDNA haplogroup to each biological sample includes constructing an mtDNA sequence for each biological sample. In some embodiments, the mtDNA sequence for each biological sample is constructed using major alleles of homoplasmy and heteroplasmy. In some embodiments, HaploGrep2 is used to assign primary mtDNA haplogroups based on the constructed mtDNA sequence (Weissensteiner et al., Nuc. Acids Res., 2016, 44, W58-W63). HaploGrep2 is an algorithm in which haplogroups are classified based on pre-calculated phylogenetic weights that correspond to the occurrence of mutations per position in the Phylotree. Similar tools for assigning primary mtDNA haplogroups include mthap (World Wide Web "dna.jameslick.com / mthap / ") and haplofind (World Wide Web "haplofind.unibo.it / ").

[0025] A biological sample having an assigned primary mtDNA haplogroup that differs from the assigned primary mtDNA haplogroup of the majority of biological samples from the same individual is an unreliable biological sample. In some embodiments, the unreliable biological sample is mislabeled. In some embodiments, the unreliable biological sample is replaced with another biological sample. In some embodiments, one or more mislabeled samples are correctly relabeled. In some embodiments, one or more mislabeled samples are discarded.

[0026] In some embodiments, the method further comprises determining the heteroplasmy number for each biological sample. In some embodiments, the heteroplasmy frequency is determined for each mutation identified in the mtDNA sequence for each biological sample. If the heteroplasmy number of the biological sample is high, the biological sample is assigned to a secondary mtDNA haplogroup. In some embodiments, the threshold for having a high heteroplasmy number is ≧10 heteroplasmies.

[0027] In some embodiments, assigning the secondary mtDNA haplogroup comprises constructing a secondary mtDNA sequence using homoplasmic and heteroplasmic minor alleles. A biological sample having an assigned secondary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroup is a contaminated and unreliable sample. The primary haplogroup selection is the same as the secondary haplogroup selection.

[0028] In some embodiments, the method further comprises determining the contamination level of the biological sample. In some embodiments, the contamination level is indicated by determining the median heteroplasmy frequency of all heteroplasmies in the contaminated sample. The higher the median heteroplasmy frequency, the higher the contamination level. There is a strong correlation between the actual contamination percentage and the median / average heteroplasmy frequency. For example, if the median heteroplasmy frequency is 6%, the contamination level is also about 6%.

[0029] In some embodiments, the method further comprises processing the mtDNA sequencing reads obtained from the nucleic acid sequencing assay for quality control and adapter sequence removal before identifying heteroplasmy and homoplasmy. In such embodiments, the mtDNA sequencing reads obtained from the nucleic acid sequencing assay are mtDNA raw sequencing reads. When performing the processing of the mtDNA raw sequencing reads, quality-controlled mtDNA sequencing reads are generated. In some embodiments, the processing of the mtDNA sequencing reads obtained from the nucleic acid sequencing assay for quality control and adapter sequence removal can be performed by using "Trimmomatic" (Bolger et al., Bioinformatics, 2014, 30, 2114-2120). This processing step improves the accuracy of subsequent mtDNA variant identification. Another tool that can be used for processing is cutadpt (World Wide Web "cutadapt.readthedocs.io / en / stable / ").

[0030] In some embodiments, the method further comprises a two-stage mapping process before identifying heteroplasmy and homoplasmy. In some embodiments, the mtDNA sequencing reads obtained from the nucleic acid sequencing assay can be used in the two-stage mapping process. In some embodiments, the quality-controlled mtDNA sequencing reads obtained from the quality control and adapter sequence removal process can be used in the two-stage mapping process. In these embodiments, the mtDNA sequencing reads (obtained from the nucleic acid sequencing assay) or the quality-controlled mtDNA sequencing reads (obtained from the quality control and adapter sequence removal process) are mapped to the mitochondrial reference genome to generate candidate mtDNA sequencing reads. In some embodiments, the mitochondrial reference genome is the revised Cambridge Reference Sequence for Mitochondrial Genome (rCRS). In some embodiments, the mapping step can be performed using "bowtie2" (Langmead et al., Nature Methods, 2012, 9, 357-359) or bwa. The candidate mtDNA sequencing reads obtained from the first mapping step are remapped to the entire human reference genome. In some embodiments, the human reference genome is GRCh38 for the nuclear genome. Additionally, GRCh37 can also be used. In some embodiments, the mapping step can be performed using "bowtie2".

[0031] After performing the two-stage mapping process, a candidate mtDNA sequencing read is retained under two circumstances: 1) if the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome, the candidate mtDNA sequencing read is retained, and 2) if the alignment mismatch count of the candidate mtDNA sequencing read is less than five mismatched bases.

[0032] In some embodiments, the method further comprises processing the mtDNA sequencing reads (obtained from the nucleic acid sequencing assay) for sorting and duplicate removal. In some embodiments, the method further comprises processing the quality-controlled mtDNA sequencing reads (obtained from the quality control and adapter sequence removal process) for sorting and duplicate removal. In some embodiments, the method further comprises performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and duplicate removal. In some embodiments, the processing for sorting and duplicate removal can be performed by using the "samtools toolkit" (Li et al., Bioinformatics, 2009, 25, 2078-2079). These processing steps are standard next-generation sequencing (NGS) data processing steps. The GATK toolkit can also be used.

[0033] In some embodiments, the method further comprises downsampling the mtDNA sequencing reads obtained from the nucleic acid sequencing assay to a desired depth before identifying heteroplasmy and homoplasmy. In some embodiments, the method further comprises downsampling the mtDNA sequencing reads obtained from the nucleic acid sequencing assay to a desired depth before processing the mtDNA raw sequencing reads for quality control and adapter sequence removal. In some embodiments, the mtDNA raw sequencing reads from the whole transcriptome dataset can be downsampled to 10 million reads. In some embodiments, the downsampling can be performed by using "seqtk" (World Wide Web "github.com / lh3 / seqtk"). The mtDNA content of RNA-seq data is usually very high. Therefore, not all sequences are required to perform the methodology described herein, as greater mtDNA coverage will require longer computation times. In some embodiments, the desired depth is about 1000, but could be as low as about 200. Additional tools that can be used include FASTQ-SAMPLE (available on the World Wide Web at "homes.cs.washington.edu / ~dcjones / fastq-tools / fastq-sample.html").

[0034] In some embodiments, the methods described herein further comprise obtaining a plurality of biological samples from the individual prior to performing a nucleic acid sequencing assay on the plurality of samples. In some embodiments, the biological samples are blood, tissue, or tumor biopsies. In some embodiments, the methods described herein further comprise amplifying nucleic acid molecules in the biological samples prior to performing a nucleic acid sequencing assay on the plurality of samples.

[0035] The present disclosure also provides a method for identifying unreliable biological samples, the method comprising: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain DNA raw sequencing reads for each biological sample; b) processing the DNA raw sequencing reads for quality control and adapter sequence removal to generate quality controlled DNA sequencing reads; c) mapping the quality controlled DNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; and d) determining whether i) the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome, and ii) determining whether the candidate mtDNA sequencing read is an alpha-nucleotide sequence. The method includes: remapping the candidate mtDNA sequencing reads to a human reference genome and retaining the candidate mtDNA sequencing reads if the mismatch count is less than 5; e) performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and removing duplicates; f) identifying heteroplasmy and homoplasmy in the retained candidate mtDNA sequencing reads for each of the biological samples; and g) assigning a primary mtDNA haplogroup to each biological sample, where any biological sample having an assigned primary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroup of the majority of biological samples from the same individual is an unreliable biological sample that is a mislabeled biological sample. The steps of the method can be performed by the process described herein.

[0036] The present disclosure also provides a method for identifying an unreliable biological sample, the method comprising: a) performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain DNA raw sequencing reads for each biological sample; b) processing the DNA raw sequencing reads for quality control and adapter sequence removal to generate quality controlled DNA sequencing reads; c) mapping the quality controlled DNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; and d) determining whether or not the candidate mtDNA sequencing read is a candidate mtDNA sequencing read if i) the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome, and ii) the alignment mismatch count of the candidate mtDNA sequencing read is less than 5. e) performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and de-duplication; f) identifying heteroplasmy and homoplasmy in the retained candidate mtDNA sequencing reads for each of the biological samples; g) assigning a primary mtDNA haplogroup to each biological sample; and h) determining a total heteroplasmy count for each biological sample, wherein if the biological sample has a high heteroplasmy count, the biological sample is assigned a secondary mtDNA haplogroup and a biological sample having an assigned secondary mtDNA haplogroup that is different from the assigned primary mtDNA haplogroup is a contaminated unreliable sample. The steps of the method can be performed by the process described herein.

[0037] In some embodiments, the methods described herein can be implemented as a workflow. For example, a number of workflow management tools, such as Pyflow (see the world wide web at github.com / Illumina / pyflow), can be used to streamline the steps together.

[0038] The method described herein has several advantages. First, the method does not require any nuclear DNA (nDNA) variant information for each sample or at the allele frequency level of the population - this kind of information is often not available for many studies, especially RNA-seq studies. Second, the method does not require intensively pre-processed sequencing data as input, such as whole genome mapped bam files, whole genome variant VCF files, etc. The method can directly take fastq files as input. Third, the method can be applied to low coverage sequencing data. Usually, nDNA variant based methods require high coverage (>50x) to detect low levels of contamination. Due to the multi-copy nature of mtDNA, even with low coverage data, e.g., 2-4x in the 1000 Genomes Project, the mtDNA coverage can still be as high as 1000-2000x, which is enough to detect contamination levels as low as 1%. Fourth, the method does not require high computational power, and a typical sample with 1000x mtDNA coverage can be processed in 10-20 min using a single processor and 4 Gb memory. Samples with high mtDNA content may take longer to process, but can be downsampled to reduce processing time. The method can be easily incorporated into standard NGS data processing pipelines and serves as a critical quality control step by identifying problematic samples and further improving the accuracy of downstream data analysis.

[0039] In order that the subject matter disclosed herein may be more efficiently understood, examples are provided below, it being understood that these examples are for illustrative purposes and should not be construed as limiting the claimed subject matter in any manner. EXAMPLES

[0040] Example 1: Identification of mtDNA Variations and Haplogroup Assignment General methodology mtDNA mutations (both homoplasmic and heteroplasmic) are identified from next generation sequencing data, for example, by performing the analysis shown in Figure 1 (see Box 1). When performing nucleic acid sequencing assays on multiple biological samples, the raw sequencing reads can be downsampled to a desired depth to reduce the computational burden, for example, using "seqtk" which can be found on the world wide web at "github.com / lh3 / seqtk" (see Figure 1, Box 1, "Step 0"). This step is optional and does not need to be performed.

[0041] The raw mtDNA sequencing reads obtained from the nucleic acid sequencing assay (optionally from a previous downsampling step) are processed for quality control and adapter sequence removal by using "Trimmomatic" (Bolger et al., Bioinformatics, 2014, 30, 2114-2120) (see Figure 1, Box 1, "Step 1").

[0042] To obtain candidate mtDNA sequencing reads, quality-controlled sequencing reads are mapped to the mitochondrial reference genome using "bowtie2" (Langmead et al., Nature Methods, 2012, 9, 357-359) (see Figure 1, Box 1, "Step 2"). Nuclear mitochondrial DNA segments (NUMTs) in the nuclear genome may mismap to the mitochondrial genome and be counted as mtDNA reads. To minimize the effects of NUMTs, a second round of mapping can be performed, whereby the mapped reads from the first round are remapped to the entire human reference genome, GRCh38 for the nuclear genome, and the revised Cambridge Reference Sequence (rCRS) for the mitochondrial genome. A read (or read pair) is retained if a) the read (or read pair) is uniquely mapped to the mitochondrial genome or has fewer mismatches to the mitochondrial genome than to the nuclear genome, and b) the alignment mismatch count is less than 5.

[0043] The retained candidate mtDNA sequencing reads were further processed by the “samtools toolkit” (Li et al., Bioinformatics, 2009, 25, 2078-2079), including sam-to-bam conversion, sorting, and duplicate removal (see Figure 1, Box 1, “Step 3”).

[0044] The "samtools mpileup function" (Li et al., Bioinformatics, 2009, 25, 2078-2079) is used to compile the retained candidate mtDNA sequencing reads for each mtDNA site, further filtering bases by sequencing quality (>=20) and identifying heteroplasmy and homoplasmy (see Figure 1, Box 1, "Step 4"). Heteroplasmy is identified by the following criteria: a) sequencing coverage >=50, b) minor allele frequency >=1%, and c) for DNA data, the minor allele must be observed at least twice from each strand and for RNA data, the minor allele must be observed at least three times. Homoplasmy is identified by the following criteria: a) sequencing coverage >10, and b1) only one allele is observed at a given site and differs from the reference allele, or b2) multiple alleles are observed and the major allele differs from the reference, but the site does not meet the heteroplasmy criteria.

[0045] The mtDNA sequence for each sample is assembled using the homoplasmy information and the major allele at the heteroplasmic site, and "HaploGrep2" (Weissensteiner et al., Nuc. Acids Res., 2016, 44, W58-W63) is used to assign haplogroups based on the assembled sequence (see Figure 1, Box 1, "Step 5"). The haplogroups assigned in this step are called the primary haplogroup of each sample.

[0046] If a particular sample has an abnormally high heteroplasmy count, the homoplasmy information and the minor allele at the heteroplasmic site are used to construct a secondary mtDNA sequence and a secondary haplogroup is assigned based on this secondary mtDNA sequence (see Figure 1, Box 1, "Step 6").

[0047] Detection of sample mislabeling / swapping (see Figure 1, Box 2) For multiple samples, each sample can be assigned a primary haplogroup as described herein. If all samples are processed correctly, all samples from the same individual will be assigned to the same haplogroup. Conversely, if more than one haplogroup is assigned between these samples, the sample(s) with the minor haplogroup assignment will be considered to be haplogroup unmatched (i.e., mislabeled or swapped with another sample). For example, in Table 1 below, sample 001 will be considered to be swapped with sample 008.

[0048] [Table 1]

[0049] Detection and quantification of sample contamination (see Figure 1, Box 3) If an abnormally high number of heteroplasmy is observed in a particular sample, the sample may be contaminated. Primary and secondary haplogroups are assigned to the suspect sample based on the major and minor alleles at the heteroplasmic sites, respectively (see Figure 1, Box 1, "Step 5" and "Step 6"). If the primary and secondary haplogroups are different, the sample is considered to be a contaminated sample.

[0050] If a sample is determined to be contaminated, the median frequency of all heteroplasmies in the sample is used to represent the level of contamination. Example 2: Using mtDNA haplogroups to detect mislabeled samples Samples are collected from multiple individuals, and each individual has multiple samples. After obtaining RNA-seq data for a batch of clinical samples, an mtDNA haplogroup is assigned to each sample (see Table 2). Samples collected from the same individual should belong to the same mtDNA haplogroup. Mismatched mtDNA haplogroups suggest that the sample may be mislabeled. Samples with haplogroup L3h1a1 should be considered mislabeled samples.

[0051] [Table 2]

[0052] Example 3: Use of mtDNA heteroplasmy to detect sample contamination Preparation of virtual contaminant samples Whole genome sequencing fastq files for two individuals, HG00290 and NA19086, were downloaded from the 1000 Genomes Project (see the ftp site at ftp.1000genomes.ebi.ac.uk / vol1 / ftp / ). Sequencing reads were sampled from the two individuals, and NA19086 reads were mixed into HG00290 in various ratios (0.1%, 0.5%, 1%, 2%, 5%, 10%, 20%, 30%, and 40%) to create virtual contamination samples.

[0053] Real-world datasets DNA sequencing data were downloaded from the 1000 Genomes Project (see the ftp site at "ftp.1000genomes.ebi.ac.uk / vol1 / ftp / "). For each individual, reads mapped to the mitochondrial genome were extracted from the bam files by samtools (Li et al., Bioinformatics, 2009, 25, 2078-2079) and then converted to paired-end fastq files. The fastq files were used as input for the methods described herein.

[0054] Fastq files of two RNA-seq studies were downloaded from GEO (GSE81266 and GSE127165). GSE81266 contained whole transcriptome data for 77 ileum and prepouch ileum samples, including 61 paired-end (2 × 75 bp) samples and 16 single-end (50 bp) samples. GSE127165 contained whole transcriptome data from 57 laryngeal squamous cell carcinoma patients, each patient had a tumor sample and an adjacent normal sample. All samples were paired-end with a read length of 150 bp.

[0055] Analytical performance A hypothetical contaminant sample was analyzed. Whole genome sequencing data of two individuals, HG00290 and NA19086, were downloaded from the 1000 Genomes Project (Auton et al., Nature, 2015, 526, 68-74). HG00290 belongs to haplogroup U5a2a1a, and one heteroplasmy was identified in the mtDNA genome of this individual (2610T>C 1.4%), while NA19086 belongs to haplogroup D4b1a1, and two heteroplasmies were identified (1646T>C 2.1%, 12785T>T 21.3%). The two individuals had 45 nucleotide differences in their mitochondrial genomes.

[0056] Virtual contaminated samples were created by mixing sequencing reads from two samples in a series ratio ranging from 0.1% to 40%. HG00290 was treated as the original sample and NA19086 was treated as the contaminant. Each contaminated sample contained 50 million read pairs with a read length of 100 bp. Virtual contaminated samples were processed by the method described herein for contamination analysis, and the results are summarized in Figure 3. When the contamination level exceeded 2%, 45 to 46 heteroplasmies were identified in the samples, which was far above the normal range (1 to 2 heteroplasmies in one individual). These heteroplasmies covered almost all expected sites (45 segregating sites between the two individuals plus the original heteroplasmy 2610 T>C in HG00290). Only one expected site was missing in the 2% sample. The primary haplogroup was U5a2a1a (same as the original sample HG00290) for all six samples, and the secondary haplogroup was D4b1a1 (same as the contaminant). When the contamination level was 1%, 29 heteroplasmies were detected. The 17 missing sites were manually checked and confirmed to all show some heteroplasmy signal, but these sites were not cut off because the heteroplasmy frequency identification cutoff was set at 1%. The secondary haplogroup of the 1% sample was correctly assigned to D4b1a1. When the contamination level was 0.1% and 0.5%, only 1 and 11 heteroplasmies were detected, and the secondary haplogroup for these two samples was still U5a2a1a, so contamination could not be reliably detected in these low contamination level samples. These results indicate that contamination as low as 1% could be detected by combining the heteroplasmy count and secondary haplogroup assignment.

[0057] The heteroplasmy frequency in artificially contaminated samples was further evaluated. Although there was some variation in the heteroplasmy frequency in each sample, the mean and median frequency were significantly correlated with the theoretical contamination level (see Figure 2 and Figure 3, for mean and median, Pearson correlation=0.996781, 0.9979935, P value=6.212e-09, 1.189e-09, respectively). Therefore, when a given sample is detected as contaminated by the method described herein, the contamination level can be relatively quantified by the mean / median heteroplasmy frequency in the sample.

[0058] Real-world data application: RNA-seq data There are several factors that may make low frequency (<5%) heteroplasmy identification more difficult in RNA-seq data than in DNA-seq data: 1) errors introduced during the reverse transcription step, 2) RNA editing / modification, and / or 3) uneven coverage across the mtDNA genome due to variations in gene expression levels. Therefore, to reduce false positive heteroplasmy, only heteroplasmies with frequencies >5% were considered to be reliable heteroplasmies in RNA data. Furthermore, three well-defined mtDNA editing sites were excluded: 295, 2617 and 13710 (Bar-Yaacov et al., Genome Res., 2013, 23, 1789-1796, and Hodgkinson et al., Science, 2014, 344, 413-415).

[0059] The methods described herein were applied to two bulk RNA-seq datasets to evaluate the context of different diseases or tissue types. First, the methods described herein were applied to a dataset containing 77 samples from 25 subjects (Huang et al., Inflamm. Bowel Dis., 2017, 23, 366-378). Most subjects in this study had samples from different tissues (ileum and anterior pouch ileum) and / or at different biopsy time points (4 months, 8 months, 12 months, etc.). Sixteen samples in this dataset were single-end samples with a 50 bp read length and 61 samples were paired-end samples with a 75 read length. For each sample, 10 million reads (pairs) were randomly sampled and tested. First, the assignment of primary haplogroups to samples was evaluated. In this dataset, all samples from the same subject were assigned to the same mtDNA haplogroup (see Figure 4), indicating no sample shuffling. Potential contamination in these samples was then evaluated. Except for sample SRR3493833, at the 5% heteroplasmy frequency cutoff, all other samples had at most 6 heteroplasmies, and the secondary haplogroup assignment was the same as the primary haplogroup (see Figure 4). In sample SRR3493833, 29 heteroplasmies were identified, which was much higher than the normal range, and the median heteroplasmy frequency was 14.8%. The secondary haplogroup of this sample was J1c8a, which was also different from the primary haplogroup U5b2a1a. These results indicated that sample SRR3493833 was potentially contaminated by another sample from the J1 haplogroup, with a contamination level of approximately 14.8%.

[0060] The methods described herein were also applied to a dataset containing tumor samples (Wu et al., Molec. Cancer, 2020, 19, 99). This dataset included samples from 57 patients with laryngeal squamous cell carcinoma, each with a tumor sample and a paired adjacent normal mucosa sample. In this dataset, paired tumor samples and adjacent normal samples from the same patient were all assigned to the same haplogroup (see Figure 5), and no sample shuffling was detected. All samples had low heteroplasmy counts and the same primary and secondary haplogroup assignments, and therefore no detectable contamination. This dataset demonstrated that the methods described herein can identify the identity of tumor samples.

[0061] Early and accurate sample shuffling and contamination detection is a key quality control step for large-scale NGS data, as it can filter out suspicious samples and improve their quality for subsequent analysis. In these examples, an efficient method is presented for detecting sample shuffling and inter-individual contamination by using mtDNA mutations identified from NGS data. The methodology can take demultiplexed fastq files as input without any data preprocessing. First, it detects any sample shuffling for individuals with multiple samples. It further detects and quantifies potential contamination, and then suggests the source sample of the contaminant. Although whole genome DNA sequencing data from the 1000 Genomes Project and two bulk RNA-seq datasets were used as working examples for these examples, the method described herein can be generalized to any NGS dataset that contains mtDNA reads, such as whole exome sequencing data including off-site mtDNA reads, single-cell RNA-seq, ATAC-seq data, etc. Simulation results described herein show that the method described herein effectively detected contamination as low as 1%.

[0062] In addition to those described herein, various modifications of the described subject matter will become apparent to those skilled in the art from the foregoing description. Such modifications are also intended to fall within the scope of the appended claims. Each reference cited in this application (including, but not limited to, journal articles, U.S. and non-U.S. patents, patent application publications, international patent application publications, gene bank accession numbers, etc.) is incorporated herein by reference in its entirety.

Claims

1. 1. A method for identifying one or more contaminated biological samples, the method comprising: performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) sequencing reads for each biological sample; For each of the biological samples, identifying heteroplasmy and homoplasmy in the mtDNA sequencing reads from the previous step; assigning a primary mtDNA haplogroup to each biological sample; determining the total heteroplasmy count for each biological sample; assigning secondary mtDNA haplogroups to biological samples having high heteroplasmy counts; detecting biological samples having an assigned secondary mtDNA haplogroup that differs from the assigned primary mtDNA haplogroup of a majority of said biological samples from the same individual; determining that a biological sample having an assigned secondary mtDNA haplogroup that differs from the primary mtDNA haplogroup assigned to a majority of said biological samples from the same individual is a contaminated biological sample; discarding the one or more contaminated biological samples; and A method comprising:

2. 2. The method of claim 1, wherein identifying the heteroplasmy includes determining sequencing coverage, minor allele frequency, and presence of a minor allele, and heteroplasmy exists if i) the sequencing coverage is > 50, ii) the minor allele frequency is > 1%, and iii) for DNA data, the minor allele is observed at least twice from each strand, or for RNA data, the minor allele is observed at least three times.

3. 2. The method of claim 1, wherein identifying homoplasmy comprises determining sequencing coverage and the presence of one or more alleles, and homoplasmy exists if i) the sequencing coverage is > 10, and ii) only one allele is observed at a particular nucleic acid site and is different from the corresponding reference allele, or multiple alleles are observed at a particular nucleic acid site and the major allele is different from the corresponding reference allele, and the particular nucleic acid site does not meet heteroplasmy criteria.

4. 2. The method of claim 1, wherein assigning the primary mtDNA haplogroup to each biological sample comprises constructing an mtDNA sequence for each biological sample using the major alleles of the homoplasmy and the heteroplasmy.

5. 2. The method of claim 1, wherein assigning the secondary mtDNA haplogroup comprises constructing a secondary mtDNA sequence using the homoplasmic and heteroplasmic minor alleles, and identifying biological samples having assigned secondary mtDNA haplogroups that differ from the assigned primary mtDNA haplogroup as contaminated and unreliable biological samples.

6. 6. The method of claim 5, further comprising determining a level of contamination of the biological sample by determining a median heteroplasmy frequency of all heteroplasmies in the contaminated biological sample, the greater the median heteroplasmy frequency, the greater the level of contamination.

7. 2. The method of claim 1, further comprising processing the mtDNA sequencing reads obtained from the nucleic acid sequencing assay for quality control and adapter sequence removal prior to identifying heteroplasmy and homoplasmy to generate quality-controlled mtDNA sequencing reads.

8. The method further comprising: mapping the quality-controlled mtDNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; mapping the candidate mtDNA sequencing reads to a human reference genome; i) the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome; and 8. The method of claim 7, further comprising: ii) retaining the candidate mtDNA sequencing read if the alignment mismatch count of the candidate mtDNA sequencing read is less than 5.

9. 9. The method of claim 8, further comprising performing post-mapping processing of the retained candidate mtDNA sequencing reads for sorting and removing duplicates.

10. 2. The method of claim 1, further comprising downsampling the mtDNA sequencing reads obtained from the nucleic acid sequencing assay to a desired depth prior to identifying the heteroplasmy and the homoplasmy and / or prior to processing the mtDNA sequencing reads for quality control and adapter sequence removal.

11. The method of claim 1 , wherein the biological sample is blood or tissue.

12. 10. The method of claim 1, wherein the method further comprises amplifying nucleic acid molecules in the biological samples prior to performing the nucleic acid sequencing assay on the plurality of biological samples.

13. 10. The method of claim 1, wherein the sequencing assay comprises next generation sequencing (NGS).

14. 14. The method of claim 13, wherein the NGS comprises whole genome sequencing.

15. 14. The method of claim 13, wherein the NGS comprises whole exome sequencing.

16. 14. The method of claim 13, wherein the NGS comprises RNA sequencing.

17. 14. The method of claim 13, wherein the NGS comprises bisulfite sequencing.

18. 1. A method for identifying one or more contaminated biological samples, the method comprising: performing a nucleic acid sequencing assay on each biological sample of a plurality of biological samples obtained from a single individual to obtain mitochondrial DNA (mtDNA) raw sequencing reads for each biological sample; processing the mtDNA raw sequencing reads for quality control and adapter sequence removal to generate quality-controlled mtDNA sequencing reads; mapping the quality-controlled mtDNA sequencing reads to a mitochondrial reference genome to generate candidate mtDNA sequencing reads; mapping the candidate mtDNA sequencing reads to a human reference genome; i) the candidate mtDNA sequencing read is uniquely mapped to the mitochondrial reference genome or has fewer mismatches to the mitochondrial reference genome than to the human reference genome; and ii) retaining the candidate mtDNA sequencing read if the alignment mismatch count of the candidate mtDNA sequencing read is less than 5; performing post-mapping processing of the retained candidate mtDNA sequencing reads for filtering and de-duplication; identifying heteroplasmy and homoplasmy in the retained candidate mtDNA sequencing reads for each of the biological samples; assigning a primary mtDNA haplogroup to each biological sample; determining the total heteroplasmy count for each biological sample; assigning secondary mtDNA haplogroups to biological samples having high heteroplasmy counts; detecting biological samples having an assigned secondary mtDNA haplogroup that differs from the assigned primary mtDNA haplogroup of a majority of said biological samples from the same individual; determining that a biological sample having an assigned secondary mtDNA haplogroup that differs from the assigned primary mtDNA haplogroup is a contaminated biological sample; discarding the one or more contaminated biological samples; and A method comprising:

Images