Automatic method for 16SrRNA amplicon analysis

Abstract

The invention provides a comprehensive 16S automatic analysis method and process based on Web, and from upstream sequencing data processing to downstream data analysis, the most common analysis items at present and 16S advanced analysis are included. And meanwhile, the analysis content is modularized, the data mining method and the like are diversified and deep, and different analysis modules can be selected according to individual needs. The method not only can help researchers without any biological information background to carry out 16S research on flora micro-ecology, but also can overcome the defects of some existing 16S analysis processes, can obtain an accurate and comprehensive analysis result, and meets various analysis requirements.

Description

Technical Field

[0001] The present invention relates to the field of intelligent automatic analysis, and in particular to an automated method for 16S rRNA amplicon analysis. Background Art

[0002] Currently, high-throughput sequencing technology is widely used in industry and scientific research and has become a routine research and detection method. The large-scale application of sequencing technology has led to a geometric increase in the amount of biological data, and how to analyze this sequencing data has become an urgent problem that needs to be solved.

[0003] To address these challenges, numerous high-throughput sequencing data analysis tools have been developed, and tools for analyzing microbial sequencing data have gradually matured and improved. Many tools are combined to complete microbial microecological analysis. These tools also require proficiency in other technical skills, such as Linux, Python, and R. This places high demands on researchers, requiring not only biological knowledge but also computer science. Some researchers with a biology background may also need to utilize sequencing technology to solve problems. Consequently, several web-based analysis tools have been developed. These tools do not require researchers to master sequencing data analysis skills; they simply upload the data, select the analysis requirements, and obtain the results.

[0004] The existing 16S rRNA amplicon analysis process is not comprehensive and does not update analysis content in real time. In addition, the original web-based analysis process only includes downstream analysis and does not include upstream sequencing data processing. Summary of the Invention

[0005] Based on the above problems in the prior art, the present invention provides an automated method for 16S rDNA amplicon analysis. The 16S rDNA sequencing data analysis process provided by the present invention uses the provided raw sequencing data as input, calls industry-standard analysis tools, and ultimately visualizes the data to obtain an easy-to-interpret and better analysis result.

[0006] The present invention relates to an automated method for 16S rDNA amplicon analysis, comprising the following steps: data import, generating a manifest file after a user uploads sequencing data, and then importing the data into a qiime2 object using QIIME2 software; after the sequencing data of the sample is converted into a qiime2 object, denoising is performed using the DADA2 method; by comparing representative sequences with reference sequences in corresponding databases, the taxonomic information corresponding to each ASV (Amplicon Sequence Variant) is obtained to achieve species annotation; performing multiple sequence alignment on the representative sequences and filtering the aligned highly variable regions, and then generating a phylogenetic tree and a rooted tree based on the filtered alignment results; Through the above steps, ASV taxonomic information, phylogenetic evolutionary trees of representative sequences, phylum or genus level feature tables, and other information, plus user-provided grouping and environmental factor information, can be obtained, which can be used as input files for subsequent data analysis steps.

[0007] The data analysis step includes at least one of the following analysis steps: species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, system evolution analysis, association statistical analysis, function prediction analysis, etc. The order does not need to be specifically limited.

[0008] In order to display the composition and structure of species, preferably, it can include species composition analysis, Alpha diversity analysis, Beta diversity analysis, and difference analysis; and can further include at least one of systematic evolution analysis, association statistical analysis, and functional prediction analysis; preferably, it can include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, systematic evolution analysis, association statistical analysis, and functional prediction analysis.

[0009] Preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, and phylogenetic analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, and association statistical analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, and function prediction analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, phylogenetic analysis, and association statistical analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, phylogenetic analysis, and function prediction analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, phylogenetic analysis, and association statistical analysis; preferably, it may include species composition analysis, Alpha diversity analysis, Beta diversity analysis, difference analysis, phylogenetic analysis, and function prediction analysis.

[0010] The above analysis methods can be performed sequentially or non-sequentially, and the selection of the sequence is mainly determined by the data source; However, in a specific embodiment, species composition analysis, Alpha diversity analysis, Beta diversity analysis, and difference analysis are performed sequentially; In a specific embodiment, the system evolution analysis, the association statistical analysis, and the function prediction analysis can be performed one by one, simultaneously, or sequentially.

