A method for constructing a metagenomic functional annotation correction model

Abstract

The invention discloses a method for constructing a metagenomic functional annotation correction model, relates to the technical field of metagenomic functional annotation, and discloses the following steps: step 1: dividing a reference amino acid sequence of known function according to the types of front-end deletion, back-end deletion, and double-end deletion; step 2: using an HMMsearch tool to import an incomplete amino acid sequence set into an HMM functional annotation model to be corrected for annotation, and obtaining scores, domain coverage, and sequence coverage results corresponding to the sequences; step 3: inputting four factors, namely, scores, domain coverage, sequence coverage results, and amino acid sequence deletion types, in the HMMsearch annotation results into a machine learning model as features, and taking whether the function of the reference amino acid is the same as that of the model as a category, performing model training, and finally obtaining a correction model of the HMM model for the incomplete amino acid sequence annotation results.

Description

Technical Field

[0001] The present invention relates to the technical field of metagenomic functional annotation, and more specifically, to a method for constructing a metagenomic functional annotation correction model. Background Art

[0002] Second-generation metagenomic sequencing generates paired reads of 150 bp in length. Metagenomic bioinformatics analysis first involves quality control of these sequences. Afterward, these reads are assembled into contigs using a splicing tool. Next, tools such as Prodigal (https: / / github.com / hyattpd / Prodigal) are used to predict genes and simultaneously obtain the corresponding amino acid sequences (often referred to as protein sequences). For functional annotation of these amino acid sequences, commonly used databases include eggNOG, GO, CAZy, MetaCyc, and KEGG. Common functional annotation software or methods fall into two main categories: sequence alignment, such as analysis software like diamond, blast, and mmseqs2; and HMM-based scoring. For example, KofamScan (https: / / doi.org / 10.1093 / bioinformatics / btz859) is often used for KEGG database annotation of amino acid sequences.

[0003] KofamScan uses the HMMER / HMMsearch tool to import sequence data into the Kofam (a custom HMM database for KEGG Orthologs) model for scoring, thereby completing the functional annotation of the KEGG database. This process scores the sequence against the HMM model of known KO (KEGG Ortholog) in the database. After obtaining the score, the score is compared with the preset confidence score threshold of the HMM model. If the sequence score exceeds the threshold, the sequence is annotated as the corresponding KO functional unit and assigned a KO number.

[0004] Compared to sequence alignment tools such as Diamond, BLAST, and mmseqs2, HMMsearch is more reliable for functional annotation. This is because these tools require a threshold, such as an E_value, to be set when performing sequence alignments. Setting a threshold is subjective and unvalidated, which can lead to uncertainty or bias in the alignment results. Furthermore, the E_value itself can vary depending on the size of the database. In contrast, HMMsearch uses an HMM model to annotate sequences by comparing the sequence score to a model confidence score threshold. If the score exceeds the threshold, the sequence is annotated. Because the model confidence score threshold is determined based on statistical information from large-scale validation data, the annotation results are more consistent and reliable. However, it is important to note that many amino acid sequences predicted from metagenomics are incomplete. Using incomplete sequences for HMM annotation can result in low model scores, which in turn fall below the annotation threshold, making functional annotation impossible for the sequence. Currently, there is no effective solution to this problem, which limits the HMMsearch annotation method. Summary of the Invention

[0005] To address this issue, we introduced a machine learning model as a corrector to calibrate the HMMsearch annotation results. Taking a certain HMM functional annotation model as an example, we artificially constructed an incomplete reference amino acid sequence set and imported four features from the HMMsearch annotation results of these sequences: score, domain_coverage, query_coverage, and amino acid sequence integrity information from Prodigal (named seq_miss_type, including the following four types: incomplete front end and complete back end, complete front end and incomplete back end, incomplete front end and incomplete back end, and complete front end and complete back end). We then trained the model based on whether the function of the reference amino acid was the same as that of the model, and finally obtained a correction model for the HMM model's incomplete amino acid sequence annotation results. Creating corresponding correction models for all HMM models that need to be annotated can not only improve the overall effect of metagenomic functional annotation, but also make the annotation results more comprehensive and reliable.

[0006] To achieve the above object, the present invention provides the following technical solutions:

[0007] A method for constructing a metagenomic functional annotation correction model comprises the following steps:

[0008] Step 1: The reference amino acid sequence with known function is segmented into the types of front-end deletion, back-end deletion, and double-end deletion to generate incomplete amino acid sequences;

[0009] Step 2: Use the HMMsearch tool to import the incomplete amino acid sequence set into the HMM functional annotation model to be corrected for annotation, and obtain the score, domain_coverage, and query_coverage results corresponding to the sequence;

[0010] Step 3: Input the four factors of score (score), domain coverage (domain_coverage), sequence coverage (query_coverage) and amino acid sequence missing type (seq_miss_type) in the HMMsearch annotation results into the machine learning model as features, and take whether the function of the reference amino acid is the same as the function of the model as the category (target) for model training, and finally obtain the correction model of the HMM model for incomplete amino acid sequence annotation results.

[0011] Furthermore, the segmentation method in step 1 is as follows: first, the length of the complete amino acid sequence is determined, and then the length of the sequence to be removed is calculated based on the proportion of the retained length of 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90%. Next, the sequence is truncated according to different segmentation types and excision lengths, and the same sequence produces an incomplete subsequence in each mode.

[0012] Compared with the prior art, the present invention has the following beneficial effects:

[0013] By importing four features from the HMMsearch annotation results—score, domain_coverage, query_coverage, and amino acid sequence completeness information from Prodigal (named seq_miss_type)—into the constructed machine learning model, a correction model for the HMM model was obtained for incomplete amino acid sequence annotation results. Creating a corresponding correction model for each HMM model that requires annotation not only improves the overall effectiveness of metagenomic functional annotation, but also significantly enhances the accuracy and sensitivity of the annotation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 A flow chart for constructing a correction model of the present invention;

[0015] Figure 2 This is a flow chart of the use of the calibration model of the present invention. DETAILED DESCRIPTION

[0016] Reference Figure 1 、 Figure 2 ,Hidden Markov model (HMM) is a statistical model used to model and analyze sequence data with potential hidden states.

[0017] HMMsearch is a tool used to search for Hidden Markov Models (HMMs) within a set of sequences. Its primary goal is to find sequences that match a specified HMM model within a given sequence database. HMM models are often used to represent specific protein families, functions, or domains. HMMsearch compares the input HMM model with sequences in the sequence database, calculates a similarity score between them, and selects the most relevant sequences.

[0018] KOfam (KEGG Orthology From Alignment) is a customized hidden Markov model (HMM) database used to annotate gene or protein sequence data. It is based on the organization and annotation standards of the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. KEGG is a widely used bioinformatics database that covers genes, proteins, metabolic pathways, and functional annotations across a wide range of biological fields. KOs are a set of functional annotations defined in the KEGG database that describe the same gene or protein family.

[0019] KofamScan uses the HMMER / HMMsearch tool to compare sequence data with Kofam (a custom HMM database for KEGG Orthologs) and assigns a KO number to it. The KO number is a unique identifier used to represent specific functional annotations in the KEGG database. These functional annotations include gene families, reactions, metabolic pathways, etc. KofamScan determines whether the sequence matches specific KOs based on the score of the HMM comparison and a predefined threshold. If the sequence comparison score is higher than the threshold, the sequence is considered to have the function represented by the KO, and the corresponding KO number is assigned to the sequence.

[0020] Random Forest is an ensemble learning machine learning model that combines multiple decision trees for prediction and classification tasks. Each decision tree is constructed independently and has no correlation with each other. This independence is achieved by generating a training set using random sampling. For each decision tree, Random Forest randomly selects a portion of the training samples and features for training. This reduces correlation between samples and increases model diversity. Building on the decision tree model, Random Forest employs two types of randomness: sample randomness and feature randomness. Sample randomness means that the training samples used by each decision tree are randomly sampled with replacement from the original data, meaning that a sample may appear repeatedly in the training sets of multiple decision trees. Feature randomness means that each decision tree randomly selects a portion of the available features for node splitting. This randomness reduces the risk of overfitting and improves the model's generalization ability. During the prediction phase, for classification tasks, Random Forest integrates the predictions of each decision tree using majority voting to determine the predicted class for the entire Random Forest. For regression tasks, Random Forest takes the average of the predictions of the decision trees as the final prediction.

[0021] The ROC curve (Receiver Operating Characteristic curve) is a commonly used tool for evaluating the performance of classification models. It is plotted with the True Positive Rate (TPR) on the vertical axis and the False Positive Rate (FPR) on the horizontal axis. Its graphical characteristic is that each point on the curve corresponds to the performance of the classification model at different thresholds. Ideally, a classifier should achieve a high TPR and a low FPR at all thresholds. Therefore, the closer the ROC curve is to the upper left corner, the better the model performance.

[0022] Area Under the ROC Curve (AUC). The AUC value represents the area under the ROC curve and ranges from 0.5 to 1, with values closer to 1 indicating better performance. AUC can be used to compare the performance of different classification models. It is generally believed that models with an AUC greater than 0.5 have better classification capabilities.

[0023] True Positive (TP): the model correctly predicts a positive sample as positive.

[0024] True Negative (TN): the model correctly predicts negative samples as negative.

[0025] False Positive (FP): The model mistakenly predicts a negative sample as positive.

[0026] False Negative (FN): The model mistakenly predicts a positive sample as negative.

