Method for annotating and classifying euryascomycete transposons

Abstract

The application provides a method for annotating and classifying transposons of Erysiphaceae fungi, comprising the following steps: based on the annotation of TE domains of powdery mildew in multiple public databases, a transposon domain integration database is constructed; based on the transposon domain integration database, the domain characteristics of a target sequence are detected; based on the domain characteristics, a quantitative scoring system of a preset dimension is constructed, each dimension score is calculated, and a feature vector of each TE sequence is extracted; the random forest classifier is used to verify the dimension scores, and a transposon candidate sequence is obtained; the transposon candidate sequence is verified again to confirm the corresponding classification; based on the confirmed classification, a non-redundant transposon classification database for a specific species is constructed, and the random forest classifier is fed back for training iteration. The application realizes end-to-end automatic processing from the input of the original genomic sequence to the construction of the classification tree, and can improve the efficiency of large-scale genomic data analysis and the accuracy of transposon classification.

Description

Technical Field

[0001] This invention relates to the field of transposon annotation and classification technology, and in particular to a method for annotating and classifying transposons of Erythrocytesaceae fungi. Background Technology

[0002] Erythritolaceae ( Erysiphaceae Powdery mildew is an important class of obligate plant parasitic pathogens that can cause powdery mildew in various crops, resulting in serious economic losses to agricultural production. In recent years, with the increasing awareness of powdery mildew species (such as...) Blumeria graminis , Podosphaera xanthii , Erysiphe necator The completion of high-throughput genome sequencing revealed that transposons constitute an extremely high proportion of the powdery mildew genome, reaching 60%-85%, far exceeding that of typical fungal genomes. In the powdery mildew genome, some transposons have lost their characteristic sequences during long-term evolution, making them difficult to identify using traditional methods; a large number of transposons overlap or tandem, forming nested structures, further increasing the difficulty of identification; due to the high degree of genomic structural variation in powdery mildew, a large number of transposons are fragmented, with a low proportion of complete elements. These characteristics place higher demands on the annotation and classification methods for transposons in powdery mildew.

[0003] Transposable elements (TEs) are DNA sequences in the genome that can change their position, and they have a profound impact on genome evolution, gene expression regulation, and species diversity. Based on different transposition mechanisms, transposons can be mainly divided into type I transposons (reverse transposons) and type II transposons (DNA transposons).

[0004] Type I transposons employ a "copy-and-paste" mechanism: they transcribe DNA into mRNA, then reverse transcribe the mRNA into cDNA, and finally integrate the cDNA into a new location in the genome under the action of integrase. Based on the presence of long terminal repeats (LTRs) at both ends, type I transposons can be divided into LTR retrotransposons and non-LTR retrotransposons. LTR retrotransposons are further divided into Copia and Gypsy types, with differences in the enzyme coding sequence in their pol genes. Non-LTR retrotransposons mainly contain long stray repeat elements (LINEs) and short stray repeat elements (SINEs).

[0005] Type II transposons employ a "cut-and-paste" mechanism, meaning they are directly excised from one location and inserted into another. They encode inverted repeat sequences at both ends of specific transposase recognition elements. Common superfamilies include hAT, CACTA, and Mutator.

[0006] Existing transposon annotation and classification techniques are based on the Wicker six-level classification system (class, subclass, order, superfamily, family, subfamily). For example, transposon de novo prediction methods (EDTA, RepeatModeler2) predict based on sequence repetition and characteristic structure, but have a low recognition rate for fragmented or highly differentiated transposons that have lost their structural features. Homology alignment methods (Dfam, Pfam databases) detect homologous transposon sequences and conserved domains based on hidden Markov models, but if a single database is not comprehensive, a large number of sequences will not be identified or classified. Machine learning methods (TEclass, DeepTE) currently have limited training data for specific species, making it impossible to accurately identify and classify transposons in the genomes of specific species. Summary of the Invention

[0007] The purpose of this invention is to provide an annotation and classification method for transposons of powdery mildew fungi. By integrating multi-source public databases and a custom database of typical structural domains of powdery mildew transposons, a classification rule system is established, and a random forest model is introduced for verification. The classification results are automatically generated into a classification tree, and the final classification and sorting results are output.

[0008] To achieve the above objectives, the present invention provides a method for annotation and classification of transposons of Erythrinae fungi, the method comprising: S11. Based on multi-source public databases and annotation of transposon domains of powdery mildew, construct an integrated database of transposon domains; S12. Based on the integrated database of transposon subdomains, detect the domain features of the target sequence; S13. Based on the structural domain features, construct a quantitative scoring system with preset dimensions, and calculate the scores of each dimension of the structural domain features. S14. Based on the structural domain features, extract the feature vector of each transposon sequence; S15. The scores and feature vectors of each dimension are validated once by a random forest classifier to obtain transposon candidate sequences. S16. Perform secondary verification on the transposon candidate sequences to confirm the classification corresponding to the transposon candidate sequences; S17. Based on the confirmed classification, construct a non-redundant transposon subclassification database for specific species and feed it back to the random forest classifier for training iteration.

[0009] Furthermore, a transposon subdomain integrated database is constructed, specifically including: S21. From at least two different public structural domain databases, according to the preset lower limit threshold of bit fraction and the preset upper limit threshold of E value corresponding to each database, the structural domains related to transposons are selected to obtain the set of transposon-related structural domains. S22. Annotate the transposon domains of powdery mildew and extract candidate transposon sequences one by one to obtain the set of transposon domains of powdery mildew. S23. Based on the transposon-related domain set and the powdery mildew transposon domain set, construct a transposon domain integrated database, which includes a reverse transcriptase-specific library, a Gag multiprotein domain library, and a DNA transposase library.

[0010] Furthermore, the structural domain features of the target sequence are detected, specifically including: S31. Based on the transposon subdomain integrated database, the domain features of the target sequence are searched using a hidden Markov model to obtain the first interval file and candidate sequences, wherein the candidate sequences are sequences for which no domain was detected. S32. Based on the transposon subdomain integrated database, the domain features of the candidate order are searched using the BLAST model to obtain the second interval file; S33. Detect and label the structural domain features that overlap between the first interval file and the second interval file to obtain the overlapping region; S34. Merge the first interval file and the second interval file, and perform deduplication on the overlapping areas according to preset rules to obtain the structural domain features of the target sequence.

[0011] Furthermore, the target sequence is searched using Hidden Markov Models and BLAST Models, specifically including: S41. Set the sequence significance threshold, the single domain significance threshold, and the inclusion threshold containing domains. Search for the domain features of the target sequence using a hidden Markov model to generate the first interval file and candidate sequences. S42. Select candidate sequences whose length is less than a preset length threshold from the candidate sequences; S43. Set BLAST search parameters, including E-value threshold, word size threshold, gap opening penalty, and gap extension penalty. When the candidate sequence length is less than or equal to the first length threshold, the E-value threshold is set to the first baseline value. When the candidate sequence length is greater than the first length threshold, the E-value threshold is calculated using the following formula: E=L1×(1 + ln(L s / L0)) In the formula, E represents the E-value threshold; L1 represents the first reference value; L s L0 represents the candidate sequence length; L0 represents the baseline length. S44. The selected candidate sequences are searched using the BLAST model to generate a second interval file.

[0012] Furthermore, a quantitative scoring system with pre-defined dimensions is constructed, specifically including: S51. The core structural domain matching score is calculated by multiplying the ratio of the target sequence bit score to the reference maximum bit score by the first weight coefficient, and at least one type of conservative motif region is given an additional weight of a preset proportion during the scoring process. S52. The detected domain order is compared with a predetermined standard topology, and the domain order conformity score is calculated by multiplying the number of misalignments by a second weighting coefficient. The number of misalignments allows for positional shifts within a predetermined deviation range. S53. The domain feature matching score is calculated by multiplying the number of perfectly matched motifs and the number of mismatches by the third weighting coefficient. S54. The coverage and integrity scores are calculated by multiplying the ratio of the detected domain length to the reference model length by a fourth weighting coefficient, wherein the fourth weighting coefficient is assigned different weight values ​​according to the predetermined interval to which the coverage belongs.

[0013] Furthermore, the feature vector of each transposon subsequence is extracted, specifically including: S61. Statistically analyze the sequence length characteristics of each transposon subsequence; S62. Based on the sequence length characteristics, calculate the supplementary features, which include the GC content characteristics, AT content characteristics, CpG dinucleotide frequency characteristics, Shannon information entropy characteristics, and the normalized occurrence frequency of each 5-mer in each TE sequence. S63. Integrate all the supplementary features calculated in step S62 into a feature vector.

[0014] Furthermore, a random forest classifier is used to intelligently validate the scores and feature vectors for each dimension, specifically including: S71. Obtain positive transposon sequences of the target fungal group from a public transposon database, obtain negative non-transposon sequences from a public fungal genome database, and set the ratio of positive to negative samples to be basically balanced. S72. Input the scores and feature vectors of each dimension to be verified into the random forest model for verification, and output the classification probability distribution; S73. Set a probability threshold, and determine the category of the sequence to be verified based on the relationship between the classification probability distribution and the threshold to obtain transposon candidate sequences.

[0015] Furthermore, the candidate transposon sequences are subjected to secondary verification, specifically including: S81. Perform multiple sequence alignment on transposon candidate sequences and extract conserved functional region sequences from the alignment results. S82. Using a phylogenetic tree construction algorithm, the extracted conserved region sequence is used as input to automatically select an alternative model and perform a predetermined number of bootstrapping tests to construct a phylogenetic tree. S83. Verify the consistency of the classification results based on the clustering position or branching relationship of each transposon candidate sequence on the phylogenetic tree; S84. If a transposon candidate sequence clusters with known transposon candidate sequences in the same evolutionary branch, the classification corresponding to the transposon candidate sequence is confirmed. If it is isolated into a branch or clustered with non-transposon candidate sequences, it is marked as pending verification or a suspected new group.

[0016] Furthermore, a non-redundant transposon classification database for specific species is constructed, specifically including: Clustering and redundancy removal are performed on the set of transposon sequences that have been confirmed and classified, generating a non-redundant transposon classification database for a specific species or group. The clustering redundancy removal process adopts a preset clustering standard, which includes a sequence similarity threshold, a sequence coverage threshold, and a minimum sequence length threshold. Transposon sequences are classified into the same family only when the similarity between two transposon sequences to be clustered reaches or exceeds the sequence similarity threshold, the alignment coverage of the two transposon sequences both reaches or exceeds the sequence coverage threshold, and the length of the transposon sequence reaches or exceeds the minimum sequence length threshold.

[0017] Compared with the prior art, the beneficial effects of the present invention are: This invention provides an annotation and classification method for transposons in Erythrocytesaceae fungi. By integrating public databases and a custom core functional domain database, it significantly improves the detection rate of highly differentiated transposons in the Erythrocytesaceae genome and reduces the Unknown rate of traditional classification methods. Hidden Markov Models (HMMs) and BLAST models are commonly used search tools in bioinformatics analysis and are relatively easy to use. By setting genome-specific parameters for Erythrocytesaceae and dynamically adjusting the E-value threshold based on sequence length, the accuracy of HMM and BLAST model search and comparison is improved. By defining classification rules to comprehensively characterize transposon features, misjudgments and omissions caused by single indicators are avoided, especially for nested transposon structures and fragmented transposons commonly found in the Erythrocytesaceae genome, improving the sensitivity and accuracy of classification. The feature vector of each TE sequence is extracted, and it still maintains an excellent detection rate for fragmented and incomplete TE fragments caused by assembly breaks in next-generation sequencing. The introduction of a random forest mechanism can adaptively correct the bias of rule matching, improving classification accuracy and robustness.

[0018] Building a tree from the classification results visually reveals the evolutionary origins and differentiation history of transposons in different categories, helping users verify the reliability of the classification results and correct errors, and facilitating the discovery of new types of transposon families. Constructing a library from the classification results helps analyze the transposon characteristics in the genome of a specific species, or mask genomic transposon elements in the gene annotation process. Inputting the transposon library from the classification results into a random forest model for iterative training continuously improves the model's classification sensitivity and accuracy. This invention achieves end-to-end automated processing from inputting the original genome sequence to constructing the classification tree, and then fully automating the process of feeding the classification results back into the random forest model for iterative training, improving the efficiency of large-scale genomic data analysis and the accuracy of transposon classification. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0020] Figure 1 A schematic flowchart illustrating a method for annotating and classifying transposons of Erythrocytesae (powdery mildew fungi) provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the process for constructing a transposon subdomain integrated database provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the process for detecting the structural domain features of a target sequence according to an embodiment of the present invention; Figure 4 A schematic diagram illustrating the process of searching for a target sequence using a hidden Markov model and a BLAST model, provided in an embodiment of the present invention. Figure 5 This is a schematic diagram of the process for constructing a quantitative scoring system with preset dimensions provided in an embodiment of the present invention; Figure 6 This is a schematic diagram illustrating the process of extracting the feature vector of each TE sequence according to an embodiment of the present invention; Figure 7 This is a schematic diagram illustrating the intelligent verification process of scoring and feature vectors for each dimension using a random forest classifier, as provided in an embodiment of the present invention. Figure 8 This is a schematic diagram of the secondary verification process for transposon candidate sequences provided in an embodiment of the present invention. Detailed Implementation

[0021] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0022] Reference Figure 1 This embodiment provides a method for annotation and classification of transposons of Erythrinae fungi, the method comprising: S11. Based on multi-source public databases and annotation of the TE domain of powdery mildew, construct an integrated database of transposon domains.