[0011] In a specific embodiment of the present invention, during sequencing data import, the sequencing data is paired-end, and the manifest file can be named manifest.txt, which can include but is not limited to three columns of information, such as sample identity information (ID), the absolute path (R1, R2) of the sample sequencing data, etc. Generally, sequencing data can be obtained using a sequencer, including commercial sequencers known in the art for high-throughput sequencing, including first-generation sequencers and second-generation sequencers, such as the Illumina Miseq sequencer, the Illumina PE250 sequencer, and the Illumina NovaSeq X sequencer.

[0012] In a specific embodiment of the present invention, in the denoising step, denoising can include quality filtering, removal of adapter sequences, double-end merging, and removal of chimeras, among other denoising methods, to generate a feature table and representative sequences; DADA2 aims to obtain biological sequences free of amplification and sequencing errors and chimeras through denoising; In a specific embodiment of the present invention, in the step of constructing the evolutionary tree, the phylogeny plug-in of the QIIME2 software is used to perform a multiple sequence alignment of the representative sequences using the MAFFT tool, and then the aligned highly variable regions are filtered, and then the FastTree tool is called to generate a phylogenetic tree based on the filtered alignment results. Finally, the root midpoint method is applied to place the root of the tree at the midpoint of the longest end-to-end distance in the unrooted tree, thereby forming a rooted tree. The file is saved in Newick format.

[0013] In a specific embodiment, the visual mapping data analysis can be performed as follows.

[0014] 1 Species composition analysis

[0015] The species composition analysis is used to study the structure of microbial communities and their dynamic changes in specific environments. It can provide basic information such as how many species there are in the community, how many individuals there are in each species, which species are dominant species, and which are rare species. It plays an important role in optimizing community structure, regulating community function, and exploring new microbial functional groups.

[0016] Species composition analysis can be performed using at least one of the following methods: heat map, Circos species relationship diagram, species common and unique analysis, and community composition histogram; Based on the obtained ASV abundance table, the abundance of each ASV in all samples was calculated, and some ASVs were selected for cluster analysis and heat map drawing, and / or sample-species abundance association analysis was performed and Circos string diagram was drawn; preferably, the top 20 ASVs with the highest abundance were selected for cluster analysis and heat map drawing, and preferably, the top 10 ASVs with the highest abundance were selected for sample-species abundance association analysis and Circos string diagram drawing.

[0017] According to the obtained ASV abundance table, the number of common and unique ASVs in each sample group was calculated, and the number of common and unique ASVs in each sample group was visually presented through the UpSet graph; According to the obtained ASV abundance table, the abundance of each ASV in all samples was calculated, and then some of the ASVs with the highest abundance were selected and the proportion of each ASV in each sample was calculated. The samples were clustered according to the number of these ASVs in each sample and cluster stacked histograms, species (facet) stacking or impact histograms, and grouped species stacking or impact histograms were drawn; preferably, the top 20 ASVs with the highest abundance were selected for processing.

[0018] 2 Alpha diversity analysis

[0019] Alpha diversity refers to the degree of species diversity within a habitat, primarily referring to the differences in the number and frequency distribution of species within a habitat or microbial community, including species richness and species evenness. Alpha diversity analysis can help understand the number and relative abundance of microbial species within each sample, revealing the health status or degree of disturbance of the ecosystem. Alpha diversity analysis can be performed using at least one of the following: alpha diversity index and box plot, dilution curve, abundance rank curve, or species accumulation curve. Calculate six diversity indices for each sample in the ASV abundance table and draw a box plot of each index, including but not limited to Chao1, ACE, Pielou, goods_coverage, Shannon, and Simpson diversity indices, preferably Chao1, ACE, Shannon, and Simpson diversity indices; The total number of sequences in each sample of the ASV abundance table was randomly sampled at different sequencing depths, the number of ASVs corresponding to the total number of sequences in each random sampling was calculated, and a dilution curve was drawn; According to the ASV abundance table, by randomly sampling and gradually increasing the sampling amount, the number of species corresponding to each sample was calculated and the species accumulation curve was drawn; Each ASV in each sample in the ASV abundance table is arranged in descending order of abundance, and the abundance level curve of each sample is drawn with the ASV level as the horizontal axis and the number of sequences contained in each ASV after log10 logarithmic transformation.

[0020] 3 Beta diversity analysis

