rANOMALY: AmplicoN wOrkflow for Microbial community AnaLYsis

Raw sequence processing

Samples must have been previously demultiplexed into one file per sample with the file name following this syntax: {sampleid}_R[12].fastq. The denoising process is handled using the dada2 R package³ which produces amplicon sequence variants (ASV) as a taxonomic unit. This improves resolution of the potential presence of microbial organisms by using a prediction model to correct sequencing errors before aggregating similar sequences. rANOMALY handles processing any region of 16S (V1 to V9) and ITS amplicon (ITS1, ITS2) sequences, in which an additional step of cutadapt⁴ removes ITS probes left in some short sequences. For 16S amplicon, primers are trimmed based on the primer sequence length. ASV identifiers are sequences translated into MD5 hashes which are unique identifiers based on the DNA sequences offering the possibility to be compared between projects. This step results in an object with representing sequences, a raw ASV counts table and a text file containing statistics from the denoising process.

Taxonomic assignment

Taxonomic assignments of ASVs are carried out by IDTAXA (part of the DECIPHER package), an algorithm based on a machine learning method⁵. We implemented a functionality to compute assignment with two reference databases. For example, 16S amplicon sequences can be assigned with an environment specific database (DAIRYdb⁶, HITdb⁷, MIDAS⁸) and a general database like SILVA⁹ or GreenGenes¹⁰. The assignment with the best confidence and the lowest rank is kept, thus increasing assignment depth and accuracy. As taxonomy can differ between reference databases, we implemented a taxonomy validation step in rANOMALY to unravel taxonomy inconsistency like taxa with multiple ancestors or empty ranks. As supplementary features, rANOMALY functions allow users to create their own IDTAXA formatted specific reference database. The first step consists of filling the taxonomic table empty fields with the last known rank, and checking for taxonomy incongruencies as in the assignment function. A taxid file as used with RDP classifier¹¹ is constructed and then a last function takes as input, the corrected taxonomy table, the fasta file and the taxid file to generate the IDTAXA formatted database.

Phylogenetic tree

Phylogenetic tree is generated in three steps. First, sequences are aligned with AlignSeqs function from DECIPHER package¹² by the guided-tree method. Then, distance matrix is calculated with the dist.ml function from the phangorn package. And finally, neighbour joining and pml function computes the likelihood of the phylogenetic tree.

The abundance, metadata, taxonomic table, reference sequences and phylogenetic tree are merged into a phyloseq object¹³.

Decontamination

ASVs have the advantage of enabling the distinction between contaminants and the real community. We have integrated the decontam package¹⁴ into the workflow. Indeed using OTU based clustering methods can agglomerate contaminant sequences with real sample sequences, the whole cluster could hence be considered as contaminant by mistake. Working with ASVs allows the use of R package decontam which will sensitively exclude contaminant ASVs. It integrates two main methods, one based on the prevalence of the contaminant in the control samples, and another one based on the DNA concentration of the samples. Moreover, our decontamination step allows users to apply various filters such as low ASVs frequency, low ASVs prevalence in real samples, and the minimum number of reads per sample.

Graphical plots and statistical analyses

Statistical analyses are key features of rANOMALY workflow. Main descriptive analyses are integrated thanks to phyloseq functionalities. In addition, we automatized graphical representations and advanced statistical tests for alpha, beta diversity and composition plots. Above all, we have included the four most up-to-date differential analyses to assess differentially abundant features.

Community composition plots, alpha and beta diversity analyses

rANOMALY allows users to explore microbial community composition with three different types of plots : classical interactive bar plot¹⁵ of raw and relative taxa abundances, rarefaction curves to check sampling effort, and Krona interactive pie charts¹⁶.

