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

Metagenomic sequencing is an increasingly common tool in environmental and biomedical sciences. While software for detailing the composition of microbial communities using 16S rRNA marker genes is relatively mature, increasingly researchers are interested in identifying changes exhibited within microbial communities under differing environmental conditions. In order to gain maximum value from metagenomic sequence data we must improve the existing analysis environment by providing accessible and scalable computational workflows able to generate reproducible results. Here we describe a complete end-to-end open-source metagenomics workflow running within Galaxy for 16S differential abundance analysis. The workflow accepts 454 or Illumina sequence data (either overlapping or non-overlapping paired end reads) and outputs lists of the operational taxonomic unit (OTUs) exhibiting the greatest change under differing conditions. A range of analysis steps and graphing options are available giving users a high-level of control over their data and analyses. Additionally, users are able to input complex sample-specific metadata information which can be incorporated into differential analysis and used for grouping / colouring within graphs. Detailed tutorials containing sample data and existing workflows are available for three different input types: overlapping and non-overlapping read pairs as well as for pre-generated Biological Observation Matrix (BIOM) files. Using the Galaxy platform we developed MetaDEGalaxy, a complete metagenomics differential abundance analysis workflow. MetaDEGalaxy is designed for bench scientists working with 16S data who are interested in comparative metagenomics. MetaDEGalaxy builds on momentum within the wider Galaxy metagenomics community with the hope that more tools will be added as existing methods mature.


Amendments from Version 1
The most significant change to the manuscript is the dramatic expansion of the software discussion and comparison sections to include more web based and Galaxy based metagenomics offerings. Additions include MG-RAST, MetaPipe, MOCAT2, FROGS, GmT, A-Game, and ANASTASIA to name a few. This is reflected in Table 3. I added a broader discussion about where MetaDEGalaxy fits in relative to the ever expanding metagenomic software environment. I expanded the manuscript to include more details on tools with differential abundance options including calypso and mothur methods metastats and lefse.

Introduction
It is now recognized that there is a strong link between microbial communities in the human body and human health 1 . 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 2 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 3 or the Genomics Virtual Lab 4 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 this, challenges remain in fast moving research areas such as metagenomics with only a handful of complete metagenomic offerings currently available within the popular Galaxy framework. Currently, existing metagenomics options in Galaxy include ASaiM 5 , FROGS 6 , GmT 7 , A-Game 8 , and ANASTASIA 9 with QIIME2 recently becoming available in the Galaxy Toolshed. While there is overlap between their workflows, MetaDEGalaxy differs in its focus on differential abundance by incorporating the capabilities of phyloseq 10 and DESeq2 11 for complex differential analysis. DESeq2 contains tests specifically developed to detect differences between groups in abundances for counts data. While DESeq2 is most commonly utilised for differential gene expression in RNASeq, recent studies have shown RNA-Seq algorithms methods perform similarly or better than metagenomic specific algorithms 12 . MetaDEGalaxy also offers extensive graphing capabilities by wrapping the comprehensive metagenomics R-package phyloseq 10 . 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 13 within the popular QIIME package 2 . Finally, MetaDEGalaxy is designed within the popular Genomic Virtual Lab 4 leveraging the functionality of this robust infrastructure.

Input
MetaDEGalaxy accepts either 454 or Illumina paired end sequence FASTQ files that can be overlapping or non-overlapping. Users may alternatively input a pre-computed BIOM file if they do not require BIOM file generation. Additional functionality requires a sample specific tab-delimited metadata file formatted according to QIIME map file standards. This metadata information can be utilised for determining the model to employ within DESeq2 and to generate graphs grouped by various metadata attributes.

Implementation
In total, there are four workflows in MetaDEGalaxy (Table 1) which utilise a combination of external software and custom code. 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) 19 . 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 15 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.

Operation
The Galaxy environment is available for testing purposes at http://203.101.224.202/galaxy/ and will be available on Galaxy Australia server by the end of 2019 (https://usegalaxy.org. au/). The minimum system requirements for installing the MetaDEGalaxy are a 64-bit unix environment at 4Gb of memory.

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).

Software comparison
MetaDEGalaxy is compared to existing software in Table 3. There are comparable web and/or GUI based tools such as QIIME/ QIIME2 2 , MetaPipe 20 , MG-RAST 21 , MOCAT2 22 , Calypso 23 , Explicet 24 , and Megan 25 , however none of these tools except QIIME2 are currently available within the popular Galaxy framework. Within Galaxy there are several metagenomics offerings including ASaiM 5 , GmT 7 , A-Game 8 , and ANASTASIA 9 .
While many of the features of the tools overlap, MetaDEGalaxy is the only option within Galaxy combining DESeq2 11 with the full graphing capability of phyloseq 10 . MetaDEGalaxy is

Quality control and predetermination of 16S workflow utilisation
To detect percentage of paired-end reads that overlap each other by 10bp. This workflow randomly selected 1000 reads from each sample to perform the detection. If over 50% of the PE reads overlap each other by at least 10bp, it is recommended to use workflow 2. If less than 50% of PE reads overlap by at least 10bp, it is recommended to use workflow 3.

16S_DE_for_overlapPE
For use with datasets that are sequenced using overlapping paired-end reads

16S_DE_for_nonoverlapPE
For use with datasets that are sequenced using non-overlapping paired-end reads.

