A method for 16s and metagenomic sequencing data analysis

Abstract

The application provides a 16S and metagenome sequencing data analysis method, relates to the technical field of gene sequencing data analysis, and comprises the following steps: obtaining 16S sequencing original data and metagenome sequencing original data of a sample to be analyzed, performing quality filtering on the 16S sequencing original data to obtain 16S high-quality sequences, and performing quality filtering on the metagenome sequencing original data to obtain metagenome high-quality sequences; the method is characterized in that gene sequences and amplicon sequence variant representative sequences are compared and sequence homologous mapping relationships are established, phylogenetic distance is transferred from 16S data to metagenome genes, species evolution information and functional gene information can be aligned in the same phylogenetic coordinate system, and the problem that two data evolution dimensions cannot be effectively unified in the prior art is solved.

Description

Technical Field

[0001] This invention relates to the field of gene sequencing data analysis technology, specifically a method for analyzing 16S and metagenomic sequencing data. Background Technology

[0002] In the process of omics bioinformatics analysis accompanying clinical diagnosis of human microbiota, 16S rRNA gene sequencing and metagenomic sequencing have become core technologies for revealing the structure and function of microbial communities. Early on, microbiota analysis mainly relied on single 16S sequencing for species composition and diversity assessment, or metagenomic sequencing alone for gene and functional annotation, with each forming a relatively independent data processing workflow. With the increasing clinical demand for precise microbiota typing, functional prediction, and disease association analysis, researchers have begun to explore the joint interpretation of 16S species information and metagenomic functional information to obtain more comprehensive community characteristics.

[0003] Although existing technologies can perform 16S amplicon sequence variant clustering, phylogenetic tree construction, metagenomic assembly, gene prediction and abundance calculation, and achieve preliminary integration through simple correlation analysis or matrix splicing in some processes, the following problems still exist in practical applications:

[0004] Existing technologies often employ separate pipelines when processing 16S and metagenomic data, lacking technical means to transfer 16S phylogenetic distances to metagenomic genes through sequence homology mapping. This results in the two types of data being unable to be effectively aligned in the evolutionary dimension, making it difficult to uniformly express species and functional information within the same phylogenetic framework.

[0005] Existing technologies for processing metagenomic gene abundance typically remain at the level of raw coverage depth or simple standardization. They do not perform weighted interpolation reconstruction of adjacent gene abundance based on phylogenetic distance, and cannot utilize the evolutionary structural constraints provided by 16S to smooth and project functional abundance. As a result, functional characteristics are greatly affected by sequencing depth fluctuations and gene redundancy.

[0006] Existing technologies often directly perform simple splicing or correlation calculations when integrating 16S abundance matrices and metagenomic abundance matrices, without first performing central log ratio transformations on both before splicing them according to sample dimensions. This makes it difficult to effectively eliminate the inherent closure effect of component data, and the integrated matrix lacks stability in downstream statistical and machine learning analyses.

[0007] To address the aforementioned issues, there is an urgent need for a method for analyzing 16S and metagenomic sequencing data. This method would enable cross-data transfer of phylogenetic distances based on sequence homology mapping, gene abundance interpolation reconstruction based on phylogenetic distances, and sample dimensionality integration after center-log ratio transformation. This would improve the unified expression capability and reliability of integrated analysis of 16S and metagenomic data within a phylogenetic framework. Summary of the Invention

[0008] Technical problems to be solved

[0009] To address the shortcomings of existing technologies, this invention provides a method for analyzing 16S and metagenomic sequencing data, which solves the problems mentioned in the background section.

[0010] To achieve the above objectives, the present invention provides the following technical solution: a method for analyzing 16S and metagenomic sequencing data, comprising the following steps: S1. Obtain the 16S sequencing raw data and metagenomic sequencing raw data of the sample to be analyzed. Perform quality filtering on the 16S sequencing raw data to obtain high-quality 16S sequences and on the metagenomic sequencing raw data to obtain high-quality metagenomic sequences. S2. Cluster the 16S high-quality sequences into amplicon sequence variants and extract representative sequences of the amplicon sequence variants. Annotate the amplicon sequence variants based on a reference classification database. Simultaneously, construct a phylogenetic tree based on the representative sequences of the amplicon sequence variants and calculate the phylogenetic distance to obtain the amplicon sequence variant abundance matrix and the phylogenetic distance matrix. S3. Assemble the high-quality metagenomic sequences into contigs and predict open reading frames, construct a gene catalog and generate a gene abundance matrix; S4. Align the gene sequences in the gene catalog with the representative sequences of the amplicon sequence variants to establish a sequence homology mapping relationship, and transfer the phylogenetic distance in the phylogenetic distance matrix to the corresponding gene according to the sequence homology mapping relationship to obtain the gene phylogenetic coordinates; S5. Based on the gene phylogenetic coordinates, the abundance values ​​in the gene abundance matrix are reconstructed by interpolation based on phylogenetic distance to obtain the phylogenetic projected gene abundance matrix. S6. Perform central logarithmic transformation on the amplicon sequence variant abundance matrix and the phylogenetic projection gene abundance matrix respectively, and then splice them according to the sample dimensions to obtain the integrated abundance matrix.

[0011] Preferably, step S1 specifically includes the following steps: S1.1 Read the raw 16S sequencing data and metagenomic sequencing data of the sample to be analyzed, and record the sequencing platform, sequencing depth and sample identification information; S1.2 Perform base quality filtering, adapter sequence removal, and low-complexity sequence removal on the raw 16S sequencing data to obtain high-quality 16S sequences; S1.3 Perform base quality filtering, adapter sequence removal, and host genome alignment removal on the raw metagenomic sequencing data to obtain high-quality metagenomic sequences.

[0012] Preferably, step S2 specifically includes the following steps: S2.1 After correcting errors and removing chimeras from the 16S high-quality sequences, the corrected sequences are clustered into amplicon sequence variants according to sequence similarity, and representative sequences of the amplicon sequence variants are extracted. S2.2. Compare the representative sequences of the amplicon sequence variants with the reference taxonomic database to complete the hierarchical species annotation; S2.3. Construct a phylogenetic tree based on the representative sequences of the amplicon sequence variants, calculate the phylogenetic distance between amplicon sequence variants to obtain the phylogenetic distance matrix, and obtain the amplicon sequence variant abundance matrix by counting the sequence counts of each amplicon sequence variant in each sample and then performing depth normalization.