Alpha diversity indices representing the specific richness are calculated (Richness, Simpson, Shannon...) with the vegan R package¹⁷. We added statistical tests such as multi-factors analysis of variance and pairwise Wilcoxon tests to assess significant differences between tested categories. We included repeated measures ANOVA to handle within-subjects variation. For example, it can be used when there are measures taken on the same individuals at different time points. This step outputs graphical representations and tables with results of statistical tests which are both saved in files and returned as list objects for markdown reports.

Beta diversity analysis allows users to estimate the community differences between two samples, it is also based on the vegan package¹⁷. rANOMALY can calculate all different distances such as the one based on ASVs abundance (BrayCurtis), rank based Jaccard indices, and phylogenetic distances as UniFrac, weighted Unifrac. Graphical representations PCoA, NMDS and more are available. This analysis can be processed at different taxonomic levels and categorical factors chosen by the user. The additional statistical test of PERMANOVA uses the distance matrix to determine if microbial communities of sample groups are significantly different from each other. We added a pairwise PERMANOVA¹⁸ test to confirm significant differences between specific group of samples.

Differential analysis

Differential analyses are meant to assess potential differentially abundant taxas between tested conditions chosen by users. rANOMALY wraps three methods:

The use of multiple methods for differential analysis allows the user to investigate which feature can be considered as differentially abundant between conditions. rANOMALY function using DESeq2 and metagenomeSeq outputs tables and plots with significant features. Metacoder outputs heat tree plots allowing users to infer differentially abundant features at each taxonomy rank and their position on the phylogenetic tree. With basic data management, and according to the identification of all significant differentially features (i.e. ASVs, genus, species) from the three methods, we generate a single table to find significant features in one or more methods and ease the interpretation. Additional information are added to the final table, like mean relative abundance for each feature and condition. A column in which condition features are significantly more abundant, features taxonomies and related sequences are added to the aggregated table. To complete differential analysis, we have included the well recognized PLS-DA from the mixOmics package²². It is a supervised classification method allowing users to identify features discriminating the sample groups.

Additional features and export

Functions and procedures are available to help the user to generate additional figures or to export the data to third party software. For instance, shared taxa between conditions are useful to explore. We use a function that can generate Venn diagrams, or for more complex visualisation, we can produce files readable by Cytoscape²³ to produce shared taxa networks. Krona diagrams can be displayed to explore sample microbial composition, where samples can be merged by a specified factor. rANOMALY allows users to generate inputs for STAMP²⁴, which is a graphical software that provides statistical hypothesis tests and exploratory plots for analysing taxonomic profiles.

Operation

rANOMALY requires R 3.6.3 or upper and can be run on any operating system with common specifications (1Go disk space, 4Go RAM, multicore CPU is recommended).

Install

Up to date code is hosted by the INRAE gitlab, users can simply download and install the package in R console with following command lines.

# Install ranomaly
install.packages("devtools")
devtools::install_git("https://forgemia.inra.fr/umrf/ranomaly")

# Call ranomaly package
library(ranomaly)
setwd("PATHTOWORKINGDIRECTORY")

# Create the analysis directory
dir.create("./analysis")

Processing of raw sequences

ASV definition with DADA2. The first step is to define ASVs thanks to the dada2 package. In rANOMALY, only one function is needed to compute all the different steps require from this package. Here sample names are extracted from the file name, thus be sure that files name match samples name in the metadata file.

# ASV resolution step, raw ASV table and representative sequence generation
dada_res = dada2_fun(amplicon = "16S", path = "./fastq_files", outpath = "./analysis/01_results_dada/",
    plot = FALSE, compress = FALSE, verbose = 1, paired = FALSE, torrent_single = FALSE,
    returnval = TRUE)

Main outputs of this function are:

# Adjust sample names to match with metadata
colnames(dada_res$otu.table) = tools::file_path_sans_ext(colnames(dada_res$otu.table))