[0023] S12. Based on the transposon subdomain integrated database, the domain features of the target sequence are detected.

[0024] S13. Based on the structural domain features, construct a quantitative scoring system with preset dimensions, and calculate the scores of each dimension of the structural domain features.

[0025] S14. Based on the structural domain features, extract the feature vector of each TE sequence.

[0026] S15. The scores and feature vectors of each dimension are validated using a random forest classifier to obtain candidate transposon sequences.

[0027] S16. Perform secondary verification on the transposon candidate sequences to confirm the classification corresponding to the transposon candidate sequences.

[0028] S17. Based on the confirmed classification, construct a non-redundant transposon subclassification database for specific species and feed it back to the random forest classifier for training iteration.

[0029] In this embodiment, firstly, transposon-related domains from multiple public databases such as Pfam, CDD, and TIGRFAM are integrated. Then, the TE domains of high-quality powdery mildew sequence sets are annotated, and candidate TE sequences are extracted one by one. Based on the transposon-related domains and candidate TE sequences, a dedicated transposon domain integrated database is constructed, expanding the database to include core functional domains specific to powdery mildew fungi. This fully utilizes public data resources and compensates for the insufficient coverage of powdery mildew transposons in general databases, significantly improving the sensitivity and specificity of subsequent domain detection. Next, based on the transposon domain integrated database, a hidden Markov model is used to scan the predicted protein sequences from the target genome for domain analysis. Short sequences are supplemented with BLAST for further detection, thereby comprehensively and sensitively capturing the characteristics of various transposon domains and avoiding missed detections due to excessive sequence length or differentiation.