16S_BIOM
Handles Biological Observation Matrix (BIOM) file from workflows 2 and 3 to generate 5 plots (e.g. sample correlation network plot, symmetric plot and 3 abundance bar plots.  similar in features to GmT 7 however the differential abundance options are limited with GmT as it lacks symmetric plots and the ability to construct highly customisable graphs grouped by sample metadata attributes.
Differential abundance tables generated by MetaDEGalaxy and Calypso both use the phyloseq_to_deseq2 function in phyloseq which converts phyloseq formatted BIOM files into a DESeq ready object containing dispersion estimates and an experimental design formula based on a combination of metadata attributes. Mothur differs from these two methods in offering metagenomic specific algorithms including metastat 26 and lefse 27 . Metastats uses a t-test with p-values derived from an empiric null distribution calculated by sample permutation while lefse applies the non-parametric Kruskal-Wallis and Wilcoxon-Mann-Whitney tests to identify differences in gene abundance between metagenomic groups. Not surprisingly, results from MetaDEGalaxy and Calypso were identical while the results from lefse and metastats were quite different as has been shown by previous studies 28 .

Use cases
To demonstrate how to use MetaDEGalaxy we offer four indepth 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,    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
Source data Data used for the tutorials are available from Zenodo:  This submission introduces MetaDEGalaxy, which is a workflow intended for 16S differential abundance analysis in the open source Galaxy platform. The workflow incorporates various currently popular open source algorithms, the proposed workflow support the application of such methods by Galaxy users. In particular, the workflow supports differential OTU abundance testing for common measurement platforms (454 and Illumina).
Step-by-step tutorials are provided to support the use.
The overall work is sound and clearly written. Appropriate references are provided, and the work is based on commonly used methodologies and open source resources. Data and software are openly available with a unique DOI and permanent archiving through Zenodo.
The work does not contribute to methods criticism, validation, or benchmarking. This work is a technical contribution that provide new software plugin for the broader Galaxy platform. This is relevant for the limited community of researchers who use Galaxy for 16S microbiome analysis. The contribution is a contribution to scientific software, rather than scientific discussion. This, in my understanding, fits the F1000Research scope. 1.

Yes
No competing interests were disclosed. 1 3.
-GmT : Mothur SOP 16S end-to-end pipelines have been made available as Galaxy workflows previously.
-FROGS: Metagenomics pipelines in Galaxy has been previously described.
-Other Galaxy environments and workflows such as A-Game, ANASTASIA, or this functional annotation workflow, and several others have also been previously described. While these examples have a different focus than MetaDEGalaxy (i.e. functional metagenomics rather than 16S), the authors make the very broad claim that ASaiM is the only other "existing metagenomics workflow on offering", which is inaccurate.
Please consider expanding the discussion of existing work and Table 3 to include some or all of the above.
In the results section, please discuss how the differential abundance results obtained with the MetaDEGalaxy pipelines compare to the results described in the Mothur SOP ( , e.g. under subsection "population level analysis" of https://www.mothur.org/wiki/MiSeq_SOP section "OTU-based analysis") Do you determine the same OTUs to be statistically significantly different between the two groups? Explain any differences in results, as well as the added value of your approach over the statistical methods used in the SOP.
The Phyloseq wrappers the authors have created do not appear to be available from the Galaxy toolshed currently While I appreciate that the authors have developed Ansible playbooks for the installation of the wrappers, such a custom approach is not recommended practice, and adding the tools to the tool shed will greatly increase their accessibility.
One option for this would be to submit the wrappers to the IUC tool repository on github ( ) where they will be reviewed and automatically uploaded https://github.com/galaxyproject/tools-iuc to the toolshed upon acceptance.

Minor Remarks:
Training materials for the use of the MetaDEGalaxy workflows are available in the form of PDF files. I would strongly urge the authors to consider contributing these materials to the Galaxy Training Network ( ) so that they are more readily available for use https://training.galaxyproject.org by the community. I think that MetaDEGalaxy tutorial would be happily accepted there, and the GTN community can provide support to transform the tutorials into the right format.
Since the 4 workflows offered in this manuscript are designed to be run in succession (e.g. workflow 1,2,4 or 1,3,4), have the authors considered creating some full end-to-end workflows, using Galaxy's concept of sub-workflows?

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes
Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Partly No competing interests were disclosed.

Competing Interests:
Reviewer Expertise: Bioinformatics, Galaxy, 16S metagenomics, training I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above. The most significant change to the manuscript is the dramatic expansion of the software discussion and comparison sections to include more web based and Galaxy based metagenomics offerings. Additions include MG-RAST, MetaPipe, MOCAT2, FROGS, GmT, A-Game, and ANASTASIA to name a few. I added a broader discussion about where MetaDEGalaxy fits in relative to the ever expanding metagenomic software environment.
2) Expansion of differential abundance comparison I expanded the manuscript to include more details on tools with differential abundance options including calypso and mothur methods metastats and lefse. I performed a small side-by-side comparison however MetaDEGalaxy had results identical to calypso (which also uses phyloseq_to_deseq) and very different to both mothur techniques so I simply state this fact in the manuscript and cite previous work that has shown this already (Jonsson V, et al. Statistical evaluation of methods for identification of differentially abundant genes in comparative metagenomics. BMC Genomics. 2016).

3) Availability of wrapper scripts and tutorials (Major remark 4 and minor remark 1)
Thank you for the suggestions regarding better ways to make the software more widely available, I wasn't aware of these resources. The tutorials will indeed be made widely available to the training network once the installation is completed within Galaxy Australia which should be done by the end of 2019. Currently, the code is installed on a demo server however it will be given a permanent home within Galaxy Australia very shortly. We will also make the wrapper scripts available as per your suggestion to the IUC tool repository on github.
Also, for testing please note the temporary IP address for the demo server changed to which is now reflected in the new manuscript. http://203.101.224.202/galaxy/ No competing interests were disclosed. Competing Interests: