Method for screening cross-species HGT

Abstract

The invention belongs to the technical field of bioinformatics, and particularly relates to a method for screening cross-species HGT. The method comprises the steps of constructing a background database, inferring a homologous group, constructing a gene phylogenetic tree, inferring HGT based on the gene tree, inferring verification of candidate HGTs based on the gene tree, distributing the HGTs to a timeline and the like. According to the method, high-throughput and automatic screening of HGT events in a large number of gene datasets is realized by integrating a plurality of steps of phylogenetic analysis, sequence alignment, gene function annotation, statistical verification and the like. The method disclosed by the invention is particularly suitable for the HGT research of an ultra-long time length and an ultra-long evolution distance crossing a biological boundary, and an important scientific tool can be provided for understanding gene flow and function enhancement in a biological evolution process.

Description

Technical Field

[0001] The present invention belongs to the technical field of bioinformatics, and in particular 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 phenomenon of exchanging genetic material between two species with non-vertical genetic relationships. This process is common in prokaryotes and is achieved through mechanisms such as conjugation, transduction and transformation. Unlike vertical gene transfer, which is the inheritance from parent to offspring, HGT transcends the boundaries of kinship and greatly enriches the dynamics of gene flow.

[0003] The history of HGT observation can be traced back to 1959, when it was documented that high-frequency transduction (Hfr) Escherichia coli could laterally transfer genetic information to a specific mutant of Salmonella typhimurium. In the same year, Tomochiro Akiba and Kunitaro Ochiai discovered resistance plasmids in pathogenic bacteria and subsequently confirmed that these plasmids could be transferred between different bacterial strains. However, the concept of HGT had not yet been formed at the time. It was not until the 1990s, with the emergence of genetically modified organisms (GEOs), especially genetically modified microorganisms (GEMs), and the emergence of many drug-resistant pathogens, whose origins could no longer be simply attributed to gene 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 genes shared in their genomes. The absence of these genes in closely related lineages may be caused by multiple independent gene loss events or by 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 protection and stability of the eukaryotic tree of life, HGT is an important force that promotes the diversification of eukaryotic species and helps them adapt to diverse environments. Although HGT has been widely recognized as a major evolutionary force in prokaryotes, its role in eukaryotic evolution remains controversial, mainly due to the complex evolutionary history, complex genome structure and frequent sequence contamination problems.

[0005] In order to accurately identify HGT events, it is urgent to develop efficient and accurate screening methods. Currently, common HGT analysis methods include evolutionary tree analysis, base composition analysis, selection pressure analysis, intron analysis, specific sequence analysis, and nucleotide composition bias analysis. However, existing HGT screening methods mostly rely on phylogenetic analysis, which has limitations such as incomplete database coverage, immature evolutionary tree automatic selection algorithms, and lack of unified analysis standards, and urgently need improvement and innovation. Summary of the invention

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

[0007] To achieve the above object, the technical solution adopted by the present invention is: a method for screening cross-species HGT, comprising the following steps: Step 1: constructing a background database, wherein the background database 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 use an unsupervised hidden Markov (HMM) clustering algorithm to cluster the homology group association network 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 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; Step 4: Retrieve the topological structure of the query sequence and the distant sequences in genetic distance in the single gene phylogenetic tree constructed with the query sequence as the center; screen out potential cross-species HGT by defining the nested position and setting the minimum support node threshold; at the same time, record the most similar matching sequence of the query sequence in the recipient group and the donor group; 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 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.

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

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

[0010] Preferably, the step 5 further comprises analyzing potential cross-species HGT using homology comparison features of genes on both sides of the HGT gene locus, specifically comprising: 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 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 were physically closely linked and belonged to the same gene family, they would be considered as a single cross-species HGT event, and all these adjacent potential cross-species HGTs would be retained.

[0011] Compared with the prior art, the method of the present invention has the following beneficial effects: (1) Efficiency: Through automated and standardized processes, high-throughput screening of HGT events in large genetic data sets is achieved, greatly improving analysis efficiency; (2) Accuracy: The accuracy of HGT screening is ensured through a multi-step process such as building a large background database, inferring homologous groups, and constructing a gene phylogenetic tree. At the same time, the introduction of multiple verification indicators and taxonomic analysis of physical flanking genes further reduces the false positive rate. (3) Wide applicability: This method is particularly suitable for HGT research 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