[0027] Accuracy is a common metric used in model evaluation, measuring the accuracy of a classification model or retrieval system across all predictions. Accuracy describes the ratio of the number of examples correctly predicted by the model to the total number of examples, and ranges from 0 to 1. A higher accuracy indicates a model with high accuracy across all predictions, while a lower accuracy indicates a model with more errors in its predictions.

[0028] The accuracy calculation formula is: Accuracy = (TP + TN) / (TP + TN + FP + FN)

[0029] Precision is an important metric in model evaluation, used to measure the accuracy of a classification model or retrieval system in making positive predictions. Precision describes the model's ability to correctly predict positive results, and ranges from 0 to 1. A higher precision indicates a higher accuracy in the model's positive predictions.

[0030] Precision calculation formula: Precision = TP / (TP+FP)

[0031] Sensitivity is a metric used to evaluate the performance of a binary classification model. It measures the model's coverage of positive samples, that is, the proportion of true positive samples correctly detected by the model. A higher sensitivity indicates that the model is better able to capture true positive samples.

[0032] The sensitivity calculation formula is: Sensitivity = TP / (TP+FN)

[0033] Specificity is an indicator commonly used to evaluate the performance of binary classification models. It measures the model's coverage ability when predicting negative samples, that is, the proportion of true negative samples correctly detected by the model.

[0034] The specificity calculation formula is: specificity = TN / (TN+FP)

[0035] The F1 score is a commonly used metric for evaluating the performance of binary classification models. It combines the weighted average of precision and sensitivity. The F1 score ranges from 0 to 1, with values closer to 1 indicating better model performance.

[0036] F1 index calculation formula: F1 index = 2*(accuracy*sensitivity) / (accuracy+sensitivity)

[0037] Youden Index is an indicator used to evaluate the performance of a binary classification model. Its value ranges from -1 to 1. The closer the value is to 1, the better the model performance.

[0038] Youden index calculation formula: Youden index = sensitivity + specificity - 1

[0039] Metagenomics is a field that studies the genetic characteristics and functions of environmental microbial communities by utilizing high-throughput sequencing and related analytical techniques. It utilizes mixed DNA extracted from environmental samples, including DNA from bacteria, archaea, fungi, viruses, and other microorganisms, captured through sequencing. Then, through genome assembly and annotation, genetic information is obtained for individual microorganisms within the entire microbial community. The goal of metagenomics is to conduct in-depth studies of the composition, structure, functional potential, biodiversity, and interactions of microbial communities to deepen our understanding of the microbial world and reveal the effects of microorganisms on the environment and their hosts. Metagenomics has broad applications in environmental science, environmental ecology, microbial ecology, agriculture, medicine, and other fields.

[0040] To address the shortcomings of HMMsearch in the annotation of incomplete amino acid sequences, we proposed an improved method for machine learning models. This method imports four features from the HMMsearch annotation results, namely score, domain_coverage, query_coverage, and the missing type seq_miss_type of the incomplete amino acid sequence, into the machine learning model for model training and prediction to correct the annotation results of the HMMER / HMMsearch tool.

[0041] The following table shows the statistics of 20 functional annotations of the artificially constructed incomplete amino acid reference sequence before and after correction using the random forest model. HMM represents the model before correction, and ML represents the model after correction. Each row is the statistical data of the corresponding KO functional unit.

[0042] Table 1 Statistics of the number of true and false positives and negatives before and after the machine learning model correction

[0043]

[0044] It can be seen from Table 1 that the average number of true positives before correction of the machine learning model is 17440.1, and the average number of true positives after correction is 92858.1, which is significantly higher than before correction; the average number of false positives before correction is 84406.5, and the average number of false positives after correction is 8988.5, which is much lower than before correction.

[0045] Table 2 Statistics of accuracy, sensitivity and specificity of machine learning models before and after correction

[0046]

[0047]

[0048] Table 2 shows that the average accuracy of the machine learning model before correction was 88.60%, and the average accuracy after correction was 98.72%, which was significantly higher than before correction. The average sensitivity of the machine learning model before correction was 18.85%, and the average sensitivity after correction was 90.44%, which was much higher than before correction.

[0049] Table 3 Statistics of F1 index and Youden index before and after machine learning model correction

[0050]

[0051] It can be seen from Table 3 that the average F1 index of the machine learning model before correction was 30.18%, and the average F1 index after correction was 92.26%. The average F1 index after correction was significantly higher than before correction; the average Youden index before correction was 18.79%, and the average Youden index after correction was 89.71%. The average Youden index after correction was much higher than before correction.

[0052] Table 4 AUC values before and after machine learning model correction

[0053]

[0054] It can be seen from Table 4 that the AUC value of the machine learning model before correction is 92.72%, and the average AUC value after correction is 99.44%. The average AUC value after correction is significantly higher than that before correction.

