1. Introduction

Metagenomics is the study of genetic material recovered directly from environmental samples. The broad field may also be referred to as environmental genomics, ecogenomics or community genomics. While traditional microbiology and microbial genome sequencing and genomics rely upon cultivated clonal cultures, early environmental gene sequencing cloned specific genes (often the 16S rRNA gene) to produce a profile of diversity in a natural sample. Such work revealed that the vast majority of microbial biodiversity had been missed by cultivation-based methods. Recent studies use either "shotgun" or PCR directed sequencing to get largely unbiased samples of all genes from all the members of the sampled communities. Because of its ability to reveal the previously hidden diversity of microscopic life, metagenomics offers a powerful lens for viewing the microbial world that has the potential to revolutionize understanding of the entire living world. As the price of DNA sequencing continues to fall, metagenomics now allows microbial ecology to be investigated at a much greater scale and detail than before.

The first metagenomic studies conducted using high-throughput sequencing used massively parallel 454 pyrosequencing. Three other technologies commonly applied to environmental sampling are the Ion Torrent Personal Genome Machine, the Illumina MiSeq or HiSeq and the Applied Biosystems SOLiD system. These techniques for sequencing DNA generate shorter fragments than Sanger sequencing; Ion Torrent PGM System and 454 pyrosequencing typically produces ~400 bp reads, Illumina MiSeq produces 400-700bp reads (depending on whether paired end options are used), and SOLiD produce 25-75 bp reads. Historically, these read lengths were significantly shorter than the typical Sanger sequencing read length of ~750 bp, however the Illumina technology is quickly coming close to this benchmark. However, this limitation is compensated for by the much larger number of sequence reads. In 2009, pyrosequenced metagenomes generate 200–500 megabases, and Illumina platforms generate around 20–50 gigabases, but these outputs have increased by orders of magnitude in recent years. An additional advantage to high throughput sequencing is that this technique does not require cloning the DNA before sequencing, removing one of the main biases and bottlenecks in environmental sampling.

2. Materials and Methods

2.1 Library construction and sequencing

DNA extractions:

DNA from different samples was extracted using The E.Z.N.A.® Stool DNA Kit (D4015-02, Omega, Inc., USA) according to manufacturer’s instructions. The reagent which was designed to uncover DNA from trace amounts of sample has been shown to be effective for the preparation of DNA of most bacteria. Sample blanks consisted of unused swabs processed through DNA extraction and tested to contain no DNA amplicons. The total DNA was eluted in 50 µl of Elution buffer by a modification of the procedure de-scribed by manufacturer (QIAGEN) and stored at -80°C until measurement in the PCR by LC Sciences.

DNA Library Construction and sequencing:

DNA library was constructed by TruSeq Nano DNA LT Library Preparation Kit (FC-121-4001) . DNA was fragmented by dsDNA Fragmentase(NEB, M0348S) by incubate at 37°C for 30min. Library construction begins with fragmented cDNA. Blunt-end DNA fragments are generated using a combination of fill-in reactions and exonuclease activity, and size selection is performed with provided sample purification beads. An A-base is then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters.Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. These adapters contain the full complement of sequencing primer hybridization sites for single, paired-end, and indexed reads. Single- or dual-index adapters are ligated to the fragments and the ligated products are amplified with PCR by the following conditions: initial denaturation at 95℃ for 3 min; 8cycles of denaturation at 98℃ for 15 sec, annealing at 60℃ for 15 sec, and extension at 72℃ for 30 sec; and then final extension at 72℃ for 5 min. And then we performed the paired-end sequencing on an Illumina Hiseq 4000 at LC Sciences following the vendor's recommended protocol.

2.2 Bioinformatics analysis

Raw sequencing reads were processed to obtain valid reads for further analysis. First, sequencing adapters were removed from sequencing reads using cutadapt v1.9 [1]. Secondly, low quality reads were trimmed by fqtrim v0.94 using a sliding-window algorithm. Thirdly, reads were aligned to the host genome using bowtie2 [2] to remove host contamination. Once quality-filtered reads were obtained, they were de novo assembled to construct the metagenome for each sample by SPAdes v3.10.0 [3]. All coding regions (CDS) of metagenomic contigs were predicted by MetaGeneMark v3.26.