[0012] 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; Figure 2 This is a schematic diagram of the method for screening cross-species HGT in an embodiment of the present invention; Figure 3 The topological structure of the evolutionary nodes of 113 organisms and the number of HGTs obtained at each node; Figure 4The codon preference of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources in 113 organisms is shown along the time axis. The codon preference index of each species is normalized relative to its respective core gene, and the x-axis represents the time order of HGT integration into the plant lineage. Among them, (a)-(b) are the amino acid lengths of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources, respectively; (c)-(d) are the amino acid lengths of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources, respectively. Amino acid aromaticity; (e)-(f) are the GC contents of the third position of synonymous codons in HGT from prokaryotes (Pr) and eukaryotes (Eu), respectively; (g)-(h) are the codon adaptation indexes in HGT from prokaryotes (Pr) and eukaryotes (Eu), respectively; (i)-(j) are the optimal codon frequencies in HGT from prokaryotes (Pr) and eukaryotes (Eu), respectively; (m)-(n) are the effective codon numbers in HGT from prokaryotes (Pr) and eukaryotes (Eu), respectively; Figure 5 Comparison of codon bias of core genes (CORE) and horizontal gene transfer (HGT); (a)-(b) are the number of effective codons (Nc) of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources of Rhodophyta organisms; (c)-(d) are the codon bias index (CBI) of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources of Rhodophyta organisms; (e)-(f) are the codon bias index (CBI) of HGT from prokaryotic (Pr) and eukaryotic (Eu) sources of Chlorophyta organisms. (g)-(h) are the HGT codon bias index (CBI) of prokaryotic (Pr) and eukaryotic (Eu) sources of Chlorophyta organisms, respectively; (i)-(j) are the HGT codon number (Nc) of prokaryotic (Pr) and eukaryotic (Eu) sources of Streptophyta organisms, respectively; (k)-(m) are the HGT codon bias index (CBI) of prokaryotic (Pr) and eukaryotic (Eu) sources of Streptophyta organisms, respectively. DETAILED DESCRIPTION

[0013] The present invention is further described below in conjunction with 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 is obviously capable of being implemented in a variety of other ways different from the description herein. For those skilled in the art, any substitution, improvement or change made to the embodiments of the present invention is within the protection scope of the present invention, and the protection scope of the present invention should not be limited by the content of this specific embodiment.

[0014] The present invention provides a method for screening cross-species HGT to facilitate large-scale, high-throughput HGT detection and explore the breadth of biological evolution driven by HGT in the Tree of Life (TOL). The present invention 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 contributions to the complexity and diversification of biological systems are inferred, providing research with the possibility of verifying many evolutionary theories. The system (HGTStart) used by this method is as follows Figure 1 As shown, the process is as Figure 2 As shown, the following steps are included:

[0015] 1. Build background database

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

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

[0018] An important issue affecting HGT inference is uneven sampling due to unbalanced data collection between different taxa. The present invention incorporates proteins from MMETSP into the database to make up for the problem of less red algae genome data. These extensive genomic data form a protein database containing 17,250,679 protein sequences from 1157 genomes, which have reasonable coverage in 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 algae lineage). Each of these 1157 complete genomes represents a representative species of its genus. The database is named "GNM1157".

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

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

[0021] To facilitate the function prediction of the query sequence, 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.

[0022] 2. Homology group inference and benchmarking

