A method for screening cross-species HGT

Abstract

The present invention belongs to the field of bioinformatics technology and specifically relates to a method for screening cross-species HGT. The method includes the steps of constructing a background database, inferring homology groups, constructing a gene phylogenetic tree, inferring HGT based on the gene tree, verifying candidate HGT inferred based on the gene tree, and assigning HGT to a timeline. By integrating a multi-step process such as phylogenetic analysis, sequence alignment, gene function annotation, and statistical verification, the method achieves high-throughput, automated screening of HGT events in large genetic datasets. The method of the present invention is particularly suitable for studying HGT across extremely long time periods and extremely long evolutionary distances across the biological kingdom, and can provide an important scientific tool for understanding gene flow and functional enhancement in the course of biological evolution.

Description

Technical Field

[0001] The present invention belongs to the technical field of bioinformatics, and particularly relates to a method for screening cross-species HGT. Background Art

[0002] Horizontal gene transfer (HGT), also known as lateral gene transfer (LGT), refers to the exchange of genetic material between two species that are not vertically related. This process is common in prokaryotes and occurs through mechanisms such as conjugation, transduction, and transformation. Unlike vertical gene transfer, which occurs from parent to offspring, HGT transcends the boundaries of kinship, significantly enriching the dynamics of gene flow.

[0003] The history of observations of HGT dates back to 1959, when it was documented that high-frequency transduction (Hfr) Escherichia coli could laterally transfer genetic information to a specific mutant form of Salmonella typhimurium. That same year, Tomochiro Akiba and Kunitaro Ochiai discovered resistance plasmids in pathogenic bacteria and subsequently demonstrated that these plasmids could be transferred between different bacterial strains. However, the concept of HGT had not yet been established. It was not until the 1990s, with the emergence of genetically modified organisms (GEOs), particularly genetically modified microorganisms (GEMs), and the emergence of numerous drug-resistant pathogens, whose origins could no longer be attributed solely to genetic mutations, that the concept of HGT gradually gained attention and became a research hotspot.

[0004] Similar phenotypes observed in genetically distant organisms are often attributed to shared genes in their genomes. Loss of these genes in closely related lineages may result from multiple independent gene loss events or from HGT between different lineages. In theory, HGT is possible between any two organisms with a DNA genome. Unlike vertical gene transmission, which is the cornerstone of biological heritage conservation and the stability of the eukaryotic tree of life, HGT is a major force driving the diversification of eukaryotic species and their adaptation to diverse environments. Although HGT is widely recognized as a major evolutionary force in prokaryotes, its role in eukaryotic evolution remains controversial, primarily due to the complex evolutionary history, complex genome structure, and frequent sequence contamination.

[0005] To accurately identify HGT events, there is an urgent need to develop efficient and accurate screening methods. Currently, common HGT analysis methods include phylogenetic tree analysis, base composition analysis, selection pressure analysis, intron analysis, specific sequence analysis, and nucleotide composition bias analysis. However, existing HGT screening methods, which mostly rely on phylogenetic analysis, suffer from limitations such as incomplete database coverage, immature automatic phylogenetic tree selection algorithms, and a lack of unified analysis standards. Improvement and innovation are urgently needed. Summary of the Invention

[0006] The purpose of the present invention is to remedy the deficiencies of the prior art and provide a method for efficiently and accurately screening cross-species HGT.

[0007] To achieve the above objectives, the present invention adopts a technical solution: a method for screening cross-species HGT, comprising the following steps:

[0008] Step 1: constructing a background database, which contains known protein sequences covering the phylogenetic lineages of prokaryotes and eukaryotes, and removing highly homologous and redundant sequences;

[0009] Step 2: Perform BLAST comparison on any two protein sequences in the background database, calculate the homology group association network based on the best mutual alignment in each pair of gene combinations, and cluster the homology group association network using an unsupervised hidden Markov model (HMM) clustering algorithm to generate a database containing multiple homology groups;

[0010] Step 3: For each query protein sequence to be analyzed, search for homologous sequences in the background database, perform strict sequence alignment and pruning, construct an accurate single-gene phylogenetic tree, and evaluate the support of each branch of the evolutionary tree through at least one of the ultra-fast bootstrap test, similar likelihood ratio approximation test, approximate Bayesian test, and fast local bootstrap probability test;

[0011] Step 4: Search for nested topological structures between the query sequence and genetically distant sequences in the single-gene phylogenetic tree constructed with the query sequence as the center. Screen for potential cross-species HGT by defining nested positions and setting a minimum support node threshold. Simultaneously, record the most similar monophyletic matching sequences of the query sequence in both the recipient and donor groups.

[0012] Step 5: Evaluate the confidence of potential cross-species HGT by calculating the outlier index, HGT score support index, HGT branch length support index and its consensus hit support index;

[0013] Step 6: Estimate the time of occurrence of cross-species horizontal transfer events by tracing the minimum taxonomic boundaries of cross-species HGT offspring and the minimum taxonomic boundaries of donor offspring; combine molecular clock methods with archaeological and fossil evidence to provide a temporal context for cross-species HGT.