[0021] Beta diversity, also known as interhabitat diversity, measures the dissimilarity of species composition between communities along an environmental gradient, or the rate of species turnover along that gradient. Beta diversity analysis is used to compare microbial community structure across samples, identify factors influencing community composition (such as environmental conditions, geographic location, and host characteristics), and explore similarities and specificities between communities.

[0022] Furthermore, Beta diversity analysis can be performed by at least one of UPGMA cluster analysis, PCA analysis, PCoA analysis, NMDS analysis, and PLS-DA analysis; Based on the obtained ASV abundance table and phylogenetic tree of ASV representative sequences, the distance matrix was calculated based on Bray-Curtis, weighed unifrac or unweighed unifrac distance, and the distance matrix was subjected to UPGMA cluster analysis respectively; Based on the obtained ASV abundance table, PCA analysis was performed on the community composition structure of the samples, and the natural distribution characteristics between samples were displayed in a two-dimensional image; Based on the obtained ASV abundance table, principal coordinate analysis (PCoA) was performed based on the Bray Curtis distance. The community distribution characteristics between samples were displayed in a two-dimensional image. The Bray Curtis, unweighed unifrac, or weighed unifrac distance matrices were subjected to Adonis tests, and 999 Monte-Carlo permutation tests were used to determine whether the differences between groups were statistically significant. Based on the obtained ASV abundance table, a distance matrix was calculated based on the Bray Curtis distance and non-metric multidimensional scaling (NMDS) analysis was performed to display the community distribution characteristics between samples in a two-dimensional image. The Bray Curtis distance matrix was tested with the anosim test, and 999 Monte Carlo permutation tests were performed to determine whether the differences between groups were statistically significant. Based on the obtained ASV abundance table, a predictive model is constructed, leveraging the sample's microbial community composition to predict the structure of the sample data. A PLS-DA partial least squares analysis is then performed, and a PLS-DA constrained ordination plot is generated based on the first two dimensions of the analyzed data. The PLS-DA analysis process typically incorporates internal validation mechanisms to evaluate model performance. These internal validation mechanisms provide a certain degree of assessment of the model's reliability and predictive power during the model construction and optimization process.

[0023] 4. Difference Analysis

[0024] Differential analysis involves comparing the composition or abundance of microbial communities between different sample groups (e.g., control and experimental groups) to identify statistically significant differences in microbial community characteristics. This analysis helps identify microbial community differences between sample groups and find microbial community components associated with specific conditions or factors.

[0025] Furthermore, the differential analysis included LEfSe analysis and random forest analysis; Based on the obtained ASV abundance table and taxonomic information, a LEFSE analysis was performed to screen key biomarkers. This analysis first used the non-parametric Kruskal rank sum test to detect species with significantly different abundances between groups. Then, the Wilcoxon rank sum test was used to test the consistency of the differences between the species detected in the previous step. Finally, linear regression (LDA) was used to estimate the contribution of each species' abundance to the differences and identify the species most likely to explain the differences between groups. Based on the obtained ASV abundance table, the random forest algorithm was used to pick out ASVs whose abundance distribution showed significant differences among different groups.

[0026] 5. Systematic evolution analysis

[0027] Phylogenetic analysis is used to reveal the evolutionary relationships and evolutionary history of microbial communities. This analysis can help researchers understand the genetic relationships between different bacterial groups in a microbial community and their evolutionary relevance.

[0028] Phylogenetic analysis includes combined analysis of evolutionary trees and heat maps; for example, software such as pheatmap and ggplot2 can be used to draw heat maps to display the abundance distribution of ASVs.

[0029] Based on the obtained ASV abundance table, ASV taxonomic information and phylogenetic tree, representative ASV sequences of interest were selected for phylogenetic analysis (one ASV with the highest abundance was selected as the representative ASV for each genus, and then the top 50 genera with the highest abundance were selected) to draw an evolutionary relationship tree diagram. At the same time, the absolute abundance of ASV in each group was combined to perform a heat map visualization.

[0030] For example, genus-level classification information for each ASV can be extracted from ASV taxonomic information to ensure that each ASV is correctly annotated to the genus level. All ASVs are grouped by genus to generate an ASV list for each genus. Within each genus, the total abundance of each ASV in all samples is calculated, and the ASV with the highest total abundance is selected as the representative ASV for that genus. If a genus has only one ASV, that ASV is directly selected. For each genus, the total abundance of all its ASVs in all samples is calculated. All genera are sorted from high to low by total abundance, and the top 50 genera with the highest total abundance are selected. From these top 50 genera, a representative ASV for each genus is obtained to ensure that these representative ASVs have high abundance and biological significance in the samples.

[0031] In order to ensure the reliability of the selected representative ASV, the following verification can be further performed: Abundance distribution check: Check the abundance distribution of each representative ASV in different samples. Ensure that the representative ASV has high abundance in multiple samples.

[0032] Taxonomic annotation consistency: Check whether the taxonomic annotation of each representative ASV is consistent. Ensure that the taxonomic annotation of the representative ASV is accurate.

[0033] Sequence quality check: Check the quality of each sequence representing the ASV to ensure that the sequence representing the ASV has no obvious errors or contamination.

[0034] Through evolutionary trees and heat maps, we can analyze kinship, correlation analysis between abundance and evolutionary relationships, and infer evolutionary history.

[0035] MAFFT and MUSCLE were used for multiple sequence alignment; RAxML, IQ-TREE, and FastTree were used for phylogenetic tree construction; FigTree and iTOL were used for phylogenetic tree visualization, and R language (ggtree, pheatmap) and Python (ETEToolkit, seaborn) were preferred for composite graph generation.

[0036] Furthermore, the above-mentioned phylogenetic analysis method screens representative ASVs based on a fixed abundance threshold (such as selecting one ASV with the highest abundance for each genus, and then selecting the top 50 genera with the highest abundance). In order to further screen ASVs with low abundance and high functional importance or ecological significance, the inventors further proposed a method for ASV representative sequence selection and evolutionary-functional joint analysis based on dynamic weight screening: not only considering the abundance of ASVs, but also combining their functional potential, ecological significance and evolutionary uniqueness, dynamically calculating the comprehensive weight of each ASV, thereby screening out more representative ASVs.

[0037] 1. Dynamic weight screening of ASV representative sequences (1) Abundance weight: Calculate the total abundance of each ASV in all samples; standardize the abundance value as the abundance weight; (2) Functional weight: predict the functional potential of each ASV based on functional prediction tools (such as PICRUSt2 and FAPROTAX); assign functional weights according to the importance of the function (such as participation in key metabolic pathways); (3) Ecological weight: Based on the distribution of ASVs in different groups (e.g., only appearing in the experimental group), ecological weights are assigned; for example, ASVs that are unique to a specific environment or condition may have a higher ecological weight; (4) Evolutionary weight: The comprehensive weight of each ASV is calculated by integrating abundance, function, ecology and evolutionary weight: Weight Total = α ⋅ Weight Abundance + β ⋅ Weight Function + γ ⋅ Weight Ecology + δ ⋅ Weight Evolution (in, α 、 β 、 γ 、 δis the weight coefficient, which can be adjusted according to the research objectives.) In each genus, the ASV with the highest comprehensive weight was selected as the representative ASV, and among all genera, the top 50 genera with the highest comprehensive weight were selected for subsequent analysis.

[0038] 2. Evolutionary-functional joint analysis (1) Constructing a phylogenetic tree: constructing a phylogenetic tree based on representative ASV sequences screened by dynamic weights; (2) Functional annotation: Functional annotation of representative ASVs was performed to predict the metabolic pathways or ecological functions in which they were involved; (3) Evolution-function association analysis: annotating the functional annotation information of each ASV on the evolutionary tree; Analyze whether ASVs with similar functions are clustered together in the evolutionary tree to reveal the relationship between function and evolution; (4) Visualization: Combine the evolutionary tree with the functional annotation information to draw an evolution-function joint diagram. Use colors or shapes to mark ASVs with different functions to intuitively show the relationship between function and evolution.

[0039] Therefore, combining ASV abundance distribution and functional annotation, dynamically grouping and displaying heatmaps reveals changes in abundance across functional groups. Dynamic weighted screening not only considers ASV abundance but also incorporates functional, ecological, and evolutionary information to select more representative ASVs. Evolutionary-functional joint analysis combines phylogenetic analysis with functional prediction to reveal the relationship between evolutionary relationships and functional potential within microbial communities. Dynamic grouped heatmaps can combine functional annotation and abundance distribution to dynamically display heatmaps in groups, providing a more flexible analysis method.