[0013] Preferably, step S3 specifically includes the following steps: S3.1. De novo assembly of the high-quality metagenomic sequence to obtain contigs, open reading frame prediction of the contigs, and output of candidate gene sequences; S3.2. Perform redundancy-free clustering on the candidate gene sequences to construct a non-redundant gene catalog; S3.3. The high-quality metagenomic sequence is back-attached to the non-redundant gene catalog, and the coverage depth of each gene is calculated and standardized to obtain the gene abundance matrix.

[0014] Preferably, step S4 specifically includes the following steps: S4.1. Use a sequence alignment algorithm to align each gene sequence in the gene catalog with a representative sequence of the amplicon sequence variant, and record the alignment coverage and sequence similarity. S4.2 Establish sequence homology mapping relationships between genes and amplicon sequence variants based on alignment coverage and sequence similarity; S4.3. Based on the sequence homology mapping relationship, the phylogenetic distance of the corresponding amplicon sequence variant in the phylogenetic distance matrix is ​​assigned to the successfully mapped gene to form gene phylogenetic coordinates.

[0015] Preferably, step S5 specifically includes the following steps: S5.1 Determine the position of each gene on the phylogenetic tree based on the gene phylogenetic coordinates; S5.2. Using phylogenetic distance as a scale, interpolate the abundance values ​​of genes at adjacent positions in the gene abundance matrix. S5.3 Replace the original abundance values ​​with the interpolated abundance values ​​to obtain the phylogenetic projected gene abundance matrix.

[0016] Preferably, step S6 specifically includes the following steps: S6.1 Perform central logarithmic transformation on the amplicon sequence variant abundance matrix and the phylogenetic projective gene abundance matrix, respectively; S6.2. The two transformed matrices are horizontally concatenated according to the sample identifiers to obtain the integrated abundance matrix.

[0017] Preferably, the phylogenetic tree in S2.3 is constructed using the maximum likelihood method, and the phylogenetic distance is calculated using the branch length path distance.

[0018] Preferably, the interpolation calculation in S5.2 involves selecting the nearest neighboring gene in phylogenetic distance for each gene, using the reciprocal of the phylogenetic distance between each neighboring gene and the target gene as the weight, and performing a weighted average of the abundance values ​​of the neighboring genes. The weighted average result is then used as the interpolated abundance value of the gene.

[0019] This invention provides a method for analyzing 16S and metagenomic sequencing data. It offers the following advantages: 1. By comparing gene sequences with representative sequences of amplicon sequence variants and establishing sequence homology mapping, the phylogenetic distance was transferred across data from 16S data to metagenomic genes, enabling species evolutionary information and functional gene information to be aligned in the same phylogenetic coordinate system, thus solving the problem that the two data evolutionary dimensions cannot be effectively unified in existing technologies.

[0020] 2. Based on gene phylogenetic coordinates, the abundance of adjacent genes is weighted and averaged by using the inverse of phylogenetic distance as the weight. This achieves gene abundance reconstruction based on phylogenetic structure, effectively utilizes 16S evolutionary information to constrain and smooth functional abundance, and reduces the interference of sequencing depth fluctuations and gene redundancy on functional characteristics.

[0021] 3. The amplicon sequence variant abundance matrix and the phylogenetic projected gene abundance matrix were transformed by the central log ratio and then spliced ​​according to the sample dimension. This eliminated the closure effect of the component data and generated a structurally stable integrated abundance matrix, which improved the reliability of subsequent statistical tests and machine learning analysis. Attached Figure Description

[0022] Figure 1 This is the main flowchart of the method of the present invention; Figure 2 This is a flowchart of step S1 of the present invention; Figure 3 This is a flowchart of step S2 of the present invention; Figure 4 This is a flowchart of step S3 of the present invention; Figure 5 This is a flowchart of step S4 of the present invention; Figure 6 This is a visual representation of the method of the present invention. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Specific Implementation Example 1: This paper presents a method for analyzing 16S and metagenomic sequencing data. The entire method consists of six main steps, such as... Figure 1 The main flowchart of the method is shown below; the specific steps are as follows: like Figure 2 As shown, a flowchart for step S1, which is the first step, is provided. The details are as follows: Step 1: Obtain the raw 16S sequencing data and raw metagenomic sequencing data of the sample to be analyzed. Perform quality filtering on the raw 16S sequencing data to obtain high-quality 16S sequences, and perform quality filtering on the raw metagenomic sequencing data to obtain high-quality metagenomic sequences. In this step, the raw 16S sequencing data and metagenomic sequencing data of the samples to be analyzed were obtained from the same batch of human microbiome clinical samples. Sample collection sites included the intestine, oral cavity, skin, or reproductive tract. After collection, total DNA was extracted according to standard operating procedures and then divided into two separate sequencing datasets. One dataset was used for 16S rRNA gene amplicon sequencing. Primers were designed targeting specific variable regions of the 16S rRNA gene for PCR amplification. The amplified products were purified and sequenced to obtain raw 16S sequencing data reflecting the microbial species composition. The other dataset was used for whole-genome shotgun metagenomic sequencing. The extracted total DNA was randomly fragmented, a sequencing library was constructed, and then sequenced to obtain raw metagenomic sequencing data containing all microbial genome fragments. Both datasets were stored in FASTQ format, with each sample corresponding to an independent file and accompanied by sample registration information.

[0025] The detailed steps for the first step are as follows: First, the raw 16S sequencing data and metagenomic sequencing data of the samples to be analyzed are read separately, and the sequencing platform, sequencing depth, and sample identification information are recorded. During reading, the corresponding FASTQ files are opened one by one according to the sample identification order, and the sequence identifiers, base sequences, and base quality values ​​in the files are all loaded into memory. Simultaneously, a sample metadata table is created, recording three key pieces of information in each entry: sequencing platform, sequencing depth, and sample identification. The sequencing platform record records the actual instrument model from which the data was generated, used for subsequent selection of adapter and primer sequences matching that platform. The sequencing depth record records the final number of raw sequences obtained for each sample, used for subsequent depth normalization calculations. A unique code is used for the sample identification to ensure that 16S data and metagenomic data from the same source correspond throughout the entire process, avoiding cross-contamination. After recording, the metadata table is saved for use in all subsequent steps.