[0014] Preferably, in step 1, the data in the background database include all available protein sequences of the reference sequences in NCBI and protein sequences from the marine microbial eukaryotic transcriptome sequencing project.

[0015] Preferably, in step 2, multiple sequence alignment and conserved domain analysis are performed on the sequences in each homology group.

[0016] Preferably, step five further comprises analyzing potential cross-species HGT using homology comparison features of genes above and below the HGT gene locus, specifically comprising:

[0017] If the potential cross-species HGT is located in a chromosome segment and 50% of the best matches of the genes in the chromosome segment are from other biological kingdoms, the HGT candidate gene is eliminated; or

[0018] If the potential cross-species HGT is located in a chromosome segment and 50% of the genes in the chromosome segment are identified as HGT genes, the HGT candidate gene is eliminated; or,

[0019] If at least one of the three closest upstream and downstream genes for potential cross-species HGT had a best hit in another kingdom, it was excluded; or

[0020] When two or more potential cross-species HGTs are physically closely linked and belong to the same gene family, they will be regarded as a single cross-species HGT event, and all these adjacent potential cross-species HGTs will be retained.

[0021] Compared with the prior art, the method of the present invention has the following beneficial effects:

[0022] (1) Efficiency: Through automated and standardized processes, high-throughput screening of HGT events in large genetic datasets is achieved, greatly improving analysis efficiency;

[0023] (2) Accuracy: The accuracy of HGT screening is ensured through a multi-step process including building a large background database, inferring homologous groups, and constructing a phylogenetic tree. At the same time, the introduction of multiple validation indicators and taxonomic analysis of physical flanking genes further reduces the false positive rate.

[0024] (3) Wide applicability: This method is particularly suitable for HGT studies in complex biological communities such as the plant kingdom, and can provide an important tool for understanding gene flow and functional enhancement in the process of biological evolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 A schematic diagram of the system architecture and module composition used in the method for screening cross-species HGT in an embodiment of the present invention;

[0026] Figure 2 Schematic diagram of the method for screening cross-species HGT in an embodiment of the present invention;

[0027] Figure 3The topological structure of the evolutionary nodes of 113 species and the number of HGTs obtained at each node;

[0028] Figure 4 Figure 3. Evolutionary trends of codon bias of HGT of prokaryotic (Pr) and eukaryotic (Eu) origin over time for 113 species. The codon bias index of each species was normalized relative to its respective core gene, and the x-axis represents the temporal order of HGT integration into plant lineages. (a)-(b) are the amino acid lengths of HGT of prokaryotic (Pr) and eukaryotic (Eu) origin, respectively; (c)-(d) are the amino acid lengths of HGT of prokaryotic (Pr) and eukaryotic (Eu) origin, respectively. Amino acid aromaticity; (e)-(f) are the GC contents of the third position of synonymous codons in HGT of prokaryotic (Pr) and eukaryotic (Eu), respectively; (g)-(h) are the codon adaptation indices in HGT of prokaryotic (Pr) and eukaryotic (Eu), respectively; (i)-(j) are the optimal codon frequencies in HGT of prokaryotic (Pr) and eukaryotic (Eu), respectively; (m)-(n) are the effective codon numbers in HGT of prokaryotic (Pr) and eukaryotic (Eu), respectively;

[0029] Figure 5 Comparison of codon bias in core genes (CORE) and horizontal gene transfer (HGT); (a)-(b) are the number of HGT effective codons (Nc) of prokaryotic (Pr) and eukaryotic (Eu) origins of Rhodophyta organisms, respectively; (c)-(d) are the codon bias index (CBI) of HGT of prokaryotic (Pr) and eukaryotic (Eu) origins of Rhodophyta organisms, respectively; (e)-(f) are the codon bias index (CBI) of HGT of prokaryotic (Pr) and eukaryotic (Eu) origins of Chlorophyta organisms, respectively. (g)-(h) are the HGT codon bias index (CBI) of prokaryotic (Pr) and eukaryotic (Eu) origins of Chlorophyta, respectively; (i)-(j) are the HGT codon bias index (Nc) of prokaryotic (Pr) and eukaryotic (Eu) origins of Streptophyta, respectively; (k)-(m) are the HGT codon bias index (CBI) of prokaryotic (Pr) and eukaryotic (Eu) origins of Streptophyta, respectively. DETAILED DESCRIPTION

[0030] The present invention is further described below with reference to specific embodiments and the accompanying drawings. More details are set forth in the following description to facilitate a full understanding of the present invention. However, the present invention can obviously be implemented in a variety of other ways different from the description. For those skilled in the art, any replacement, improvement or transformation made to the embodiments of the present invention is within the scope of protection of the present invention, and the scope of protection of the present invention should not be limited by the content of this specific embodiment.

[0031] This paper provides a method for screening cross-species HGT to facilitate large-scale, high-throughput HGT detection and explore the breadth of HGT-driven biological evolution in the Tree of Life (TOL). This paper applies this method to all species, generates a list of candidate HGT genes, and reveals the distribution patterns of HGT within and between kingdoms. By arranging HGT in chronological order, their contribution to the complexity and diversification of biological systems is inferred, providing research with the possibility of verifying numerous evolutionary theories. The system used in this method (HGTStart) is as follows Figure 1 As shown, the process is as Figure 2 As shown, the following steps are included:

[0032] 1. Build a background database

[0033] To conduct efficient HGT surveys, we included as many high-quality protein sequences as possible from whole-genome data across a wide range of prokaryotic and eukaryotic lineages, as well as individual sequences from major protein databases.

[0034] We examined the genome assembly summary information from RefSeq (ftp: / / ftp.ncbi.nlm.nih.gov / genomes / README_assembly_summary.txt) (updated as of May 2021) to download all available proteins from genomes completed before the update date. We grouped genomes from species within the same genus and downloaded only the genome representing the least fragmented assembly within that group (if multiple versions of the genome were available, the latest version was selected). We also searched for genomes from the JGI (https: / / genome.jgi.doe.gov / portal / ) and other databases.

[0035] A significant issue affecting HGT inference is uneven sampling due to imbalanced data collection across different taxa. This study incorporates proteins from the MMETSP database to compensate for the limited availability of red algal genomic data. This extensive genomic data results in a protein database containing 17,250,679 protein sequences from 1,157 genomes, with reasonable coverage across most lineages in the tree of life (540 bacteria, 45 archaea, 431 metaflagellates, 15 red algae, 83 green algae, and 43 genomes from the vesicular lineage). Each of these 1,157 complete genomes represents a representative species from its genus. This database has been designated "GNM1157."

[0036] In addition to GNM1157, we constructed a larger database encompassing as many known protein sequences as possible. We constructed this background database by downloading the NCBI RefSeq database (release 82, ftp: / / ftp.ncbi.nlm.nih.gov / refseq / ). Given the underrepresentation of algal lineages in RefSeq, we also included enriched algal protein data from the MMETSP project and other public sources to expand the taxonomic span of the background database.

[0037] The collected protein sequences (including GNM1157, Refseq, and MMETSP data) were integrated into the main database REFAL, and highly similar sequences (sequence identity ≥ 90%) within each taxonomic group (such as Brassicaceae or Primates) were removed using CD-HIT version 4.5.4. Ultimately, a protein database containing 39.9 million sequences from more than 7786 taxa was constructed, ensuring reasonable coverage of most lineages in the tree of life.

[0038] To facilitate function prediction of query sequences, the present invention predicted the protein sequence of GNM1157 by accessing large databases such as Interproscan, EGGNOG, PANTHER, Pfam, and SUPERFAMILY. The query sequence was associated with its gene hits in GNM1157 to retrieve the corresponding predicted functions.

[0039] 2. Homology group inference and benchmarking

[0040] Homology between genes, including orthologs and paralogs, is a cornerstone of comparative biology studies, such as horizontal gene transfer analysis. Tree-based databases such as PhylomeDB, Ensembl-Compara, EggNOG, and TreeFam have been developed in this area. However, we found that these existing databases do not fully cover the entire gene set encoded by GNM1157. Therefore, the construction of a customized homology database is crucial. We clustered the 17,250,679 proteins encoded by GNM1157 and combined them with protein-based topological analysis. All protein sequences were exhaustively analyzed using the genome-wide all-against-all BLAST workflow (version 2.2.28, with an e-value threshold of 1e-10 and a local identity threshold of 20%). Subsequently, an association network of homology groups was calculated using OrthoFinder V2.3.7 software based on the reciprocal best-length-normalized hits (RBNH) for each genome pair. OrthoFinder excels at correcting for the dependence of gene similarity on gene length and phylogenetic distance, thereby improving the accuracy of orthologous group delineation. Furthermore, an unsupervised Markov clustering algorithm (MCL3.0) was independently run to cluster this orthologous group map. The MCL algorithm was run multiple times using a set of gradient inflation parameters (1.2, 1.4, 1.6, 1.8, 2.0, and 2.2), and the best run was selected. Local downloads of more than a dozen popular online orthologous databases (e.g., EGGNOG, PANTHERN, and SUPERFAMILY) were used as benchmarks to calibrate the orthologous inference parameters. Precision, recall, and F0 values ​​for the gradient set were calculated using the NumPy library in Python 3.0 to determine the optimal value for the inflation factor.

[0041] Where TP is the number of true positive homology group assignments (i.e., correct assignments), FP is the number of false positive homology group assignments (i.e., incorrect assignments), and FN is the number of false negative homology group assignments (i.e., missed assignments). The F-score is the harmonic mean of precision and recall. It balances these two metrics and is more sensitive to lower values, making it useful for comprehensively evaluating the performance of homology group assignment.