Taxonomic assignment. assign_taxo_fun function uses IDTAXA function from DECIPHER package, and allows to use two different databases. It keeps the best assignment on two criteria, resolution (depth in taxonomy assignment) and confidence (value givenby IDTAXA). The final taxonomy is validated by checking for multiple ancestors (i.e. same species assigned to different genus) and incongruities (i.e. empty fields or incomplete lineage) correction step.

We share the latest databases we use in the IDTAXA format in this link. Users can also generate your own IDTAXA formatted database following those instructions and scripts we provide at this page.

# Taxonomy assignment step
tax.table = assign_taxo_fun(dada_res = dada_res, id_db = c("PATHTO/SILVA_SSU_r138_2019.RData",
    "PATHTO/DAIRYdb_v1.2.4_20200604_IDTAXA.rdata"), verbose = 3, output = "./analysis/02_idtaxa_silva_ddb")

Main file outputs:

Phylogenetic tree. The phylogenetic tree from the representative sequences is generated using phangorn and DECIPHER packages.

tree = generate_tree_fun(dada_res)

Phyloseq object. To create a phyloseq object, we need to merge four objects and one file:

# Phyloseq object generation
data = generate_phyloseq_fun(dada_res = dada_res, tax.table = tax.table, tree = tree,
    metadata = "./sample-metadata.csv", output = "./analysis/03_phyloseq/")

Decontamination. The decontam_fun function uses decontam R package along with control samples (PCR control) to filter out contaminants. The decontam package offers two main methods, frequency and prevalence (users can also combine those methods). For frequency method, it is mandatory to have the DNA concentration of each sample in the phyloseq object (and hence in the metadata.csv). The prevalence method does not need DNA quantification, this method allows to compare presence/absence of ASV between real samples and control samples and then identify contaminants.

Tips: sequencing plateforms often quantify the DNA before sequencing, but do not usually give the information. Just ask for it ;).

Our function integrates the basic ASV frequency $\left(Freq=\frac{nb\_reads\_ASV}{nb\_total\_reads}\right)$ and minimum prevalence in overall samples filtering. We have also included an option to filter out ASV based on their taxa names for known laboratory recurrent contaminants.

Main outputs:

Here we are going to filter out ASVs representing less than 0.1% (freq) of the reads and that are present in more than 4 samples (prev). Moreover, we are excluding "unassigned" taxas for this use case. Our sample dataset do not contain control samples, this step will be skipped.

data_filtered = decontam_fun(data = data, domain = "Bacteria", skip = TRUE, number = 4000,
    freq = 0.001, prev = 4, unassigned = TRUE, krona = TRUE, output = "./analysis/04_decontam/")

We obtain the final phyloseq object used for downstream analysis:

> data_filtered
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 59 taxa and 48 samples ]
sample_data() Sample Data:       [ 48 samples by 34 sample variables ]
tax_table()   Taxonomy Table:    [ 59 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 59 tips and 57 internal nodes ]
refseq()      DNAStringSet:      [ 59 reference sequences ]

Plots, diversity and statistic analyses

Rarefaction curves. In order to observe the sampling depth of each sample we start by plotting rarefaction curves. Those plots are generated by plotly which makes them interactive. (Figure 2)

rareplot = rarefaction(data_filtered, "source_location", 100)

Figure 2. Rarefaction plot for all samples, separated by source_location factor.

Community composition plot. The bars_fun function allows user to generate interactive community composition plot. Figure 3 presents the composition plot with relative abundances for the top 20 genera existing in our samples. The function allows to plot at different taxonomy rank and to modify the number of taxa to show.

Figure 3. Composition plot of relative abundance of 20 most abundant genus.

Here two arguments controlling the composition plot aesthetics:

source_location factor shows very different bacterial community between milk, cheese and cow teats environments.

# Relative abondance
prel = bars_fun(data_filtered, "Genus", top = 20, Ord1 = "source_location", Fact1 = "source_location",
    relative = TRUE, outfile = "plot_compo_rel.html")