[0030] Then, a quantitative scoring system is constructed based on four dimensions: core domain matching degree, domain arrangement order conformity, feature motif matching degree, and coverage and integrity. The detected domain features are weighted in multiple dimensions to obtain a feature system specific to powdery mildew transposon structures. This system can fully utilize the unique domain composition and arrangement patterns of transposons to effectively distinguish complete transposons, pseudogenes, and randomly similar sequences in the powdery mildew genome, reducing false positives. Subsequently, the four-dimensional scores, along with supplementary features such as local GC content bias and codon usage preference, are input into a random forest classifier. The trained model performs an intelligent verification of candidate sequences and outputs a "transposon," "pending," or "non-transposon" conclusion based on a probability threshold. This system can automatically learn the complex nonlinear relationships between multi-dimensional features, exhibiting stronger generalization ability and a lower false positive rate compared to fixed threshold rules. The "pending" mechanism also provides a buffer for manual review, ensuring the reliability of the results.

[0031] For transposon candidate sequences that have been validated in the first round, conserved regions (such as reverse transcriptase regions) are further extracted through multiple sequence alignment. A phylogenetic tree is constructed using IQ-TREE and subjected to bootstrapping testing. The classification results are then validated a second time based on the clustering positions of the sequences on the tree. This independent verification of classification rationality from an evolutionary perspective effectively eliminates misclassifications caused by horizontal transfer, gene duplication, or convergent evolution, playing a crucial role in identifying new families or reclassifications. Finally, the double-validated transposon sequences are clustered according to the "80-80-80" principle (≥80% similarity, ≥80% coverage, ≥80 bp length) to remove redundancy, generating a non-redundant transposon classification database for Erythrocyte Erythrocytesaceae fungi. This database is then fed back into a random forest classifier for training iterations. This process not only creates a species-specific, high-quality, non-redundant resource library but also continuously optimizes the classification model, allowing subsequent annotation work to improve with data accumulation. This achieves a closed-loop process from data integration, detection, scoring, intelligent validation, evolutionary confirmation to resource accumulation and model iteration.