[0042] To identify biologically significant homology clusters, we tested the performance of group-specific conserved domains within each orthologous group (OG) and eliminated noisy and weakly associated sequences. First, we performed multiple alignments of sequences within each OG using MAFFT v7.455. Next, we constructed an HMM profile for conserved domains for each OG using HMMER 3.0. We then used the Hmmsearch program to search all protein sequences within each OG against their corresponding HMM profiles, using an E-value threshold of 0.00001. Sequences below the threshold were removed from the group and pooled for a second round of sequence homology construction. After this round of screening, each sequence was accurately assigned to an orthologous group. Ultimately, this clustering process generated 27,631 orthologous groups, each containing at least two members. For orthologous groups with five or more members, we further examined their phylogenetic relationships by constructing maximum likelihood trees. Specifically, the sequences were aligned using MUSCLE version 3.8.31 with default settings, and the alignments were then trimmed using TrimAI version 1.232 in automatic mode (-automated1), retaining an alignment length of at least 50 amino acids.

[0043] These trimmed alignments (≥50 amino acids) were used to construct phylogenetic trees using the "WAG+CAT" model in FastTree version 2.1.7. Four rounds of minimal evolutionary SPR (subtree pruning and rejoining) moves (-spr 4) and exhaustive maximum likelihood nearest neighbor exchange (-mlacc 2 -slownni) were performed. Branch support was estimated using the Shimodaira-Hasegawa (SH) test. The successful reconstruction of phylogenetic trees using all members of the OG demonstrates that the strategy of clustering homologous groups based on sequence similarity is strongly supported by the phylogenetic strategy. These trees clearly distinguish nested positions formed by taxonomically distant sequences at variable evolutionary rates, thus elucidating gene relationships through horizontal gene transfer (HGT).

[0044] 3. Constructing a phylogenetic tree

[0045] The present invention designed a "RoutineTree" model for unsupervised homology search, sequence alignment, tree construction and tree screening to predict HGT for each given protein sequence. In short, RoutineTree searches for homologous sequences in REFAL and constructs a single-gene phylogenetic tree for each custom query protein sequence. To speed up the search process, the HMM model and probabilistic inference method integrated in HMMER3.0 are used to construct database segmentation based on functional conserved domains. 159,424 profiles from the PFAM database were scanned using hmscan, and the REFAL database was divided into families, and then sequences were retrieved using esl-fetch. Each family is named with its corresponding PFAM number to facilitate subsequent retrieval. Excluded REFAL sequences are collected into a FASTA file.

[0046] Initially, the query sequence is split into individual sequence files in the first step and scanned for PFAM HMM profiles using hmmscan. A temporary BLAST / diamond database is constructed based on the PFAM hit number for each split query sequence, from which the corresponding FASTA files in REFAL are retrieved and merged. After creating a search index in BLAST / diamond format, the temporary database is searched using a default e-value threshold of 1e-05. For each query, the top 10,000 significant hits are recorded in descending order of positional score (by default). Sequences corresponding to the hits are retrieved from the temporary database, with no more than three sequences per genus and no more than ten sequences per phylum. Significant hits in query-hit alignments with a length of at least 120 amino acids are then re-ranked in descending order of query-hit identity. The homologous sequences plus the query are assembled and aligned using MUSCLE version 3.8.31 with default settings. The resulting alignment is trimmed using TrimAI version 1.2 in automatic mode (-automated1). We discarded queries with significantly different amino acid compositions compared to the remaining sequences in the alignment (P < 0.05). These pruned alignments (≥50 amino acids) were used to construct phylogenetic trees using the best model calculated using ModelFinder in IQtree, using IQtree Multicore version 1.6.12. Branch support was estimated using the ultra-fast bootstrap (UFboot, -bb 1500) test, the Shimodaira-Hasegawa-like likelihood ratio approximation test (SH-aLRT, -alrt1200), the approximate Bayesian test (-abayes), and the fast local bootstrap probability test (-bb 1500).

[0047] 4. Tree-based HGT inference

[0048] In phylogenetic trees, the present invention is dedicated to searching for heterologous nested topological structures of query sequences. To this end, the present invention integrates NestedIn, a phylogenetic tree scanning tool implemented in Java, into the Routinetree model to screen for heterologous nested positions in gene trees. The so-called heterologous nested positions refer to two or more monophyletic groups formed by the query sequence and its genetically distant sequences in the tree, and these monophyletic groups are supported by different nodes. This essentially reflects the conflict between the gene tree and the species tree. In order to enhance the applicability and flexibility of the tool, NestedIn specifically provides a user-friendly parameter "--donor" that allows users to input any taxonomic node they wish to test for horizontal gene transfer (HGT) relationships with the query sequence.

[0049] On the query (recipient) side, in order to ensure that all descendant genes in the monophyletic group after horizontal gene transfer (HGT) occur, the NestedIn tool provides a "--optional" parameter that allows users to enter the ancestral level of the query sequence. Specifically, the NestedIn workflow is as follows:

[0050] First, it reads the tree file in Newick format and parses the tree topology. Next, starting from the query sequence, it traverses upwards to its parent nodes. During this process, NestedIn records each node that matches the user-set nesting criteria until it reaches the first leaf node that is neither a donor nor an optional node. Multiple nodes may meet the filtering criteria. The final node is determined according to the following principles: 1) To minimize the interference of contamination and recent HGT events, NestedIn removes singletons from both the donor and recipient genes and allows the user to define the minimum number requirement; 2) Only nested positions that are multiply supported by a custom threshold (default setting: SH test ≥ 0.70 and aByes test ≥ 0.70) are retained in the supporting nodes; 3) Among all remaining nodes, the node with the highest taxonomic level is selected as the final candidate HGT node.

[0051] For a given nested position, NestedIn provides the most similar hits (MMSHs) of the query sequence in the recipient group (MMSH_IN) and the donor group (MMSH_OUT). This is done to obtain predicted functional information (MMSH_IN) and list a representative donor gene (MMSH_OUT) to facilitate analysis of evolutionary diversity between the recipient and donor after HGT.

[0052] 5. Verification of Tree-Inferred HGT Candidates

[0053] To construct a systematic method for rapid, comprehensive and reliable HGT detection, the present invention utilizes the outlier index (AI), HGT score support index (hU), HGT branch length support index (hBL) and its consensus hit support index to test the confidence of HGT candidates.

[0054] First, the present invention defines some concepts for querying sequences in monophyletic groups in gene trees:

[0055] INGROUP: All sequences within the user-specified hierarchy level that contain the query sequence.

[0056] OUTGROUP: All sequences outside the user-defined hierarchy level (that is, sequences that do not belong to INGROUP).

[0057] SKIPGROUP: The query sequence itself and the INGROUP boundaries are user-specified lower-level sequences. Sequences in SKIPGROUP are assumed to be homologous sequences and may have originated after an HGT event.

[0058] (1) Anomaly Index (AI)

[0059] The present invention uses the e-value of BLAST index to calculate the AI ​​score of each query gene:

[0060] AI = (bbhO / bbhG);

[0061] Where BbhG is the e-value of the best hit in the INGROUP lineage, and bbhO is the e-value of the best BLAST hit in the OUTGROUP lineage. The e-values ​​in SKIPGROUP are skipped because they are assumed to be homologous sequences after HGT. When no significant BLAST hit is detected, the corresponding bbhG or bbhO is set to 1. The AI ​​score is an indicator of the degree of similarity between the query sequence and its homologous sequence in the INGROUP compared to its homologous sequence in the OUTGROUP. Since all BLASTp searches are performed using the same database, it is reasonable to apply a uniform threshold for all query taxa.

[0062] During the initial screening phase, we selected a relatively loose AI score threshold (AI>0). In subsequent screening phases, this threshold can be adjusted based on specific analytical needs. As a practical option, in our study, we selected a moderately less stringent threshold (AI>10) compared to previous studies, as it produced the desired results.

[0063] (2) HGT score support index (hU)

[0064] The present invention calculates the hU score of each query gene based on the best alignment score of INGROUP and OUTGROUP:

[0065] hU = (OUTGROUP's best hit bit score) - (INGROUP's best hit bit score).

[0066] Bit scores in SKIPGROUP were skipped because they were assumed to be homologous sequences after HGT. When no significant BLAST hits were detected, the best hit bit scores for OUTGROUP and INGROUP were set to 0, respectively. The hU score reflects the degree of similarity between the query and its homologous sequences in INGROUP compared to the homologous sequences in OUTGROUP. Since all BLASTp searches are performed using the same REFAL database, it is reasonable to apply a uniform threshold to all query taxa. Considering its combination with other criteria in HGT inference, a relatively loose hU score threshold (hU>0) was selected in the initial screening stage.

[0067] (3) HGT branch length support index (hBL)

[0068] The present invention develops a new index hBL (HGT Branch Length Support Index) for each query gene, based on the minimum branch length of the query INGROUP and OUTGROUP:

[0069] hBL = (minimum branch length from INGROUP to query) - (minimum branch length from OUTGROUP to query).

[0070] Branch length values ​​in SKIPGROUP are skipped because they are assumed to be homologous sequences after HGT. When no genes are detected, the minimum branch length to the query is set to 100 for both INGROUP and OUTGROUP. The hBL score serves as an indicator of the degree of similarity between the query and its homologous sequences in the INGROUP compared to its homologous sequences in the OUTGROUP. For each leaf node in the tree, the branch length to the query is determined by summing the lengths of all branches connecting the leaf node to the query. Since all branch lengths are specific to the same tree, a direct comparison between the INGROUP and OUTGROUP can be made by summing these lengths. During the initial screening phase, the present invention applies a relatively relaxed hBL score threshold (hBL>0), considering its integration with other criteria for HGT inference. In subsequent screening phases, users retain the flexibility to adjust this threshold based on their specific analytical needs.

[0071] (4) Consensus hit support

[0072] Considering the possibility of accidentally introducing sequence contamination into an INGROUP or OUTGROUP, we calculated consensus hit support for AI, hU, and hBL. Consensus hit support measures the degree of support provided by all genes in the OUTGROUP (not just the best hit gene). The following are the specific metrics:

[0073] Consensus hit support - E-value (CHE): This metric represents the ratio of the number of genes in the OUTGROUP with an E-value smaller than bbhG (the best hit in the INGROUP) to the total number of genes in the OUTGROUP. CHE serves as a confidence indicator for the AI>0 condition.

