1 Project Overview

1.1 Project Information

Project Background：Metagenomic sequencing of 6 samples was completed, yielding 6Gb of data per sample and a total of 36Gb,The data volume per sample shall not be less than 90% of the target data volume ，and complete standard metagenomic bioinformatics analysis.

1.2 Basic Sample Information

Table 1 Basic Sample Information Table

The first column is the name of the analyzed sample; the second column is the grouping information

2 Analysis Overview

Metagenomics is a scientific method for studying the collective genetic material of microbial communities. It does not rely on the isolation and cultivation of individual microorganisms, but rather directly extracts all DNA from environmental samples for study. The process includes sample preparation, nucleic acid extraction quality control, library construction, sequencing, data quality control, and data analysis. Metagenomics; Provides powerful tools for understanding the complexity of microbial communities, which is of significant importance in fields such as biodiversity conservation, disease prevention and control, environmental management, and biotechnological applications.

3 Technical route

Metagenomic sequencing is a high-throughput sequencing technology that studies the genetic information of all microorganisms in environmental samples. The process includes steps such as sample preparation, nucleic acid extraction quality control, library construction, sequencing, data quality control, and data analysis. First, sample testing is conducted to ensure sample quality, followed by DNA extraction and nucleic acid testing to obtain high-quality DNA samples; Subsequently, library construction is performed to convert DNA samples into libraries suitable for sequencing. Next, sequencing is carried out to obtain a large amount of sequence data; thereafter, data quality control is performed to remove low-quality data and ensure data accuracy. Finally, data analysis is conducted, including data quality control, metagenome assembly, non-redundant gene set construction, community functional annotation, differential gene analysis, species composition, and diversity analysis.

Figure 1 Technical roadmap

4 Analysis process

Metagenomic data analysis includes the following steps: data quality control, metagenome assembly, gene prediction and construction of a non-redundant gene set, gene functional annotation, gene abundance and differential analysis, species composition and diversity analysis, and species differences. Through these steps, researchers can gain insights into the diversity, composition, and functions of microorganisms in the samples, thereby providing important information for research in fields such as ecology, medicine, and agriculture.

1、Processing and statistics of sequencing data, including filtering low-quality sequences, data output, and quality control statistics;

2、Metagenome assembly;

3、Genome component analysis, including gene prediction and non-redundant gene set construction;

4、Functional annotation, including annotations from general databases such as KEGG, Nr, UniProt, GO, COG, KOG, Pfam, RefSeq, Tigerfam, and special databases such as CAZY, ARDB, VFDB, CARD, PHI, CYPED, TCDB;

5、Gene abundance estimation, differential gene statistics, and functional enrichment analysis;

6、Species composition and diversity analysis, including species annotation, species abundance charts, α diversity, β diversity, etc.;

5 Results Presentation

5.1 Conclusion Preview

A total of 7.20 Gb of raw data was obtained from 6 samples using second-generation Illumina sequencing, with an average of 1.20 Gb of data per sample.

5.2 Data quality control

5.2.1 Statistics of Second-Generation Sequencing Data

Due to the possibility of low-quality sequences and adapter sequences in the raw sequencing data, the raw sequencing data was processed using fastp software to remove adapters and low-quality sequences, resulting in clean reads for subsequent analysis.

Table 2 Second-Generation Sequencing Data Summary Table

Note: SampleId refers to the sample identifier; RawReads refers to the number of raw sequences; RawBases refers to the amount of raw data; CleanReads refers to the number of filtered sequences; CleanBases refers to the amount of filtered data; Q20 refers to the ratio of Q20 reads after filtering; Q30 refers to the ratio of Q30 reads after filtering; GC_content(%) is the GC ratio of filtered reads.

5.3 Metagenome assembly

MEGAHIT software was used to perform metagenome assembly on second-generation data, resulting in the final assembly for subsequent analysis.

Table 3 Assembly Results Data Statistics

Note: Sample: sample name; Contig Num: number of contigs after assembly; Total Len: total number of bases in all contigs; Largest Len: number of bases in the longest contig; N50: when contigs are sorted from longest to shortest and the cumulative length reaches half of the total length of all contigs upon adding a specific contig, the length of that contig is defined as N50; GC: percentage calculated as G+C/A+T+G+C.