Alpha diversity. The alpha_diversity_fun function can computes various alpha diversity indexes. It uses the estimate_richness function from phyloseq (Available measures : Observed, Chao1, ACE, Shannon, Simpson, InvSimpson, Fisher). Here we calculate ASV richness and Shannon index and carry out an analysis of variance on the source_location factor. Sequencing depth is automatically taken into account in this test. A pairwise wilcoxon test is added to ANOVA to define which group might be significantly different from others. Figure 4 shows boxplots of diversity indices, cow teats environment has much more ASV than other environments. Shannon index reveals more differences between cheese rinds and cheese cores. Cheese rinds show a higher Shannon index highlighting a more balanced bacterial community.

divAlpha = diversity_alpha_fun(data = data_filtered, output = "./analysis/07_all_res/plot_div_alpha/",
     column1 = "source_location", measures = c("Observed", "Shannon"))

Figure 4. Boxplot of diversity alpha indexes.

Results of the analysis of variance and pairwise wilcoxon test on Shannon index:

# > divAlpha$Shannon
# $anova
#                  Df Sum Sq Mean Sq F value  Pr(>F)
# Depth             1   0.83    0.83   9.203 0.00409 **
# source_location  3  98.37   32.79 365.661 < 2e-16 ***
# Residuals        43   3.86    0.09
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# $wilcox_col1
#           milk past_rind raw_rind
# past_rind    0        NA       NA
# raw_rind     0         0       NA
# teat         0         0        0

The alpha_diversity_fun function returns a list which contains:

And for each of the computed indices :

Beta diversity. The diversity_beta_light function allows to generate specific tests and figures ready to publish in rmarkdown report as in the example below. It is based on the vegan package function vegdist for the distance calculation and phyloseq-extended in addition to ordinate funtion for the ordination plot.

We include statistical tests to ease the interpretation of results. A permutational ANOVA is carried out on matrix distance to compare groups by testing if centroids and dispersion are equivalent for all groups. User have to inform col argument and optionally cov (covariable) to assess PERMANOVA to determine significant differences between groups. A pairwise-PERMANOVA is processed to determine which condition is significantly different from another (based on p-value).

As a return, you will get a list that contains:

Here we present results of beta diversity analysis on source_location factor, Figure 5 show ordination plot and it confirms the big differences between community of cheese, milk and cow teats.

divBeta = diversity_beta_light(psobj = data_filtered, output = "./analysis/07_all_res/plot_div_beta/",
     col = "source_location", cov = NULL, dist0 = "bray", ord0 = "MDS")

Figure 5. Ordination plot: MDS on Braycurtis distance.

The permanova tests on BrayCurtis distance shows significant p-value for the source_location factor. Pairwise permanova test is used to define which level of the factor tested is significantly different from others.

divBeta$permanova
# Permutation: free
# Number of permutations: 1000
#
# Terms added sequentially (first to last)
#
#                 Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)
# Depth            1    2.4044  2.4044  27.072 0.14822 0.000999 ***
# source_location  3    9.9981  3.3327  37.523 0.61634 0.000999 ***
# Residuals       43    3.8191  0.0888         0.23544
# Total           47   16.2217                 1.00000
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

divBeta$pairwisepermanova
#                        pairs Df SumsOfSqs   F.Model        R2 p.value p.adjusted sig
# 1        milk vs cheese_core  1 5.2215038 85.356493 0.7950753   0.001     0.0012   *
# 2        milk vs cheese_rind  1 4.8836384 53.310404 0.7078757   0.001     0.0012   *
# 3               milk vs teat  1 3.9377389 26.676358 0.5480352   0.001     0.0012   *
# 4 cheese_core vs cheese_rind  1 0.2640387  4.573476 0.1721068   0.008     0.0080   *
# 5        cheese_core vs teat  1 4.7206284 41.504931 0.6535702   0.001     0.0012   *
# 6        cheese_rind vs teat  1 4.3806123 30.384779 0.5800307   0.001     0.0012   *