[0074] Consensus hit support score (CHS): CHS represents the ratio of the number of genes in the OUTGROUP with a score greater than bbhG to the total number of genes in the OUTGROUP. This metric serves as a confidence indicator for the hU>0 condition.

[0075] Consensus hit support - branch length (CHBL): CHBL is calculated as the ratio of the number of genes in the OUTGROUP with branch length less than bbhG to the total number of genes in the OUTGROUP. It serves as a confidence indicator for the condition hBL>0.

[0076] These consensus hit support metrics provide additional insights into the reliability of AI, hU, and hBL conditions by considering a broader range of genes in OUTGROUP.

[0077] In addition to the above strategies, the present invention also implemented the following criteria to evaluate the taxonomy of potential candidate HGT physically flanking genes to eliminate the possibility of contamination:

[0078] (1) If an HGT candidate gene is located in a contig (a set of continuous DNA fragments, i.e., a contig) in which 50% of the best hits are in other kingdoms, the candidate gene will be excluded;

[0079] (2) If the HGT candidate gene was located in a contig in which 50% of the genes were primarily identified as HGT genes, it was excluded;

[0080] (3) If the best hit of at least one of the three closest flanking genes (upstream and downstream) of the HGT candidate gene is in other kingdoms, it is excluded;

[0081] (4) When two or more HGT genes are physically closely linked and belong to the same gene family, they will be considered as a single HGT event, and all of these HGTs will be retained for further analysis.

[0082] 6. Assign HGT to timeline

[0083] To better understand the impact of HGT on Earth's evolutionary history and geological change, and to place these events in a temporal context, this paper aims to estimate the timing of these HGT events. However, tracking the exact timing of gene divergence over time based on the accumulation of nucleotide mutations has become increasingly challenging. Over time, the probability of multiple nucleotide substitutions at a single site increases, as do unexpected events such as gene loss and gene duplication. This paper aims to provide a broad overview of HGT during the life of an organism, rather than delving into the complexities of evolutionary algorithms.

[0084] Nevertheless, the present invention can establish general upper and lower bounds on the timing of HGT events based on the following consensus principle: vertical inheritance is unidirectional. For example, if a successful prokaryotic-to-eukaryotic HGT event occurred at a specific evolutionary node (e.g., a single-celled organism that is the direct ancestor of modern brown algae), then the homologous genes generated by HGT could be transmitted via vertical inheritance to descendants of the brown algae node, but not to the ancestral node, the anisoflagellate.

[0085] For a given query protein sequence identified as potentially HGT, if the present invention can locate all descendants derived from the initial HGT event, the present invention can trace the sequence back to a common ancestor. The timing of this common ancestor can be inferred using molecular clock methods combined with archaeological and fossil evidence. Technically, the minimum taxonomic boundary of HGT descendants can be traced using all gene members in the INGROUP, while the minimum taxonomic boundary of donor descendants can be traced using all gene members in the OUTGROUP.

[0086] Accurate taxonomy is crucial for representing the hierarchical structure and relationships of evolutionary nodes. To this end, we used the NCBI classification system, which is widely accepted and a reasonable classification framework for classifying all life forms on Earth. We manually introduced a kingdom-level node, Chromalveolate, into the NCBI classification system and proposed reclassifying Cryptophytes, Rhizopods, Apicomplexa, Anisochorates, and Dinoflagellates into this kingdom based on recent research. Timeline data for the interval nodes was primarily sourced from the TimeTree database.

[0087] 7. Output and Visualization

[0088] For each individual query protein, the present invention provides its homologous genes, corresponding alignments, and phylogenetic trees. For the ultimately identified HGT genes, detailed information is provided in a tab-delimited text file, including the donor and recipient nodes and their occurrence time, the MMSH genes in the INGROUP and OUTGROUP, AI, hU, and hBL scores and their support index, as well as predicted functional information from multiple databases. The accession numbers listed in the table are hyperlinked to corresponding external databases such as GO and KEGG. Table 1 shows the output information for some query proteins.

[0089] Table 1 Output information of protein query

[0090] .

[0091] The feasibility of the method of the present invention was verified by the following method:

[0092] 1. Identify core genes and conduct in-depth comparison with HGT genes

[0093] Core genes are typically present in the genomes of nearly all members of a species or taxonomic group and constitute the most conserved set of genes. Identifying these core genes and comparing them with HGT is a fundamental pursuit of evolutionary genomics. This effort allows researchers to assess adaptation and evolution through HGT using host core genes as a reference point, revealing the genetic mechanisms underlying adaptation and diversity within a domain of organisms. To achieve these goals, we compared GC ratios, codon usage, and gene structure to identify differences between HGT and CORE genes.

[0094] (1) Determination of core genes

[0095] The basic criterion is that if a gene family can encompass at least 70% of the species within a phylum, it is considered the core gene family of the phylum and its members are designated as core genes.

[0096] (2) Codon usage

[0097] Codon usage is an important factor in determining HGT fate because it is compatible with the host's transcriptional machinery and tRNA repertoire. We calculated codon usage and GC content using CodonW version 1.4.4 (http: / / codonw.sourceforge.net). Correlations between CAI and gene expression were tested using the Spearman rank correlation analysis tool (P. Wessa, Free Statistical Software, Office of Research and Education Development, version 1.1.23-r7, https: / / www.wessa.net / ).