[0023] Homology between genes, including orthologs and paralogs, is the cornerstone of comparative biology studies (e.g., horizontal gene transfer analysis). Currently, tree-based databases such as PhylomeDB, Ensembl-Compara, EggNOG, and TreeFam have been developed in this field. However, the present invention found that these existing databases do not fully cover the entire gene set encoded in GNM1157, so it is particularly important to build a customized homology database. The present invention clustered the 17,250,679 proteins encoded in GNM1157 and combined them with protein-based topological structure analysis. All protein sequences were analyzed in detail using the all-to-all BLAST workflow for whole genome pairs (version 2.2.28, with an e-value threshold of 1e-10 and a local consistency threshold of 20%). Subsequently, the association network of homology groups was calculated based on the mutual best length normalized hit (RBNH) information in each pair of genome combinations using OrthoFinder V2.3.7 software. OrthoFinder performs well in correcting the dependence of gene similarity on gene length and phylogenetic distance, thereby improving the accuracy of homology group definition. In addition, the present invention independently runs an unsupervised Markov clustering algorithm (MCL3.0) to cluster this homology group map. By adopting a set of gradient expansion parameters (1.2, 1.4, 1.6, 1.8, 2.0, 2.2), the present invention runs the MCL algorithm multiple times and selects the best running result. The present invention locally downloads more than a dozen of the most popular online homology databases (e.g., EGGNOG, PANTHERN, and SUPERFAMILY) as benchmarks to calibrate the homology inference parameters of the present invention. The precision, recall, and F0 values ​​of the gradient set are calculated using the NumPy library in Python 3.0 to determine the optimal value of the expansion coefficient.

[0024] 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. wrong 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, which balances these two metrics and is more sensitive to lower metrics, which is conducive to comprehensively evaluating the performance of homology group assignments.

[0025] To identify biologically significant homology clusters, the sequences in each homology group (OG) were tested for group-specific conserved domain performance and noise and weakly associated sequences were removed. First, MAFFT v7.455 was used to perform multiple alignments of sequences within each OG. Then, HMM profiles of conserved domains were constructed for the alignment results of each OG using HMMER 3.0. Subsequently, the Hmmsearch program was used to search all protein sequences in each OG against their corresponding HMM profiles with an E-value threshold of 0.00001. Sequences below the threshold were removed from the group and collected for a second round of sequence homology construction. After this round of screening, each sequence was accurately assigned to a homology group. Ultimately, this clustering process generated 27,631 homology groups, each containing at least two members. For homology groups with five or more members, the present invention 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.

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

[0027] 3. Constructing a phylogenetic tree

[0028] The present invention designs 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 builds a single gene phylogenetic tree for each custom query protein sequence. In order 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 hmmscan, 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 for subsequent retrieval. The excluded REFAL sequences are collected into a FASTA file.

[0029] At the beginning, the query sequence will be split into individual sequence files in the first step and scanned for PFAM HMM profiles using hmmscan separately. A temporary BLAST / diamond database will be built based on the hit PFAM number of each split query sequence, from which the corresponding FASTA files in REFAL will be retrieved and merged. After making the search index in BLAST / diamond format, the temporary database will be searched using the default e-value threshold = 1e-05. For each query, the top 10,000 significant hits in descending order of position score (by default) are recorded. 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 of the query-hit alignment with a length of at least 120 amino acids are subsequently re-ranked in descending order based on query-hit agreement. 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 by ModelFinder in IQtree under IQtree Multicore version 1.6.12. Branch support was estimated by ultra-fast bootstrap (UFboot, -bb 1500) test, Shimodaira-Hasegawa-like likelihood ratio approximation test (SH-aLRT, -alrt1200), approximate Bayesian test (-abayes), and fast local bootstrap probability test (-bb 1500).

[0030] 4. Tree-based HGT inference

[0031] In the phylogenetic tree, 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 heterologous nested positions in the gene tree. The so-called heterologous nested position refers 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 enter any classification node they wish to test for horizontal gene transfer (HGT) relationship with the query sequence.

[0032] 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:

[0033] First, it reads the tree file in Newick format and parses the topology of the tree. Then, starting from the query sequence, it traverses its parent nodes step by step. In this process, NestedIn records each node that matches the nesting criteria set by the user until it reaches the first leaf node that is neither a donor nor an optional node. There may be multiple nodes that meet the filtering criteria. The determination of the final node follows the following principles: 1) In order to exclude the interference of contamination and recent HGT events as much as possible, NestedIn will remove singletons in donor and recipient genes and allow users to customize the minimum number requirements; 2) Only nested positions that are multiply supported by custom thresholds (default settings: SH test ≥ 0.70 and aByes test ≥ 0.70) in supporting nodes are retained; 3) Among all the remaining nodes, the node with the highest classification level is selected as the final candidate HGT node.

[0034] 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 analyze the evolutionary diversity between the recipient and the donor after HGT.

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

[0036] In order to construct a systematic method to detect HGT quickly, comprehensively and reliably, 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.

[0037] First, the present invention defines some concepts for querying sequences in monophyletic groups in gene trees: INGROUP: All sequences within the user-specified hierarchy level that contain the query sequence.

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

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

[0040] (1) Anomaly Index (AI) The present invention uses the e-value of the BLAST index to calculate the AI ​​score for each query gene: AI = (bbhO / bbhG); Where BbhG is the e-value of the best hit in the INGROUP pedigree, and bbhO is the e-value of the best BLAST hit in the OUTGROUP pedigree. The e-values ​​in SKIPGROUP are skipped because they are assumed to be homologous sequences after HGT. When no significant BLAST hits were detected, the corresponding bbhG or bbhO was set to 1. The AI ​​score is an indicator of the degree of similarity of the query sequence to its homologous sequences in INGROUP compared to its homologous sequences in OUTGROUP. Since all BLASTp searches are performed using the same database, it is reasonable to apply a uniform threshold to all query taxa.

[0041] In the initial screening stage, the present invention selected a relatively loose AI score threshold (AI>0). In the subsequent screening stage, this threshold can be adjusted according to its specific analysis needs. As a practical choice, in the present invention's study, the present invention selected a moderately less stringent threshold (AI>10) compared to previous studies because it produced the desired results.

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

[0043] The present invention calculates the hU score of each query gene based on the best comparison score of INGROUP and OUTGROUP: hU = (best hit bit score of OUTGROUP) - (best hit bit score of INGROUP).

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

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

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

[0047] The branch length values ​​in SKIPGROUP are skipped because they are assumed to be homologous sequences after HGT. When no gene is detected, the minimum branch length to the query is set to 100 for INGROUP and OUTGROUP, respectively. The hBL score serves as an indicator of how similar the query is to its homologous sequences in INGROUP compared to its homologous sequences in OUTGROUP. For each leaf node in the tree, the branch length to the query is determined by summing up all branch lengths connected to the leaf node to the query. Since all branch lengths are specific to the same tree, a direct comparison can be made between INGROUP and OUTGROUP by summing these lengths. In the preliminary screening stage, the present invention applies a relatively loose hBL score threshold (hBL>0), considering that it is combined with other criteria in HGT inference. In the subsequent screening stage, users retain the flexibility to adjust this threshold according to their specific analysis needs.

[0048] (4) Consensus hit support

[0049] Considering the possibility of accidentally introducing sequence contamination into INGROUP or OUTGROUP, the present invention calculates the consensus hit support of AI, hU and hBL. Consensus hit support measures the support provided by all genes in OUTGROUP (not just the best hit gene). The following are the specific indicators: Consensus hit support - E-value (CHE): This metric represents the ratio of the number of genes in OUTGROUP with an E-value less than bbhG (the best hit in INGROUP) to the total number of genes in OUTGROUP. CHE is used as a confidence metric for the AI>0 condition.

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

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

[0052] 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 the OUTGROUP.