[0026] Then, base quality filtering, adapter sequence removal, and low-complexity sequence removal were performed on the raw 16S sequencing data to obtain high-quality 16S sequences. The specific steps are as follows: (1) Based on the sequencing platform information recorded in the sample metadata table, retrieve the standard adapter sequence and amplification primer sequence corresponding to that platform. Perform local alignment of both ends of each original sequence with the retrieved adapter and primer sequences, respectively. During alignment, scan base by base from the end of the sequence inward. When the number of consecutive matching bases reaches the set matching length requirement, determine that the segment is an adapter or primer remnant and completely remove it from the sequence. After removal, only the unmatched insert fragment in the middle is retained to form the adapter-free sequence.

[0027] (2) Perform a base-by-base quality check on each sequence after the adapter removal. Read the quality value character corresponding to each base in the FASTQ file and convert the character into the corresponding numerical quality score according to the Phred quality encoding rules. Starting from the first base on the left end of the sequence, set a sliding window of fixed length, move one base to the right each time, and calculate the arithmetic mean of the quality scores of all bases in the current window. If the average quality score of a certain window is lower than the preset window quality threshold, then all bases on the right side of the sequence are truncated from the beginning of that window. At the same time, count the number of bases in the entire sequence whose quality score is lower than the single base quality threshold, calculate the proportion of this number to the total length of the sequence, and if the proportion exceeds the preset upper limit of the proportion of low-quality bases, then the entire sequence is discarded directly. The sequence retained after the above truncation and discarding operations is called the quality filtered sequence.

[0028] (3) Calculate the complexity index for the quality filtering sequences and remove low-complexity sequences. Two parallel determination methods are used for complexity calculation. The first method is the highest base frequency method: The frequency of each of the four bases (adenine, thymine, cytosine, and guanine) in the sequence is counted, and the base with the highest frequency is identified. This frequency is divided by the total sequence length to obtain the highest base frequency percentage. If this percentage is higher than a pre-set frequency threshold, the sequence is determined to be a low-complexity sequence. The second method is the Shannon entropy method: The frequency of each of the four bases is counted, and the negative logarithm of each frequency value is calculated. The negative logarithm is then multiplied by the corresponding frequency, and the four products are summed to obtain the Shannon entropy value. If the Shannon entropy value is lower than a pre-set entropy threshold, the sequence is also determined to be a low-complexity sequence. Any sequence determined to be low-complexity by either method is removed.

[0029] (4) Length screening of the remaining sequences, check the total number of bases in each sequence, and retain only the sequences whose length is greater than or equal to the lower limit and less than or equal to the upper limit, and finally obtain the 16S high-quality sequence.

[0030] Finally, base quality filtering, adapter sequence removal, and host genome alignment were performed on the raw metagenomic sequencing data to obtain high-quality metagenomic sequences. The specific steps are as follows: First, based on the sequencing platform information recorded in the sample metadata table, the standard adapter sequence corresponding to that platform is retrieved. Both ends of each original sequence are locally aligned with the retrieved adapter sequence. During alignment, the sequence is scanned base by base from the very end inwards. When the number of consecutive matching bases reaches the set matching length requirement, that segment is determined to be adapter remnant and is completely removed from the sequence. After removal, only the unmatched insert fragments in the middle are retained, forming the adapter-free sequence.

[0031] Second, each sequence after connector removal undergoes a base-by-base quality check. The quality value character corresponding to each base in the FASTQ file is read and converted into a numerical quality score according to Phred quality encoding rules. Starting from the first base on the left of the sequence, a fixed-length sliding window is set, moving one base to the right each time, and the arithmetic mean of the quality scores of all bases within the current window is calculated. If the average quality score of a window is lower than a pre-set window quality threshold, all bases on the right side of the sequence are truncated from the beginning of that window. Simultaneously, the number of bases in the entire sequence with quality scores lower than the individual base quality threshold is counted, and the proportion of this number to the total sequence length is calculated. If the proportion exceeds a pre-set upper limit for the proportion of low-quality bases, the entire sequence is discarded. The sequence retained after the above truncation and discarding operations is called the quality-filtered sequence.

[0032] Third, the quality-filtered sequences are aligned with the pre-prepared host reference genome. First, the host reference genome is constructed as an alignment index. Then, each quality-filtered sequence is used as a query sequence, and a seed-expansion alignment strategy is employed for alignment with the index: a short, fixed-length seed fragment is extracted from the query sequence, and possible matching positions are quickly located in the index. Alignment is then extended outwards from these positions, calculating the number of matching bases, mismatched bases, and gaps in the alignment region. After alignment, the sequence is determined to have successfully aligned to the host reference genome based on the proportion of matching bases to the total length and the number of mismatches and gaps. Sequences deemed successfully aligned are marked as host-origin sequences and removed, retaining only sequences that failed to align to the host reference genome.

[0033] Fourth, the complexity index of the retained non-host sequences is calculated, and low-complexity sequences are eliminated. Two parallel determination methods are used for complexity calculation. The first method is the highest base frequency method: the frequency of each of the four bases (adenine, thymine, cytosine, and guanine) in the sequence is counted, the most frequent base is identified, and this frequency is divided by the total sequence length to obtain the highest base frequency percentage. If this percentage is higher than a pre-set frequency threshold, the sequence is determined to be a low-complexity sequence. The second method is the Shannon entropy method: the frequency of each of the four bases is counted, the negative logarithm of each frequency value is calculated, the negative logarithm is multiplied by the corresponding frequency, and the four products are summed to obtain the Shannon entropy value. If the Shannon entropy value is lower than a pre-set entropy threshold, the sequence is also determined to be a low-complexity sequence. Sequences determined to be low-complexity by either method are eliminated.

[0034] Fifth, the remaining sequences are filtered by length. The total number of bases in each sequence is checked, and only sequences with a length greater than or equal to the lower limit and less than or equal to the upper limit are retained to obtain high-quality metagenomic sequences.

[0035] In a practical application scenario, taking a microecological testing project for inflammatory bowel disease at a top-tier hospital as an example, researchers collected fecal samples from patients and healthy controls. After DNA extraction, they performed 16S amplicon sequencing and metagenomic sequencing, respectively. After inputting the obtained raw FASTQ file and sample registration form into the first step of this method, the system automatically completed data reading, platform depth and labeling recording, and output high-quality 16S sequences and high-quality metagenomic sequences, respectively, after adapter removal, quality filtering, low-complexity removal, and host removal. Subsequent analysts directly used these two sets of high-quality sequences for species annotation and gene function mining, significantly reducing false positive results caused by adapter contamination, low-quality bases, and host DNA interference.

[0036] The first step involved acquiring the raw 16S and metagenomic data, recording key information, and performing systematic quality filtering. This step serves as the data entry point and quality control stage for the entire methodology. Its significance lies in strictly linking the two sets of data through a unified sample identifier and removing adapter residues, low-quality bases, low-complexity sequences, and host sequences. This provides clean and reliable input data for subsequent amplicon sequence variant clustering, phylogenetic tree construction, and gene catalog generation, ensuring that the entire integration analysis is built upon a high-quality sequence foundation.

[0037] like Figure 3 As shown, a flowchart for step S2, which is the second step, is provided. The details are as follows: Step 2: Cluster the high-quality 16S sequences into amplicon sequence variants and extract representative sequences of the amplicon sequence variants. Annotate the amplicon sequence variants based on a reference classification database. Simultaneously, construct a phylogenetic tree based on the representative sequences of the amplicon sequence variants and calculate the phylogenetic distance to obtain the amplicon sequence variant abundance matrix and phylogenetic distance matrix. The detailed steps for the second step are as follows: First, after error correction and chimera removal of the high-quality 16S sequences, the corrected sequences were clustered into amplicon sequence variants based on sequence similarity, and representative sequences of the amplicon sequence variants were extracted. The specific steps are as follows: (1) Error Correction Process: After arranging all high-quality 16S sequences in the same direction, each base position is statistically analyzed independently. During the analysis, all sequences are traversed, and the frequency of adenine, thymine, cytosine, and guanine at each position is recorded and their respective frequencies are calculated. If the frequency of a certain base at a certain position is significantly higher than that of the other three, and this frequency exceeds a pre-set dominant frequency threshold, then the correct base at that position is determined to be the dominant base. Subsequently, the actual base at that position in each sequence is checked. If it is inconsistent with the determined correct base, then the base at that position in that sequence is replaced with the correct base. After completing the above statistics and replacement for all positions, the corrected sequence set is obtained.

[0038] (2) Chimera Removal Process: For each corrected sequence, it is split into a left and right half at the middle site. The left half is locally aligned with the left half of all other sequences, and the top few sequences with the highest number of matching bases are recorded; similarly, the right half is locally aligned with the right half of all other sequences, and the top few sequences with the highest number of matching bases are recorded. If the sequence with the highest matching degree in the left half and the sequence with the highest matching degree in the right half do not originate from the same original sequence, the current sequence is determined to be a chimera and is removed from the set. After all sequences have been detected, non-chimera sequences are retained.

[0039] (3) Clustering process: Pairwise alignment of the retained non-chimeric sequences is performed. A dynamic programming method is used to construct a score matrix during alignment. Matching bases receive positive scores, mismatched bases receive negative scores, and gap introduction and extension are penalized differently. After backtracking the score matrix to obtain the optimal alignment path, the number of matching bases on the path and the total length of the alignment region are counted to calculate the sequence similarity. The calculation formula is as follows: ; in Indicates the similarity between two sequences; This indicates the number of bases that completely match after the two sequences are aligned. The formula represents the total length of the alignment region. It directly reflects the degree of consistency between two sequences after alignment by the ratio of the number of matched bases to the total length of the alignment region. The higher the ratio, the more identical bases the two sequences have at corresponding positions. Therefore, it can objectively measure the similarity between sequences and thus serve as the basis for classifying them into the same amplicon sequence variant.

[0040] When the similarity between two sequences exceeds a preset clustering threshold, they are grouped into the same amplicon sequence variant. Clustering employs a hierarchical aggregation method. Initially, each sequence forms an independent cluster. Each time, the two clusters with the highest current similarity are merged until the sequence similarity within all clusters is no less than the threshold. After clustering, the frequency of each sequence within each amplicon sequence variant is counted, and the sequence with the highest frequency is selected as the representative sequence of that variant. If the frequencies are the same, the longest sequence is selected.

[0041] Then, representative sequences of amplicon sequence variants are compared with reference taxonomic databases to complete hierarchical species annotation; The specific operation process is as follows: A representative sequence from each amplicon sequence variant is used as the query sequence, and a global alignment is performed with each known sequence in the reference classification database. The global alignment also employs a dynamic programming method to construct a complete score matrix, allowing gaps at both ends and within the sequence. After alignment, the number of matching bases, mismatched bases, and gaps are counted, and the similarity is calculated using the following formula: ; in Similarity; To match the number of bases; The number of mismatched bases; The formula uses the number of matching bases as the numerator and the sum of matches, mismatches, and gaps as the denominator. This allows it to consider both consistency and differences on similarity, avoiding bias caused by using only the number of matches, and thus more accurately reflecting the overall similarity between the query sequence and the database sequence.

[0042] For each query sequence, the sequence with the highest similarity is selected from all database sequences. The complete taxonomic path from kingdom to species, pre-annotated in the database sequence, is read and assigned to the corresponding amplicon sequence variant, thus completing the hierarchical species annotation.

[0043] Finally, a phylogenetic tree was constructed based on representative sequences of amplicon variants. The phylogenetic distance between amplicon variants was calculated to obtain the phylogenetic distance matrix, and the sequence counts of each amplicon variant in each sample were depth-normalized to obtain the amplicon variant abundance matrix. The phylogenetic tree was constructed using the maximum likelihood method, and the phylogenetic distance was calculated using branch length path distance. The specific process is as follows: (1) Phylogenetic tree construction: First, multiple sequence alignment was performed on representative sequences of all amplicon sequence variants. A progressive alignment strategy was adopted, first aligning the most similar sequence pairs, and then gradually adding the remaining sequences to obtain a site alignment matrix. Based on the alignment matrix, the frequency of base substitution at each site was statistically analyzed, and a suitable substitution model was selected. Under this model, a heuristic search method was used to traverse possible tree topologies. For each topology, the likelihood value of the current observation data was calculated, and the topology with the largest likelihood value and its corresponding branch length were retained to obtain the final phylogenetic tree.

[0044] (2) Phylogenetic distance calculation: For any two leaf nodes in the phylogenetic tree, find the unique path connecting the two leaf nodes, and add up the length values ​​of each branch on the path. The sum is the phylogenetic distance between the two amplicon sequence variants. Arrange all pairwise distances in the order of the row and column corresponding amplicon sequence variants to form a phylogenetic distance matrix.

[0045] (3) Abundance matrix generation: Count the number of sequences belonging to each amplicon sequence variant in each sample to obtain the original counting matrix. The depth normalization calculation formula is as follows: ; in, For the sample amplicon sequence variants Normalized abundance; For the original count; For the sample The total number of sequences, To achieve a uniform scaling factor, the formula first divides the original count by the total depth of the sample to eliminate the influence of differences in sequencing volume among different samples, and then multiplies it by a uniform scaling factor to make the abundance values ​​of different samples comparable. Therefore, it can obtain normalized abundance that can be directly used for subsequent analysis.

[0046] The normalized values ​​form the amplicon sequence variant abundance matrix.

[0047] In a practical application scenario, the gastroenterology department of a top-tier hospital conducted a comparative study of gut microbiota between patients with inflammatory bowel disease and healthy individuals. After collecting fecal samples and performing 16S sequencing, the research team directly imported the high-quality 16S sequences obtained in the first step into this process. The system automatically performed error correction, chimera removal, amplicon sequence variant clustering, species annotation, and phylogenetic tree construction, outputting an amplicon sequence variant abundance matrix with complete taxonomic hierarchy and quantitative evolutionary distance. Clinical researchers used this matrix to quickly identify significantly elevated or decreased gut microbiota groups in the patient group and provided standardized species-side input for subsequent integration with metagenomic functional data, significantly shortening the time from raw sequences to interpretable species results.

[0048] The second step, which completes the transformation from high-quality 16S sequences to amplicon sequence variant abundance and phylogenetic distance matrices, is a core step in obtaining species composition and evolutionary relationships. Its significance lies in transforming the original sequences into structured data with clear taxonomic positions and quantitative evolutionary distances through meticulous error correction, chimera removal, similarity clustering, taxonomic annotation, and maximum likelihood phylogenetic reconstruction. This provides the necessary foundation for subsequent phylogenetic information transfer to metagenomic genes.

[0049] like Figure 4 As shown, a flowchart for step S3, which is the third step, is provided. The details are as follows: Step 3: Assemble high-quality metagenomic sequences into contigs and predict open reading frames, construct a gene catalog and generate a gene abundance matrix; The detailed steps for the third step are as follows: First, high-quality metagenomic sequences are de novo assembled to obtain contigs. Open reading frames are then predicted from these contigs to output candidate gene sequences. The specific steps are as follows: Assembly process: The entire metagenomic high-quality sequence is decomposed into fixed-length short fragments. Each short fragment is treated as a node, and the overlap relationships between short fragments are used as directed edges to construct a de Bruin graph. Starting from a node with an in-degree of zero or an out-degree greater than its in-degree, an Eulerian path is sought along the edge direction. Short fragments along the path are connected sequentially to form longer, continuous sequences, i.e., contigs. After traversing all possible paths, all contigs are output.

[0050] Open reading frame prediction process: For each contiguous group, scan both the forward and reverse strands. During scanning, check each possible start codon position sequentially, and read downstream from that position until a stop codon is encountered. If the base length between the start and stop codons exceeds a pre-defined minimum gene length, the base sequence of that interval is extracted as a candidate gene sequence. After scanning all contiguous groups, output all candidate gene sequences.

[0051] Then, the candidate gene sequences are clustered to remove redundancy, constructing a non-redundant gene catalog. The specific steps are as follows: ① Perform pairwise alignments on all candidate gene sequences using the same dynamic programming method as the sequence similarity calculation in step two. Construct a score matrix, awarding positive points for matching bases and negative points for mismatched bases. Different penalties are applied for gap introduction and extension. After backtracking the score matrix to obtain the optimal alignment path, count the number of matching bases on the path and the total length of the alignment region. The similarity calculation formula is as follows: ; in Similarity; To match the number of bases; The length of the alignment region is denoted by 1. This formula directly reflects the degree of similarity between two sequences after alignment by the ratio of the number of matched bases to the total length of the alignment region. The higher the ratio, the more identical bases the two sequences have at corresponding positions, thus objectively measuring the level of similarity between sequences.

[0052] When the similarity between two candidate gene sequences is higher than a preset redundancy removal threshold, they are grouped into the same gene cluster.

[0053] ② Within each gene cluster, the longest sequence is selected as the representative sequence, and the remaining sequences are marked as redundant and deleted. All representative sequences are arranged according to a uniform number to form a non-redundant gene catalog.

[0054] Finally, the high-quality metagenomic sequences were back-linked to a non-redundant gene catalog, and the coverage depth of each gene was calculated and normalized to obtain a gene abundance matrix. The specific steps are as follows: (1) Alignment process: Each high-quality metagenomic sequence is locally aligned with all gene sequences in the non-redundant gene catalog. The alignment adopts a seed-expansion strategy, first extracting short seeds from the query sequence and locating them in the gene sequence, then extending them to both sides and allowing for limited mismatches. After a successful alignment, the starting and ending coordinates of the query sequence on the gene are recorded.

[0055] (2) Coverage depth calculation: For each gene, a position count array is established over its entire length, with all initial values ​​set to zero. Whenever a sequence aligns to a region of the gene, the count values ​​of all positions within that region are incremented by one. After all sequences have been processed, the sum of all values ​​in the position count array is divided by the gene length to obtain the average coverage depth of the gene.

[0056] (3) Calculate the standardized abundance, where the calculation formula is as follows: ; in For the sample Zhonggen Standardized abundance, The original average coverage depth, This represents the total number of genes in the non-redundant gene catalog. To unify the scaling factor, the formula first divides the coverage depth of each gene by the sum of the coverage depths of all genes in the sample, so that the sum of abundance within the same group is a fixed value, and then multiplies it by the unified scaling factor, thereby eliminating the influence of differences in the total sequencing volume of different samples and making gene abundance comparable across samples.

[0057] (4) The standardized values ​​are arranged according to the sample and the gene to obtain the gene abundance matrix.

[0058] In a practical application scenario, the Department of Infectious Diseases at a tertiary hospital conducted a study on changes in intestinal functional genes in patients with Clostridium difficile infection. After completing metagenomic sequencing and obtaining high-quality sequences, the research team input the data into this step. The system automatically performed contig assembly, open reading frame prediction, gene redundancy removal, and coverage depth normalization, outputting a gene abundance matrix. Clinical researchers directly used this matrix to compare the abundance differences of key drug resistance genes and toxin-related genes between patients and healthy controls, providing a quantitative basis for subsequent integration of functional genes with 16S phylogenetic information.

[0059] Summary of this step: Step three completes the transformation from high-quality metagenomic sequences to a gene abundance matrix, which is the core step in obtaining quantitative information on functional genes. Its significance lies in transforming the raw metagenomic sequence into structured gene abundance data through de novo assembly, open reading frame identification, redundancy removal clustering, and coverage depth normalization, providing functional input for subsequent integration with the 16S phylogenetic distance.

[0060] like Figure 5 As shown, a flowchart for step S4, which is the fourth step, is provided. The details are as follows: Step 4: Align the gene sequences in the gene catalog with the representative sequences of amplicon sequence variants to establish sequence homology mapping relationships. Based on the sequence homology mapping relationships, transfer the phylogenetic distance in the phylogenetic distance matrix to the corresponding genes to obtain the gene phylogenetic coordinates. The detailed steps for step four are as follows: First, a sequence alignment algorithm is used to align each gene sequence in the gene catalog with a representative sequence of amplicon variants, and the alignment coverage and sequence similarity are recorded. The specific steps are as follows: For each gene sequence, it is used as the query sequence and locally compared with representative sequences of all amplicon sequence variants one by one. The alignment employs a seed-expansion strategy: first, a fixed-length short seed fragment is extracted from the query sequence, and all perfectly matching seed positions are quickly located in representative sequences of the target amplicon sequence variant; then, extension is performed from each seed position outwards. During extension, for each base or gap added, a score is accumulated based on the matching situation: a matching base adds a fixed positive score, a mismatch deducts a fixed penalty score, introducing a gap deducts a higher penalty score, and gap extension deducts a lower penalty score. When the accumulated score falls below a pre-set stopping threshold, extension in that direction is stopped. After extension, the alignment interval between the query sequence and the target sequence is obtained. Within the alignment interval, the number of perfectly matching bases and the total number of bases involved in the alignment are counted. The formula is: ; Calculate the alignment coverage; where For coverage; This is the number of bases in the query sequence that participated in the alignment; The total length of the query sequence is given by the formula, which divides the number of bases involved in the alignment by the total length of the query sequence. This directly reflects the proportion of the query sequence that was successfully aligned. The higher the proportion, the more complete the alignment region.

[0061] According to the formula: ; Calculate sequence similarity, where Similarity; This represents the number of completely matching bases within the alignment interval. This represents the number of bases involved in the alignment. The formula divides the number of matching bases by the number of bases involved in the alignment, reflecting the degree of homology within the aligned region; a higher value indicates stronger homology.

[0062] Record the coverage and similarity of each gene with each amplicon sequence variant.

[0063] Then, sequence homology mapping relationships between genes and amplicon sequence variants are established based on alignment coverage and sequence similarity. The specific steps are as follows: For each gene, its coverage and similarity records with all amplicon sequence variants are iterated. If the coverage of a gene with an amplicon sequence variant is higher than a preset coverage threshold and the similarity is higher than a preset similarity threshold, the gene is considered to have sequence homology with that amplicon sequence variant. If a gene satisfies both conditions with multiple amplicon sequence variants, the similarity values ​​of these variants are compared, and the one with the highest similarity is selected to establish a mapping. All successfully established mapping relationships are summarized into a sequence homology mapping table, where each record contains a gene identifier, the corresponding amplicon sequence variant identifier, coverage, and similarity.

[0064] Finally, based on sequence homology mapping, the phylogenetic distances of the corresponding amplicon sequence variants in the phylogenetic distance matrix are assigned to the successfully mapped genes, forming gene phylogenetic coordinates. The specific steps are as follows: Read the sequence homology mapping table and locate the corresponding amplicon variant for each successfully mapped gene. Extract all distance values ​​from the phylogenetic distance matrix for the row containing the amplicon variant; these values ​​represent the phylogenetic distances between the amplicon variant and all other amplicon variants. Completely copy this set of distance values ​​to the gene as its phylogenetic coordinates. After assigning values ​​to all successfully mapped genes, a set of gene phylogenetic coordinates is obtained. Genes that were not successfully mapped are not assigned coordinates.

[0065] In a practical application scenario, the translational medicine center of a top-tier hospital conducted a study on the synergistic changes in the gut microbiota between species and function. After obtaining representative sequences of amplicon variants, a phylogenetic distance matrix, and a gene catalog, the research team simultaneously input these three sets of data into this step. The system automatically performed local alignment between genes and amplicon variants, two-condition homology determination, and phylogenetic distance transfer, outputting the phylogenetic coordinates of each successfully mapped gene. Based on this, the researchers placed previously independent functional genes into their corresponding positions on the 16S phylogenetic tree, achieving direct alignment between functional abundance and species evolutionary relationships, providing the necessary coordinate foundation for subsequent interpolation and reconstruction.

[0066] The fourth step completes the cross-data transfer of phylogenetic distance from amplicon sequence variants to genes, serving as a crucial bridge for aligning 16S and metagenomic information. Its significance lies in assigning evolutionary distance to functional genes through precise local alignment, dual-threshold homology determination, and complete distance replication, enabling species and functional information to be uniformly expressed within the same phylogenetic coordinate system.

[0067] Step 5: Based on the gene phylogenetic coordinates, the abundance values ​​in the gene abundance matrix are reconstructed by interpolation based on phylogenetic distance to obtain the phylogenetic projected gene abundance matrix; The detailed steps for step five are as follows: First, the position of each gene on the phylogenetic tree is determined based on gene phylogenetic coordinates; The specific work steps are as follows: Read the gene phylogenetic coordinates output from step four. Each gene's coordinates are a vector representing its distances to all amplicon sequence variants. In the constructed phylogenetic tree, each leaf node also corresponds to a distance vector to all other leaf nodes. For the target gene's coordinate vector, compare it element-wise with the distance vectors of each leaf node: subtract the values ​​at corresponding positions in the two vectors and take the absolute value, then sum all the absolute values ​​to obtain the total difference between the gene and that leaf node. After traversing all leaf nodes, select the leaf node with the smallest total difference value and record the branch position of that leaf node as the gene's phylogenetic position.

[0068] Then, using phylogenetic distance as a scale, interpolation is performed on the abundance values ​​of genes at adjacent positions in the gene abundance matrix. The interpolation calculation involves selecting the nearest phylogenetic neighbor for each gene, using the reciprocal of the phylogenetic distance between each neighboring gene and the target gene as the weight, and then averaging the abundance values ​​of the neighboring genes. The weighted average result is used as the interpolated abundance value for that gene. The specific steps are as follows: For each target gene, calculate its phylogenetic distance to all other genes in the gene abundance matrix, sort them by distance from smallest to largest, and select the top few genes with the smallest distance as neighboring genes. Calculate the weight for each neighboring gene using the following formula: ; in For the first The weight of each adjacent gene; For the target gene and the first The formula takes the inverse of the distance as the weight, and the closer the distance between adjacent genes, the greater the weight. This is consistent with the biological hypothesis that "genes that are phylogenetically closer have similar abundances". Therefore, it can reasonably allocate the contribution of different adjacent genes to the abundance of the target gene.

[0069] Then, a weighted average is calculated, as follows: ; in The interpolated abundance value of the target gene; The number of adjacent genes; As weight; For the first The original abundance values ​​of each neighboring gene in the gene abundance matrix are used. This formula first amplifies or reduces the abundance of neighboring genes using weights, then divides by the sum of the weights to obtain a weighted average abundance, thus incorporating phylogenetic neighbor relationships while preserving the original abundance information. The calculated... Record the interpolated abundance value of the target gene.

[0070] Finally, the interpolated abundance values ​​are used to replace the original abundance values ​​to obtain the phylogenetic projected gene abundance matrix. Details are as follows: Iterate through each gene position in the gene abundance matrix, replacing the original abundance value at that position with the interpolated abundance value calculated in the previous step. After all replacements are completed, maintain the original order of sample rows and gene columns to form a phylogenetic projective gene abundance matrix.

[0071] In a practical application scenario, taking a gut microbiota functional evolution analysis project at a top-tier hospital as an example, researchers input the gene phylogenetic coordinates output from step four and the gene abundance matrix output from step three into this step. The system automatically completes location positioning and distance-weighted interpolation, outputting a phylogenetic projection gene abundance matrix. The functional abundance in this matrix is ​​constrained by the 16S evolutionary structure, making it more suitable for joint analysis with species abundance.

[0072] The fifth step, gene abundance interpolation reconstruction based on phylogenetic distance, is the core step in optimizing functional abundance using evolutionary information. Its significance lies in the fact that by precisely selecting neighboring genes and using a distance-inverse weighted average, abundance information of phylogenetically adjacent genes is mutually transferred, reducing the impact of sequencing depth fluctuations and gene redundancy on the abundance of individual genes, thus making functional abundance smoother and more stable within the phylogenetic framework.

[0073] Step 6: Perform central log ratio transformation on the amplicon sequence variant abundance matrix and the phylogenetic projected gene abundance matrix respectively, and then splice them according to the sample dimensions to obtain the integrated abundance matrix; The detailed steps for step six are as follows: First, the amplicon sequence variant abundance matrix and the phylogenetic projected gene abundance matrix were subjected to central log ratio transformation. The specific steps are as follows: For each sample in the amplicon sequence variant abundance matrix, first extract all amplicon sequence variant abundance values ​​corresponding to that sample. Check if any of these abundance values ​​are zero; if so, replace them with a pre-defined, extremely small positive number to avoid mathematical errors when taking the logarithm later. Then calculate the geometric mean of these abundance values: multiply all abundance values ​​together, and then take the square root of the product with the same power as the number of abundance values ​​to obtain the geometric mean of that sample.

[0074] Next, a transformation is performed on each abundance value: the abundance value is divided by the geometric mean obtained earlier, and the natural logarithm of the quotient is taken to obtain the transformed value. For each sample in the phylogenetic projected gene abundance matrix, the exact same zero-value replacement, geometric mean calculation, and logarithmic transformation process is performed to obtain the transformed gene abundance value.

[0075] The central logarithmic ratio transformation serves to convert the original component abundance data into the real number space. The original abundance data suffers from a closure effect, meaning the sum of the abundances of all features is fixed, and direct analysis introduces spurious correlations. By dividing by the geometric mean and taking the logarithm, this closure constraint is eliminated, allowing the transformed values ​​to be used more freely for subsequent statistical and machine learning analyses, while preserving the key information about relative abundance from the original data.

[0076] Next, the two transformed matrices are horizontally concatenated according to the sample identifiers to obtain the integrated abundance matrix. The specific steps are as follows: Check the sample identifier lists of the two transformed matrices to ensure that the identifier content is completely identical and arranged in the same order. Use all transformed values ​​corresponding to amplicon sequence variants as the left column block of the new matrix, and all transformed values ​​corresponding to phylogenetic projection genes as the right column block, then directly concatenate them column-wise. The rows of the concatenated matrix correspond to samples, and the columns correspond to all amplicon sequence variant features and all gene features, resulting in a coherent abundance matrix.

[0077] In a practical application scenario, a clinical research center of a top-tier hospital constructed a multi-omics joint diagnostic model for patients with inflammatory bowel disease. After obtaining the amplicon sequence variant abundance matrix and the phylogenetic projected gene abundance matrix, the research team input both matrices into this step. The system automatically performed central logarithmic ratio transformation and sample dimension concatenation, outputting an integrated abundance matrix. This matrix simultaneously contains evolutionarily constrained species abundance and functional abundance information. Researchers directly used this as a feature input to a machine learning classifier to construct a model that distinguishes between patients and healthy controls. The model performance was significantly improved compared to using only a single omics dataset, providing a more comprehensive data foundation for clinical auxiliary diagnosis.

[0078] Step six, which completes the joint transformation and integration of amplicon sequence variant abundance and phylogenetic projected gene abundance, is the final output of the entire method. Its significance lies in eliminating the inherent closure effect of component data through central logarithmic ratio transformation, and then unifying species and functional information into a single matrix through rigorous sample label alignment and lateral splicing. This provides downstream statistical analysis and machine learning with structurally stable and informationally complete integrated data.

[0079] It should be noted that this method can be integrated and deployed on existing online analysis platforms. A dedicated visual integration console can be developed within the platform to host and run the processing flow of this method, such as... Figure 6 As shown in the diagram, a visual interface is provided. The diagram shows the data input status visualization interface. After dragging the two sets of raw data into the corresponding file fields, click "Start Analysis," and the platform can analyze the data according to the steps described in this method.