[0032] As a preferred embodiment, constructing a transposon subdomain integrated database specifically includes: S21. From at least two different public structural domain databases, according to the preset lower limit threshold of bit fraction and the preset upper limit threshold of E value corresponding to each database, the structural domains related to transposons are selected to obtain the set of transposon-related structural domains. S22. Annotate the TE domains of powdery mildew and extract candidate TE sequences one by one to obtain the set of transposon domains of powdery mildew. S23. Based on the transposon-related domain set and the powdery mildew transposon domain set, construct a transposon domain integrated database, which includes a reverse transcriptase-specific library, a Gag multiprotein domain library, and a DNA transposase library.

[0033] In this embodiment, refer to Figure 2 Based on the integration of multiple public databases such as Pfam, CDD, and TIGRFAM, this study focuses on expanding the core functional domains of transposons specific to fungi and powdery mildew families. Transposon-related domains (such as PF00078 transposase and PF00595 reverse transcriptase) were extracted from the Pfam database, with the selection criteria being a bit score greater than or equal to 50 and an E value threshold less than or equal to 1e-10. Transposon-related domains (such as cd00236 reverse transcriptase) were extracted from the CDD database, with the selection criteria being a bit score greater than or equal to 40 and an E value threshold less than or equal to 1e-5. Transposon-related domains were extracted from the TIGRFAM database, with the selection criteria being a bit score greater than or equal to 55 and an E value threshold less than or equal to 1e-5. Domain annotation was performed on a high-quality TE sequence set of powdery mildew, and candidate TE sequences were extracted one by one to obtain a sample set of typical transposon domains in the powdery mildew genome. The transposon domains from the multiple databases and the sample set of typical transposon domains of powdery mildew were then integrated into a dedicated domain integration database.

[0034] Based on the extracted transposon-related domains, a custom core functional domain library was constructed, including a reverse transcriptase-specific library, a Gag multi-protein domain library, and a DNA transposase library. Specifically, the reverse transcriptase-specific library includes frequently occurring Copia and Gypsy-type reverse transcriptase-specific motifs (such as the YMDD motif) in the Ascomycetes or Powdery mildew family genomes; the Gag multi-protein domain library includes Gag precursor proteins and their protease cleavage sites (such as cysteine / histidine motifs) unique to Ascomycetes or Powdery mildew family genomes; and the DNA transposase library includes DDE catalytic core domains of DNA transposons active in the Ascomycetes or Powdery mildew family genomes, such as hAT, Mutator, and CACTA. Redundant data was removed during database construction based on 90% sequence similarity and 80% sequence coverage, and a unified naming convention was adopted, such as using the format "DB_Superfamily_Domain Type_Version"; simultaneously, a mapping table between domain families and transposon superfamilies / families was established.

[0035] As a preferred embodiment, the detection of structural domain features of the target sequence specifically includes: S31. Based on the transposon domain integration database, the domain features of the target sequence are searched using a hidden Markov model to obtain the first interval file and candidate sequences, wherein the candidate sequences are sequences for which no domain was detected.

[0036] S32. Based on the transposon subdomain integrated database, the domain features of candidate sequences are searched using the BLAST model to obtain the second interval file; The target sequence is searched using Hidden Markov Models and BLAST models, specifically including: S41. Set the sequence significance threshold, the single domain significance threshold, and the inclusion threshold for domains. Search the target sequence using a hidden Markov model to generate the first interval file and candidate sequences.

[0037] S42. Select candidate sequences whose length is less than a preset length threshold from the candidate sequences; S43. Set BLAST search parameters, including E-value threshold, word size threshold, gap opening penalty, and gap extension penalty. When the candidate sequence length is less than or equal to the first length threshold, the E-value threshold is set to the first baseline value. When the candidate sequence length is greater than the first length threshold, the E-value threshold is calculated using the following formula: E=L1×(1 + ln(L s / L0)) In the formula, E represents the E-value threshold; L1 represents the first reference value; L s L0 represents the candidate sequence length; L0 represents the baseline length.

[0038] The first length threshold is determined based on the length distribution characteristics of validated transposon sequences from powdery mildew fungi. Specifically, it is determined by statistically analyzing the coverage and length distribution of core structural domains in high-confidence transposon sequences, selecting a threshold that preserves the complete transposon structural features while reducing non-specific matches in short fragments. The first benchmark value is determined based on the structural domain coverage distribution characteristics of complete transposon sequences in the high-confidence transposon dataset of powdery mildew fungi, and is used to screen candidate transposon sequences with high structural integrity. The benchmark length is predetermined during the method establishment phase based on the statistical analysis results of the high-confidence transposon dataset of powdery mildew fungi. When the length of a candidate sequence is lower than the benchmark length, it may not be possible to accurately determine its transposon type due to insufficient number of structural domains or incomplete structural information, and therefore it will not be considered a complete candidate transposon sequence for subsequent classification.

[0039] S44. The selected candidate sequences are searched using the BLAST model to generate a second interval file.

[0040] S33. Detect and label the structural domain features that overlap between the first interval file and the second interval file to obtain the overlapping region; S34. Merge the first interval file and the second interval file, and perform deduplication on the overlapping areas according to preset rules to obtain the structural domain features of the target sequence.

[0041] In this embodiment, refer to Figure 3 and Figure 4To scan the protein sequences predicted from the powdery mildew genome using a Hidden Markov Model (HMM), the following commands were executed: A transposon domain integrated database was used as the query to search the input sequences. The E-value thresholds for sequences, domains, and containing domains were set to 0.01. For short sequences (e.g., <300 bp) not detected by the HMM, a BLAST model was used for supplementary detection. This involved searching the input sequences using a transposon domain integrated database, with the E-value threshold set to 1e-5, word size to 3, gap opening penalty to 8, and gap extension penalty to 2.

