MetaDEGalaxy: Galaxy workflow for differential abundance analysis of 16s metagenomic data

Introduction

It is now recognized that there is a strong link between microbial communities in the human body and human health¹. While the importance of such communities is understood, the composition and function of the human micro-biome largely remains a mystery. Uncovering how the composition and function of the micro-biome impacts human health represents a significant area of growth. Another important area of research growth is the study of environmental microbial communities in fields such as agriculture, marine science, and ecology. By identifying the composition of microbial communities, researchers are able to link microbes to specific environments and using comparative metagenomics identify how microbial communities’ changes under altered environmental conditions.

Central to elucidating the link between the metagenomic data and human health or altered environmental conditions is sequencing; however, obtaining useful research outcomes from large volumes of unprocessed sequence data represents a challenge for many bench scientists. The major bottleneck in obtaining value from such data is the huge computational and logistic task required for analysing the large volumes of sequencing data routinely generated in a single sequencing run.

The sequencing of entire microbial communities requires metagenomic analysis tools. These tools rely on the ability to analyse unbroken sequence reads covering the 16S variable regions. Due to limitations of short read sequencing platforms such as IIlumina, the longest fragment of variable regions of a 16S gene that can be sequenced is shorter than the ideal full 600 bp. Illumina paired-end sequencing of 300 bp on forward read and reverse read produces only 550 bp to allow for stitching the forward end and reverse end together. With 550 bp fragment length, the reads can cover both variable region 3 (V3) and variable region 4 (V4). The length of V3 and V4 are 393bp and 440bp respectively.

A major challenge for bench scientists working with metagenomic data is that many popular software programs requires a 64-bit Linux environment, an environment often unavailable and unfamiliar to researchers. Furthermore, even when such an environment is available, the complexity of the rapidly changing metagenomic algorithms means no gold standard methodologies exist. As such, there are currently over 100 metagenomic analysis tools available, making it challenging to select the appropriate software. For example, the popular metagenomic tool QIIME² consists of more than 150 python scripts, many of which are wrappers to external programs.

An increasingly common alternative for the growing number of non-bioinformaticians working with NGS data is the availability of user-friendly interfaces. These interfaces are typically attached to significant compute resources with pre-installed software packages readily available. Interfaces such as Galaxy³ or the Genomics Virtual Lab⁴ are examples of powerful platforms that grant non-bioinformaticians access to the latest NGS methodologies. The Galaxy platform enables scientists to use bioinformatics tools in an easy to use graphical user interface (GUI) environment, where tool resource management is handled by the administrators of each Galaxy service. The platform’s functionality power comes from the ability to chain tools into workflows, and share the data and workflows. Further, the flexibility of Galaxy platform allows developers to integrate new tools and workflows into the platform. Galaxy maintains a single tool shed repository of pre-wrapped tools that cover an abundance of next generation sequence analyses.

Despite these solutions however, challenges remain in fast moving research areas such as metagenomics with limited metagenomic workflows currently available within the popular Galaxy framework. Currently, there is one end-to-end existing metagenomics workflow offering, the recently published ASaiM⁵ for both 16S and shotgun metagenomic analysis. While there is overlap between their workflows, MetaDEGalaxy differs in its focus on differential abundance by incorporating the capabilities of phyloseq⁶ and DESeq2⁷ for complex differential bacterial analysis. DESeq2 contains tests specifically developed to detect significant differences in abundances for counts data. While DESeq2 is most commonly utilised for differential gene expression in RNASeq, using the phyloseq API a Biological Observation Matrix (BIOM) file is formatted for use within DESeq2 and differences in taxa abundances discovered. MetaDEGalaxy also offers extensive graphing capabilities by wrapping the comprehensive metagenomics R-package phyloseq⁶. Extensive graphing options are available within MetaDEGalaxy wrapping most functions offered within phyloseq which offer the user a high level of control. Additionally, user supplied metadata files can be input to DESeq2 for model generation and to phyloseq for enhanced graphing capabilities allowing for grouping, clustering, and colouring of all graph types based on metadata information. All software wrapped within the workflow is open-source software, a current limitation of existing workflows such as usearch⁸ within the popular QIIME package ². Finally, MetaDEGalaxy is designed within the popular Genomic Virtual Lab⁴ leveraging the functionality of this robust infrastructure.

Methods

Workflows

Four comprehensive MetaDEGalaxy tutorial are currently available in github which demonstrate how to work with both overlapping and non-overlapping 16S paired end Illumina reads.

Tutorial #1 details the workflow for data QC and the detection of paired end overlap in sequencing data and preparing FastQ files for metagenomic analysis (Figure 1). Tutorial #2 details the entire workflow for overlapping paired end Illumina reads (Figure 2) using the same data set employed by the Mothur_SOP run with the popular Mothur software (v1.35.1)¹⁴. This workflow inputs a group of paired-end MiSeq files and a metadata map file and generates overlapping FASTQ files, an annotated BIOM file, a DESeq2 table of differentially expressed microbes, and a variety of phyloseq graphs. Tutorial #3 details the entire workflow for non-overlapping paired end Illumina reads and is similar to tutorial #2 with the exception of pre-processing steps transforming FASTQ files into a Fasta file where PEAR¹⁰ software is not run. Finally, tutorial #4 details a workflow for BIOM file processing and analyses detailing how to utilise the platform for analyses starting from an input BIOM file.