[0053] In addition to the above strategies, the present invention also implemented the following criteria to evaluate the taxonomy of potential candidate HGT physical flanking genes to eliminate the possibility of contamination: (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 of the genes are in other kingdoms, the candidate gene will be excluded; (2) If the HGT candidate gene was located in a contig in which 50% of the genes were mainly identified as HGT genes, it was excluded; (3) If at least one of the three closest flanking genes (upstream and downstream) of the HGT candidate gene had a best hit in other kingdoms, it was excluded; (4) When two or more HGT genes are physically closely linked and belong to the same gene family, they will be regarded as a single HGT event, and all of these HGTs will be retained for further analysis.

[0054] 6. Assign HGT to timeline

[0055] To better understand the impact of HGT events on Earth's evolutionary history and geological changes, 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 based on the accumulation of nucleotide mutations has become increasingly challenging over time. Over time, the probability of multiple nucleotide substitutions at a single site increases, as well as the occurrence of unexpected events such as gene loss and gene duplication. The goal of this paper is to provide a broad overview of HGTs during the existence of an organism, rather than delving into the complexity of evolutionary algorithms.

[0056] Nevertheless, the present invention can establish general upper and lower limits on the time of HGT events based on the following consensus principles: vertical inheritance is unidirectional. For example, if a successful prokaryotic to eukaryotic HGT event occurred at a specific evolutionary node (such as the direct ancestor of modern brown algae, a unicellular organism), then the homologous genes generated by HGT can be transmitted to the descendants of the brown algae node through vertical inheritance, but cannot be transmitted to the ancestral node, the anisochoroidallum.

[0057] For a given query protein sequence determined as a potential HGT, if the present invention can find all descendants derived from the initial HGT event, the present invention can be traced back to a common ancestor node. The time of occurrence of this common ancestor can be inferred by molecular clock methods combined with archaeological and fossil evidence. Technically, the minimum classification boundary of HGT descendants can be tracked using all gene members in INGROUP, while the minimum classification boundary of donor descendants can be tracked using all gene members in OUTGROUP.

[0058] Accurate taxonomy is essential for representing the hierarchical structure and organism relationships of evolutionary nodes. To this end, the present invention uses the NCBI classification system, which is used to classify all life forms on Earth due to its widely accepted and reasonable classification framework. The present invention manually introduced a kingdom-level node in the NCBI classification system, Chromalveolate, and proposed to reclassify cryptophytes, rhizopods, apicomplexans, anisoflagellates, and dinoflagellates to this kingdom based on recent research. The timeline data of the interval nodes mainly comes from the TimeTree database.

[0059] 7. Output and Visualization

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

[0061] Table 1 Output information of protein query .

[0062] The feasibility of the method of the present invention is verified by the following method:

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

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

[0065] (1) Determination of core genes

[0066] 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.

[0067] (2) Codon usage

[0068] Codon usage is an important factor in determining the fate of HGT because it is compatible with the gene transcription machinery and tRNA pool in the host. Codon usage and GC content index were calculated using CodonW version 1.4.4 (http: / / codonw.sourceforge.net). The correlation test between CAI and gene expression was performed 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 / ).

[0069] The present invention firstly analyzed various gene parameters (number of synonymous codons (L_sym), total number of amino acids (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 index (codon adaptation index (CAI), optimal codon frequency (Fop), number of effective codons (Nc), codon bias index (CBI)) and amino acid index (hydrophilicity (gravy score), protein chromaticity (Aromo)) of each species using CodonW. The codon usage preference of each gene in each species was analyzed. Then, the present invention determined the most conserved gene in each species as the core gene of the species. For each species, the present invention classified the HGT genes according to the transfer node (red plant, green algae, green algae, chain plant, embryo plant, tracheophyte, mesoangiosperm, rose plant and plant) and donor category (prokaryotes and eukaryotes).

[0070] In order to determine the test method of codon usage preference, the normality and variance homogeneity tests were first performed on the core genes (CORE) and HGT codon preference indexes of each species using SPSS. The KS test was performed on 14 codon usage preference indices (T3s, C3s, A3s, G3s, CAI, CBI, Fop, Nc, GC3s, GC, L_sym, L_aa, Gravy and Aromo) of 113 species. If the asymptotic significance is > 0.05, it means that the index conforms to the normal distribution, otherwise it does not conform. Taking 9 species such as Cyanopa as an example, the codon usage preference index of each species is less than 0.05, which means that they do not conform to the normal distribution. Since the parametric test requires that each group of data must conform to the normal distribution, a non-parametric test method is required to test the differences between CORE and HGT genes and adjacent node transfer genes between different phyla. KS test was performed on 14 codon usage preference indicators (T3s, C3s, A3s, G3s, CAI, CBI, Fop, Nc, GC3s, GC, L_sym, L_aa, Gravy and Aromo) of 113 species. If the asymptotic significance is > 0.05, it means that the indicator conforms to the normal distribution, otherwise it does not conform. Taking 9 species such as Cyanopa as an example, the codon usage preference indicators of each species are less than 0.05, which means that they do not conform to the normal distribution. Since the parametric test requires that each group of data must conform to the normal distribution, non-parametric test methods are required to test the differences between CORE and HGT genes and adjacent node transfer genes between different phyla. The results are shown in Tables 2 and 3: Table 2 Statistical description a ; Table 3 One-sample Kolmogorov-Smirnov test .

[0071] Since the data of each species cannot satisfy normality and homogeneity of variance at the same time, the present invention uses the Mann-WhitneyU 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 shown in Figure 5 As shown in the figure, the red dots represent core genes, the blue dots represent HGT from prokaryotes (Pr), and the green dots represent HGT from eukaryotes (Eu). Each dot on the same horizontal axis represents a species. The horizontal axes represent the core genes and the transfer nodes arranged from ancient to recent in time from left to right. The connecting lines represent the differences between adjacent horizontal axes of the species. The red line indicates that the codon preference index of the recent node HGT of the species is significantly increased compared with the ancient node HGT, the blue line indicates a significant decrease, and the gray line indicates no significant difference.

[0072] 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 the 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 gene, that is, the codon preference index of HGT is divided by the codon preference index of the core gene. A violin plot is used to present the codon usage preferences of different species at all nodes, and the 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 in the figure, each point on the same horizontal axis represents the normalized codon usage bias index of a species. The HGT genes obtained at different evolutionary nodes are compared according to their time sequence. The parameters compared include L_aa, L_sym, GC3s, CAI, Fop and Nc, and the detailed explanations of each parameter are as follows: Before 1000MYA, the L_aa and L_sym of genes transferred from eukaryotes and prokaryotes to plants were similar, significantly shorter than core genes, and about 75% of core genes. It is widely reported that long fragments of foreign DNA undergo faster inactivation and are less likely to be fixed in the genome of the recipient organism than short insertions. However, a significant difference is that the L_aa and L_sym of genes transferred from eukaryotes to plants have steadily increased over time, and in modern times, they have been very close to core genes. This difference may be due to the closeness of plants to the donor. Obviously, plants are closer to other eukaryotes than prokaryotes, have a more similar genome composition, and are more conducive to the fixation of long fragments of foreign genes.

[0073] GC3s show completely different trends in genes of eukaryotic and prokaryotic origin. The GC3s of prokaryotic genes are always higher than those of the host core genes, and they are getting higher and higher over time. The GC3s of eukaryotic genes are always lower than those of the host core genes, and they are getting lower and lower over time. This well demonstrates the process of "domestication", that is, the earlier acquired genes evolve with the host genome, and the codon usage is getting closer and closer to the core genes of the host itself.

[0074] CAI refers to the degree of consistency between the usage frequency of synonymous codons in the coding region and the optimal codon usage. The CAI of a specific gene can be determined by comparing the codon usage frequency of a specific gene with a reference set of highly expressed genes of the species. The CAI score of a gene is calculated based on the usage frequency of all codons in the gene. CAI is often used to evaluate 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 genes of prokaryotic origin is lower than that of the host core genes, but it is getting closer and closer to the CAI value of the host core genes over time. The CAI of genes of eukaryotic origin is always very close to that of the host core genes.

[0075] Fop refers to the ratio of the optimal codon to its synonymous codon. The higher the Fop, the more frequently the optimal codon is used. Over time, the Fop of prokaryotic genes has become higher and higher, while the Fop of eukaryotic genes has always been very close to that of the host core genes.

[0076] 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 equal synonymous codons in 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, the value of NC ranges 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. The Nc values ​​of genes transferred from eukaryotes and prokaryotes to plants both show a decreasing trend over time and become increasingly closer to core genes.

[0077] (4) Gene function annotation and enrichment

[0078] The present invention collects annotation information from 9 major databases, namely OG1157, eggNOG, PANTHER, SuperFamily, Interproscan, GO, KEGG, Pfam and KO. The present invention predicts the function of the query sequence based on the sequence with the smallest e-value of the query sequence in the BLAST comparison. In addition, the present invention also uses the eggNOG Mapper tool (European Molecular Biology Laboratory; http: / / eggnog-mapper.embl.de) 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, wherein the background database 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 use an unsupervised hidden Markov clustering algorithm to cluster the homology group association network 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: Retrieve the topological structure of the query sequence and the distant sequences in genetic distance in the single gene phylogenetic tree constructed with the query sequence as the center; screen out potential cross-species HGT by defining the nested position and setting the minimum support node threshold; at the same time, record the most similar matching sequence of the query sequence in the recipient group and the donor group; 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 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 step 2, 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, characterized in that: The step 5 also includes analyzing potential cross-species HGT using homology comparison features of genes on both sides of 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 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 were physically closely linked and belonged to the same gene family, they would be considered as a single cross-species HGT event, and all these adjacent potential cross-species HGTs would be retained.

Images