[0042] In steps S31 and S32, the E-value threshold can be dynamically adjusted according to the sequence length to improve matching accuracy. The adjustment algorithm is as follows: if the sequence length is less than or equal to the first length threshold of 1000, then the E-value threshold is set to 1e-5; otherwise, the E-value threshold is set to 1e-5 * (1 + ln(L s / 1000), where 1000 is the baseline length set in this embodiment, used to dynamically adjust the E value threshold, L s The candidate sequence length is specified. Simultaneously, in step S33, a custom script or bedtools can be used to detect and label overlapping structural domains: using the Hidden Markov Model (HMM) search results as the first interval file and the BLAST model search results as the second interval file, interval overlap comparisons are performed. An overlap ratio threshold of 0.5 is set, meaning that at least 50% of the interval length in the HMM result must overlap with the interval length in the BLAST result to be considered a match.

[0043] Output all records in the Hidden Markov Model (HMM) results (regardless of overlap), along with the overlapping BLAST results (the corresponding field will be empty if there is no overlap) and the actual number of overlapping bases. Merge the HMM and BLAST search results, sort them based on the highest bit score, and for overlapping regions, use a strategy of "retaining the highest bit score + highest coverage" to remove duplicates.

[0044] As a preferred embodiment, a quantitative scoring system with preset dimensions is constructed, specifically including: S51. Calculate the core structural domain matching score by multiplying the ratio of the target sequence's bit score to the reference maximum bit score by the first weighting coefficient, and give an additional weight to at least one type of conservative motif region by a preset proportion during the scoring process.

[0045] S52. The detected domain order is compared with a predetermined standard topology, and the domain order conformity score is calculated by multiplying the number of misalignments by a second weighting coefficient. The number of misalignments allows for positional shifts within a predetermined deviation range.

[0046] S53. The domain feature matching score is calculated by multiplying the number of perfectly matched motifs and the number of mismatches by the third weighting coefficient.

[0047] S54. The coverage and integrity scores are calculated by multiplying the ratio of the detected domain length to the reference model length by a fourth weighting coefficient, wherein the fourth weighting coefficient is assigned different weight values ​​according to the predetermined interval to which the coverage belongs.

[0048] In this embodiment, refer to Figure 5 Based on the detected domain features, a four-dimensional quantitative scoring system is constructed: 1. Core structural domain matching degree The matching degree algorithm is defined as follows: Score_match (match score) = (bitscore_i / max_bitscore) × 0.3 Where bitscore_i is the bit score of sequence i, and max_bitscore is the maximum bit score of the corresponding domain in the training set. 0.3 represents the weight coefficient of the domain matching feature in the quantitative scoring system, used to measure the importance of domain matching information for transposon classification. This weight coefficient is determined by comprehensively considering the contribution of different features in the validated transposon dataset to classification during the construction of the powdery mildew fungal transposon classification model. In calculating the matching score using the matching algorithm, a weighted increase of 1.5 can be applied to the reverse transcriptase YMDD motif and the DNA transposase DDE motif regions. Finally, the Z-score normalization method is used to normalize the original scores.

[0049] 2. Consistency of the order of structural domains The standard topology of the LTR retrotransistor of the white powder family is defined as follows: Copia class: 5'-LTR → Gag → Protease → RT → RNase H → 3'-LTR Gypsy class: 5'-LTR → Gag → Protease → RT → RNase H → Integrase →3'-LTR Among them, LTR (long terminal repeat); Gag (Gag protein); Protease; RT (reverse transcriptase); RNase H (RNA endonuclease); and Integrase (integrase).

[0050] Define the standard topology of DNA transposons, which is the relative position of the terminal inverted repeat (TIR) ​​sequence and the internal transposase.

[0051] Define the scoring algorithm: Score_arrangement = 1 - (number of misalignments × 0.1) × weight Here, "1" represents the maximum sequence matching score when the domain arrangement in the candidate transposon sequence completely conforms to the standard structural pattern of the corresponding transposon type, i.e., the upper limit of the normalized score under ideal conditions. The misalignment refers to the number of amino acids misaligned from the standard topological structure, allowing for a spacing deviation of ±20%. "0.1" represents the base penalty coefficient for domain misalignment, used to measure the degree of influence of each domain position deviation on the sequence score. The weighting is defined as follows: for the reverse transcriptase YMDD motif region, weight = 1.0; for the Gag protease cleavage site region, weight = 0.8; for other regions, weight = 0.5.

[0052] 3. Structural domain feature matching degree Define the domain ordering matching algorithm: Score_feature = (Number of perfect matching radix - Number of mismatches × 0.2) × 0.2 In the formula, the number of perfectly matched motifs is the number of motifs in the candidate sequence that are consistent with the reference motif and whose arrangement order and relative position are normal; the number of misaligned motifs is the number of motifs in the candidate sequence whose position or arrangement order is offset relative to the reference structural pattern; the number of misaligned motifs × 0.2 is a penalty term, where 0.2 is the penalty coefficient; the other 0.2 is the weight of the structural feature scoring term in the quantitative scoring system.

[0053] The feature score parameters of the domain motif matching algorithm were optimized specifically for the YMDD motif. A perfect match of the YMDD motif earns 1 point, a mismatch deducts 0.1 points, and a maximum of two mismatches deduct 0.2 points. A perfect match of the Gag protease cleavage site earns 1 point, and a mismatch of one amino acid deducts 0.3 points. If the LTR terminal TSD sequence is detected, 1 point is earned, and if not detected, 0 points are earned.

[0054] 4. Coverage and completeness Define coverage and integrity algorithms: Score_coverage (coverage and integrity score) = (0.5 × (detection length / model length) + 0.5 × (number of detected core domains / number of predefined core domains for this transposon subtype)) × weight In the formula, the detection length represents the coverage of the core structural domains actually detected in the candidate transposon sequence, obtained through the structural domain detection step. The model length represents the sum of the theoretical core structural domain lengths in the complete structural model of the corresponding transposon type. The two "0.5"s represent the weights of coverage and completeness, respectively.

[0055] The weighting is as follows: complete components (coverage ≥ 80%) have a weight of 0.3; partially incomplete components (30% ≤ coverage < 80%) have a weight of 0.2; and extremely incomplete components (coverage < 30%) have a weight of 0.1. Complete components represent core structural domains that are basically complete and can represent typical TE structures; partially incomplete components represent components with some missing structural domains but still have clear TE characteristics; and extremely incomplete components represent components that retain only a small number of structural fragments and have severely degraded structures.

[0056] As a preferred embodiment, the feature vector of each TE sequence is extracted, specifically including: S61. Statistically analyze the sequence length characteristics of each TE sequence.

[0057] S62. Based on the sequence length characteristics, calculate the supplementary features, which include the GC content characteristics, AT content characteristics, CpG dinucleotide frequency characteristics, Shannon information entropy characteristics, and the normalized occurrence frequency of each 5-mer in each TE sequence.

[0058] S63. Integrate all the supplementary features calculated in step S62 into a feature vector.

[0059] In this embodiment, the specific method for extracting the feature vector of each TE sequence is shown in Table 1: Table 1. Calculation method for extracting the feature vector of each TE sequence.

[0060] Specifically, refer to Figure 6The sequence length feature L of each TE sequence is calculated, which directly counts the total number of A, T, C, and G bases in each TE sequence. The GC content feature is calculated as GC = (G+C) / (A+T+G+C). The overall GC content of the powdery mildew genome is 37-43%, and different TE superfamilies exhibit characteristic differences in GC preference. The AT content feature is calculated as AT = (A+T) / (A+T+G+C), which complements the GC content to enhance the distinguishing ability of transposons of different taxa. The CpG dinucleotide frequency feature is calculated, reflecting the base combination preference within the transposon sequence. The Shannon information entropy feature is calculated as H = -Σ(pi × log2(pi)), where pi is the frequency of each A, T, C, and G base in the sequence. This feature comprehensively reflects the repetition degree of the sequence and the uniformity of base distribution. Constructing a 1024-dimensional high-dimensional feature vector of 5-mer frequencies: The system enumerates all 4^5=1024 5-mer combinations (AAAAA to TTTTT), and counts the normalized frequency of each 5-mer in the sequence. The above features are then combined to construct a 1029-dimensional feature vector.

[0061] As a preferred embodiment, a random forest classifier is used to intelligently validate the scores and feature vectors for each dimension, specifically including: S71. Obtain positive transposon sequences of the target fungal group from a public transposon database, obtain negative non-transposon sequences from a public fungal genome database, and set the ratio of positive to negative samples to a basic balance.

[0062] S72. Input the scores and feature vectors of each dimension to be verified into the random forest model for verification, and output the classification probability distribution.

[0063] S73. Set a probability threshold, and determine the category of the sequence to be verified based on the relationship between the classification probability distribution and the threshold to obtain transposon candidate sequences.

[0064] In this embodiment, refer to Figure 7 A random forest classifier is constructed and trained to intelligently validate the four-dimensional rule matching results and feature vectors. Transposon sequences of Ascomycetes or Erythrocytesaceae fungi are extracted from the Dfam database (positive samples, transposons); fungal housekeeping genes (such as act1, tub2, etc.) are downloaded from the MycoCosm database (negative samples, non-transposons); the ratio of positive to negative samples is set to approximately 1:1 to ensure class balance. The four-dimensional rule feature vectors Score_match (matching score), Score_arrangement (arrangement score), Score_feature (feature score), and Score_coverage (coverage score) and the feature vectors are input into the random forest model.

[0065] The random forest model parameters are set as follows: n_estimators=500, meaning 500 decision trees are built to balance accuracy and computational efficiency; max_features='sqrt', meaning the number of features considered when splitting a tree is the square root of the total number of features; max_depth=30, meaning the maximum tree depth is set to 30 to prevent overfitting; min_samples_split=2, meaning a node must contain at least 2 samples before splitting; Criterion='gini', meaning the Gini coefficient is used as the node splitting criterion. StratifiedKFold(n_splits=5) is used to divide the sample dataset into 5 stratified subsets, with each stratified subset serving as a fold for cross-validation. Because the stratified partitioning maintains consistency between the proportions of different classes in each fold and the original dataset, each fold has a class distribution similar to the overall dataset, making the model evaluation results more stable and reliable. n_splits=5 specifies the number of folds to split. Here it is set to 5, which means that the data is divided into 5 parts. Each time, 4 parts are used for training and 1 part is used for validation. A total of 5 training and validation cycles are performed.

[0066] The four-dimensional classification rule weights are defined as follows: core structural domain matching degree = 0.3; structural domain arrangement order matching degree = 0.3; structural feature matching degree = 0.25; coverage and completeness = 0.15. The four-dimensional features and feature vectors are input into the trained random forest model, which outputs a classification probability distribution. A threshold of 0.85 is set. If the maximum probability > 0.85, the classification is adopted; if it is between 0.6 and 0.85, it is marked as "pending" and manual review is requested; if it < 0.6, it is determined as "non-transposon" or "suspected new family". See Table 2 for details. Table 2. Detailed Parameter Table of Transposon Classification Rules

[0067] As a preferred embodiment, secondary verification of the transposon candidate sequence is performed, specifically including: S81. Perform multiple sequence alignment on the candidate transposon sequences and extract the conserved functional region sequences from the alignment results.

[0068] S82. Using a phylogenetic tree construction algorithm, the extracted conserved region sequence is used as input to automatically select an alternative model and perform a predetermined number of bootstrapping tests to construct a phylogenetic tree.

[0069] S83. Verify the consistency of the classification results based on the clustering position or branching relationship of each transposon candidate sequence on the phylogenetic tree.

[0070] S84. If a transposon candidate sequence clusters with known transposon candidate sequences in the same evolutionary branch, the classification corresponding to the transposon candidate sequence is confirmed. If it is isolated into a branch or clustered with non-transposon candidate sequences, it is marked as pending verification or a suspected new group.

[0071] In this embodiment, refer to Figure 8 The candidate transposon sequences were subjected to multiple sequence alignment to extract the conserved reverse transcriptase region sequence. The IQ-TREE algorithm was used for phylogenetic tree construction, with the conserved reverse transcriptase region sequence as the alignment file. ModelFinder was enabled to automatically select the best alternative model, and 1000 ultra-fast bootstrapping tests were performed. The classification results were validated based on the position of each transposon classification result on the phylogenetic tree.

[0072] After confirming the classification of the transposon candidate sequences, adjacent TE prediction fragments belonging to the same category are merged to restore the genomic boundaries of the complete TE. All high-confidence TE prediction fragments are sorted according to their genomic coordinates (sequence ID and start position). The sorted fragments are traversed; if adjacent fragments have the same predicted TE category and adjacent genomic coordinates, they are merged into the same TE merge group. The TE interval of the merge group is the minimum start coordinate and maximum end coordinate of all fragments within the group, and the TE confidence of the merge group is the coverage length-weighted average of the predicted probabilities of each fragment within the group. If adjacent fragments predict different TE categories, they are retained as independent TE sites and not merged. Finally, a TE annotation result file in standard bioinformatics format is output.

[0073] As a preferred embodiment, a non-redundant transposon classification database for a specific species is constructed, specifically including: Clustering and redundancy removal are performed on the set of transposon sequences that have been confirmed for classification, generating a non-redundant transposon classification database for specific species or taxa.

[0074] The clustering redundancy removal process uses preset clustering criteria, which include sequence similarity threshold, sequence coverage threshold, and minimum sequence length threshold.

[0075] Transposon sequences are classified into the same family only when the similarity between two transposon sequences to be clustered reaches or exceeds the sequence similarity threshold, the alignment coverage of the two transposon sequences both reaches or exceeds the sequence coverage threshold, and the length of the transposon sequence reaches or exceeds the minimum sequence length threshold.

[0076] In this embodiment, the classified sequence set is clustered to remove redundancy, generating a non-redundant transposon subclassification database for a specific species. The clustering and redundancy removal process uses preset clustering criteria, defined as follows: classified sequences are clustered into family levels based on ≥80% sequence similarity (identity), ≥80% length coverage, and a minimum sequence length of ≥80 bp (the 80-80-80 principle). The clustering method sets a sequence similarity threshold of 0.8, meaning that sequences are only grouped into one class when their similarity reaches 80% or more. Simultaneously, both longer and shorter sequences are required to have an alignment coverage of at least 0.8, ensuring that the alignment length of the sequences meets the coverage requirements during clustering.

[0077] Based on the transposon annotation and classification method for powdery mildew fungi constructed above, this study was conducted using three representative powdery mildew fungi (Gracilaria spp.). Erysiphe necator Powdery mildew Blumeria graminis f. sp tritici Powdery mildew of cucurbits Podosphaera xanthii The method was tested on the genome and compared with the current mainstream transposon annotation tools EDTA and RepeatModeler. The proportion of the total length of transposon sequences predicted by each method to the total length of the corresponding species' genome was calculated as an indicator of transposon genome coverage. To further evaluate the statistical significance of the differences between different methods, the genomes of each species were divided into fixed-length windows. Whether different methods detected transposon regions within each window was counted, and paired significance analysis was performed using the McNemar test based on the window-level detection results. P < 0.05 indicated a significant difference between the two methods. The results showed that the proposed method was significantly higher than existing tools in transposon detection rate (Table 3), while the proportion of unknown sequences was significantly lower than existing tools (Table 4). This indicates that the proposed method can not only identify more transposon sequences, but also more accurately classify them into known transposon families, thereby significantly reducing the proportion of unclassified sequences. Specific data are as follows: Table 3 Comparison of transposon detection rates

[0078] Table 4 Comparison of Unknown Rates for Categories

[0079] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for annotation and classification of transposons of Erythrinae fungi, characterized in that, The method includes: S11. Based on multi-source public databases and annotation of transposon domains of powdery mildew, construct an integrated database of transposon domains; S12. Based on the integrated database of transposon subdomains, detect the domain features of the target sequence; S13. Based on the structural domain features, construct a quantitative scoring system with preset dimensions, and calculate the scores of each dimension of the structural domain features. S14. Based on the structural domain features, extract the feature vector of each transposon sequence; S15. The scores and feature vectors of each dimension are validated once by a random forest classifier to obtain transposon candidate sequences. S16. Perform secondary verification on the transposon candidate sequences to confirm the classification corresponding to the transposon candidate sequences; S17. Based on the confirmed classification, construct a non-redundant transposon classification database for powdery mildew and feed it back to the random forest classifier for training iteration; Construct a quantitative scoring system with predefined dimensions, specifically including: S51. The core structural domain matching score is calculated by multiplying the ratio of the target sequence bit score to the reference maximum bit score by the first weight coefficient, and at least one type of conservative motif region is given an additional weight of a preset proportion during the scoring process. S52. The detected domain order is compared with a predetermined standard topology, and the domain order conformity score is calculated by multiplying the number of misalignments by a second weighting coefficient. The number of misalignments allows for positional shifts within a predetermined deviation range. S53. The domain feature matching score is calculated by multiplying the number of perfectly matched motifs and the number of mismatches by the third weighting coefficient. S54. The coverage and integrity scores are calculated by multiplying the ratio of the detected domain length to the reference model length by the fourth weighting coefficient. The fourth weighting coefficient is assigned different weight values ​​according to the predetermined interval to which the coverage belongs. Secondary verification of transposon candidate sequences is performed, specifically including: S81. Perform multiple sequence alignment on transposon candidate sequences and extract conserved functional region sequences from the alignment results. S82. Using a phylogenetic tree construction algorithm, the extracted conserved region sequence is used as input to automatically select an alternative model and perform a predetermined number of bootstrapping tests to construct a phylogenetic tree. S83. Verify the consistency of the classification results based on the clustering position or branching relationship of each transposon candidate sequence on the phylogenetic tree; S84. If a transposon candidate sequence clusters with known transposon candidate sequences in the same evolutionary branch, the classification corresponding to the transposon candidate sequence is confirmed. If it is isolated into a branch or clustered with non-transposon candidate sequences, it is marked as pending verification or a suspected new group.

2. The annotation and classification method for transposons of Erythrinae fungi according to claim 1, characterized in that, Constructing a transposon subdomain integrated database specifically includes: S21. From at least two different public structural domain databases, according to the preset lower limit threshold of bit fraction and the preset upper limit threshold of E value corresponding to each database, the structural domains related to transposons are selected to obtain the set of transposon-related structural domains. S22. Annotate the transposon domains of powdery mildew and extract candidate transposon sequences one by one to obtain the set of transposon domains of powdery mildew. S23. Based on the transposon-related domain set and the powdery mildew transposon domain set, construct a transposon domain integrated database, which includes a reverse transcriptase-specific library, a Gag multiprotein domain library, and a DNA transposase library.

3. The annotation and classification method for transposons of Erythrinae fungi according to claim 1, characterized in that, Detection of structural domain features of the target sequence, specifically including: S31. Based on the transposon subdomain integrated database, the domain features of the target sequence are searched using a hidden Markov model to obtain the first interval file and candidate sequences, wherein the candidate sequences are sequences for which no domain was detected. S32. Based on the transposon subdomain integrated database, the domain features of candidate sequences are searched using the BLAST model to obtain the second interval file; S33. Detect and label the structural domain features that overlap between the first interval file and the second interval file to obtain the overlapping region; S34. Merge the first interval file and the second interval file, and perform deduplication on the overlapping areas according to preset rules to obtain the structural domain features of the target sequence.

4. The annotation and classification method for transposons of Erythrinae fungi according to claim 3, characterized in that, The structural domain features of the target sequence are searched using Hidden Markov Models and BLAST models, specifically including: S41. Set the sequence significance threshold, the single domain significance threshold, and the inclusion threshold containing domains. Search for the domain features of the target sequence using a hidden Markov model to generate the first interval file and candidate sequences. S42. Select candidate sequences whose length is less than a preset length threshold from the candidate sequences; S43. Set BLAST search parameters, including E-value threshold, word size threshold, gap opening penalty, and gap extension penalty. When the candidate sequence length is less than or equal to the first length threshold, the E-value threshold is set to the first baseline value. When the candidate sequence length is greater than the first length threshold, the E-value threshold is calculated using the following formula: E=L1×(1 + ln(L s / L0)) In the formula, E represents the E-value threshold; L1 represents the first reference value; L s L0 represents the candidate sequence length; L0 represents the baseline length. S44. The selected candidate sequences are searched using the BLAST model to generate a second interval file.

5. The annotation and classification method for transposons of Erythrinae fungi according to claim 1, characterized in that, Extracting the feature vector of each transposon subsequence, specifically including: S61. Statistically analyze the sequence length characteristics of each transposon subsequence; S62. Based on the sequence length characteristics, calculate the supplementary features, which include the GC content characteristics, AT content characteristics, CpG dinucleotide frequency characteristics, Shannon information entropy characteristics, and the normalized occurrence frequency of each 5-mer in each TE sequence. S63. Integrate all the supplementary features calculated in step S62 into a feature vector.

6. The annotation and classification method for transposons of Erythrinae fungi according to claim 1, characterized in that, Intelligent validation of scores and feature vectors across various dimensions is performed using a random forest classifier, specifically including: S71. Obtain positive transposon sequences of the target fungal group from a public transposon database, obtain negative non-transposon sequences from a public fungal genome database, and set the ratio of positive to negative samples to be basically balanced. S72. Input the scores and feature vectors of each dimension to be verified into the random forest model for verification, and output the classification probability distribution; S73. Set a probability threshold, and determine the category of the sequence to be verified based on the relationship between the classification probability distribution and the threshold to obtain transposon candidate sequences.

7. The annotation and classification method for transposons of Erythraceae fungi according to claim 1, characterized in that, Constructing a non-redundant transposon classification database for powdery mildew, specifically including: Clustering and redundancy removal are performed on the set of transposon sequences that have been confirmed and classified to generate a non-redundant transposon classification database for powdery mildew. The clustering redundancy removal process adopts a preset clustering standard, which includes a sequence similarity threshold, a sequence coverage threshold, and a minimum sequence length threshold. Transposon sequences are classified into the same family only when the similarity between two transposon sequences to be clustered reaches or exceeds the sequence similarity threshold, the alignment coverage of the two transposon sequences both reaches or exceeds the sequence coverage threshold, and the length of the transposon sequence reaches or exceeds the minimum sequence length threshold.

Images