Methods and kits for screening mitochondrial genetic markers

Abstract

The present application relates to a method and a kit for screening mitochondrial genetic markers. The method comprises: (a) obtaining a mitochondrial sequencing dataset of training samples, the training samples comprising a plurality of abnormal group samples and a plurality of normal group samples; (b) determining mitochondrial genetic mutation information of the training samples based on the mitochondrial sequencing dataset, the mitochondrial genetic mutation information comprising a gene mutation heterogeneity level and at least one selected from a gene mutation site, a base type and a gene mutation type, and optionally, the gene mutation type comprises that the gene mutation is a homogeneous mutation or a heterogeneous mutation; and (c) selecting, based on the mitochondrial genetic mutation information determined in step (b), a mitochondrial genetic mutation capable of being used to distinguish the abnormal group samples from the normal group samples as the mitochondrial genetic marker. The method described in the present application can be used to screen markers for detecting mitochondrial diseases.

Description

Technical Field

[0001] This application relates to the field of biotechnology, specifically to methods and kits for screening mitochondrial gene markers, and more specifically to methods, mitochondrial gene markers, kits, and uses for determining mitochondrial gene markers. Background Technology

[0002] Mitochondria are organelles found in most eukaryotic cells, comprising approximately 10% of the cellular proteome and performing a variety of functions related to cellular metabolism and homeostasis. Cellular energy production via oxidative phosphorylation (OXPHOS) is a hallmark of mitochondria. Furthermore, mitochondria also play roles in calcium homeostasis, the initiation of caspase-dependent apoptosis, cellular stress responses, heme biosynthesis, sulfur metabolism, and cytoplasmic protein degradation.

[0003] Human mitochondrial DNA (mtDNA) is a double-stranded circular DNA molecule, 16,569 bp in length, comprising 13 protein-coding genes, 22 transfer RNA (tRNA) genes, 2 ribosomal RNA (rRNA) genes, and a D-loop region. Unlike the double copies of nuclear DNA, each human cell contains hundreds or thousands of mtDNA copies. Considering the multiple mtDNA copies within each cell, mutations can affect all mtDNA molecules (called homogeneous sites) or a subset of mtDNA molecules (called heterogeneous sites). The proportion of these damaged or mutated mtDNA molecules determines whether an individual will exhibit a clinical phenotype.

[0004] Methods for detecting mitochondrial heterogeneity sites can be categorized into three types: Sanger sequencing, microarray-based methods, and next-generation sequencing (NGS)-based methods. Sanger sequencing, also known as dideoxy chain termination, utilizes the DNA replication principle. This method can detect a minimum heterogeneity level of approximately 15% and is only suitable for qualitative studies of small sample sizes of heterogeneous sites. Microarray-based methods, compared to Sanger sequencing, offer advantages such as shorter detection cycles, lower cost, and higher throughput, but can only detect known mutations and cannot discover new heterogeneous sites. Next-generation sequencing methods offer even higher throughput and provide the possibility of detecting heterogeneous sites with even lower mutation levels.

[0005] Mitochondrial heterogeneous sites are widespread in healthy individuals, but their pathogenic potential is not well described. Studying the pathogenic potential of heterogeneous sites will help further understand the role of mtDNA in processes such as aging, epilepsy, and neurodegeneration. However, there is still a lack of unified standards for the detection and pathogenicity definition of mitochondrial heterogeneous sites, resulting in a lack of effective analytical methods and procedures for studying the association between mitochondrial mutations and diseases.

[0006] Therefore, research on mitochondrial gene markers still has great potential for development. Summary of the Invention

[0007] This application is made by the inventors based on their findings regarding the following problems and facts:

[0008] Currently, the main methods for detecting mitochondrial heterogeneity sites include Sanger sequencing, microarray-based methods, and next-generation sequencing (NGS) methods. Sanger sequencing is only suitable for small sample sizes and can only qualitatively characterize heterogeneous sites with a heterogeneity level greater than 15%, not quantitatively. Microarray-based methods are only applicable to known heterogeneous sites. NGS methods, which rely on variant detection software for heterogeneity site detection, cannot detect the level of heterogeneity and cannot consider the different impacts of different heterogeneity levels on characterization. Methods that use the ratio of mutant to wild-type reads as the level of heterogeneity fail to consider the interference of sequencing quality and alignment quality on the level of heterogeneity, resulting in a high false positive rate.

[0009] This application aims to address at least one of the technical problems existing in the prior art. To this end, one objective of this application is to propose a method for screening mitochondrial gene markers. Based on existing next-generation sequencing data, heterogeneous sites in mitochondria from a large-scale sample are detected and analyzed.

[0010] In a first aspect of this application, a method for screening mitochondrial gene biomarkers is proposed. According to an embodiment of this application, the method includes: (a) acquiring a mitochondrial sequencing dataset of training samples, the training samples including multiple abnormal group samples and multiple normal group samples; (b) based on the mitochondrial sequencing dataset, determining mitochondrial gene mutation information of the training samples, the mitochondrial gene mutation information including gene mutation heterogeneity level and at least one selected from gene mutation site, base type, and gene mutation type, optionally, the gene mutation type including the gene mutation being a homogeneous mutation or a heterogeneous mutation; and (c) based on the mitochondrial gene mutation information determined in step (b), selecting mitochondrial gene mutations capable of distinguishing the abnormal group samples and the normal group samples as the mitochondrial gene biomarkers.

[0011] According to embodiments of this application, mitochondrial sequencing datasets from multiple abnormal and normal groups are obtained using the above method, establishing a relatively comprehensive database of mitochondrial gene mutation information. This makes the screened mitochondrial gene biomarkers more representative and reliable. Furthermore, by determining mitochondrial gene mutation information, even when certain gene mutation types cannot be accurately identified, sites with high levels of change and high mutation rates can be selected as mitochondrial gene biomarkers. This helps improve the sensitivity and specificity of the biomarkers.

[0012] According to embodiments of this application, the above-described method for screening mitochondrial gene markers may also have at least one of the following additional technical features:

[0013] According to embodiments of this application, the level of gene mutation heterogeneity is determined through the following steps:

[0014] (1) The mitochondrial sequencing data is compared with the reference genome to determine the mitochondrial DNA site and base type corresponding to the sequencing read;

[0015] (2) Based on the results of the alignment, determine the major base type of at least one of the mitochondrial DNA sites, wherein the major base type has the most supporting sequencing reads;

[0016] (3) For a given base type at a given site, the heterogeneity level of the base type is determined by maximum likelihood estimation based on the relevant information of multiple sequencing reads corresponding to the given site.

[0017] It should be noted that the sequencing depth in the mitochondrial sequencing data is no less than 100X.

[0018] It should be noted that the “reference genome” refers to the genome sequence of the species corresponding to the known sample. It can be a genome sequence obtained through public channels or obtained through sequencing assembly. It can be the entire genome sequence or a part of the genome of interest. For example, when analyzing human samples, multiple versions of human genome sequences provided by public databases can be used, such as hg19.

[0019] For example, the mitochondrial DNA site corresponding to the sequencing read can be determined. For instance, at position 11719 of the mtDNA, if the reference genome sequence is NC_012920.1, and a sequencing read is mapped to this position after alignment with the reference genome, then it can be determined that the read corresponds to position 11719 in the MT-ND4L gene, and the base type at position 11719 in the MT-ND4L gene is the base type of the sequencing read.

[0020] According to embodiments of this application, determining the heterogeneity level of a base type through maximum likelihood estimation means determining the heterogeneity level of that base using maximum likelihood estimation. Specifically, assuming there are two base types, such as A and G, at a certain site, the heterogeneity level of these two types at this position is determined using sequencing data. Since the heterogeneity level of each base is a value between 0 and 1, a parameter p can be defined to represent the heterogeneity level p of the G base. Using the maximum likelihood estimation method, based on the statistical information of the sequencing data, the p value that maximizes the probability of obtaining observed data is found, i.e., maximizing the probability of observing data. This p value is the optimal estimate of the heterogeneity level.

[0021] According to embodiments of this application, the relevant information includes: sequencing quality of sequencing reads, number of supporting sequencing reads for the given base type, and number of supporting sequencing reads for the major base type.

[0022] According to an embodiment of this application, in step (4), the heterogeneity level of the given base at the given site is determined by maximum likelihood estimation using the following likelihood function:

[0023] ;

[0024] in,

[0025] f is the level of heterogeneity for the given base type at the given site;

[0026] l is the number of supporting sequencing reads of the given base type at the given site;

[0027] i represents the i-th supporting sequencing read of the given base type at the given site;

[0028] This represents the sequencing error rate of the i-th supporting sequencing read of the given base type at the given site;

[0029] k represents the number of supporting sequencing reads for the major base type at the given site;

[0030] j represents the j-th supporting sequencing read of the major base type at the given site;

[0031] This represents the sequencing error rate of the j-th supporting sequencing read of the major base type at the given site.

[0032] According to an embodiment of this application, after determining the heterogeneity level of the given base at the given site, a log-likelihood ratio test is further performed on the heterogeneity level.

[0033] According to an embodiment of this application, the log-likelihood ratio test is performed using the following formula:

[0034]

[0035] in,

[0036] f ^m1 represents the level of heterogeneity of the given base at the given site, as determined by the maximum likelihood estimation;

[0037] f ^m0 represents the frequency when the given base at the given site is assumed to be a homogeneous mutation;

[0038] f mtDB This indicates the frequency of the non-dominant allele at the given locus in the mtDB database (Ingman, M. & Gyllensten, U. mtDB: Human Mitochondrial Genome Database, a resource for population genetics and medical sciences. Nucleic Acids Res. 34, D749-751 (2006).).

[0039] It should be noted that the non-dominant alleles mentioned are also called minor alleles. Taking a biseleural locus as an example, if one allele occurs more frequently than another, the allele with the higher frequency is usually called the major allele, while the allele with the lower frequency is called the minor allele. For example, in a population with a certain SNP at a certain gene locus, 178 individuals have the genotype of uppercase letter A, while 22 individuals have the genotype of uppercase letter G. In this case, uppercase letter A is the dominant allele, and uppercase letter G is the minor allele.

[0040] According to embodiments of this application, a heterogeneity level not greater than 0.9 and an LLR score not less than 5 for the heterogeneous site are indicators that the given gene site belongs to a heterogeneous site; a heterogeneity level greater than 0.9 is an indicator that the given gene site belongs to a homogeneous site.

[0041] According to embodiments of this application, the gene mutation sites are obtained through conservation fraction annotation of bases at each site, mutation harmfulness annotation of heterogeneous sites of tRNA-encoding genes in the mitochondrial genome, and / or functional annotation information.

[0042] According to an embodiment of this application, the conservation score annotation information is obtained using site conservation database analysis.

[0043] It should be noted that the conservation database can be constructed by a method including the following steps: performing multiple sequence alignment based on mitochondrial genome sequences of multiple species, calculating the conservation score of each site, and obtaining the conservation database based on the conservation score of each site.

[0044] According to an embodiment of this application, the mutation harmfulness annotation information is obtained using tRNA database analysis.

[0045] It should be noted that the tRNA database can be constructed using a method including the following steps: based on the effect of mutations in the secondary structure of tRNA downloaded by tRNAscan-SE 2.0 on base complementary pairing, the sites are annotated to obtain site annotation results; based on the site annotation results, the position of each tRNA base is statistically analyzed to determine whether it is on the "stem" or "loop" of the tRNA "cloverleaf" structure, and the position of each tRNA base is statistically analyzed to obtain the tRNA database based on the statistical results.

[0046] According to embodiments of this application, the functional annotation information is obtained using known mitochondrial data analysis.

[0047] According to embodiments of this application, the mitochondrial gene mutations that distinguish between the abnormal group samples and the normal group samples are determined by comparing the heterogeneity levels of the gene mutation sites in different groups.

[0048] Specifically, by comparing the frequency or heterogeneity of mitochondrial mutations in different groups, sites showing significant differences were identified, and gene biomarkers were screened based on site annotation results. The identification of these significantly different sites was calculated using tools in the R language. For continuous variables, the rank-sum test was used; for discontinuous variables, Fisher's exact test or logistic regression was used.

[0049] In a second aspect of this application, a mitochondrial gene marker is proposed. According to embodiments of this application, the mitochondrial gene marker is obtained based on the method described in the first aspect. The method described in the first aspect can be used to identify mitochondrial gene markers associated with mitochondrial diseases.

[0050] In a second aspect, this application proposes a mitochondrial gene marker. According to embodiments of this application, the mitochondrial gene marker includes C12906A and C15401T. This gene marker can be used to detect mitochondrial-related diseases.

[0051] In a third aspect of this application, the use of the mitochondrial gene markers described in the second aspect in the preparation of a kit is proposed. According to embodiments of this application, the kit is used to detect age-related diseases, epilepsy, or neurodegenerative diseases.

[0052] According to embodiments of this application, the age-related diseases include coronary heart disease.

[0053] In a fourth aspect, this application provides a kit. According to embodiments of this application, the kit includes reagents for detecting mitochondrial gene markers, said markers including C12906A and C15401T. According to embodiments of this application, the kit has the advantage of being simple and convenient to use in detecting age-related diseases, epilepsy, or neurodegenerative diseases.

[0054] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0055] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0056] Figure 1 This is a diagram of the tRNA secondary structure (tryptophan). A mutation that changes a previously complementary (WC) base pair on the stem to a non-complementary pair is defined as "abolished" (e.g., an A mutation at position 5512 to T / C / G); a mutation that changes a non-WC pair on the stem to a WC pair is defined as "reconstructed" (e.g., a T mutation at position 5524 to C).

[0057] Figure 2 This is for Sanger site validation (C12906A site, C15401T site). The top left is the Sanger plot corresponding to 12906bp of the mutated sample, and the bottom left is the Sanger plot corresponding to 12906bp of the non-mutated sample; the top right is the Sanger plot corresponding to 15401bp of the mutated sample, and the bottom right is the Sanger plot corresponding to 15401bp of the non-mutated sample. Detailed Implementation

[0058] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0059] Definitions and Explanations

[0060] Unless otherwise specified in this article, the term "sequencing read" is used interchangeably with the term "read" or "segment". It refers to a nucleic acid sequence obtained during sequencing. This sequence is referred to as "sequencing read" or "read" in this article. Typically, the length of sequencing reads in first-generation sequencing and third-generation sequencing is between several thousand and tens of thousands of bp, while the sequencing reads in second-generation sequencing are relatively shorter, averaging tens to hundreds of bp.

[0061] Prior art related to this application:

[0062] Regarding methods for detecting heterogeneous sites:

[0063] (1) Sanger sequencing

[0064] The Sanger sequencing reaction system includes the target DNA fragment, deoxyribonucleotide triphosphates (dNTPs), dideoxynucleotide triphosphates (ddNTPs), sequencing primers, and DNA polymerase. The core of the sequencing reaction is the ddNTPs used: lacking a 3'-OH group, they cannot form a phosphodiester bond with another dNTP; these ddNTPs can be used to terminate DNA chain elongation. Furthermore, these ddNTPs are attached with radioactive isotopes or fluorescent labeling groups, allowing them to be detected by automated instruments or gel imaging systems.

[0065] The Sanger sequencing method is only suitable for small sample sizes and can only qualitatively identify heterogeneous sites with a heterogeneity level >15%, detecting the specific numerical value of the site heterogeneity level.

[0066] (2) Microarray chip-based method

[0067] Microarray-based methods refer to methods that involve attaching oligonucleotide probes to the solid surface of a microarray chip, and then hybridizing these probes with the fragments to be detected, thereby detecting mutations.

[0068] The microarray chip-based method is only applicable to known heterogeneous or homogeneous sites.

[0069] (3) Methods based on next-generation sequencing

[0070] Regarding methods for detecting heterogeneous sites, one approach is based on mutation detection software, such as GATK and Samtools. Another approach is based on the ratio of mutant to wild-type sequencing reads as the level of heterogeneity at a site.

[0071] Methods based on next-generation sequencing. Methods that detect heterogeneous or homogeneous sites using variant detection software cannot account for the different impacts of sites with varying levels of heterogeneity on characterization. Methods that use the ratio of mutant to wild-type sequencing reads as a measure of heterogeneity fail to consider the interference of sequencing quality and alignment quality on heterogeneity levels, resulting in a high false positive rate.

[0072] The method described in this application

[0073] In one aspect of this application, a method for screening mitochondrial gene markers is proposed. According to embodiments of this application, the method includes:

[0074] (a) Obtain the mitochondrial sequencing dataset of the training samples, which includes multiple abnormal group samples and multiple normal group samples;

[0075] (b) Based on the mitochondrial sequencing dataset, determine the mitochondrial gene mutation information of the training samples, wherein the mitochondrial gene mutation information includes the level of gene mutation heterogeneity and at least one selected from gene mutation sites, base types, and gene mutation types.

[0076] Optionally, the gene mutation type includes whether the gene mutation is a homogeneous mutation or a heterogeneous mutation; and

[0077] (c) Based on the mitochondrial gene mutation information determined in step (b), select a mitochondrial gene mutation that can be used to distinguish the abnormal group sample from the normal group sample as the mitochondrial gene marker.

[0078] Beneficial effects

[0079] This application, based on next-generation sequencing data, provides a method for detecting heterogeneous and homogeneous sites in mitochondria from large-scale samples, along with subsequent analysis methods for these sites. By acquiring mitochondrial sequencing datasets from multiple abnormal and normal groups using this method, a relatively comprehensive database of mitochondrial gene mutation information is established, making the selected mitochondrial gene biomarkers more representative and reliable. Furthermore, by determining mitochondrial gene mutation information, even when certain gene mutation types cannot be accurately identified, sites with high levels of change and high mutation rates can be selected as mitochondrial gene biomarkers. This helps improve the sensitivity and specificity of the biomarkers. In addition, based on the above-mentioned mitochondrial gene biomarker screening method, this application also proposes a set of mitochondrial gene biomarkers, including C12906A and C15401T. The mitochondrial gene biomarkers discovered in this application can be used to detect age-related diseases, epilepsy, or neurodegenerative diseases.

[0080] It should be noted that the features and technical effects described in this article for different aspects can be mutually referenced, and will not be elaborated further here.

[0081] The following examples illustrate this application, but should not be construed as limiting the scope of the subject matter of this application to the following examples. All technologies implemented based on the above content of this application fall within the scope of this application. The compounds or reagents used in the following examples are commercially available or prepared using conventional methods known to those skilled in the art; the experimental instruments used are commercially available.

[0082] Example 1: Detecting mitochondrial heterogeneity sites to obtain gene markers

[0083] 1. Data Sources and Processing

[0084] 1.1 Obtaining Sample Data

[0085] The samples in this embodiment of the invention include samples from the discovery phase and samples from the validation phase. The discovery phase samples include 950 patients with early-onset coronary artery disease and 1000 geographically matched healthy individuals; the validation phase samples include 1005 samples of early-onset coronary artery disease and 2704 samples from the public database mtDB. Mitochondrial genome next-generation sequencing data were obtained from the discovery phase and validation phase samples: except for the mtDB database which contained site frequency information, the next-generation sequencing data for the remaining samples were exon capture sequencing data.

[0086] 1.2 Obtaining Reference Genome Data

[0087] The human mitochondrial genome sequence rCRS (NCBI Reference Sequence: NC_012920.1 (02-SEP-2020), which consists of 16,569 bases) and the human genome sequence hg19 were downloaded from GenBank. The mitochondrial sequence carried on the hg19 sequence was removed, and the hg19 sequence with the mitochondrial sequence removed was merged with the rCRS. An index was then created to obtain the human reference genome sequence for subsequent analysis.

[0088] 2. Data comparison to obtain comparison files

[0089] 2.1 Data Comparison

[0090] The sequencing reads from the sequencing data of the samples in the discovery and validation phases were aligned to the human reference genome obtained in step 1.2 using the bwa mem software. Only the sequencing reads that were aligned to the mitochondrial genome were retained, and a bam format alignment file was generated.

[0091] 2.2 Filtering Samples and Detection Areas for Heterogeneous Sites

[0092] 2.2.1 Recalculate the base pairing quality for each site.

[0093] The mpileup command –E parameter in the Samtools software was used to recalculate the base alignment quality (BAQ) of each base site in the BAM format alignment file, and an mpileup format file was generated. This file includes the alignment information of each base site, the position information of the base site, the minor or major allele information, the sequencing error rate of the base site, and the alignment quality value of the sequencing read length (related literature: Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N. et al. The Sequence Alignment / Map format and SAMtools. Bioinformatics. 25, 2078-2079 (2009)).

[0094] 2.2.2 Filtering sample alignment data and detection areas of heterogeneous sites

[0095] To reduce false positives in subsequent analyses, it is necessary to filter out abnormal samples (including samples with abnormal sequencing depth, closely related samples, etc.) and remove sites with high false positive rates (including sites with excessively low sequencing depth, sites with unbalanced forward and reverse strand sequencing read length ratios, homopolymer sites, etc.). The sequencing depth and base site depth of the samples are statistically analyzed based on the mpileup format file.

[0096] After filtering, at the sample level, 832 individuals with early-onset coronary heart disease and 850 control individuals remained in the phase sample, while 894 individuals with early-onset coronary heart disease and 2704 samples from the mtDB database remained in the validation phase sample. At the locus level, 14,585 bp regions in the phase sample and 10,712 bp regions in the validation phase sample met the criteria and can be used for subsequent analysis.

[0097] 3. Detection of heterogeneous sites, homogeneous sites, and variant sites.

[0098] 3.1 Detection of heterogeneous and homogeneous sites

[0099] Using base quality values ​​of sequencing reads supporting mutant / wild-type sequences, alignment quality values, and the frequency of the site in the mtDB2 database, the heterogeneity level of each site in the mpileup format file is calculated based on the maximum probability algorithm.

[0100] Specifically, the heterogeneity level and LLR score are calculated based on the maximum probability method.

[0101] To further avoid false positives, bases with a base quality <20 or an alignment quality <50 in the sequencing read were excluded from the heterogeneity level calculation. The formula for calculating the heterogeneity level is shown in Equation I (related literature: Ye, K., Lu, J., Ma, F., Keinan, A. & Gu, Z. Extensive pathogenicity of mitochondrial heteroplasmy in healthy human individuals. Proc Natl Acad Sci US A. 111, 10654-10659(2014)). The formula for calculating the LLR score is shown in Equation II.

[0102] (Formula I);

[0103] in,

[0104] f is the level of heterogeneity for the given base type at the given site;

[0105] l is the number of supporting sequencing reads of the given base type at the given site;

[0106] i represents the i-th supporting sequencing read of the given base type at the given site;

[0107] This represents the sequencing error rate of the i-th supporting sequencing read of the given base type at the given site;

[0108] k represents the number of supporting sequencing reads for the major base type at the given site;

[0109] j represents the j-th supporting sequencing read of the major base type at the given site;

[0110] This represents the sequencing error rate of the j-th supporting sequencing read of the major base type at the given site.

[0111]

[0112] (Formula II);

[0113] in,

[0114] f ^m1 represents the level of heterogeneity of the given base at the given site, as determined by the maximum likelihood estimation;

[0115] f^m0 represents the frequency when the given base at the given site is assumed to be a homogeneous mutation;

[0116] f mtDB This indicates the frequency of the non-dominant allele at the given locus in the mtDB database.

[0117] The LLR threshold for high-quality heterogeneous sites can be selected based on sequencing depth. Higher sequencing depth allows for a higher threshold, while lower sequencing depth allows for a lower threshold. In this example, an LLR ≥ 5 is used as the threshold for high-quality heterogeneous sites. A heterogeneity level (L(f)) less than or equal to 0.9 is considered a heterogeneous site; a heterogeneity level greater than 0.9 is considered a homogeneous site. In the discovery phase samples, 822 heterogeneous sites (L(f) ≤ 0.9 and LLR ≥ 5) were found at 292 mitochondrial genomic locations, and 56,180 homogeneous sites (L(f) > 0.9 and LLR ≥ 5) were found at 2,496 mitochondrial genomic locations.

[0118] To verify the accuracy of heterogeneous and homogeneous loci, 38 samples corresponding to 7 loci were randomly validated, with a consistency rate of 100% (Sanger plots of some loci are shown below). Figure 2 ).

[0119] 3.2 Detection of variant sites

[0120] Varscan4 software was used for mutation detection (related literature: Koboldt, DC, Zhang, Q., Larson, DE, Shen, D., McLellan, MD, Lin, L. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568-576 (2012)). The frequency threshold for Varscan can be found in the aforementioned literature on this disease and set according to the threshold for pathogenicity heterogeneity level corresponding to the disease type.

[0121] 4. Annotation database construction and site annotation

[0122] 4.1 Construction of a database of conserved mitochondrial genome sites

[0123] Multiple sequence alignments were performed on mitochondrial genome sequences from multiple species to calculate the conservation of each site and obtain a conservation score for each site.

[0124] Mitochondrial genome sequences from 14 vertebrates were downloaded from GenBank (see Table 1). Multiple sequence alignment was performed using Muscle software. Using human (Homo sapiens) bases as wild-type bases, the conservation score of each base site relative to the wild-type bases in the 14 species was calculated, ultimately obtaining a mitochondrial genome site conservation database.

[0125] Table 1: Mitochondrial genome sequence information of 14 vertebrates

[0126]

[0127] 4.2 tRNA Database Construction

[0128] The effects of mutations in the secondary structure of tRNA downloaded using tRNAscan-SE 2.0 on base pairing were analyzed, and the sites were annotated (related literature: Chan, PP, Lin, BY, Mak, AJ&Lowe, TM tRNAscan-SE2.0: improved detection and functional classification of transfer RNA genes. Nucleic Acids Res. 49, 9077-9096 (2021)).

[0129] Specifically, the secondary structures of 22 tRNAs were downloaded from tRNAscan-SE 2.0. Based on the secondary structure image of each tRNA (e.g., ...), ... Figure 1 As shown in the diagram, the position of each base in the tRNA is determined by analyzing whether it is located on the "stem" or the "loop" of the tRNA's "cloverleaf" structure to obtain a tRNA database. If a base on the "stem" mutates, it affects the hydrogen bonds formed between complementary bases, leading to a change in the three-dimensional structure at that location. During the tRNA's function, the three-dimensional structure influences its functionality.

[0130] Harmful mutations on tRNA fall into two categories: "abolition" and "reconstruction." "Abolition" occurs when a site mutation causes a non-complementary base pairing (WC) structure in the tRNA secondary structure stem to become a non-complementary pairing structure. "Reconstruction" occurs when a site mutation causes a non-WC pair in the tRNA secondary structure stem to become a WC pair structure. Figure 1 ).

[0131] 4.3 Downloading existing database data and literature data

[0132] Download data from the MT_ensGene database (including the gene at each site) and the Mitimpact database (including annotations for OXPHOS_complex, Ensembl_gene_id, Uniprot_name, PolyPhen2_score, FatHmm_score, etc.); and download the Mutpred analysis results (related literature: Pereira, L., Soares, P., Radivojac, P., Li, B. & Samuels, DC Comparing phylogeny and the predicted pathogenicity of protein variations reveals equal purifying selection across the global humanmtDNA diversity. Am J Hum Genet. 88, 433-439 (2011)).

[0133] 4.4 Site annotation was performed for both heterogeneous and homogeneous sites.

[0134] Site annotation was performed on heterogeneous sites (L(f)≤0.9 and LLR≥5) and homogeneous sites in the mitochondrial genome using a constructed site conservation database, tRNA database, and existing database and literature data (MT_ensGene, Mitimpact, and Mutpred). This facilitates the screening of mutations that are more likely to be associated with the disease and is used for subsequent analysis.

[0135] The site annotation process includes: annotating the conservation score of each heterogeneous and homogeneous site using the conservation database constructed in step 4.1; annotating the mutational harmfulness of heterogeneous and homogeneous sites on tRNA-encoding genes using the tRNA database constructed in step 4.2; and annotating the function of each heterogeneous and homogeneous site using existing database data and literature data (MT_ensGene, Mitimpact, and Mutpred).

[0136] 5. Gene marker analysis

[0137] 5.1 Single-point correlation analysis

[0138] By comparing the frequencies of heterogeneous or homogeneous loci in disease (early-onset coronary heart disease) samples with those in control normal individuals, loci with different frequencies are identified, and heterogeneous or homogeneous loci with higher frequencies in disease samples are used as mitochondrial gene markers.

[0139] Heterogeneous or homogeneous sites that are more frequent in control samples are used as irrelevant gene markers.

[0140] 5.1.1 Differential loci analysis between two groups of samples with no difference in sequencing depth

[0141] For association analysis of populations with high sequencing depth (greater than or equal to 150X) and no significant difference in sequencing depth between the disease (early-onset coronary heart disease samples) and the control (normal individual samples), two types of analysis are performed:

[0142] 1) Based on logistic regression or rank-sum test, calculate the loci (homogeneous or heterogeneous loci) where the heterogeneity level differs significantly between the case and control sample populations. These loci are candidate gene markers.

[0143] 2) Samples (individuals) with heterogeneity levels above a certain threshold (0.6 for early-onset coronary heart disease in this example) are considered to carry mutation sites. Sites with different carrying rates between cases and the control group are calculated, and these sites are candidate gene markers.

[0144] 5.1.2 Differential loci analysis between two groups of samples with different sequencing depths

[0145] When performing association analysis on populations with low sequencing depth (less than 150X) or significant differences in sequencing depth between cases and controls, two types of analysis are performed:

[0146] 1) Based on logistic regression and using the number of mitochondrial heterogeneity sites as a covariate for correction, sites with significant differences in heterogeneity levels between cases and controls (homogeneous or heterogeneous sites) were calculated as candidate gene markers; then, gene markers were determined based on the annotation results of the candidate gene markers.

[0147] 2) Based on variant detection methods, variant sites that differ between cases and controls are identified as candidate gene biomarkers. The gene biomarkers are then determined based on their annotation results. Differences in sequencing depth have a relatively small impact on variant detection, but the number of variant sites still needs to be corrected for as a covariate when performing association analyses.

[0148] Since the average sequencing depth of samples in the discovery phase was 183X and the average sequencing depth of samples in the validation phase was 121X, with a large difference in depth, only the data from the discovery phase were analyzed for heterogeneity sites, while the data from both phases were analyzed for homogeneity sites (Table 2).

[0149] In single-point association analysis, association analysis is performed only on functional sites within heterogeneous, homogeneous, or variable sites. Functional sites include those with a carrier proportion of less than 0.05 (a threshold for low-frequency carrier rate commonly used in genomics research) and located in mitochondrial genome protein-coding genes or in other regions (non-protein-coding gene regions) with a conservation score >0.6.

[0150] In the discovery phase, the association between heterogeneous loci and early-onset coronary artery disease (OCD) was investigated. Based on the rank-sum test and logistic regression, heterogeneous loci showing significant differences in heterogeneity levels between disease samples and normal control samples were calculated. The annotation results of these loci were analyzed to determine if they were gene biomarkers. However, due to the small sample size and low depth of analysis, no significantly different loci associated with early-onset COD were found (P < 0.05). Based on thresholds for different heterogeneity levels (heterogeneity level > 0.6, i.e., the proportion of mutated DNA greater than 60%), heterogeneous loci showing differences in carrier rates between case samples and normal control samples were calculated. The annotation results of these loci (including mutation type, harmfulness of software annotations, frequency in existing databases, and conservation status) were analyzed to determine if they were gene biomarkers. Fisher's exact test was performed, but no significantly different loci associated with early-onset COD were found (Fisher's exact test P < 0.05).

[0151] This study investigated the association between homogeneous loci and early-onset coronary artery disease (CAD). Based on the rank-sum test and logistic regression, in the discovery phase samples, homogeneous loci showing significant differences in heterogeneity between case samples and normal control samples were calculated. The annotation results of these loci (including mutation type, harmfulness of software annotations, frequency in existing databases, and conservation status) were used to analyze whether these homogeneous loci were gene biomarkers. The results showed that 9 homogeneous loci reached the significance threshold (Fisher's exact test, P < 0.05). In the validation phase samples, 5 homogeneous loci reached the significance threshold (Fisher's exact test, P < 0.05). Combining the data from both phases, three homogeneous loci (T16136C, A2363G, C15402T) reached the significant difference threshold in both phase samples (Fisher's exact test, P < 0.05 / 9) (Table 2), and can be considered as gene biomarkers associated with early-onset CAD.

[0152] Table 2: Results of Single-Point Association Analysis of Homogeneous Sites

[0153]

[0154] Note: inf represents infinity.

[0155] 5.2 Gene-level association analysis

[0156] For populations with high sequencing depth and no significant difference in sequencing depth between cases and controls, gene-level analysis is performed only on rare pathogenic heterogeneous or homogeneous sites. This involves calculating pathogenic heterogeneous or homogeneous sites where the heterogeneity level differs significantly between case samples and normal control samples. Pathogenic heterogeneous sites must meet the following criteria: mutations on protein-coding genes must be missense mutations with multiple software predictions of harmful or nonsense mutations, frameshift mutations, and splicing mutations; mutations on tRNA genes must be conserved (conservation score > 0.8) and alter base pairing (WC). Since cells only exhibit corresponding characteristics when the mutant / wild-type ratio exceeds a certain threshold, gene-level analysis is performed on heterogeneous sites that meet a certain threshold (0.6 for early-onset coronary heart disease in this example).

[0157] For protein-coding genes in the mitochondrial genome, loss-of-function site mutations must meet any one of the following conditions: 1) nonsense mutation, frameshift mutation, or splicing mutation; 2) missense mutations predicted as pathogenic by five software programs (CADD, PolyPhen, MutPred, MutationAssessor, SIFT). For tRNA genes in the mitochondrial genome, loss-of-function site mutations must simultaneously meet the following two conditions: 1) site annotation conservation score > 0.8; 2) site annotation is a mutation on the tRNA that causes the WC bond to be abolished or reconstructed. For rare pathogenic mutations (heterogeneity level of the variant site > 0.6), gene-level analysis was performed to identify susceptibility genes associated with early-onset coronary artery disease (ECD). The results showed that no significant genes associated with early-onset ECD were found (Fisher exact test, P < 0.05).

[0158] 5.3 Association Analysis Based on Complexes

[0159] Mitochondrial genes function in a tightly packed manner, with several genes forming a complex to perform their functions. Generally, mitochondria have four complexes, each encoded by a different gene, forming a protein complex that performs its function. A harmful mutation in any gene encoding this complex will affect its function. Therefore, association analysis is performed based on these complexes.

[0160] Mitochondrial protein-coding genes can be divided into four categories based on the proteins they encode: NADH reductases (ND1, ND2, ND3, ND4, ND4L, ND5, ND6), cytochrome b (CYB), cytochrome C (COX1, COX2, COX3), and ATPases (ATP6, ATP8). Different categories perform different functions.

[0161] For populations with high depth and no significant difference in depth between cases and controls, when performing complex / gene set-level analysis, association analysis is only performed on rare pathogenic heterogeneous sites. This involves calculating pathogenic heterogeneous or homogeneous sites where the heterogeneity level differs significantly between case samples and normal individual control samples. Pathogenic heterogeneous or homogeneous sites must meet the following criteria: nonsense mutations, frameshift mutations, or multiple software-predicted harmful missense mutations; and the tRNA gene must be conserved and have altered base pairing. Since cells only exhibit corresponding characteristics when the mutant / wild-type ratio exceeds a certain threshold, complex-level analysis is required for potential pathogenic heterogeneous / homogeneous sites that meet a certain threshold (0.6 for early-onset coronary artery disease in this example).

[0162] For protein-coding genes, loss-of-function site mutations must meet any one of the following conditions: 1) nonsense mutation, frameshift mutation, or splicing mutation; 2) mutations predicted as pathogenic by five software programs (CADD, PolyPhen, MutPred, MutationAssessor, SIFT). For tRNA genes, loss-of-function site mutations must simultaneously meet the following two conditions: 1) site annotation conservation score > 0.8; 2) site annotation is a mutation on tRNA that causes the WC bond to be abolished or reconstructed. For rare pathogenic mutations (heterogeneity level > 0.6), a complex-level analysis was performed to identify susceptibility complexes associated with coronary artery disease. The results showed that no significant complexes associated with early-onset coronary artery disease were found (Fisher's exact test, P < 0.05).

[0163] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0164] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A method of screening for mitochondrial genetic markers, characterized by, include: (a) Obtain the mitochondrial sequencing dataset of the training samples, which includes multiple abnormal group samples and multiple normal group samples; (b) Based on the mitochondrial sequencing dataset, determine the mitochondrial gene mutation information of the training samples, wherein the mitochondrial gene mutation information includes the level of gene mutation heterogeneity and at least one selected from gene mutation sites, base types, and gene mutation types; and (c) Based on the mitochondrial gene mutation information determined in step (b), select mitochondrial gene mutations that can be used to distinguish between the abnormal group samples and the normal group samples as the mitochondrial gene markers. The level of gene mutation heterogeneity was determined through the following steps: (1) The mitochondrial sequencing data is compared with the reference genome to determine the mitochondrial DNA sites and base types corresponding to the sequencing reads; (2) Based on the results of the alignment, determine the major base type of at least one of the mitochondrial DNA sites, wherein the major base type has the most supporting sequencing reads; (3) For a given base type at a given site, the heterogeneity level of the base type is determined by maximum likelihood estimation based on the relevant information of multiple sequencing reads corresponding to the given site.

2. The method according to claim 1, characterized in that, The relevant information includes: sequencing quality of sequencing reads, number of supporting sequencing reads for the given base type, and number of supporting sequencing reads for the major base type.

3. The method according to claim 2, characterized in that, In step (3), the level of heterogeneity of the given base at the given site is determined by maximum likelihood estimation using the following likelihood function: ； in, f is the level of heterogeneity for the given base type at the given site; l is the number of supporting sequencing reads of the given base type at the given site; i represents the i-th supporting sequencing read of the given base type at the given site; This represents the sequencing error rate of the i-th supporting sequencing read of the given base type at the given site; k represents the number of supporting sequencing reads for the major base type at the given site; j represents the j-th supporting sequencing read of the major base type at the given site; This represents the sequencing error rate of the j-th supporting sequencing read of the major base type at the given site.

4. The method according to claim 3, characterized in that, After determining the level of heterogeneity of the given base at the given site, a log-likelihood test is further performed on the level of heterogeneity.

5. The method according to claim 4, characterized in that, The log-likelihood ratio test is performed using the following formula: ； in, This represents the level of heterogeneity of the given base at the given site, as determined by the maximum likelihood estimation. This indicates the frequency when the given base at the given site is assumed to be a homogeneous mutation; This indicates the frequency of the non-dominant allele at the given locus in the mtDB database.

6. The method according to any one of claims 1 to 5, characterized in that, The heterogeneity level being no greater than 0.9 and the LLR score of the heterogeneous site being no less than 5 are indicators that the given site belongs to a heterogeneous site. A heterogeneity level greater than 0.9 indicates that the given site belongs to a homogeneous site.

7. The method according to claim 1, characterized in that, The gene mutation sites were obtained through conserved fractional base annotation at each site, mutation harmfulness annotation at heterogeneous sites of tRNA-encoding genes in the mitochondrial genome, and / or functional annotation information.

8. The method according to claim 7, characterized in that, The conservation score annotation information was obtained using site conservation database analysis.

9. The method according to claim 7, characterized in that, The mutation harmfulness annotation information was obtained using tRNA database analysis.

10. The method according to claim 7, characterized in that, The functional annotation information was obtained using known mitochondrial data analysis.

11. The method according to claim 1, characterized in that, The mitochondrial gene mutations that distinguish between the abnormal group samples and the normal group samples are determined by comparing the heterogeneity level of the gene mutation sites in different groups.

12. The method according to claim 1, characterized in that, The gene mutation type includes whether the gene mutation is a homogeneous mutation or a heterogeneous mutation.

13. A group of mitochondrial gene markers, characterized in that, The mitochondrial gene markers are obtained based on the method described in any one of claims 1 to 9.

Images