[0098] The present invention first used CodonW to analyze various gene parameters (synonymous codon number (L_sym), total amino acid number (L_aa), G+C content (GC), G+C content GC3 position (GC3s), A / T / C / G content synonymous codon third (A3s, T3, G3, C3s)), codon usage indices (codon adaptation index (CAI), optimal codon frequency (Fop), number of effective codons (Nc), codon bias index (CBI)), and amino acid indices (hydropathy (GRAVY score), protein color (Aromo)) for each species. The codon usage bias of each gene in each species was analyzed. The present invention then identified the most conserved genes in each species as the core genes of that species. For each species, the present invention classified HGT genes according to the transfer node (red plant, green alga, green alga, catenella plant, embryophyte, tracheophyte, mesoangiosperm, rose plant, and plant) and the donor class (prokaryotic and eukaryotic).

[0099] To determine the codon usage bias testing method, we first used SPSS to test normality and homogeneity of variance for the core genes (CORE) and HGT codon bias indices of each species. KS tests were performed on 14 codon usage bias indices (T3s, C3s, A3s, G3s, CAI, CBI, Fop, Nc, GC3s, GC, L_sym, L_aa, Gravy, and Aromo) across 113 species. Asymptotic significance > 0.05 indicated that the indices conformed to a normal distribution; otherwise, they did not. For nine species, including Cyanopa, the codon usage bias indices for each species were all less than 0.05, indicating that none conformed to a normal distribution. Because parametric tests require that each data set conform to a normal distribution, nonparametric tests were used to test differences between CORE and HGT genes, as well as between adjacent node transfer genes between different phyla. A KS test was performed on 14 codon usage bias indicators (T3s, C3s, A3s, G3s, CAI, CBI, Fop, Nc, GC3s, GC, L_sym, L_aa, Gravy, and Aromo) for 113 species. If the asymptotic significance is > 0.05, it indicates that the indicator conforms to a normal distribution; otherwise, it does not. For nine species, including Cyanopa, the codon usage bias indicators for each species were all less than 0.05, indicating that none conformed to a normal distribution. Because parametric tests require that each data set must conform to a normal distribution, nonparametric tests are required to test differences between CORE and HGT genes, as well as genes transferred between adjacent nodes in different phyla. The results are shown in Tables 2 and 3.

[0100] Table 2 Statistical description a

[0101] ;

[0102] Table 3 One-sample Kolmogorov-Smirnov test

[0103] .

[0104] Since the data of each species cannot satisfy normality and homogeneity of variance at the same time, the present invention uses the Mann-Whitney U test in non-parametric tests to detect the differences in CORE and HGT genes and adjacent node transfer genes between different phyla. The results are as follows Figure 5 As shown in the figure, red dots represent core genes, blue dots represent HGT from prokaryotes (Pr), and green dots represent HGT from eukaryotes (Eu). Each dot on the same horizontal axis represents a species. The horizontal axis, from left to right, represents core genes and the incoming nodes, arranged from ancient to recent. Connecting lines indicate the differences between adjacent horizontal axes within the species. A red line indicates a significant increase in the codon bias index for recent HGT compared to ancient nodes, a blue line indicates a significant decrease, and a gray line indicates no significant difference.

[0105] In order to better compare HGT genes with core genes, the present invention uses core genes to normalize HGT genes, uses violin plots to compare the codon usage preferences of all nodes, and establishes a correlation between codon usage preference and transfer time. Since the evolutionary rate of each species is different, if a horizontal comparison is to be made, the present invention uses the core genes of each species as a reference standard. The codon preference index of HGT of each species is normalized relative to the codon preference index of its respective core genes, that is, the codon preference index of HGT is divided by the codon preference index of the core genes. Violin plots are used to present the codon usage preferences of different species on all nodes, and a correlation between codon usage preference and transfer time is established. Figure 3 and Figure 4 The temporal trend of codon preference in horizontal gene transfer (HGT) is shown. Figure 4As shown, each point on the same horizontal axis represents the normalized codon usage bias index for a species. HGT genes acquired at different evolutionary nodes were compared according to their timing. The parameters compared included L_aa, L_sym, GC3s, CAI, Fop, and Nc. Detailed explanations of each parameter are as follows: Before 1000 MYA, genes transferred from eukaryotes and prokaryotes to plants had similar L_aa and L_sym values, significantly shorter than those of core genes, accounting for approximately 75% of the core. It has been widely reported that long fragments of foreign DNA undergo more rapid inactivation and are less likely to be fixed in the recipient genome than short insertions. However, a notable difference is that L_aa and L_sym values ​​for genes transferred from eukaryotes to plants have steadily increased over time, and in recent times, have become very close to those of core genes. This difference may stem from the proximity of plants to the donor organism. Clearly, plants are more closely related to other eukaryotes than prokaryotes, sharing a more similar genomic organization, which is more conducive to the fixation of long foreign genes.

