Model for calculating age based on excrement DNA and construction method

Abstract

The invention discloses a model for calculating age based on excrement DNA and a construction method, and belongs to the technical field of bioinformatics. According to the method, the problem that in the prior art, required materials are not suitable to obtain is solved by taking excrement as a material, and the method for obtaining the mtDNA heterogeneity site is disclosed for the first time, so that the target mtDNA heterogeneity site can be more accurately captured. A set of methodological process for screening age-related detection sites is constructed based on target mtDNA heterogeneity sites and machine learning, a model for calculating age based on mtDNA heterogeneity sites in excrement DNA is constructed, and the model is mainly applied to calculation of age of large animals or wild animals. The model solves the problem of high cost in the prior art.

Description

Technical Field

[0001] The invention belongs to the technical field of bioinformatics, and in particular relates to an age estimation method based on mtDNA heterogeneity. Background Art

[0002] Age is one of the basic parameters for population research and an important factor in animal ecology research. Age is of great indicative significance for determining life history characteristics such as birth rate and mortality rate of animal populations. Age usually includes actual age and biological age. The actual age refers to the total time experienced since birth, while the biological age is the age inferred based on the anatomical and physiological developmental status, reflecting the actual state of the body's tissue structure and physiological function. The degree of conformity between biological age and actual age reflects the interaction relationship and interaction mechanism between the body and the environment, and also plays a key role in population monitoring, development trend prediction, life history research, and protection and management decision-making. Therefore, it is particularly important to accurately determine the age of species.

[0003] However, there are no age records for animals in the wild. Traditionally, the age of wild animals is estimated mainly based on appearance morphology, anatomical indicators (growth rings and wear of hard tissues, bone morphology measurements, size of food cuts, growth curves), physiological and biochemical indicators (sex hormone levels, chromosome telomere length, DNA methylation levels) and omics indicators (transcriptome, proteome, metabolome, methylome and metagenome). The appearance morphology observation method is often limited to long-lived animals, whose body shape does not change much in adulthood, making it impossible to accurately determine their age. However, the other three methods are limited by the selection of species materials and usually rely on animal tissues, organs, tissues or blood. This type of material is difficult to collect and is restricted by many factors. It can only roughly estimate the age of animals (juvenile, sub-adult, adult and old). Therefore, traditional methods have obvious shortcomings in calculating wild animals or large animals.

[0004] Usually, feces is the most easily accessible material for wild animals. The continuous division and differentiation of stem cells and the shedding of intestinal epithelial cells constitute a highly dynamic intestinal epithelium. Intestinal epithelial cells are one of the cells with the highest cell turnover rate (they are completely renewed approximately every 3 to 5 days). After completing their mission, these daughter cells fall off into the intestinal cavity and are excreted from the body with feces. The DNA obtained from feces mainly comes from the shed intestinal epithelial cells, and the huge number of intestinal epithelial cells in feces (the number of cells shed by mice every 24 hours can reach 2×10 8 , while humans can reach 10 11 The existing method of using fecal DNA for methylation site detection is costly, so a method with lower cost needs to be developed.

[0005] Mitochondrial DNA (mtDNA) is the genetic material in mitochondria, which can generate energy (ATP) for cells. It is a special form of deoxyribonucleic acid found in the mitochondria of cells. Mitochondria are organelles that provide energy (ATP) for cells. There are generally multiple DNA molecules in a mitochondria. The mtDNA in feces mainly comes from shed intestinal epithelial cells. The intestinal epithelium not only has a huge surface area (up to 111 times (humans) to 119 times (mice) of the body surface area), but also has a very fast cell renewal rate (it will be completely renewed every 3 to 5 days). Taking humans as an example, each intestinal crypt contains 4 to 6 stem cells, and 300 daughter cells must be generated every day to maintain the integrity of the intestinal epithelium. Such a massive cell proliferation is also accompanied by a massive proliferation of mitochondria and their genomes, and the resulting mutations are accumulated in the daughter epithelial cells. After completing their mission, these daughter cells fall off into the intestinal cavity and are excreted from the body with feces. Therefore, using the huge number of intestinal epithelial cells in feces (the number of cells shed by mice every 24 hours can reach 2×10 8 , while humans can reach 10 11 The heteroplasmy of mtDNA can be detected.

[0006] The changes in mtDNA heteroplasmy and the accumulation of heteroplasmy are closely related to the age of individuals. In previous studies on humans, it was found that mtDNA heteroplasmy showed a downward trend with increasing age, with an average decrease of 0.4 copies per year. At the same time, heteroplasmic variants are usually more numerous and more frequent in the elderly. The mtDNA heteroplasmy of individuals over 70 years old is 58.5% higher than that of individuals under 40 years old on average. Moreover, this phenomenon of heteroplasmy being associated with increased age is not limited to the human species, but has also been verified in many other mammals. However, in these studies, the relationship between mtDNA heteroplasmy and age was not further explored, and then used to infer age. Therefore, in order to calculate the age of large animals (including wild animals) more conveniently, accurately and at a low detection cost, it is necessary to construct an age estimation method based on mtDNA heteroplasmy. Summary of the invention

[0007] The present invention solves the problem in the prior art that the required materials are difficult to obtain and the cost is high by constructing a computational model based on the correlation between mtDNA heteroplasmy sites and age in DNA in feces.

[0008] A method for constructing a model for calculating age based on fecal DNA, wherein the method is to construct a model for calculating age based on mtDNA heteroplasmy sites in fecal DNA.

[0009] The model is an elastic regression model, a support vector machine model, an XGBoost regression model or a random forest model.

[0010] The specific steps of the method are as follows:

[0011] Step 1: Extract genomic DNA from feces;

[0012] Step 2: construct a library based on genomic DNA and sequence it, then use trimmomatic to process the sequence in the library to obtain the processed library sequence;

[0013] Step 3: Obtain a complete report of mtDNA heterogeneity for each sample through the processed library sequence;

[0014] Step 4: Data processing of the complete report of mtDNA heteroplasmy for each sample;

[0015] Step 5: Screening target heterogeneous sites;

[0016] Step 6: Construct a machine learning model of age and target heterogeneous sites based on the xgboost package, randomForest package, e1071 package and glmnet package.

[0017] The nuclear CT value of the DNA in step 1 is less than or equal to 27.

[0018] The bacterial CT value of the DNA in step 1 is greater than 7.

[0019] The sequencing depth in step 2 is greater than or equal to 5000X.

[0020] The condition for screening the target heteroplasmic sites in step five is that the correlation between age and the target mtDNA heteroplasmic sites is limited to be greater than or equal to 0.5, or less than or equal to -0.5.

[0021] The data processing described in step 4 is to use the dplyr package to set the outliers to null values ​​according to the 3sigma principle, and use the KNN model to interpolate the missing values.

[0022] The specific steps for obtaining a complete report of mtDNA heterogeneity for each sample are:

[0023] Step 31: Use bwa to align the processed library data with the mtgenome genome, and then process it through samtools to obtain the mitochondrial related reading file mt1.bam;

[0024] Step 32: Convert mt1.bam to mt1.fq sequence and use trimmomatic to process the sequence, and then perform a second alignment with the Numt sequence to obtain the secondary alignment file nunt-sorted.bam, which is then deduplicated using samtools software.

[0025] Step 33: Convert the nunt-sorted.bam after deduplication to the nunt.fq sequence, and then perform trimmomatic processing and compare it with the mitochondrial genome for the third time to obtain numttomt.bam;

[0026] Step 3-4: Integrate the file mt1.bam in step 1 and the numttomt.bam in step 3, and remove the numt-sorted.bam in step 2 to obtain all mitochondrial readings;

[0027] Step 35: Use gatk to extract mtDNA heteroplasmy information from all mitochondrial reads, and use mutserve v2.0.1 to extract mtDNA heteroplasmy information from all mitochondrial reads. Integrate the heteroplasmy information obtained from the two to obtain the mtDNA heteroplasmy information report for each sample.

[0028] Step 36: Use bcftools to establish a consensus sequence, construct a circular sequence through Python processing, then build a library and sequence based on the circular sequence, then use trimmomatic to process the sequence in the library to obtain the processed library sequence, repeat steps 1 to 5 to obtain a new heterogeneity report for each sample, and finally integrate the new heterogeneity information of each sample with the mtDNA heterogeneity information report of each sample obtained in step 5 to obtain a complete report on the mtDNA heterogeneity of each sample.

[0029] A model for calculating age based on mtDNA heteroplasmy sites in fecal DNA, wherein the model is constructed by the above method.

[0030] A model for calculating age based on mtDNA heteroplasmy sites in fecal DNA and its application in calculating the age of large animals or wild animals.

[0031] Beneficial Effects

[0032] Compared with the prior art, the invention constructs a model for calculating age based on mtDNA heteroplasmic sites in fecal DNA, and the required materials will not cause damage to animals and are relatively easy to obtain.

[0033] The present invention discloses for the first time a method for obtaining mtDNA heteroplasmic sites, which can more accurately capture age-related mtDNA heteroplasmic sites.

[0034] The invention constructs a model for calculating age based on mtDNA heteroplasmic sites in fecal DNA, which has higher prediction accuracy than the prior art.

[0035] The present invention discloses for the first time a methodological process for screening age-related detection sites based on fecal mtDNA heterogeneity and machine learning. The detection cost is low and the accuracy is high. It can be applied to the screening of age-related heterogeneous sites in more animals, and also provides a reference for constructing age identification methods for other animals. DETAILED DESCRIPTION

[0036] The present invention is further described below in conjunction with implementation examples, but the present invention is not limited to the following implementation examples.

[0037] Implementation condition 1: The nuclear CT value of the feces DNA selected in the present invention does not exceed 27, and the bacterial CT value is greater than 7;

[0038] Implementation condition 2: Screening of target mtDNA heteroplasmic sites must be carried out according to the relevant steps described in the detailed method;

[0039] Implementation condition three: The correlation between the screening age and the target mtDNA heteroplasmic site in the present invention is limited to greater than or equal to 0.5 or less than or equal to -0.5.

[0040] Implementation condition four: The present invention constructs a computational model of the correlation between target mtDNA heteroplasmy sites and age based on the elastic regression model, support vector machine model, random forest model and eXtreme Gradient Boosting (XGBoost) regression model, and determines the optimal parameters and evaluation indicators of each model through leave-one-out cross-validation.

[0041] Implementation condition five: The average mitochondrial sequencing depth of the samples selected by the present invention must be greater than or equal to 5000X.

[0042] The present invention is further described in detail below in conjunction with specific implementation modes.

[0043] 1. Sample DNA extraction and library construction.

[0044] 1. Sample preparation.

[0045] The samples were 33 pieces of Siberian tiger feces from the Siberian Tiger Park and Hengdaohezi in Heilongjiang Province. They were numbered and recorded with the corresponding date of birth, sampling date, and gender, and stored in a -80℃ refrigerator.

[0046] 2. Sample DNA extraction, library construction and WGS sequencing.

[0047] The genomic DNA was extracted from 33 stool samples according to the method established by our research group (patent document: method recorded in CN113186185A), and the extracted DNA was annotated and stored at -20°C for future use. The qualified DNA in the 33 samples (DNA nuclear CT value not exceeding 27, bacterial CT value greater than 7) was constructed and the constructed library was sequenced. The sequencing depth was required to be greater than or equal to 5000X. The library construction and sequencing were commissioned to Shenzhen BGI. The detection cost of mtDNA heterogeneity for each sample was 648 yuan (the detection cost of methylation sites for each sample in the traditional method was 4224 yuan).

[0048] 2. Extract heterogeneous information.

[0049] 1. Find the Numts sequence.

[0050] The nuclear genome and mitochondrial genome sequences were obtained from NCBI, and the NUMTFinder V0.5.4 tool was used to obtain the Numts sequence of the sample through the nuclear genome and mitochondrial genome sequences.

[0051] 2. Library data processing.

[0052] Trimmomatic (version: 0.36) was used to process the sequenced data, including adapter trimming, quality filtering, and trimming low-quality bases, to improve the quality and reliability of the sequence data. The standard for processing sequencing data is to remove adapter sequences, with a maximum of 2 mismatches allowed and a threshold of 15 for adapter removal; use a sliding window for quality filtering with a window size of 5 bases, and perform trimming when the average quality value is less than 20; retain sequences with a length greater than or equal to 35 bases and an average quality greater than or equal to 20; do not trim the front and back ends of the sequence. The entire process is processed in parallel using two threads, unpaired sequences are discarded, and paired sequences are saved in the specified file.

[0053] 3. Three comparisons.

[0054] The purpose of the three alignments is to ensure that the Numts reads are removed while obtaining the mitochondrial reads, and to recover the mitochondrial reads that were misidentified as Numts.

[0055] 3.1. Obtaining preliminary mitochondrial readings.

[0056] Use bwa (version: 0.7.17-r1188) to align the processed library data with the mtgenome genome (from NCBI), and then process the aligned data through samtools (Version: 1.15.1) to obtain the mitochondrial related reading file mt1.bam.

[0057] 3.2. Obtain Numts readings in mitochondria.

[0058] The file mt1.bam (file of mitochondrial related readings) obtained in 3.1 was converted into the mt1.fq sequence, and then the mt1.fq sequence was trimmed by trimmomatic to remove low-quality sequences (trimming was performed when the average quality value was less than 20), and then aligned with the Numt sequence for the second time to obtain the secondary alignment file nunt-sorted.bam, which is the Numts readings contained in the first alignment, and was sorted and deduplicated by the samtools software.

[0059] 3.3. Recovering mitochondrial reads misidentified as Numts.

[0060] The nunt-sorted.bam (i.e., the Numts reads obtained from the second alignment) after deduplication was converted into the nunt.fq sequence. The nunt.fq sequence was trimmed by trimmomatic to remove low-quality sequences (trimming was performed when the average quality value was less than 20) and aligned with the mitochondrial genome for the third time to obtain numttomt.bam, and the mitochondrial reads that were misidentified as Numts were restored.

[0061] 3.4. Integrate bam files.

[0062] Integrate the mitochondrial readings mt1.bam obtained in step 3.1 with the recovered mitochondrial readings numttomt.bam obtained in step 3.3, and remove the Numts readings numt-sorted.bam in the second step. All mitochondrial readings can be obtained.

[0063] 4. Extract heterogeneous information.

[0064] The Mutect2 module in gatk (Version: 4.3.0.0) was used to extract the mtDNA heteroplasmy information of all mitochondrial reads obtained in 3.4. At the same time, mutserve v2.0.1 was used to extract the mtDNA heteroplasmy information of all mitochondrial reads obtained in 3.4. The heteroplasmy information obtained by the two was integrated to obtain the mtDNA heteroplasmy information report of each sample.

[0065] 5. Create a consensus sequence for each sample.

[0066] Use bcftools to establish a consensus sequence based on the heteroplasmy information of each sample, and construct a circular sequence through Python processing. This circular sequence can be used as a reference genome. Repeat steps 2 and 3 to obtain a new heteroplasmy report for each sample. Finally, the new heteroplasmy information of each sample and the heteroplasmy information generated in step 4 are reintegrated to obtain a complete report on the heteroplasmy of mtDNA for each sample.

[0067] 3. Data processing and screening of age-related mtDNA heteroplasmy.

[0068] 1. Data processing.

[0069] The dplyr package (R language) was used to integrate the mtDNA heterogeneity information of 33 samples, outliers were detected according to the 3sigma principle, outliers were set as null values, and the KNN model was used to interpolate missing values.

[0070] 2. Screening for age-related heterogeneity sites.

[0071] Based on Spearman correlation analysis, the conditions for screening the correlation between heterogeneous sites and age were limited to less than or equal to -0.5 or greater than or equal to 0.5.

[0072] 4. Model, evaluate and predict heterogeneous sites and age.

[0073] Four machine learning models for age and target heterogeneity sites were constructed based on the xgboost package, randomForest package, e1071 package, and glmnet package, namely elastic regression model, support vector machine model, XGBoost regression model, and random forest model. The leave-one-out cross-validation method was used to evaluate each model, and the evaluation index of each model (R 2 , MSE, RSME) and the optimal parameters of each model. At the same time, the correlation between the predicted age and the actual age of the four different models needs to be determined to select the optimal model. 2 It represents correlation, MAE represents the mean of absolute errors (hereinafter referred to as errors), and RMSE represents the root mean square error. The evaluation criteria are that the higher the correlation, the smaller the error, and the higher the accuracy of the model.

[0074] 5. Results

[0075] A total of 15 target heteroplasmic sites with a correlation with age greater than or equal to 0.5 or less than or equal to -0.5 were screened out through Spearman correlation analysis. In order to further verify the accuracy of the estimated age, the leave-one-out cross-validation method was combined with four machine learning models: elastic regression, support vector machine, XGBoost regression, and random forest to predict the age of the 15 points. The results are shown in Table 1. The results show that all four models can be used to construct a model for calculating age based on mtDNA heteroplasmic sites in fecal DNA.

[0076] When the elastic regression model is used to calculate age, in the training set, the correlation between the calculated age and the actual age can reach 0.859, with an error of only 6.18 months. In the test set, the correlation between the predicted age and the actual age can reach 0.878, with an error of only 13.87 months. When the support vector machine regression model is selected as the analysis model, in the training set, the correlation between the predicted age and the actual age can reach 0.995, with an error of only 1.93 months. In the test set, the correlation between the predicted age and the actual age can reach 0.805, with an error of only 11.01 months. When the random forest regression model is selected as the analysis model, the correlation between the predicted age and the actual age can reach 0.956 in the training set, with an error of only 5.41 months, and the correlation between the predicted age and the actual age can reach 0.811 in the test set, with an error of only 9.93 months; when the XGBoost regression model is selected as the analysis model, the correlation between the predicted age and the actual age can reach 0.987 in the training set, with an error of only 1.75 months, and the correlation between the predicted age and the actual age can reach 0.905 in the test set, with an error of only 8.91 months. Since the XGBoost regression model has the highest correlation and the smallest error, the XGBoost regression model is a superior model in the calculation of age based on fecal mtDNA heterogeneity.

[0077] Table 1 Regression evaluation indicators of different machine learning modeling

[0078]

[0079]

Claims

1. A method for constructing a model for calculating age based on fecal DNA, characterized in that: The method is to construct a model for calculating age based on mtDNA heteroplasmy sites in fecal DNA.

2. The method according to claim 1, characterized in that The model is an elastic regression model, a support vector machine model, an XGBoost regression model or a random forest model.

3. The method according to claim 1, characterized in that: The specific steps of the method are as follows: Step 1: Extract genomic DNA from feces; Step 2: construct a library based on genomic DNA and sequence it, then use trimmomatic to process the sequences in the library; Step 3: Obtain a complete report of mtDNA heterogeneity for each sample through the processed library sequence; Step 4: Data processing of the complete report of mtDNA heteroplasmy for each sample; Step 5: Screening target heterogeneous sites; Step 6: Construct a machine learning model of age and target heterogeneous sites based on the xgboost package, randomForest package, e1071 package and glmnet package.

4. The method according to claim 3, characterized in that The nuclear CT value of the genomic DNA in step 1 is less than or equal to 27.

5. The method according to claim 3, characterized in that: The bacterial CT value of the genomic DNA in step 1 is greater than 7.

6. The method according to claim 3, characterized in that The sequencing depth in step 2 is greater than or equal to 5000X.

7. The method according to claim 3, characterized in that The condition for screening the target heteroplasmic sites in step five is that the correlation between age and the target mtDNA heteroplasmic sites is greater than or equal to 0.5, or less than or equal to -0.

5.

8. The method according to claim 3, characterized in that The data processing described in step 4 is to use the dplyr package to set the outliers to null values ​​according to the 3sigma principle, and use the KNN model to interpolate the missing values.

9. The method according to claim 3, characterized in that: The method for obtaining a complete report of mtDNA heteroplasmy for each sample as described in step 3 is: Step 31: Use bwa to align the processed library data with the mtgenome genome, and then process it through samtools to obtain the mitochondrial related reading file mt1.bam; Step 32: Convert mt1.bam to mt1.fq sequence and use trimmomatic to process the sequence, and then perform a second alignment with the Numt sequence to obtain the secondary alignment file nunt-sorted.bam, and then use samtools software to remove duplicates; Step 33: Convert the deduplicated nunt-sorted.bam to the nunt.fq sequence, and perform a third alignment with the mitochondrial genome after trimmomatic processing to obtain numttomt.bam; Step 3-4: Integrate the file mt1.bam in step 1 and numttomt.bam in step 3, and remove numt-sorted.bam in step 2 to obtain all mitochondrial readings; Step 35: Use gatk and mutserve v2.0.1 to extract mtDNA heteroplasmy information from all mitochondrial reads simultaneously, integrate the heteroplasmy information obtained from the two, and obtain the mtDNA heteroplasmy information report for each sample; Step 36: Use bcftools and Python to construct circular sequences, then build libraries and sequence based on the circular sequences, then use trimmomatic to process the sequences in the library, and then repeat steps 31 to 35 to obtain a new heterogeneity report for each sample, and finally integrate it with the mtDNA heterogeneity information report for each sample obtained in step 5 to obtain a complete report on the mtDNA heterogeneity of each sample.

10. A model for calculating age based on mtDNA heteroplasmic sites in fecal DNA, characterized in that: The model is constructed by the method described in any one of claims 1 to 9.

Images