[0040] 6. Association Statistical Analysis

[0041] Correlation analysis can investigate the relationship between microbial community structure and other factors. These factors can include environmental factors (such as temperature, humidity, and pH), disease states, and biological characteristics. Correlation analysis can explore the correlation between these factors and microbial community composition and assess whether the relationship is significant.

[0042] Further association statistical analyses included CCA / RDA analysis, correlation heat map, and network analysis; Based on the obtained ASV abundance table and ASV taxonomic information, a DCA analysis of species abundance advancement at a specified taxonomic level was performed. If the maximum value on the first axis of the DCA ordination exceeded 4, a unimodal model (CCA) was selected; if it was less than 3, a linear model (RDA) was selected; if it was between 3 and 4, both CCA and RDA were acceptable. An ANOVA analysis was then performed, and statistical significance was determined using 9999 permutation tests. A sample-by-environmental factor-by-species ordination plot was generated. Alternatively, model parameters (such as the weight of the environmental factor) were optimized through cross-validation. Model fit was assessed using Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC).

[0043] Based on the obtained ASV abundance table, a Spearman correlation analysis was performed with the environmental factor information provided by the user. Finally, 20 species related to environmental factors were selected and displayed as a heat map based on their absolute abundance; Based on the obtained ASV abundance table and ASV taxonomic information, the Spearman rank correlation coefficient between the top 150 taxa in terms of abundance was calculated, and an association network was constructed and visualized for the related dominant bacteria with a value greater than |rho|>0.6 and a P value <0.05.

[0044] 7 Functional Prediction Analysis

[0045] Functional prediction can infer the functional potential or metabolic capacity of microbial communities based on 16S rRNA gene sequence information.

[0046] Furthermore, the functional prediction analysis includes the use of EC library predicted abundance table, KO library predicted abundance table, COG library predicted abundance table, KEGG library predicted abundance table, functional difference boxplot, and functional difference heatmap.

[0047] Based on the obtained ASV abundance table and representative sequences, PICRUST2 software was used to predict the functions of the community samples based on functional categories in the KO / KEGG / EC / METACYC databases. Functional differences among the sample communities were demonstrated through PCA analysis, and heat maps were used to display the functional gene categories and abundance changes encoded in each sample. PICRUST2 maps amplicon sequence variants (ASVs) in 16S rRNA gene sequencing data to a reference phylogenetic tree and uses evolutionary relationships and genomic information from known species to predict the gene content of each ASV. These predictions are then integrated to infer the functional profile of the entire microbial community, thereby enabling the indirect inference of potential microbial community functions from marker gene sequences.

[0048] Furthermore, in order to improve database coverage and prediction accuracy, it can also include gene-level prediction, metabolic-level prediction, and ecological function-level prediction, thereby more comprehensively predicting the functional potential of the microbial community.

[0049] The gene-level prediction uses PICRUSt2 to predict KO (KEGG Orthology) and EC (Enzyme Commission) gene functions, supplemented by Tax4Fun2 or FAPROTAX to predict specific ecological functions (such as nitrogen cycle and sulfur cycle).

[0050] The metabolic level prediction is based on the KO prediction results, mapped to the KEGG metabolic pathway database to predict the metabolic potential of the microbial community. The METACYC database is used to supplement the prediction of specific metabolic pathways (such as antibiotic synthesis and secondary metabolite synthesis).

[0051] The ecological function level prediction is based on the FAPROTAX database, which predicts the ecological functions of microbial communities (such as pathogens and nitrogen-fixing bacteria). Combined with environmental factor data, the relationship between function and environment is analyzed.

[0052] Furthermore, dynamic weights are assigned to each database based on its coverage and prediction accuracy. (For example, the KEGG database, with its high coverage, receives a higher weight; the FAPROTAX database, with its focus on ecological functions, receives a medium weight.) For each functional category, the predicted abundance across the different databases is calculated. The predictions are then integrated based on the weights to generate a comprehensive functional abundance table. The integrated results are then evaluated for confidence, with high- and low-confidence function predictions labeled.

[0053] Finally, interactive visualization tools are used to dynamically display functional prediction results and differential analysis.

[0054] Beneficial effects

[0055] The present invention provides a comprehensive, web-based automated 16S analysis method that combines the most advanced analysis projects, from upstream sequencing data processing to downstream data analysis. Simultaneously, the analysis content is modularized, and the data mining methods are diverse and in-depth, allowing users to select different analysis contents based on their needs. This method not only helps researchers without any bioinformatics background conduct 16S studies of microbial microecology, but also overcomes some shortcomings of existing 16S analysis processes (such as the lack of a web page for pre-processing sequence analysis, the lack of real-time updates to the analysis software and database, and the lack of real-time updates to the analysis content), resulting in accurate and comprehensive analysis results that meet a variety of analytical needs. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1Overview of automated methods for 16S rDNA amplicon analysis.

[0057] Figure 2 CCA diagram of samples and environmental factors.

[0058] Figure 3 Predict abundance tables for functions. DETAILED DESCRIPTION

[0059] In a specific embodiment, the analysis content included in the method is as follows Figure 1 As shown: 1. Sequencing data import The FastDNA kit was used to extract DNA from samples such as human intestinal contents, water, and soil. The 16S rRNA gene of the samples was sequenced using the IIlumina Miseq sequencer, and raw paired-end sequencing data was generated. The raw data was used as input data and imported into the QIIME2 software using the QIIME2 software to generate a QIIME2 object.

[0060] 2. Denoising For the QIIME2 object generated in step 1, we used the DADA2 plugin built into QIIME2 to perform quality filtering, primer removal, double-end merging, and chimera removal to generate a signature table (ASV table) and representative sequence files in one step. The generated signature table and representative sequence files are both qiime2 objects, which we then converted to text using QIIME2 and the BIOM tool.

[0061] 3. Species annotation Using the feature-classifier plug-in of the QIIME2 software, the representative sequences generated in the previous step were compared with the reference sequences in the database using the default parameters to obtain the taxonomic information corresponding to each representative sequence; at the same time, in the QIIME2 software, the taxonomic information obtained in the previous step and the feature table generated in step 2 were used to generate the phylum (Phylum) and genus (Genus) level feature tables.

[0062] 4. Constructing an evolutionary tree The phylogeny plug-in of QIIME2 software was used to perform multiple sequence alignment of representative sequences using the MAFFT tool. The highly variable regions of the alignment were then filtered, and the FastTree tool was called to generate a phylogenetic tree based on the filtered alignment results.

[0063] 5. Species composition analysis: Using the phylum or genus level characteristic table and grouping information generated in step 3 as input files, use the R language Heatmap package to draw a heat map of the top 20 species in abundance; Using the phylum or genus level feature table generated in step 3 as the input file, use the circlize package in R language to draw a Circos chord diagram of the sample-species abundance association for the top 10 species in abundance. Using the feature table and grouping information generated in step 2 as input files, use the R language upset package to analyze the number of ASV-specific common features of each sample (group) and visualize them; Using the phylum or genus level characteristic table and grouping information generated in step 3 as input files, use the R language ggplot2 package to draw cluster stacked histograms and impact histograms for the top 20 species in abundance.

[0064] 6 Alpha Diversity Analysis Using the feature table and grouping information generated in step 2 as input files, use the estimateR package in R to calculate the observed_species, chao, and ace indices, and the diversity package to calculate the Shannon and Simpson indices and draw box plots for each diversity index; Using the feature table generated in step 2 as the input file, use the rarecurve package in R to randomly sample the total number of sequences for each sample in the feature table at different sequencing depths. Then use the ggplot2 package to plot the rarefaction curve based on the number of randomly sampled sequences at each depth and the corresponding number of ASVs. Using the feature table generated in step 2 as the input file, use the R language specaccum package to draw a species accumulation curve for the total number of ASVs corresponding to each sample in the feature table; Using the feature table generated in step 2 as the input file, the abundance value of each ASV in the feature table was log10 transformed using the R language ggplot2 package to draw the abundance level curve of each sample.