Differential analysis. We choose three different methods to process differential analysis which is a key step of the workflow. The main advantage of the use of multiple methods is to cross validate deferentially abundant taxa between tested conditions. For this use case, we choose to focus on milk and cow teat environment to compare community at genus level.

Metacoder Metacoder is the most simple differential analysis tool of the three. Counts are normalized by total sum scaling to minimize the sample sequencing depth effect and it uses a Kruskal-Wallis test to determine significant differences between groups. The metacoder_fun function allows the user to choose the taxonomic rank, which factor to the test (column1), and a specific pairwise comparison (comp) to launch the differential analysis.

It produces pretty graphical trees, representing taxas present in both groups and coloring branches depending on which group this taxa is more abundant (Figure 6). Two trees are produced, a raw one, where everything is displayed and a filtered one where only significant features are represented (p-value <= 0.05).

# Metacoder
mtc = metacoder_fun(data = data_filtered, output = "./analysis/07_all_res/DA/metacoder",
    column1 = "source_location", rank = "Genus", signif = TRUE, plottrees = TRUE,
    min = "10", comp = "milk~teat")

Figure 6. Heattree of significant features generated by metacoder package.

Main output is a list with :

DESeq2 DESeq2 is a widely used method, primarily for RNAseq applications, for assessing differentially expressed genes between controlled conditions. Its use for metabarcoding datas is sensibly the same and well documented. The deseq2_fun allows to process differential analysis as metacoder_fun, and users can choose the taxonomic rank, the factor to test and which condition to compare. DESeq2 algorithm uses negative binomial generalized linear models with VST normalization (Variance Stabilizing Transformation).

Main output is a list with:

# DESeq2
deseq = deseq2_fun(data = data_filtered, output = "./analysis/07_all_res/DA/deseq/",
    column1 = "source_location", verbose = 1, rank = "Genus", comp = "milk~teat")

MetagenomeSeq MetagenomeSeq uses a normalization method able to control for biases in measurements across taxonomic features and a mixture model that implements a zero-inflated Gaussian distribution to account for varying depths of coverage. As deseq2_fun, metagenomeseq_fun returns a table with statistics and a plot with significant features for each comparison.

# MetagenomeSeq
mgSeq = metagenomeseq_fun(data = data_filtered, output = "./analysis/07_all_res/DA/metagenomeseq/",
    column1 = "source_location", verbose = 1, rank = "Genus", comp = "milk~teat")

Results aggregation step The aggregate_fun function allows to merge the results from the three differential analyses methods computed previously to obtain one unique table with all informations of significant differentially abundant features. Figure 7 shows the most significant and differential abundant Genera between the two environments.

# Aggregate results
resF = aggregate_fun(data = data_filtered,
    metacoder = "./analysis/07_all_res/DA/metacoder/metacoder_signif_Family.csv",
    deseq = "./analysis/07_all_res/DA/deseq/",
    mgseq = "./analysis/07_all_res/DA/metagenomeseq/",
    output = "./analysis/07_all_res/DA/aggregate_diff/",
    column1 = "source_location", column2 = NULL, verbose = 1,
    rank = "Genus", comp = "milk~teat")

Figure 7. Top significant and differentially abundant genera between Milk and Teat samples.

The generated table include the following fields:

Here is an overview of the aggregate_fun table informing in which methods each feature is significant, their DESeq2 LogFoldChange value, taxonomy and representative sequences:

as.data.frame(head(apply(resF$table, 2, stringr::str_sub, 1, 10)))
#        seqid Comparaison DESeq metagenomeSeq metacoder sumMethods   DESeqLFC absDESeqLFC
# 1 a0bc0ee9c9  milk_vs_te     1             1         1          3  4.9664456   4.9664456
# 2 13b622bf4d  milk_vs_te     1             1         1          3  8.6935035   8.6935035
# 3 9c7d3af28a  milk_vs_te     1             1         1          3 -3.3346494   3.3346494
# 4 d47e07353a  milk_vs_te     1             0         1          2 -2.2647154   2.2647154
# 5 7f32bad90a  milk_vs_te     1             0         1          2 -2.7186798   2.7186798
# 6 ae592ee947  milk_vs_te     1             1         1          3  5.2998345   5.2998345
#   MeanRelAbcond1 MeanRelAbcond2 Condition     Domain  ...     Species   sequence
# 1     1.177130e-     2.996581e-      milk k__Bacteri  ...  s__Lactoco ATCTTCGGCA
# 2     9.035658e-     2.193079e-      milk k__Bacteri  ...  s__Enterob ATATTGCACA
# 3     1.310357e-     1.481454e-      teat k__Bacteri  ...  s__Peptost ATATTGCACA
# 4     1.510200e-     8.532074e-      teat k__Bacteri  ...  s__Coryneb ATATTGCACA
# 5     8.272207e-     7.662663e-      teat k__Bacteri  ...  s__Turicib ATCTTCGGCA
# 6     2.477752e-     5.291575e-      milk k__Bacteri  ...  s__Sphingo ATATTGGACA

Miscellaneous function

Heatmap. User can generate an interactive heatmap with heatmap_fun function to explore relative abundance of top taxa through samples as showed in Figure 8.

# heatmap
heatmap_plot = heatmap_fun(data = data_filtered, column1 = "source_location", top = 100,
     output = "./analysis/07_all_res/plot_heatmap/", rank = "Genus")

Figure 8. Heatmap of top genuses relative abundance.

Venn diagram of shared taxa. Figure 9 shows Venn diagram comparing the three environment of this study, this method is useful to determine shared taxa between group of samples. This function uses venn::venn function which handles up to 7 groups comparison.

# venn
outvenn = ASVenn_fun(data = data_filtered, output = "./analysis/07_all_res/ASVenn/",
     rank = "Genus", column1 = "source_location", shared = TRUE)

Figure 9. Venn plot of shared genuses.

ASVenn_fun generate a table allowing user to find which taxa are shared between conditions.

outvenn$TABf  #table with shared taxa, taxonomies and sequences

as.data.frame(head(apply(outvenn$TABf, 2, stringr::str_sub, 1, 10)))
#   teat milk cheese_rind cheese_core        ASV   taxonomy        seq
# 1    1    1           1           1 a0bc0ee9c9 k__Bacteri ATCTTCGGCA
# 2    1    1           1           1 13b622bf4d k__Bacteri ATATTGCACA
# 3    1    1           1           1 15875b7f67 k__Bacteri ATATTGCACA
# 4    0    0           1           1 fe8453c0a4 k__Bacteri ATATTGCACA
# 5    1    1           1           1 907e61bc24 k__Bacteri ATATTGCACA
# 6    0    0           1           1 fe1733ab7e k__Bacteri ATATTGCACA

Export and compatibility

phy2cyto_fun allows to generate input files (SIF format) for Cytoscape which is useful to visualise shared taxa and easily modify each nodes and arrows position and aesthetic.

phy2tsv_fun function allows user to generate tabulated format abundance table at different taxonomic rank with following commnands:

for (i in c("ASV", "Species", "Genus", "Family")) {
     phy2tsv_fun(data_prune, output = glue::glue("./analysis/08_tsv_table/{i}_table/"),
          rank = i)
}

export_to_stamp_fun function creates two files that can be imported into the STAMP software.

export_to_stamp_fun(data = data, output = "./analysis/stamp/", correc = TRUE)

csv2phyloseq_fun function allows user to import data from other bioinformatic pipeline like FROGS, Qiime2. This function needs tabulated ASV, taxonomy, metadata and DNA sequence table as inputs. It can generate phylogenetic tree if missing and output a phyloseq object ready for downstream analyses.