[0080] In summary, this method starts with raw 16S sequencing data and metagenomic sequencing data from the same batch of human microbiome samples. First, high-quality 16S sequences and high-quality metagenomic sequences are obtained through quality filtering. Then, error correction, chimera removal, similarity clustering, species annotation, and phylogenetic tree construction are performed on the high-quality 16S sequences to obtain an amplicon sequence variant abundance matrix and a phylogenetic distance matrix. Simultaneously, the high-quality metagenomic sequences are assembled, open reading frames are predicted, redundancy is removed, and coverage depth is normalized to obtain a gene abundance matrix. Subsequently, homology mapping between genes and amplicon sequence variants is established through sequence alignment, and phylogenetic distance is transferred to genes to form gene phylogenetic coordinates. Then, distance-weighted interpolation is performed on gene abundance based on these coordinates to reconstruct the phylogenetic projected gene abundance matrix. Finally, the amplicon sequence variant abundance matrix and the phylogenetic projected gene abundance matrix are subjected to central logarithmic ratio transformation and spliced ​​according to sample dimensions to output the integrated abundance matrix. The entire process achieves deep integration of 16S species evolution information and metagenomic functional abundance information within a unified phylogenetic framework, providing a structurally stable and information-complementary multi-omics data foundation for clinical diagnosis and research of human microecology.

[0081] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising a reference structure" does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes the element.

[0082] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for analyzing 16S and metagenomic sequencing data, characterized in that, Includes the following steps: S1. Obtain the 16S sequencing raw data and metagenomic sequencing raw data of the sample to be analyzed. Perform quality filtering on the 16S sequencing raw data to obtain high-quality 16S sequences and on the metagenomic sequencing raw data to obtain high-quality metagenomic sequences. S2. Cluster the 16S high-quality sequences into amplicon sequence variants and extract representative sequences of the amplicon sequence variants. Perform species annotation on the amplicon sequence variants based on a reference classification database. At the same time, construct a phylogenetic tree based on the representative sequences of the amplicon sequence variants and calculate the phylogenetic distance to obtain the amplicon sequence variant abundance matrix and the phylogenetic distance matrix. S3. Assemble the high-quality metagenomic sequences into contigs and predict open reading frames, construct a gene catalog and generate a gene abundance matrix; S4. Align the gene sequences in the gene catalog with the representative sequences of the amplicon sequence variants to establish a sequence homology mapping relationship, and transfer the phylogenetic distance in the phylogenetic distance matrix to the corresponding gene according to the sequence homology mapping relationship to obtain the gene phylogenetic coordinates; S5. Based on the gene phylogenetic coordinates, the abundance values ​​in the gene abundance matrix are reconstructed by interpolation based on phylogenetic distance to obtain the phylogenetic projected gene abundance matrix. S6. Perform central logarithmic transformation on the amplicon sequence variant abundance matrix and the phylogenetic projection gene abundance matrix respectively, and then splice them according to the sample dimensions to obtain the integrated abundance matrix.

2. The method for analyzing 16S and metagenomic sequencing data according to claim 1, characterized in that, S1 specifically includes the following steps: S1.1 Read the raw 16S sequencing data and metagenomic sequencing data of the sample to be analyzed, and record the sequencing platform, sequencing depth and sample identification information; S1.2 Perform base quality filtering, adapter sequence removal, and low-complexity sequence removal on the raw 16S sequencing data to obtain high-quality 16S sequences; S1.3 Perform base quality filtering, adapter sequence removal, and host genome alignment removal on the raw metagenomic sequencing data to obtain high-quality metagenomic sequences.

3. The method for analyzing 16S and metagenomic sequencing data according to claim 2, characterized in that, S2 specifically includes the following steps: S2.1 After correcting errors and removing chimeras from the 16S high-quality sequences, the corrected sequences are clustered into amplicon sequence variants according to sequence similarity, and representative sequences of the amplicon sequence variants are extracted. S2.

2. Compare the representative sequences of the amplicon sequence variants with the reference taxonomic database to complete the hierarchical species annotation; S2.

3. Construct a phylogenetic tree based on the representative sequences of the amplicon sequence variants, calculate the phylogenetic distance between amplicon sequence variants to obtain the phylogenetic distance matrix, and obtain the amplicon sequence variant abundance matrix by counting the sequence counts of each amplicon sequence variant in each sample and then performing depth normalization.

4. The method for analyzing 16S and metagenomic sequencing data according to claim 1, characterized in that, S3 specifically includes the following steps: S3.

1. De novo assembly of the high-quality metagenomic sequence to obtain contigs, open reading frame prediction of the contigs, and output of candidate gene sequences; S3.

2. Perform redundancy-free clustering on the candidate gene sequences to construct a non-redundant gene catalog; S3.

3. The high-quality metagenomic sequence is back-attached to the non-redundant gene catalog, and the coverage depth of each gene is calculated and standardized to obtain the gene abundance matrix.

5. The method for analyzing 16S and metagenomic sequencing data according to claim 1, characterized in that, S4 specifically includes the following steps: S4.

1. Use a sequence alignment algorithm to align each gene sequence in the gene catalog with a representative sequence of the amplicon sequence variant, and record the alignment coverage and sequence similarity. S4.2 Establish sequence homology mapping relationships between genes and amplicon sequence variants based on alignment coverage and sequence similarity; S4.

3. Based on the sequence homology mapping relationship, the phylogenetic distance of the corresponding amplicon sequence variant in the phylogenetic distance matrix is ​​assigned to the successfully mapped gene to form gene phylogenetic coordinates.

6. The method for analyzing 16S and metagenomic sequencing data according to claim 1, characterized in that, S5 specifically includes the following steps: S5.1 Determine the position of each gene on the phylogenetic tree based on the gene phylogenetic coordinates; S5.

2. Using phylogenetic distance as a scale, interpolate the abundance values ​​of genes at adjacent positions in the gene abundance matrix. S5.3 Replace the original abundance values ​​with the interpolated abundance values ​​to obtain the phylogenetic projected gene abundance matrix.

7. The method for analyzing 16S and metagenomic sequencing data according to claim 1, characterized in that, S6 specifically includes the following steps: S6.1 Perform central logarithmic transformation on the amplicon sequence variant abundance matrix and the phylogenetic projective gene abundance matrix, respectively; S6.

2. The two transformed matrices are horizontally concatenated according to the sample identifiers to obtain the integrated abundance matrix.

8. The method for analyzing 16S and metagenomic sequencing data according to claim 3, characterized in that, In S2.3, the phylogenetic tree is constructed using the maximum likelihood method, and the phylogenetic distance is calculated using branch length path distance.

9. A method for analyzing 16S and metagenomic sequencing data according to claim 6, characterized in that, The interpolation calculation in S5.2 involves selecting the nearest neighboring gene in phylogenetic distance for each gene, using the reciprocal of the phylogenetic distance between each neighboring gene and the target gene as the weight, and averaging the abundance values ​​of the neighboring genes. The weighted average result is then used as the interpolated abundance value of the gene.

Images