Figure 1. Workflow 1 in MetaDEGalaxy for data QC and detecting PE read overlap.

Figure 2. Workflow 2, 3, and 4 for differential abundance detection of operational taxonomic units (OTUs).

Both workflow 2 and 3 use all the components in the workflow, the only difference is workflow 2 takes in paired-end reads data as input and workflow 3 take single-end reads data as input. The workflow 4 is the subset of the main workflow which starts with blue boxes and ends with all plots generated.

Results

To demonstrate some of the advanced functionality of MetaDEGalaxy, we follow tutorial #2 using the Mothur_SOP data to first generate a normalised count table and a table of differentially abundant OTUs (Table 2). The differentially abundant OTU table is formatted in DESeq2 output with additional taxonomic information appended to each row.

We use this table of differentially abundant OTUs to next generate a symmetric plot. Users are able to select any taxonomic level as well as any metadata variable for comparison and further to pick two values of this variable for direct comparison (Figure 3). In this example, we pick Phylum for our taxonomy level and time as our variable of interest and group the graph according to ‘Early’ or ‘Late’. The resulting symmetric plot shows the differences in OTUs for ‘Early’ and ‘Late’ samples across different phylum (Figure 4). We are also able to generate alpha diversity abundance plots according to various sample attributes grouped here for ‘Replicate Group’ and coloured by ‘Food’ (Figure 5). As a final example, we generate a network plot where we select ‘Replicate group’ for the correlation and select ‘Food’ as the legend (Figure 6).

Figure 3. MetaDEGalaxy menu options for generating symmetric plots for differentially abundant operational taxonomic units (OTUs).

Users are able to select the taxonomic rank to examine in addition to two values within any user-defined metadata category.

Figure 4. Symmetric plot of the most differentially abundant operational taxonomic units (OTUs) grouped by ‘Time’ with ‘Early’ and ‘Late’ samples compared.

Figure 5. Alpha diversity abundance plots grouped for replicate group and coloured by ‘Food’.

Figure 6. Network plots grouped for replicate group and coloured by ‘Food’.

Use cases

To demonstrate how to use MetaDEGalaxy we offer four in-depth tutorials describing available workflows. Tutorials 1, 2 and 4 utilise the same input data as the well-documented Mothur_SOP while tutorial 3 utilises custom 300bp paired end, non-overlapping Illumina MiSeq data. In either use case, reads can be accessed and pre-processed via Galaxy Interface with the following steps:

1) click on "Operations on multiple dataset" on the top of the history panel

2) check the box for all paired-end files listed on the history panel

3) click on the "For all selected..." button the top of the history panel

4) click on "Build list of Dataset Pairs" on the drop-down menu

5) Type in a common field of the file name for both forward and reverse paired end data

6) click on the "Auto-pair"

7) Enter a name for the collection of paired datasets and click "Create list"

Apart from the paired-end reads in data collection, users are required to have loaded the metadata table and both 16S reference genome and annotation files. When the paired-end reads from a data collection is imported into a Galaxy history, an important step for the later in the workflow is the renaming of the FASTA sequence header by appending the sample ID to end at the end of each read ID using the reheader tool in Galaxy. This information will be used as the column header for OTU table generated by the workflows.

Workflow 1 (Figure 1) is designed to detect the status of overlapped paired-end reads data using PEAR. Users should proceed with workflow 2 if the percentage of overlapped paired-end reads data is high. Otherwise, workflow 3 should be used for non-overlapping reads. Both workflow 2 and 3 are fundamentally the same (Figure 2), however, workflow 3 can take single-end reads data as input when the overlapped paired-end reads are not overlapping.

Workflow 4 is designed to take a precomputed BIOM file as input. BIOM file format is designed to store OTU counts, metadata, and OTU annotation into one file. When users input a BIOM file, workflow 4 can be used to add metadata to an existing BIOM file and create abundance bar plot, network plot and symmetric plots using phyloseq R package.

More detailed tutorial documentation is available in the github repository.

Conclusion

MetaDEGalaxy is a complete end-to-end Galaxy workflow for 16S differential abundance analysis. Harnessing the power of open source algorithms such as vsearch, phyloseq, and DESeq2, MetaDEGalaxy offers users high-level of control over their data and analysis options. Focusing on discovering the most differentially abundant OTUs between samples, MetaDEGalaxy allows users to assess the impact of different environmental condition on overall microbial community composition.

Data availability

Software availability

Software available from: http://203.101.224.124/galaxy

Source code available from: https://github.com/QFAB-Bioinformatics/jcu.microgvl.ansible.playbook

Archived source code at time of publication: https://doi.org/10.5281/zenodo.2658835¹⁶

Licence: GNU General Public License v3.0 for all script/wrappers