[0106] GC3s show distinct trends in eukaryotic and prokaryotic genes. Prokaryotic genes consistently have higher GC3s than host core genes, increasing over time. Conversely, eukaryotic genes consistently have lower GC3s than host core genes, decreasing over time. This clearly illustrates the process of "domestication," whereby earlier acquired genes evolve alongside the host genome, with codon usage increasingly resembling that of the host's core genes.

[0107] CAI refers to the degree to which synonymous codon usage in the coding region matches the optimal codon usage frequency. The CAI of a specific gene can be determined by comparing its codon usage frequency with a reference set of highly expressed genes in the species. A gene's CAI score is calculated based on the usage frequency of all codons in that gene. CAI is often used to assess the expression level of exogenous genes in the host; the higher the CAI, the higher the expression level of the exogenous gene in the host. The CAI of prokaryotic genes is lower than that of the host core genes, but over time it has become closer to the CAI value of the host core genes. The CAI of eukaryotic genes remains very close to that of the host core genes.

[0108] Fop refers to the ratio of optimal codons to their synonymous codons. The higher the Fop, the more frequently the optimal codon is used. Over time, the Fop of prokaryotic genes has become increasingly higher, while the Fop of eukaryotic genes has remained very close to that of the host's core genes.

[0109] Nc is a measure of the degree of codon bias in a gene and quantifies the extent to which a gene uses the same or equivalent synonymous codons for each amino acid class. NC is the best overall estimator of absolute synonymous codon usage bias and can be easily calculated from codon usage data alone. For each gene, NC values ​​range from 20 (extreme bias when only one codon is used per amino acid) to 61 (when all codons are used uniformly). Higher NC values ​​indicate weaker codon usage bias. Nc values ​​for genes transferred from both eukaryotes and prokaryotes to plants have shown a decreasing trend over time, becoming increasingly closer to core genes.

[0110] (4) Gene function annotation and enrichment

[0111] This study collected annotation information from nine major databases: OG1157, eggNOG, PANTHER, SuperFamily, Interproscan, GO, KEGG, Pfam, and KO. The function of the query sequence was predicted based on the sequence with the lowest e-value in the BLAST alignment. Furthermore, the eggNOG Mapper tool (European Molecular Biology Laboratory; http: / / eggnog-mapper.embl.de) was used to perform annotation analysis on the complete gene set of the species.

Claims

1. A method for screening cross-species HGT, characterized in that: The following steps are involved: Step 1: constructing a background database, which contains known protein sequences covering the phylogenetic lineages of prokaryotes and eukaryotes, and removing highly homologous and redundant sequences; Step 2: Perform BLAST comparison on any two protein sequences in the background database, calculate the homology group association network based on the best mutual alignment in each pair of gene combinations, and cluster the homology group association network using an unsupervised hidden Markov clustering algorithm to generate a database containing multiple homology groups; Step 3: For each query protein sequence to be analyzed, search for homologous sequences in the background database, perform strict sequence alignment and pruning, construct a single-gene phylogenetic tree, and evaluate the support of each branch of the evolutionary tree by at least one of the ultra-fast bootstrap test, similar likelihood ratio approximation test, approximate Bayesian test, and fast local bootstrap probability test; Step 4: Search for nested topological structures between the query sequence and genetically distant sequences in the single-gene phylogenetic tree constructed with the query sequence as the center. Screen for potential cross-species HGT by defining nested positions and setting a minimum support node threshold. Simultaneously, record the most similar monophyletic matching sequences of the query sequence in both the recipient and donor groups. Step 5: Evaluate the confidence of potential cross-species HGT by calculating the outlier index, HGT score support index, HGT branch length support index and its consensus hit support index; Step 6: Estimate the time of occurrence of cross-species horizontal transfer events by tracing the minimum taxonomic boundaries of cross-species HGT offspring and the minimum taxonomic boundaries of donor offspring; combine molecular clock methods with archaeological and fossil evidence to provide a temporal context for cross-species HGT.

2. The method for screening cross-species HGT according to claim 1, characterized in that In the step 1, the data in the background database includes all available protein sequences of the reference sequences in NCBI and protein sequences from the marine microbial eukaryotic transcriptome sequencing project.

3. The method for screening cross-species HGT according to claim 1, characterized in that In the second step, multiple sequence alignment and conserved domain analysis are performed on the sequences in each homology group.

4. The method for screening cross-species HGT according to claim 1, wherein The step 5 also includes analyzing potential cross-species HGT using homology comparison features of genes above and below the HGT gene locus, specifically including: If the potential cross-species HGT is located in a chromosome segment and 50% of the best matches of the genes in the chromosome segment are from other biological kingdoms, the HGT candidate gene is eliminated; or If the potential cross-species HGT is located in a chromosome segment and 50% of the genes in the chromosome segment are identified as HGT genes, the HGT candidate gene is eliminated; or, If at least one of the three closest upstream and downstream genes for potential cross-species HGT had a best hit in another kingdom, it was excluded; or When two or more potential cross-species HGTs are physically closely linked and belong to the same gene family, they will be regarded as a single cross-species HGT event, and all these adjacent potential cross-species HGTs will be retained.

Images