CDS sequences of all samples were clustered by CD-HIT v4.6.1 [4] to obtain unigenes. Unigene abundance for a certain sample were estimated by TPM based on the number of aligned reads by bowtie2 v2.2.0 [2]. The lowest common ancestor taxonomy of unigenes were obtained by aligning them against the NCBI NR database by DIAMOND v 0.7.12 [5]. Similarly, the functional annotation(GO, KEGG, eggnog, CAZy, CARD, PHI) of unigenes were obtained. Based on the taxonomic and functional annotation of unigenes, along with the abundance profile of unigenes, the differential analysis were carried out at each taxonomic or functional or gene-wise level by Fisher's exact test(non-replicated groups) or Kruskal-Wallis test(replicated groups).

Bioinformatics pipeline for megagenomics

Reference:

1. Martin M (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. Embnet Journal 17.

2. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9: 357-359.

3. Bankevich, A., et al. (2012). SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing. Journal of Computational Biology, 19(5), 455.

4. Li, W., & Godzik, A. (2006). CD-HIT: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics, 22(13), 1658.

5. Buchfink, B., Xie, C., & Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nature Methods, 12(1), 59.

3. Project information

3.1 Sample information

Sample type: human intestinal

Objects: bacteria/archaea/fungus/virus

3.2 Database

3.3 Bioinformatics software

4. Overview of results

4.1 Sample information

document location: summary/1_rawdata/sample_map.xlsx

4.2 Quality control

document location: summary/2_CleanData/DEMO_data_stat.xlsx

4.3 Assembly

4.3.1 Results of assembly

We utilize metaSPAdes (http://bioinf.spbau.ru/en/spades3.7) to assemble reads of each sample into contigs first. Then we merged all samples into a "MIX" sample and continued to assemble "MIX" sample. The final result with FASTA format was deposited in summary/3_Assembly. Only those sequences longer than 500nt would be remained.

Example:

NODE 1 length 299711 cov 17.4623 ID 9988

AACCTAAAACTCCTGTAAGAAGAAGCGCCCCTAGAAATTGAGTCATCTTATCACCTACTTTTTATTTTATTATAACACATTTAGTACCTAGTACTAAATT

ACGGGTAGCCCGACTACCCTTATTATTTTTTAATATTTTATAGAACATACGTTCTTGCAGGAGGTATAAACATGTGGATTGAAAAATTTAAAAACAAAAAT

AACGAAACTAAATACAGATACTACGAGAAGTATAAAGAT...

Explaination: NODE 1=serial number; length 299711=length of contigs; cov 17.4623=coverage of maximum kmer

document location: summary/3_Assembly/*/contigs.min500.fasta

4.3.2 Statistics of assembled results

document location:

summary/3_Assembly/assembly_stats.html

summary/3_Assembly/assembly_stats.xlsx

4.4 Gene prediction and abundance estimation

4.4.1 Abundance estimation

Coding regions (CDS) of assembled contigs were predicted by MetaGeneMark and those contigs shorter than 100nt were removed. CD-HIT was utilized to produce a set of non-redundant and to output a cluster file (identity=95%, coverage=90%). Counts of each unigene was calculated by bowtie2 through mapping reads of each sample on unigene. Abundance estimation was calculated by TPM as shown below.

counts of unigene:

document location: summary/4_GenePredict/4_Unigenes_count.txt

TPMs of unigene:

document location: summary/4_GenePredict/4_Unigenes_abund.txt

4.4.2 Boxplot of unigene abundance

document location: summary/4_GenePredict/4_Unigenes_abund.boxplot.png

4.4.3 Density of unigene abundance

document location: summary/4_GenePredict/4_Unigenes_abund.density.png

4.4.4 Statistics of unigene CDS

document location:

summary/4_GenePredict/2_Unigenes_len.png

summary/4_GenePredict/2_Unigenes_CDS_stats.xlsx

4.4.5 Dilution curve of core/pan genes

core genes = intersection of genes from different samples

pan genes = union of genes from different samples

Dilution curve of core genes:

Dilution curve of pan genes:

document location:

summary/4_GenePredict/5_core_gene.png

summary/4_GenePredict/5_pan_gene.png

4.4.6 Venn diagram of different groups

We calculated unigenes owned by all samples and shared particularlly through venn diagram.

1111_venn_diagram_Group_stage1_vs_stage2_vs_control.pngvenn_diagram_Group_stage1_vs_stage2_vs_control.png

document location:

summary/5_TaxonomicProfiling/*/taxonomy_*_stacked_bar.png

summary/5_TaxonomicProfiling/*/taxonomy_*_stacked_bar.xlsx

5.3.3 Heatmap

We selected most abundanct classifications (TOP20) and have drawn heatmaps from SuperKingdom to species in each sample.

taxonomy_SuperKingdom_heatmap.pngtaxonomy_Phylum_heatmap.pngtaxonomy_Class_heatmap.pngtaxonomy_Order_heatmap.pngtaxonomy_Family_heatmap.pngtaxonomy_Genus_heatmap.pngtaxonomy_Species_heatmap.png

document location:

summary/5_TaxonomicProfiling/*/taxonomy_*_heatmap.png

summary/5_TaxonomicProfiling/*/taxonomy_*_heatmap.xlsx

5.3.4 Principal Component Analysis

1_SuperKingdom_PCA_Group.png2_Phylum_PCA_Group.png3_Class_PCA_Group.png4_Order_PCA_Group.png5_Family_PCA_Group.png6_Genus_PCA_Group.png7_Species_PCA_Group.png

document location: summary/5_TaxonomicProfiling/*/PCA/PCA_Group.png

5.3.5 Cluster analysis

taxonomy_SuperKingdom_cluster.pngtaxonomy_Phylum_cluster.pngtaxonomy_Class_cluster.pngtaxonomy_Order_cluster.pngtaxonomy_Family_cluster.pngtaxonomy_Genus_cluster.pngtaxonomy_Species_cluster.png

document location: summary/5_TaxonomicProfiling/*/taxonomy_*_cluster.png

5.4 Difference analysis

Group/sample differences of abundance were analyzed using Fisher's exact test (non biological repetitions) or Kruskal-Wallis test (biological repetitions and more than two groups). We took Phylum as an example.

document location: summary/5_TaxonomicProfiling/*/diff_abundance/*vs*/*vs*.diff.xlsx

PCA:

Heatmap:

document location:

summary/5_TaxonomicProfiling/*/diff_abundance/*vs*/*vs*.PCA.png

summary/5_TaxonomicProfiling/*/diff_abundance/*vs*/*vs*.heatmap.png

6. Functional annotation

6.1 GO database

6.1.1 GO annotation of unigene

Gene ontology (GO) (http://www.geneontology.org) is a major bioinformatics initiative to unify the representation of gene and gene product attributes across all species. More specifically, the project aims to: 1) maintain and develop its controlled vocabulary of gene and gene product attributes; 2) annotate genes and gene products, and assimilate and disseminate annotation data; and 3) provide tools for easy access to all aspects of the data provided by the project, and to enable functional interpretation of experimental data using the GO, for example via enrichment analysis.

GO is part of a larger classification effort, the Open Biomedical Ontologies (OBO).

Although gene nomenclature itself aims to maintain and develop controlled vocabulary of gene and gene products, the Gene Ontology extends the effort by using markup language to make the data (not only of the genes and their products but also of all their attributes) machine readable, and to do so in a way that is unified across all species (whereas gene nomenclature conventions vary by biologic taxon).

document location: summary/6_FunctionalProfiling/GO/Unigenes_GO.txt

6.1.2 Stacked bar of GO functions/IDs

GO_GO_Function_stacked_bar.pngGO_GO_ID_stacked_bar.png

document location: summary/6_FunctionalProfiling/GO/*/GO_*_stacked_bar.png

6.1.3 Heatmap of GO functions/IDs

GO_GO_Function_heatmap.pngGO_GO_ID_heatmap.png

document location: summary/6_FunctionalProfiling/GO/*/GO_*heatmap.png

6.1.4 PCA analysis

1_GO_Function_PCA_Group.png2_GO_ID_PCA_Group.png

document location: summary/6_FunctionalProfiling/GO/*/PCA/PCA_Group.png

6.1.5 Cluster analysis of GO functions/IDs

GO_GO_Function_cluster.pngGO_GO_ID_cluster.png

document location: summary/6_FunctionalProfiling/GO/*/GO_*_cluster.png

6.2 KEGG database

6.2.1 KEGG annotation of unigene

KEGG (Kyoto Encyclopedia of Genes and Genomes) (http://www.kegg.jp/) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances. KEGG is utilized for bioinformatics research and education, including data analysis in genomics, metagenomics, metabolomics and other omics studies, modeling and simulation in systems biology, and translational research in drug development.

Category of enriched KEGG pathways:

document location:

summary/6_FunctionalProfiling/KEGG/Unigenes_KEGG.txt

summary/6_FunctionalProfiling/KEGG/KEGG_category.png

6.2.2 Stacked bar of KEGG pathways

KEGG_KEGGLevel1_stacked_bar.pngKEGG_KEGGLevel2_stacked_bar.pngKEGG_PathwayEntry_stacked_bar.pngKEGG_KOEntry_stacked_bar.png

document location: summary/6_FunctionalProfiling/KEGG/*/KEGG_*_stacked_bar.png

6.2.3 Heatmap of KEGG pathways

KEGG_KEGGLevel1_heatmap.pngKEGG_KEGGLevel2_heatmap.pngKEGG_PathwayEntry_heatmap.pngKEGG_KOEntry_heatmap.png

document location: summary/6_FunctionalProfiling/KEGG/*/KEGG_*_heatmap.png

6.2.4 PCA analysis

1_KEGGLevel1_PCA_Group.png2_KEGGLevel2_PCA_Group.png3_PathwayEntry_PCA_Group.png4_KOEntry_PCA_Group.png

document location: summary/6_FunctionalProfiling/KEGG/*/PCA/PCA_Group.png

6.2.5 Cluster analysis of KEGG pathways

KEGG_KEGGLevel1_cluster.pngKEGG_KEGGLevel2_cluster.pngKEGG_PathwayEntry_cluster.pngKEGG_KOEntry_cluster.png

document location: summary/6_FunctionalProfiling/KEGG/*/PCA/PCA_Group.png

6.2.6 Difference analysis of KEGG pathways

We took KEGG Level1 as an example for difference analysis of KEGG pathways between different samples.

document location: summary/6_FunctionalProfiling/KEGG/*/diff_abundance/*vs*/*vs*.diff.xlsx

Heatmap and PCA analysis of different KEGG pathways according to their functions. We took PathwayEntry as an example.

PCA:

Heatmap:

document location:

summary/6_FunctionalProfiling/KEGG/3_PathwayEntry/diff_abundance/*vs*/*vs*.PCA.png

summary/6_FunctionalProfiling/KEGG/3_PathwayEntry/diff_abundance/*vs*/*vs*.heatmap.png

6.3 eggNOG database

6.3.1 eggNOG annotation of unigenes

The eggNOG database (http://eggnogdb.embl.de/) is a database of biological information hosted by the ensembl. It is based on the original idea of COGs (clusters of orthologous groups) and expands that idea to non-supervised orthologous groups constructed from numerous organisms. The database was created in 2007 and updated to version 4.5 in 2015. eggNOG stands for evolutionary genealogy of genes: Non-supervised Orthologous Groups.

We aligned unigenes against eggNOG database by DIAMOND (blastp mode, evalue ≤ 1e-5). Hits of highest score were annotated on unigenes with several level of eggNOG.

document location: summary/6_FunctionalProfiling/eggNOG/Unigenes_eggNOG.txt

Number of unigenes on category of eggNOG:

document location: summary/6_FunctionalProfiling/eggNOG/eggNOG_category.png

Taxonomy counts of eggNOG:

document location: summary/6_FunctionalProfiling/eggNOG/eggNOG_taxon_count.xlsx

6.3.2 Stacked bar of eggNOG

This section is similar with 6.2.2 Stacked bar of KEGG pathways

document location: summary/6_FunctionalProfiling/eggNOG/*/eggNOG_*_stacked_bar.png

6.3.3 Heatmap of eggNOG

This section is similar with 6.2.3 Heatmap of KEGG pathways

document location: summary/6_FunctionalProfiling/eggNOG/*/eggNOG_*_heatmap.png

6.3.4 PCA analysis

This section is similar with 6.2.4 PCA analysis

document location: summary/6_FunctionalProfiling/eggNOG/*/PCA/PCA_Group.png

6.3.5 Cluster analysis of eggNOG

This section is similar with 6.2.5 Cluster analysis of KEGG pathways

document location: summary/6_FunctionalProfiling/eggNOG/*/eggNOG_*_cluster.png

6.3.6 Difference analysis of eggNOG

This section is similar with 6.2.6 Difference analysis of KEGG pathways

document location:

summary/6_FunctionalProfiling/eggNOG/*/diff_abundance/*vs*/*vs*.diff.xlsx

summary/6_FunctionalProfiling/eggNOG/*/diff_abundance/*vs*/*vs*.PCA.png

summary/6_FunctionalProfiling/eggNOG/*/diff_abundance/*vs*/*vs*.heatmap.png

6.4 CAZy database

6.4.1 CAZy annotation of unigenes

CAZy (http://www.cazy.org/) is a database of Carbohydrate-Active enZYmes (CAZymes). The database contains a classification and associated information about enzymes involved in the synthesis, metabolism, and transport of carbohydrates. Included in the database are glycoside hydrolases, glycosyltransferases, polysaccharide lyases, carbohydrate esterase and carbohydrate-binding families.

CAZy was established in 1999 in order to provide online and constantly updated access to the family classification of CAZymes. New genomes are added shortly after they appear in the daily releases of GenBank. As of November 2013, CAZy contains sequence information on nearly 340,000 CAZymes.

The CAZy database is coupled with the CAZypedia online encyclopedia, which was launched in June 2009 and intended to be a wiki-based encyclopedia of CAZymes. As of 2014, CAZy is developed by the Glycogenomics group at AFMB, a research centre affiliated with the French National Centre for Scientific Research and Aix-Marseille University.

We aligned unigenes against CAZy database by DIAMOND (blastp mode, evalue ≤ 1e-5). Hits of highest score were annotated on unigenes.

document location: summary/6_FunctionalProfiling/CAZy/Unigenes_CAZy.txt

document location: summary/6_FunctionalProfiling/CAZy/CAZy_category.png

6.4.2 Stacked bar of CAZy

This section is similar with 6.2.2 Stacked bar of KEGG pathways

document location: summary/6_FunctionalProfiling/CAZy/*/CAZy_*_stacked_bar.png

6.4.3 Heatmap of CAZy

This section is similar with 6.2.3 Heatmap of KEGG pathways

document location: summary/6_FunctionalProfiling/CAZy/*/CAZy_*_heatmap.png

6.4.4 PCA analysis

This section is similar with 6.2.4 PCA analysis

document location: summary/6_FunctionalProfiling/CAZy/*/PCA/PCA_Group.png

6.4.5 Cluster analysis of CAZy

This section is similar with 6.2.5 Cluster analysis of KEGG pathways

document location: summary/6_FunctionalProfiling/CAZy/*/CAZy_*_cluster.png

6.4.6 Difference analysis of CAZy

This section is similar with 6.2.6 Difference analysis of KEGG pathways

document location:

summary/6_FunctionalProfiling/CAZy/*/diff_abundance/*vs*/*vs*.diff.xlsx

summary/6_FunctionalProfiling/CAZy/*/diff_abundance/*vs*/*vs*.PCA.png

summary/6_FunctionalProfiling/CAZy/*/diff_abundance/*vs*/*vs*.heatmap.png

6.5 CARD database

6.5.1 CAZy annotation

The Comprehensive Antibiotic Resistance Database ("CARD") (https://card.mcmaster.ca/) provides data, models, and algorithms relating to the molecular basis of antimicrobial resistance. The CARD provides curated reference sequences and SNPs organized via the Antibiotic Resistance Ontology ("ARO"). These data can be browsed on the website or downloaded in a number of formats. These data are additionally associated with detection models, in the form of curated homology cut-offs and SNP maps, for prediction of resistome from molecular sequences. These models can be downloaded or can be used for analysis of genome sequences using the Resistance Gene Identifier ("RGI"), either online or as a stand-alone tool.

The CARD was designed and developed by the laboratories of Drs. Gerry Wright and Andrew G. McArthur of McMaster University's Department of Biochemistry & Biomedical Sciences (Hamilton, Ontario, Canada) with the help of a global team of collaborators. It is built entirely using open source software and tools. This research has been supported by funds from the Canadian Foundation for Innovation, Canadian Institutes of Health Research, Natural Sciences and Engineering Research Council of Canada, Medical Research Council (UK), and Ontario Research Fund, as well as a Cisco Research Chair in Bioinformatics supported by Cisco Systems Canada, Inc. (Dr. McArthur), Canada Research Chair (Dr. Wright), and Killam Research Fellowship (Dr. Wright).

document location: summary/6_FunctionalProfiling/CARD/Unigenes_CARD.txt

6.6 PHI database

6.6.1 PHI annotation

The Pathogen-Host Interaction database (PHI-base) (http://www.phi-base.org/) contains expertly curated molecular and biological information on genes proven to affect the outcome of pathogen-host interactions. The database was created and is maintained by researchers at Rothamsted Research and external collaborators since 2005. Since December 2016 PHI-base is part of ELIXIR, the European life-science infrastructure for biological information via its ELIXIR-UK node.

document location: summary/6_FunctionalProfiling/PHI/Unigenes_PHI.txt

7. Functional analysis of differential genes

7.1 Analysis of differential genes

We took Group.stage1_vs_control as an example.

document location: summary/7_Different_expression/Group.stage1_vs_control/Group.stage1_vs_control_diff_exp.txt

7.2 Volcano of differential unigenes

Group.stage1_vs_control_diff_exp_volcano.pngGroup.stage2_vs_control_diff_exp_volcano.pngGroup.stage2_vs_stage1_diff_exp_volcano.png