[0065] 7 Beta diversity analysis Using the feature table generated in step 2, the phylogenetic tree generated in step 4, and the grouping information as input files, use the R language vegdist package to calculate the Bray-Curtis distance matrix and / or the UniFrac package to calculate the (un)weighted_unifrac distance matrix. Then, use the hclust package to perform UPGMA cluster analysis based on the distance matrix. Finally, use the ggtree package to draw a phylogenetic tree to visualize the cluster analysis results. Using the feature table and grouping information generated in step 2 as input files, use the R language prcomp package to perform PCA analysis, and use the ggplot2 package to visualize the PCA analysis results; Using the feature table and grouping information generated in step 2 as input files, use the R language vegdist package to calculate the Braycurtis distance and the UniFrac package to calculate the unweighed unifrac and weighed unifrac distances. Then, use the cmdscale package to perform PCoA analysis based on the distance matrix and use the adonis2 package for significance testing. Finally, use the ggplot2 package to visualize the PCoA analysis results. Using the feature table and grouping information generated in step 2 as input files, we used the R language metaMDS package to perform NMDS analysis and the anosim package to perform significance tests. Finally, we used ggplot2 to visualize the NMDS results. Take the feature table and grouping information generated in step 2 as input files, and use Using the R language plsda package PLS-DA analysis was performed and visualized using ggplot2.

[0066] 8. Difference Analysis Take the feature table generated in step 2, the taxonomy information and grouping information generated in step 3 as input files, and use use The R language microeco package was used to perform Lefse analysis to screen key biomarkers, and the ggplot2 package was used to visualize the Lefse analysis results; Using the feature table and grouping information generated in step 2 as input files, the R language randomForest package was used to screen ASVs with significantly different abundances among different groups, and then the ggplot2 package was used to visualize the analysis results.

[0067] 9 Phylogenetic Analysis Using the feature table generated in step 2, the taxonomic information generated in step 3, and the phylogenetic tree generated in step 4 as input files, the ggtree package in R was used to visualize the phylogenetic tree file, and the gheatmap package was used to add a heat map. Finally, phylogenetic analysis was performed in the form of a phylogenetic tree diagram + heat map.

[0068] 10. Association Statistical Analysis Using the feature table, environmental factors, and grouping information generated in step 2 as input files, we used the R language cca and rda packages to perform CCA and RDA analyses, then used the anova.cca package for significance testing. Finally, we used ggplot2 to visualize the analysis results and generate a sample-environmental factor-species ordination plot. Using the phylum or genus level feature table and environmental factors generated in step 3 as input files, use the R language corrplot package to perform correlation heat map analysis, and use the ggplot2 package to visualize the heat map analysis results; Take the feature table generated in step 2, the taxonomy information and grouping information generated in step 3 as input files, and use use The R language ggClusterNet package calculated the Spearman rank correlation coefficient between the top 150 taxa in abundance, and constructed and visualized the association network for the dominant species with |rho|>0.6 and P value <0.05; 11 Functional Prediction The feature table, representative sequences, and grouping information generated in step 2 were used as input files to perform functional prediction using PICRUST2 software. PCA analysis was used to show the functional similarity between samples, and a heat map was used to show the abundance changes of functional genes in different samples. Figure 3 shown.

[0069] All of the above steps can be integrated into a shell script in the Linux operating system via command lines, allowing for a single execution, automating the entire analysis process and improving efficiency. Alternatively, individual functions can be executed individually to complete the analysis as needed.

Claims

1. An automated method for 16S rDNA amplicon analysis, comprising the following steps: data import, generating a manifest file after the user uploads sequencing data, and then importing the data into a qiime2 object using QIIME2 software; After the sequencing data of the samples were converted into qiime2 objects, the DADA2 method was used to denoise them. By comparing the representative sequences with the reference sequences of the corresponding databases, the taxonomic information corresponding to each ASV was obtained to achieve species annotation. Multiple sequence alignment was performed on the representative sequences and the highly variable regions of the alignment were filtered. Then, phylogenetic trees and rooted trees were generated based on the filtered alignment results. The obtained ASV abundance table, ASV taxonomic information, grouping information, phylogenetic evolutionary tree of representative sequences, environmental factors, and phylum or genus level feature table information were used as input files for the data analysis step, which included species composition analysis, alpha diversity analysis, beta diversity analysis, difference analysis, phylogenetic analysis, association statistical analysis, and functional prediction analysis.

2. The method according to claim 1, wherein denoising comprises quality filtering, removal of linker sequences, double-end merging, and removal of chimeras to generate a feature table and a representative sequence.

3. The method according to claim 1, using the phylogeny plug-in of QIIME2 software, uses the MAFFT tool to perform multiple sequence alignment of representative sequences, then filters the aligned highly variable regions, and then calls the FastTree tool to generate a phylogenetic tree based on the filtered alignment results. Finally, the root midpoint method is applied to place the root of the tree at the midpoint of the longest end-to-end distance in the unrooted tree to form a rooted tree.

4. The method according to claim 1, wherein the species composition analysis is performed using at least one of a heat map, a Circos species relationship diagram, a species common-specific analysis, and a community composition histogram; the Alpha diversity analysis is used to analyze the number and relative abundance of microbial species in each sample, using at least one of an Alpha diversity index and a box plot, a dilution curve, an abundance rank curve, and a species accumulation curve; the Beta diversity analysis is used to compare the microbial community structure between different samples and identify factors affecting community composition, and is performed by at least one of UPGMA cluster analysis, PCA analysis, PCoA analysis, NMDS analysis, and PLS-DA analysis; the difference analysis can compare the composition or abundance of microbial communities to determine microbial community characteristics with statistically significant differences, including LEfSe analysis and random forest analysis; the phylogenetic analysis is used to analyze the kinship and evolutionary correlation between different bacterial communities in the microbial community, including a combined analysis of evolutionary trees and heat maps; the association statistical analysis can study the degree of correlation between microbial community structure and other factors, including RDA / CCA analysis, correlation heat maps, and network analysis; the functional prediction analysis is based on 16S rRNA gene sequence information is used to infer the functional potential or metabolic capacity of the microbial community, including the use of ES library predicted abundance table, KO library predicted abundance table, COG library predicted abundance table, KEGG library predicted abundance table, functional difference boxplot, and functional difference heatmap.

5. The method according to claim 4, wherein In the species composition analysis, based on the obtained ASV abundance table, the abundance of each ASV in all samples and the number of common and unique ASVs in each sample group were calculated, and some ASVs were selected for specific mapping analysis.

6. The method according to claim 1, wherein in the Alpha diversity analysis, a diversity index is calculated for each sample in the ASV abundance table, the total number of sequences in each sample in the ASV abundance table is randomly sampled at different sequencing depths, the number of ASVs corresponding to the total number of sequences in each random sampling is calculated, and the ASVs are arranged according to their abundance, so as to perform mapping analysis separately.

7. The method according to claim 1, wherein the phylogenetic analysis further uses ASV representative sequences based on dynamic weight screening, first performing dynamic weight screening of ASV representative sequences, selecting the ASV with the highest comprehensive weight as the representative, and then performing evolutionary-functional joint analysis, specifically including construction of an evolutionary tree, functional annotation, evolutionary-functional association analysis, and visual display.

8. The method according to claim 4, wherein the LEfSe analysis is used to screen key biomarkers, and the differences of the species detected between groups are consistent in different groups. Finally, linear regression (LDA) is used to estimate the contribution of each species abundance to the differences, and determine the species that is most likely to explain the group differences.

9. The method according to claim 1, wherein the association statistical analysis can generate a result graph of the association between samples and environmental factors, and further generate a dominant sample association network.

10. The method according to claim 1, wherein the functional prediction analysis includes realizing gene-level prediction, metabolic-level prediction, and ecological-function-level prediction for predicting the functions of community samples, performing functional prediction of community samples based on functional categories in the KO / KEGG / EC / METACYC database, and displaying the functional differences of sample communities through PCA analysis; assigning dynamic weights to each database, calculating the predicted abundance of each functional category in different databases, integrating the prediction results according to the weights, generating a comprehensive functional abundance table, performing confidence assessment on the integrated results, marking high-confidence and low-confidence functional predictions, and dynamically displaying the functional prediction results and difference analysis through interactive visualization tools.

11. The method according to claim 1, wherein the beta diversity analysis is performed by calculating the distance matrix based on the phylogenetic tree, and performing UPGMA cluster analysis on the weighted and unweighted distance matrices respectively; PCA analysis, PCoA principal coordinate analysis, and NMDS non-metric multidimensional scaling analysis are performed to display the community distribution characteristics between samples in two-dimensional images, and the validity of PCoA and NMDS analysis is tested by Adonis test and anosim test; PLS-DA analysis is used to improve the predictive ability of the model by adjusting the model parameters, and the reliability and predictive ability of the model are evaluated during the model construction and optimization process.

Images