Result file：01.assemble/assemble.stat.txt

5.4 Gene set construction

5.4.1 Gene prediction

This project utilized Prokka software to call Prodigal for ORF prediction on the assembly results. The prediction results for each sample are shown in the table below:

Table 4 Gene Prediction Results Statistics

Note: Sample: Sample Name; GeneNum: Number of Predicted Genes; geneLen: Total Base Count of All Predicted Genes; minLen: Minimum Length of Predicted Genes; maxLen: Maximum Length of Predicted Genes; meanLen: Average Length of Predicted Genes; N50: Genes are sorted from longest to shortest, and the cumulative length is calculated. When adding a certain gene causes the cumulative length to equal half of the total length of all genes, the length of that gene is defined as N50.

Result file：02.gene/geneStat.txt

5.4.2 Construction of Non-Redundant Gene Set

After gene prediction is completed, CD-HIT is used to remove redundancy from the annotated genes. The statistics of the number of genes after redundancy removal are shown in the table below:

Table 5 Non-redundant Gene Prediction Results Statistics

Note: Geneset_number: Number of predicted gene sets; Total_Len: Total number of bases of all predicted genes; Average: Average number of bases of predicted genes; Max_length: Maximum length of predicted genes; Min_length: Minimum length of predicted genes.

Result file：02.gene/

5.5 Prediction of genomic functional elements

5.5.1 CRISPR sequence prediction

CRISPR prediction of the genome was performed using the minced software, and the CRISPR prediction statistics for each sample are shown in the table below:

Table 6 CRISPR Prediction Results Statistics for Each Sample

Note: sample: Sample name; numbers: Number of predicted CRISPRs; length: Total length of predicted CRISPRs.

Result file：03.genomeElement/CRISPR/

5.5.2 Prediction of prophages

Prophages in the genome were predicted using the PhiSpy software, and the statistical results for each sample are shown in the table below.

Table 7 Statistical results of prophage predictions for each sample

Note: sample: sample name; numbers: number of predicted prophages; length: total length of predicted prophages.

Result file：03.genomeElement/Prophage/

5.5.3 Prediction of genomic islands

Genomic islands present in the genome were predicted using the IslandPath software, and the statistical results for each sample are shown in the table below.

Table 8 Prediction results of genomic islands

Note: sample: sample name; numbers: number of predicted genomic islands; length: total length of predicted genomic islands.

Result file：03.genomeElement/Island/

5.6 Universal Database Annotation

5.6.1 Annotation Proportion Statistics

Using methods based on BLAST, Diamond, or HMMER alignment, the predicted protein sequences of the genes are compared with functional databases such as UniProt, KEGG, GO, RefSeq, COG, and Pfam to obtain gene functional annotation results. The annotation statistics for each database are shown in the table below:

Table 9 Database Annotation Proportion Statistics

Note: In the first column, ‘All’ refers to all genes; ‘Annotation’ refers to genes with at least one annotation; ‘KEGG’ refers to genes annotated to the KEGG database; ‘Pathway’ refers to genes annotated to the KEGG Pathway database; ‘Nr’ refers to genes annotated to the Nr database; ‘Uniprot’ refers to genes annotated to the Uniprot database; ‘GO’ refers to genes annotated to the GO database; ‘COG’ refers to genes annotated to the COG database; ‘KOG’ refers to genes annotated to the KOG database; ‘Pfam’ refers to genes annotated to the Pfam database; ‘Refseq’ refers to genes annotated to the Refseq database; ‘Tigerfam’ refers to genes annotated to the Tigerfam database; the second column shows the number corresponding to each item; the third column shows the proportion of each item.

Result file：04.regularFunction/stat.xls

5.6.2 KEGG annotation

KEGG (Kyoto Encyclopedia of Genes and Genomes) is a major public database for systematically analyzing the metabolic pathways of gene products in cells and the functions of these gene products. Utilizing KEGG allows for further research into the complex behaviors of genes in biology. After performing KEGG annotation on the genes, they are classified based on the KEGG metabolic pathways they participate in, as shown in the figure below.

Figure 2 Classification of KEGG pathway results