[0055] By analyzing the number of true positives, false positives, accuracy, sensitivity, specificity, F1 index, Youden index, AUC value and ROC curve before and after machine learning model correction, we came to the following conclusions.

[0056] First, the number of true positives, accuracy, sensitivity, F1 index, Youden index, and AUC values after correction by the machine learning model were all higher than those before correction, demonstrating that the machine learning model effectively improved the accuracy of annotation results. By introducing four features (score, domain_coverage, query_coverage, and seq_miss_type) for model training and prediction, the functional type of amino acid sequences can be more accurately determined, thereby improving the quality of annotation.

[0057] Secondly, the number of false positives after the machine learning model correction was much lower than before correction, indicating that the machine learning model was able to identify and correctly annotate a large number of amino acid sequences that were originally misidentified after correction. This is crucial for increasing the completeness and reliability of the annotation results.

[0058] In summary, by introducing four features from the HMMsearch annotation results—score, domain_coverage, query_coverage, and amino acid sequence miss type seq_miss_type—into the machine learning model to correct the HMMsearch annotation results, we can effectively improve annotation accuracy and further enhance the completeness of the annotation results. Improvements in sensitivity, F1 index, Youden index, and AUC value further validate the effectiveness and reliability of the machine learning model in improving HMMsearch annotation results.

[0059] This paper aims to improve the accuracy and completeness of HMMsearch in annotating specific functional HMM models. We take the following steps to build a machine learning model:

[0060] Reference amino acid sequences with known functions are segmented according to leading, trailing, and paired-end deletion types. The segmentation method is as follows: First, the length of the complete amino acid sequence is determined. Then, the sequence length to be removed is calculated based on the percentage of retained length, which is 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90%. Next, sequence truncation is performed based on the segmentation type and the cutoff length. For example, using the 10% cutoff length, after calculating the sequence length to be removed, this length is removed from the front of the sequence to produce an incomplete sequence with a leading 10% missing sequence. This length is then removed from the back of the sequence to produce an incomplete sequence with a trailing 10% missing sequence. Finally, half of the sequence length is removed from the front and half of the sequence length is removed from the back, resulting in an incomplete sequence with a paired-end 10% missing sequence. This method simulates the incompleteness of amino acid sequences generated by metagenomic sequencing and simultaneously records the amino acid sequence missing type (seq_miss_type).

[0061] The reference amino acid sequence can be a sequence of corresponding function from the KEGG reference database, or an annotated amino acid sequence downloaded from NCBI or UniProt.

[0062] Use the HMMsearch tool to import the incomplete amino acid sequence set into the HMM functional annotation model that needs to be corrected for annotation, and obtain the corresponding score, domain_coverage, and query_coverage results of the sequence.

[0063] Among them, the HMMHMM functional annotation model can come from KOFAM, PFAM, etc.

[0064] The four factors of score, domain_coverage, query_coverage and amino acid sequence missing type seq_miss_type in the HMMsearch annotation results are input into the machine learning model as features, and whether the function of the reference amino acid is the same as the function of the model is used as the category (target) for model training to finally obtain the correction model.

[0065] The machine learning model may be any classification model, such as a random forest model and a support vector machine model, but is not limited to these two models.

[0066] The above description is merely a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiment. All technical solutions based on the concept of the present invention are within the scope of protection of the present invention. It should be noted that for those skilled in the art, certain improvements and modifications that do not depart from the principles of the present invention should also be considered within the scope of protection of this template.

[0067] The above is a detailed description of an embodiment of the present invention. However, the content described is only a preferred embodiment of the present invention and should not be considered to limit the scope of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the scope of the patent coverage of the present invention.

Claims

1. A method for constructing a metagenomic functional annotation correction model, characterized in that: The steps include: Step 1: The reference amino acid sequence with known function is segmented into the types of front-end deletion, back-end deletion, and double-end deletion to generate incomplete amino acid sequences; Step 2: Use the HMMsearch tool to import the incomplete amino acid sequence set into the HMM functional annotation model that needs to be corrected for annotation, and obtain the score, domain coverage domain_coverage, and sequence coverage query_coverage results corresponding to the sequence; Step 3: Input the four factors of the HMMsearch annotation results, namely score, domain coverage domain_coverage, sequence coverage query_coverage, and amino acid sequence missing type seq_miss_type, into the machine learning model as features, and use whether the function of the reference amino acid is the same as that of the model as the category target for model training, and finally obtain the correction model of the HMM model for incomplete amino acid sequence annotation results.

2. The method for constructing a metagenome functional annotation correction model according to claim 1, wherein The segmentation method in step 1 is as follows: first determine the length of the complete amino acid sequence, then calculate the length of the sequence to be removed based on the ratio of retained length to 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90%. Sequence truncation is performed according to different segmentation types and excision lengths. The same sequence produces an incomplete subsequence in each mode.

Images