Note: The vertical axis represents KEGG level 2 classification; the horizontal axis represents the number of annotated genes, with different colors representing different KEGG level 1 classifications.

Result file：04.regularFunction/kegg/

5.6.3 Gene Ontology classification

Gene Ontology (abbreviated as GO) is an internationally standardized gene function classification system that provides a dynamically updated controlled vocabulary to comprehensively describe the attributes of genes and gene products in organisms. After simplifying the GO (Gene Ontology) annotation information, the GOslim classification is obtained. The functions of genes are summarized from three aspects: cellular components, molecular functions, and biological processes. The top 20 GOslim secondary classifications with the most annotations in each category are selected for plotting, as shown in the figure below.

Figure 3 GO Annotation Result Classification Chart

Note: The horizontal coordinate represents the annotated GO classification; the vertical coordinate represents the number of annotated genes; different colors represent different major categories.

Result file：04.regularFunction/go/

5.6.4 COG Classification

Cluster of Orthologous Groups of proteins (COG) is a database for the classification of gene products based on direct homologs. Each COG protein is assumed to originate from an ancestral protein. The COG database is constructed based on the coding proteins of bacteria and algae with complete genomes and their phylogenetic relationships. By aligning the Unigene and COG databases, we predict the possible functions of genes and perform functional classification statistics, thereby gaining a macro understanding of the gene function distribution characteristics of the species. COG is divided into 26 groups. The genes annotated by COG are classified according to the COG groups, as shown in the figure below.

Figure 4 COG Results Classification Chart

Note: The horizontal coordinate represents the 26 COG categories; the vertical coordinate represents the number of annotated genes.

Result file：04.regularFunction/cog/

5.6.5 Nr Annotation

Nr (Non-Redundant Protein Sequence Database) is a non-redundant protein database that includes all non-redundant protein sequences from GenBank, EMBL, DDBJ, and PDB. For all known or potential coding sequences, the NR records provide the corresponding amino acid sequences (inferred from known or potential reading frames) as well as sequence numbers from specialized protein databases. After performing Nr annotation on the genes, annotations are made based on their homologous sequences in the alignment, and the species to which the aligned genes belong are extracted. The top ten species with the highest number of aligned genes are plotted, as shown in the figure below.

Figure 5 Species distribution map of Nr results

Note: The figure legend corresponds to the species indicated by color.

Result file：04.regularFunction/nr/

5.7 Advanced Database Annotation

5.7.1 CAZy Database Annotation

Using the hmmer software, the protein sequences of non-redundant genes were aligned with the CAZy database to identify potential carbohydrate-active enzymes present in the gene sequences. The distribution ratios of various types of carbohydrate enzymes are shown in the figure below.

Figure 6 Distribution Proportion of Carbohydrate Enzymes

Result file：05.specificFunction/cazy/

5.7.2 ARDB Database Annotation

The Antibiotic Resistance Genes Database (ARDB) is primarily constructed through literature review, database searches, comparative analysis, and organization. In addition to resistance genes, it also includes 12 drug action targets.

Table 10 ARDB Database Annotation (20 rows displayed)

table

Note: SeqID is the gene set gene ID; ProteinID is the annotated gene name; Pident is the annotation consistency; Qcovs is the annotation coverage; ResistanceType is the resistance gene type name; Description is the description of the antibiotic resistance gene; ResistanceProfile is the name of the resistance gene; Require is the common name of the gene annotated to the resistance gene.

Result file：05.specificFunction/ardb/

5.7.3 CARD database annotation

Using the tool software rgi from the CARD database, the protein sequences of non-redundant genes were aligned with the database to identify corresponding sequences in the database and obtain related resistance genes and resistance-related information.

Table 11 CARD database annotation (20 rows displayed)

table

Note: ORF_ID is the gene set gene id; Best_Hit_ARO is the best matching ARO; Best_Identities is the highest alignment rate; ARO is the ARO ID; Drug Class is the type of drug resistance; Resistance Mechanism is the mechanism of resistance; AMR Gene Family is the resistance gene family to which the resistance gene belongs.

Result file：05.specificFunction/card/

5.7.4 VFDB database annotation

The protein sequences of non-redundant genes were aligned with the VFDB database using Blast, and the alignment results are as follows:

Table 12 VFDB Annotation (Displaying 20 rows)

table

Note: SeqID refers to the gene set gene id; Pident refers to the consistency of the annotation; Qcovs refers to the coverage of the annotation; Accession refers to the annotated VFDB id; Annotation refers to the VFDB annotation information.

Result file：05.specificFunction/vfdb/

5.7.5 Pathogen-Host Interaction Analysis

Using Diamond blastp to align the protein sequences of non-redundant genes to the PHI database, the results are shown in the table below:

Table 13 PHI Annotation (Displaying 20 rows)

table

Note: SeqID refers to the gene set gene id; Pident refers to the consistency of the annotation; Qcovs refers to the coverage of the annotation; ProteinID refers to the annotated PHI sequence ID; PhiID refers to the annotated PHI id; PathogenGene refers to the pathogen gene; PathogenSpecies refers to the species of the pathogen; MutantPhenotype refers to pathogenicity.

Result file：05.specificFunction/phi/

5.7.6 Cytochrome P450 Annotation

Aligning the protein sequences of non-redundant genes to the CYPED database, the statistical results are as follows:

Table 14 CYPED Annotation (showing 20 rows)

table

Note: SeqID is the gene set gene id; Pident is the consistency of the annotation; Qcovs is the coverage of the annotation; SequenceID is the p450 number of the annotation; ProteinID is the p450 protein number of the annotation; HomologousFamilies is the homologous family number of the annotation; Superfamilies is the sub-number of the homologous family of the annotation; GI is the GI number of the annotated gene on NCBI; TaxonID is the NCBI species taxonomy number; Species is the Latin name of the annotated species.

Result file：05.specificFunction/cyped/

5.7.7 TCDB Transport Protein Classification Annotation

Aligning the protein sequences of non-redundant genes to the TCDB database, the statistical results are as follows:

Figure 7 TCDB classification bar chart

Note: The horizontal coordinate represents the type of transport protein; the vertical coordinate represents the number of genes annotated to this type; different colors represent different types.

Result file：05.specificFunction/tcdb/

5.8 Gene abundance analysis

Quantification of the non-redundant gene set was performed using the salmon software, with gene abundance standardized using the TPM (Transcripts Per Million) algorithm.

5.8.1 Overall situation of gene abundance

A density plot based on standardized gene abundance is drawn:

Figure 8 Gene abundance density plot

Note: The horizontal coordinate is log10(TPM); the vertical coordinate is density; different colors represent different samples.

Result file：06.geneExpression/sample.density.png

Based on the standardized gene abundance, the correlation between samples is calculated, and a heatmap is drawn to display the similarity and clustering relationships among the specimens.

Figure 9 Sample Correlation Clustering Heatmap

Result file：06.geneExpression/corHeatMap/

Principal Component Analysis (PCA) is a technique for analyzing and simplifying datasets by decomposing variance, reflecting the differences of multiple datasets in a two-dimensional coordinate graph. By analyzing the gene abundance of different samples, the differences and distances between the samples can be reflected; the closer the distance between two samples in the graph, the more similar their compositions are.

Figure 10 Gene Abundance PCA Graph

Note: The points represent each sample; different colors indicate different groups; the horizontal coordinate represents the first principal component, and the percentage indicates the contribution value of the first principal component to the sample differences; the vertical coordinate represents the second principal component, and the percentage indicates the contribution value of the second principal component to the sample differences; different colors represent different groups.

Result file：06.geneExpression/PCA

5.8.2 Differential gene analysis

When there are multiple groups and each group has at least one replicate, DESeq2 is used for differential analysis between pairs of samples; When there is a single group with multiple samples or multiple groups without replicates, edgeR is used for differential analysis between pairs of samples. The statistical table of the number of differential genes for each differential combination is as follows:

Table 15 Differential gene count statistical table

The first column is the name of the differential combination; the second column is the number of significantly differential genes; the third column is the number of upregulated genes; the fourth column is the number of downregulated genes.

The volcano plot of differential genes for group.A_vs_B is as follows: (Other differential combination results can be found in the result file)

Figure 11 Volcano plot of differential genes

Note: Each point in the figure represents a gene; the horizontal coordinate is the logarithmic value (base 2) of the abundance difference of a gene between two samples; the vertical coordinate represents the logarithmic value (base 10) of the p-value (padj) after multiple testing correction, multiplied by -1.

Result file：06.geneExpression/DEG/

5.8.3 Differential gene enrichment analysis

Based on the R package clusterProfiler, KEGG and GO functional enrichment analyses were conducted on differential genes to investigate whether the functional differences in abundance are concentrated in certain specific functional classifications.

This section only presents the KEGG and GO enrichment analysis results of differential genes for group.A_vs_B; results for other differential combinations can be found in the result file.

Figure 12 Bubble chart of KEGG enrichment for differential genes

Note: In the above figure, the horizontal axis represents GeneRatio (the proportion of differential genes in this Pathway classification to all differential genes), the vertical axis represents KEGG Pathway classification, the color of the points represents q-value, and the size of the points represents the number of differential genes.

Figure 13 Bubble chart of GO enrichment for differential genes

Note: In the above figure, the horizontal axis represents GeneRatio (the proportion of differential genes in this GO classification to all differential genes), the vertical axis represents Gene Ontology functional classification, the color of the points represents q-value, and the size of the points represents the number of differential genes.

Result file：06.geneExpression/DEG/

5.9 Species composition and diversity analysis

5.9.1 Species composition analysis

This project uses Kraken2 software for read-based metagenomic species annotation. This method is more comprehensive and accurate, eliminating the resource consumption of assembly and gene prediction processes, allowing for the annotation of more low-abundance species, and significantly improving the accuracy of species annotation and relative abundance.

Table 16 Statistical table of species classification levels

The first column is the sample ID, and the second to eighth columns represent the number of annotations at the kingdom, phylum, class, order, family, genus, and species levels, respectively.

5.9.1.1 Analysis of endemic species

Using Venn diagrams (suitable for 2-5 samples or groups) or petal diagrams (for sample or group numbers greater than or equal to 2) can display the number of shared and endemic species between different samples or groups, visually representing the overlap of species among different samples or groups, and identifying shared endemic microorganisms in different environments. The results are shown in the figure below.

5.9.1.1.1 Venn diagram

Figure 14 Venn diagram

Note: The numbers in the overlapping areas between multiple colors in the figure represent the number of shared species among multiple samples, while the non-overlapping areas represent the number of endemic species for each sample.

5.9.1.1.2 Petal diagrams

Figure 15 petal diagrams

Note: The number in the center circle of the figure indicates the total number of species shared by all samples, the numbers on the petals pointing outward indicate the number of species unique to each sample, and the numbers on the petals pointing inward indicate the number of species shared with other samples.

Result file：07.taxonomy/taxonomy/venn/

5.9.1.2 Species abundance at various levels

To optimize the visual effect, only the top ten species by abundance are displayed, and other species are combined and shown as ‘Others’ in the figure. ‘Unclassified’ represents species that have not been assigned a taxonomic annotation.

5.9.1.2.1 phylum

Figure 16 Bar chart of species composition at the phylum level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

5.9.1.2.2 class

Figure 17 Bar chart of species composition at the class level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

5.9.1.2.3 order

Figure 18 Bar chart of species composition at the order level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

5.9.1.2.4 family

Figure 19 Bar chart of species composition at the family level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

5.9.1.2.5 genus

Figure 20 Bar chart of species composition at the genus level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

5.9.1.2.6 species

Figure 21 Bar chart of species composition at the species level

Note: The horizontal coordinate represents samples, the vertical coordinate represents relative abundance, and different colors represent different species.

Result file：07.taxonomy/taxonomy/barplot/

5.9.1.3 Species abundance clustering heatmap

A heatmap is a graphical representation that uses a color gradient to represent the magnitude of values in a data matrix and clusters based on the similarity of species or sample abundance. High and low abundance species are clustered and aggregated, reflecting the similarity and differences in the community composition of multiple samples through color gradients and degrees of similarity. Based on the relative abundance of species in each sample, the top 50 species (excluding unclassified ones) are selected at the phylum, class, order, family, genus, and species taxonomic levels for heatmap generation.