We share this use case to show all features of rANOMALY package that allow to highlight microbial community differences between groups of samples. rANOMALY main features are :

  • ASV (Amplpicon Sequence Variant) computation thanks to dada2 algorithm

  • Taxonomic assignment thanks to IDTAXA from DECIPHER package. rANOMALY allow to assign with up to 2 complementary reference databases (eg. SILVA which is generalist, DAIRYdb which is specific to dairy environments), keep the best and more confident taxonomic assignment, check for taxonomy incongruencies like empty fields or multiple ancestors.

  • phyloseq format to handle data and downstream statistical analysis.

  • Decontamination step thanks to decontam package. rANOMALY allows to filter ASV depending on ASV overall frequency (low abundance filter), and ASV prevalence (presence in 1 or more samples), specific filters like manual suppression of specific taxa are available to.

  • Complete statistical analysis workflow: rarefaction curves, community composition plots, alpha/beta diversity analysis with associated statistical tests and multivariate analysis ; differential analysis using four different methods (metacoder, DESeq2, metagenomeSeq, plsda).

  • At each step, plots and R objects are returned in command line AND saved in folders to be able to relaunch the analysis and to export plots.

1 Help

Each function has a detailed help accessible in R via ?{function}.

2 Tests datasets

The dataset can be downloaded via this link.

This tutorial assumes that you have extracted all the read file in a folder named reads along with the sample-metadata.csv file.

We share a 24 DNA samples test dataset extracted from rats feces at two different time (t0 & t50) and in two nutrition conditions. Rats were fed with two strains of some bacteria (wild type & mutant). Also, extraction control sample (blank) are included.

sm <- read.table("sample_metadata.csv", sep="\t",header=TRUE)
DT::datatable(sm)

3 Processing of raw sequences

3.1 ASV definition with DADA2

The first step is the creation of ASVs thanks to the dada2 package. In rANOMALY, only one function is needed to compute all the different steps require from this package. The function is designed to process Illumina paired end sequences in fastq format.

Sample names are extracted from the file name, from the beginning to the first underscore (_), so files must be formatted as followed : {sample-id1}_R1.fastq.gz {sample-id1}_R2.fastq.gz etc… Those must match the sample names stored in sample-metadata.csv file.

dada_res = dada2_fun(path="./reads", compress=TRUE, extension = "_R1.fastq.gz")

Main outputs:

  • ./dada2_out/read_tracking.csv summarizes the read number after each filtering step.

  • ./dada2_out/raw_otu-table.csv the raw ASV table.

  • ./dada2_out/rep-seqs.fna: fasta file with all representative sequences for each ASV.

  • ./dada2_out/robjects.Rdata with saved dada_res list containing raw ASV table and representative sequences in objects otu.table, seqtab.export & seqtab.nochim.

Read tracking table:

  • input: raw read number.

  • filtered: after dada2 filtering step: no N’s in sequence, low quality, and phiX.

  • denoisedF & denoisedR: after denoising. Forward & Reverse.

  • merged: after merging R1 & R2.

  • nonchim: after chimeras filtering.

DT::datatable(read.table("dada2_out/read_tracking.csv",sep="\t",header=TRUE), filter = "top")

3.2 Taxonomic assignment

assign_taxo_fun function uses IDTAXA function from DECIPHER package, and allows to use 2 different databases. It keeps the best assignment on 2 criteria, resolution (depth in taxonomy assignment) and confidence. The final taxonomy is validated for multiple ancestors and incongruities correction step.

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

tax.table = assign_taxo_fun(dada_res = dada_res, id_db = c("path_to_your_banks/silva/SILVA_SSU_r132_March2018.RData","path_to_your_banks/DAIRYdb_v1.2.0_20190222_IDTAXA.RData") )

Main file outputs:

  • ./idtaxa/robjects.Rdata with taxonomy in phyloseq format in tax.table object.

  • final_tax_table.csv the final assignation table that will be use in next steps.

  • allDB_tax_table.csv raw assignations from the two databases, mainly for debugging.

3.3 Phylogenetic Tree

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

tree = generate_tree_fun(dada_res)

Main outputs: - tree_robjects.Rdata with phylogenetic tree object in phyloseq format.

3.4 Phyloseq object

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

  • the asv table dada_res and the representative sequences seqtab.nochim from dada2_robjects.Rdata

  • a taxonomy table (taxtable) taxo_robjects.Rdata from taxo_robjects.Rdata

  • the phylogenetic tree tree from tree_robjects.Rdata

  • metadata from sample-metadata.csv

data = generate_phyloseq_fun(dada_res = dada_res, tax.table = tax.table, tree = tree, metadata = "./sample_metadata.csv")

Main output:

  • ./phyloseq/robjects.Rdata with phyloseq object in data for raw counts and data_rel for relative abundance.

3.5 Decontamination

The decontam_fun function uses decontam R package with control samples to filter contaminants. The decontam package offers two main methods, frequency and prevalence (and then you can combine those methods). For frequency method, it is mandatory to have the DNA concentration of each sample in phyloseq (and hence in the sample-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 automaticaly give the information. Just ask for it ;).

Our function integrates the basic ASV frequency (nb_reads_ASV/nb_total_reads) and overall prevalence (nb_sample_ASV/nb_total_sample) filtering. We included an option to filter out ASV based on their taxa names (for known recurrent contaminants).

data_decontam = decontam_fun(data = data, domain = "Bacteria", column = "type", ctrl_identifier = "control", spl_identifier = "sample", number = 100, krona= TRUE)

Or users can simply apply filters on ASVs frequencies and their prevalence in samples (without decontam features):

data_decontam2 <- decontam_fun(data = data, skip = TRUE, number = 100, freq = 5e-5, prev = 1, krona= TRUE)
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 146 taxa and 24 samples ]
## sample_data() Sample Data:       [ 24 samples by 5 sample variables ]
## tax_table()   Taxonomy Table:    [ 146 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 146 tips and 144 internal nodes ]
## refseq()      DNAStringSet:      [ 146 reference sequences ]

Main outputs:

  • robjects.Rdata with contaminant filtered phyloseq object named data.

  • Exclu_out.csv list of filtered ASVs for each filtering step.

  • Krona plot before and after filtering.

  • raw_asv-table.csv & relative_asv-table.csv.

  • venndiag_filtering.png.

4 Plots, diversity and statistics

We are currently developping a ShinyApp to visualize your data, sub-select your samples/taxons and do all those analyses interactively !!! ExploreMetabar

4.1 Rarefaction curves

In order to observe the sampling depth of each samples we start by plotting rarefactions curves. Those plots are generated by plotly which makes the plot interactive.

rareplot = rarefaction(data_decontam, "strain_time", 100 )
htmltools::tagList(list(rareplot))

4.2 Composition plots

The bar_fun function outputs a composition plot, in this example it reveals the top 10 genus present in our samples. The function allows to plot at different rank and to modify the number of taxas to show.

  • Ord1 option order the sample along the X axis.

  • Fact1 option control labels of the X axis. Fact1="sample.id" if you don’t want the sample to be renamed.

4.2.1 Raw abundance

bars_fun(data = data_decontam, top = 10, Ord1 = "strain_time", rank="Genus", relative = FALSE, split = TRUE, verbose = FALSE)

4.2.2 Relative abundance

bars_fun(data = data_decontam, top = 10, Ord1 = "strain_time", rank="Genus", relative = TRUE, split = TRUE, verbose = FALSE)

4.3 Diversity analyses

4.3.1 Alpha diversity

This function computes various alpha diversity indexes, it uses the estimate_richness function from phyloseq.

Available measures : c(“Observed”, “Chao1”, “ACE”, “Shannon”, “Simpson”, “InvSimpson”, “Fisher”).

It returns a list which contains:

  • a boxplot comparing conditions. ($plot)

  • a table of indices values. ($alphatable)

And for each of the computed indices :

  • an ANOVA analysis. (${measure}$anova)

  • a pairwise wilcox test result comparing conditions and giving the pvalue of each comparison tested. (${measure}$wilcox_col1, ${measure}$wilcox_col2_fdr, ${measure}$wilcox_col2_collapsed)

  • a mixture model if your data include repetition in sampling, ie. column3 option. (${measure}$anovarepeat, ${measure}$mixedeffect)

All this in a single function.

alpha <- diversity_alpha_fun(data = data_decontam, output = "./plot_div_alpha/", column1 = "strain", column2 = "time", column3 = "", supcovs = "", measures = c("Observed", "Shannon") )

4.3.2 Table of diversity index values

The table of values for each indices you choose to compute.

pander(alpha$alphatable, style='rmarkdown')
Table continues below
  sample.id time type strain strain_time
SB1-Sauv0 SB1-Sauv0 t0 sample wildtype wildtype_t0
SB10-Mut0 SB10-Mut0 t0 sample mutant mutant_t0
SB11-Mut0 SB11-Mut0 t0 sample mutant mutant_t0
SB12-Mut0 SB12-Mut0 t0 sample mutant mutant_t0
SB13-Sauv50 SB13-Sauv50 t50 sample wildtype wildtype_t50
SB14-Sauv50 SB14-Sauv50 t50 sample wildtype wildtype_t50
SB15-Sauv50 SB15-Sauv50 t50 sample wildtype wildtype_t50
SB16-Sauv50 SB16-Sauv50 t50 sample wildtype wildtype_t50
SB17-Sauv50 SB17-Sauv50 t50 sample wildtype wildtype_t50
SB18-Sauv50 SB18-Sauv50 t50 sample wildtype wildtype_t50
SB19-Mut50 SB19-Mut50 t50 sample mutant mutant_t50
SB2-Sauv0 SB2-Sauv0 t0 sample wildtype wildtype_t0
SB20-Mut50 SB20-Mut50 t50 sample mutant mutant_t50
SB21-Mut50 SB21-Mut50 t50 sample mutant mutant_t50
SB22-Mut50 SB22-Mut50 t50 sample mutant mutant_t50
SB23-Mut50 SB23-Mut50 t50 sample mutant mutant_t50
SB24-Mut50 SB24-Mut50 t50 sample mutant mutant_t50
SB3-Sauv0 SB3-Sauv0 t0 sample wildtype wildtype_t0
SB4-Sauv0 SB4-Sauv0 t0 sample wildtype wildtype_t0
SB5-Sauv0 SB5-Sauv0 t0 sample wildtype wildtype_t0
SB6-Sauv0 SB6-Sauv0 t0 sample wildtype wildtype_t0
SB7-Mut0 SB7-Mut0 t0 sample mutant mutant_t0
SB8-Mut0 SB8-Mut0 t0 sample mutant mutant_t0
SB9-Mut0 SB9-Mut0 t0 sample mutant mutant_t0
  Observed Shannon Depth
SB1-Sauv0 41 1.477 29356
SB10-Mut0 40 2.073 25949
SB11-Mut0 51 2.178 19918
SB12-Mut0 38 2.116 30395
SB13-Sauv50 46 2.691 21859
SB14-Sauv50 57 2.905 29029
SB15-Sauv50 50 2.793 25273
SB16-Sauv50 52 2.8 26802
SB17-Sauv50 49 2.624 24186
SB18-Sauv50 54 2.831 26666
SB19-Mut50 66 2.638 18077
SB2-Sauv0 26 2.099 24234
SB20-Mut50 72 2.721 28898
SB21-Mut50 79 3.062 21017
SB22-Mut50 81 2.81 25344
SB23-Mut50 84 3.175 28073
SB24-Mut50 90 3.148 30773
SB3-Sauv0 19 0.1962 34506
SB4-Sauv0 41 2.52 22110
SB5-Sauv0 46 1.923 30165
SB6-Sauv0 46 1.067 27763
SB7-Mut0 33 2.256 21560
SB8-Mut0 58 2.089 29000
SB9-Mut0 50 2.237 27296

4.3.3 Boxplots

The boxplots of those values.

alpha$plot

4.3.4 Tests on Observed index

4.3.4.1 ANOVA results

For each indices, you have access to the ANOVA test. Here we present the result for the “Observed” indice.

pander(alpha$Observed$anova)
Analysis of Variance Model
  Df Sum Sq Mean Sq F value Pr(>F)
Depth 1 49.36 49.36 0.5091 0.4838
strain 1 1877 1877 19.36 0.0002764
time 1 3649 3649 37.64 5.392e-06
Residuals 20 1939 96.96 NA NA

4.3.4.2 Wilcox test

Wilcox tests are made on each factor you have entered, and the combination of the two. Here “strain” and “time”.

4.3.4.3 Wilcox test for “strain” factor

pander(alpha$Observed$wilcox_col1)
  mutant
wildtype 0.043

4.3.4.4 Wilcox test for “time” factor

pander(alpha$Observed$wilcox_col2_fdr)
  t0
t50 0.001

4.3.4.5 Wilcox test for the collapsed factors

pander(alpha$Observed$wilcox_col2_collapsed)
  mutant_t0 mutant_t50 wildtype_t0
mutant_t50 0.002 NA NA
wildtype_t0 0.377 0.005 NA
wildtype_t50 0.336 0.002 0.008

4.4 Beta diversity

The diversity_beta_light function is equivalent to the first function but allows to generate specific tests and figures ready to publish in rmarkdown as in the example below. It is based on the vegan package function vegdist for the distance calculation and ordinate for the ordination plot.

We include statistical tests to ease the interpretation of your results. A permutational ANOVA to compare groups and test that centroids and dispersion of the groups are equivalent for all groups. User can inform col and cov arguments to assess PERMANOVA to determine significant differences between groups (eg. factor “strain” here and covariable “time”). A pairwise-PERMANOVA is processed to determine which condition is significantly different from another (based on p-value).

As return, you will get an object that contains:

  • An ordination plot ($plot)

  • The permANOVA results ($permanova)

  • The pairwise permANOVA ($pairwisepermanova)

beta_strain_time = diversity_beta_light(data_decontam, col = "strain_time", dist0 = "bray", ord0 = "MDS", output="./plot_div_beta_strain_time/", tests = TRUE)

4.4.1 Ordination plot

htmltools::tagList(list(ggplotly(beta_strain_time$plot)))

4.4.2 Statistical tests

  • permanova
pander(beta_strain_time$permanova)
Permutation test for adonis under reduced model
  Df SumOfSqs R2 F Pr(>F)
Depth 1 0.5075 0.06845 4.643 0.000999
strain_time 3 4.83 0.6515 14.73 0.000999
Residual 19 2.077 0.2801 NA NA
Total 23 7.415 1 NA NA
  • pairwise permanova on concatenated factors
pander(beta_strain_time$pairwisepermanova)
Table continues below
pairs Df SumsOfSqs F.Model R2 p.value
wildtype_t0 vs mutant_t0 1 0.9528 4.64 0.317 0.005
wildtype_t0 vs wildtype_t50 1 2.021 28.97 0.7434 0.004
wildtype_t0 vs mutant_t50 1 2.197 26 0.7223 0.005
mutant_t0 vs wildtype_t50 1 1.681 11.59 0.5369 0.002
mutant_t0 vs mutant_t50 1 1.57 9.826 0.4956 0.004
wildtype_t50 vs mutant_t50 1 1.818 75.22 0.8827 0.004
p.adjusted sig
0.005 *
0.005 *
0.005 *
0.005 *
0.005 *
0.005 *

4.5 Differential analyses

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 differentially abundant taxa between tested conditions.

4.5.1 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 sample 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 also produces pretty graphical trees, representing taxas present in both groups and coloring branches depending on the group in which this taxa is more abundant. Two of those trees are produced, the raw one and the non significant features filtered one (p-value <= 0.05).

out1 = metacoder_fun(data = data_decontam, output = "./metacoder", column1 = "strain_time", rank = "Family", signif = TRUE, plottrees = TRUE, min ="10", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
  • Table
 DT::datatable(out1$table, filter = "top", options = list(scrollX = TRUE))
  • Comparaison 1
 out1$wildtype_t0_vs_mutant_t0$signif
## [[1]]

  • Comparaison 2
 out1$wildtype_t50_vs_mutant_t50$signif
## [[1]]

4.5.2 DESeq2

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

Main output is a list with :

  • list\({comparison}\)plot: plot with significant features
  • list\({comparison}\)table: table with statistics (LogFoldChange, pvalue, adjusted pvalue…)
#fig.keep='all', fig.align='left', fig.width = 15, fig.height = 10
out2 = deseq2_fun(data = data_decontam, output = "./deseq/", column1 = "strain_time", verbose = 1, rank = "Family", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(out2$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(out2$wildtype_t50_vs_mutant_t50$table, filter = "top", options = list(scrollX = TRUE))

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

out3 = metagenomeseq_fun(data = data_decontam, output = "./metagenomeseq/", column1 = "strain_time", verbose = 1, rank = "Family", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(out3$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(out3$wildtype_t50_vs_mutant_t50$table, filter = "top", options = list(scrollX = TRUE))

4.5.4 Aggregate methods

The aggregate_fun function allows to merge the results from the three differential analysis methods used before, to obtain one unique table with all informations of significant differentially abundant features.

The generated table include the following informations :

  • seqid: ASV ID

  • Comparaison: Tested comparison

  • Deseq: differentially abundant with this method (0 no or 1 yes)

  • metagenomeSeq: differentially abundant with this method (0 no or 1 yes)

  • metacoder: differentialy abundant with this method (0 no or 1 yes)

  • sumMethods: sum of methods in which feature is significant.

  • DESeqLFC: Log Fold Change value as calculated in DESeq2

  • absDESeqLFC: absolute value of Log Fold Change value as calculated in DESeq2

  • MeanRelAbcond1: Mean relative abundance in condition 1

  • MeanRelAbcond2: Mean relative abundance in condition 2

  • Condition: in which the mean feature abundance is higher.

  • Taxonomy & representative sequence.

resF = aggregate_fun(data = data_decontam, metacoder = "./metacoder/metacoder_signif_Family.csv", deseq = "./deseq/", mgseq = "./metagenomeseq/", output = "./aggregate_diff/", column1 = "strain_time", column2 = NULL, verbose = 1, rank = "Genus", comp = "wildtype_t50~mutant_t50,wildtype_t0~mutant_t0")
ggplotly(resF$wildtype_t0_vs_mutant_t0$plot)
ggplotly(resF$wildtype_t50_vs_mutant_t50$plot)
DT::datatable(resF$table, filter = "top", options = list(scrollX = TRUE))

4.5.5 PLSDA (Partial Least Squares Discriminant Analysis)

Our PLS-DA/sparse PLS-DA analysis is based on the mixOmics package. These are supervised classification methods. They allow to define most important features discriminating our groups. The plsda_res function outputs a list of graphs and tables:

  • plsda$plotIndiv ordination plot from the PLSDA analysis.

  • splsda$plotIndiv ordination plot from the sparse PLSDA analysis.

  • loadings$comp{n} all the loadings plots for each component from the sPLSDA analysis. Loadings weights are displayed and colors correspond to the group in which each feature has its maximum mean value.

  • splda.loading_table Loading of all taxas for each component.

  • splsda$plotArrow arrow plot of samples.

plsda_res <- plsda_fun(data = data, output = "./plsda_family/", column1 = "strain_time", rank = "Family")
plsda_res$splsda.plotIndiv$graph

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp1.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp2.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp3.png')

knitr::include_graphics('plsda_family/splsda_loadings_strain_time_Family_comp4.png')

5 Miscellaneous function

5.1 Heatmap

heatmap_fun allows to plot relative abundance of the top abundant taxas (20 by default). Users can choose which taxonomic rank to plot and the factor used to separate plots.

heatmap_plot = heatmap_fun(data = data_decontam, column1 = "strain_time", top = 20, output = "./plot_heatmap/", rank = "Species")
heatmap_plot$plot

5.2 Shared taxa

  • Venn Diagram

ASVenn_fun allows to define shared taxa between specific conditions (up to 5 conditions). It generates a venn diagram and a table informing presence of each ASV in each condition with its taxonomy and sequence. Counts can be verified with this table. rank argument allows to generate output with desired taxonomic level.

outvenn = ASVenn_fun(data = data_decontam,output = "./ASVenn/", rank = "ASV", column1 = "strain_time", shared = TRUE)
  • Venn diagram:
outvenn$venn_plot

  • Table:
DT::datatable(outvenn$TABf, filter = "top", options = list(scrollX = TRUE))
  • Cytoscape phy2cyto_fun allows to generate input files (SIF format) for Cytoscape which is useful to visualize shared taxa. You can see below an example of representation with Cytoscape on our test data. Each nodes and arrows position and aesthetic can be easily modified.

5.3 Supplementary tools and features

csv2phyloseq_fun allows to import the 4 files (ASV, taxonomy, tree and sequences) in tabulated format to generate phyloseq object ready to use with rANOMALY.

assign_fasta_fun is used to assign sequences from any fasta files.

export_to_stamp_fun allows you to generate input files for STAMP.

split_table_fun allows to split the phyloseq object in multiple sub-object according to one factor (column).

update_metadata_fun allows you to easily modify sample_data part of the phyloseq object.

6 R Session info

sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Paris
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ranomaly_1.0.0 pander_0.6.5   DT_0.26       
## 
## loaded via a namespace (and not attached):
##   [1] fs_1.6.2                    matrixStats_0.63.0         
##   [3] bitops_1.0-7                devtools_2.4.5             
##   [5] psadd_0.1.3                 httr_1.4.6                 
##   [7] RColorBrewer_1.1-3          doParallel_1.0.17          
##   [9] profvis_0.3.7               tools_4.3.2                
##  [11] backports_1.4.1             utf8_1.2.3                 
##  [13] R6_2.5.1                    vegan_2.6-4                
##  [15] lazyeval_0.2.2              mgcv_1.9-0                 
##  [17] permute_0.9-7               urlchecker_1.0.1           
##  [19] withr_2.5.0                 prettyunits_1.1.1          
##  [21] gridExtra_2.3               VennDiagram_1.7.3          
##  [23] textshaping_0.3.6           cli_3.6.1                  
##  [25] Biobase_2.60.0              formatR_1.14               
##  [27] labeling_0.4.2              sass_0.4.6                 
##  [29] qdapTools_1.3.5             readr_2.1.4                
##  [31] Rsamtools_2.16.0            systemfonts_1.0.4          
##  [33] decontam_1.20.0             sessioninfo_1.2.2          
##  [35] plotrix_3.8-3               GA_3.2.3                   
##  [37] limma_3.56.0                readxl_1.4.1               
##  [39] rstudioapi_0.14             RSQLite_2.3.1              
##  [41] generics_0.1.3              shape_1.4.6                
##  [43] hwriter_1.3.2.1             gtools_3.9.4               
##  [45] crosstalk_1.2.0             vroom_1.6.3                
##  [47] car_3.1-2                   dplyr_1.1.2                
##  [49] Matrix_1.6-2                interp_1.1-4               
##  [51] biomformat_1.28.0           futile.logger_1.4.3        
##  [53] fansi_1.0.4                 S4Vectors_0.38.1           
##  [55] ggfittext_0.10.0            DECIPHER_2.28.0            
##  [57] abind_1.4-5                 lifecycle_1.0.3            
##  [59] yaml_2.3.7                  carData_3.0-5              
##  [61] SummarizedExperiment_1.30.1 gplots_3.1.3               
##  [63] rhdf5_2.30.1                Rtsne_0.16                 
##  [65] grid_4.3.2                  blob_1.2.4                 
##  [67] promises_1.2.0.1            crayon_1.5.2               
##  [69] venn_1.11                   miniUI_0.1.1.1             
##  [71] lattice_0.22-5              cowplot_1.1.1              
##  [73] pillar_1.9.0                knitr_1.42                 
##  [75] GenomicRanges_1.52.0        corpcor_1.6.10             
##  [77] mixOmics_6.23.4             admisc_0.32                
##  [79] codetools_0.2-19            fastmatch_1.1-3            
##  [81] glue_1.6.2                  ShortRead_1.58.0           
##  [83] data.table_1.14.8           remotes_2.4.2              
##  [85] vctrs_0.6.2                 png_0.1-8                  
##  [87] cellranger_1.1.0            gtable_0.3.3               
##  [89] cachem_1.0.8                xfun_0.39                  
##  [91] metacoder_0.3.6             S4Arrays_1.0.4             
##  [93] mime_0.12                   metagenomeSeq_1.42.0       
##  [95] survival_3.5-7              iterators_1.0.14           
##  [97] ellipsis_0.3.2              nlme_3.1-163               
##  [99] usethis_2.1.6               phyloseq_1.44.0            
## [101] bit64_4.0.5                 GenomeInfoDb_1.36.0        
## [103] rprojroot_2.0.3             bslib_0.4.2                
## [105] KernSmooth_2.23-22          colorspace_2.1-0           
## [107] BiocGenerics_0.46.0         DBI_1.1.3                  
## [109] ade4_1.7-22                 phangorn_2.11.1            
## [111] DESeq2_1.41.1               tidyselect_1.2.0           
## [113] processx_3.8.1              bit_4.0.5                  
## [115] compiler_4.3.2              chron_2.3-61               
## [117] microbiome_1.22.0           glmnet_4.1-7               
## [119] desc_1.4.2                  DelayedArray_0.26.2        
## [121] plotly_4.10.1               scales_1.2.1               
## [123] caTools_1.18.2              quadprog_1.5-8             
## [125] callr_3.7.3                 stringr_1.5.0              
## [127] digest_0.6.31               dada2_1.28.0               
## [129] rmarkdown_2.21              XVector_0.40.0             
## [131] htmltools_0.5.5             pkgconfig_2.0.3            
## [133] jpeg_0.1-10                 MatrixGenerics_1.12.0      
## [135] highr_0.10                  fastmap_1.1.1              
## [137] rlang_1.1.1                 htmlwidgets_1.6.2          
## [139] shiny_1.7.4                 farver_2.1.1               
## [141] jquerylib_0.1.4             jsonlite_1.8.4             
## [143] BiocParallel_1.34.1         RCurl_1.98-1.12            
## [145] magrittr_2.0.3              GenomeInfoDbData_1.2.10    
## [147] Rhdf5lib_1.20.0             munsell_0.5.0              
## [149] Rcpp_1.0.10                 ape_5.7-1                  
## [151] ranacapa_0.1.0              stringi_1.7.12             
## [153] zlibbioc_1.46.0             MASS_7.3-60                
## [155] plyr_1.8.8                  pkgbuild_1.4.0             
## [157] parallel_4.3.2              ggrepel_0.9.3              
## [159] deldir_1.0-6                Biostrings_2.68.1          
## [161] splines_4.3.2               multtest_2.56.0            
## [163] hms_1.1.3                   locfit_1.5-9.7             
## [165] ps_1.7.5                    igraph_1.4.2               
## [167] ggpubr_0.6.0                ggsignif_0.6.4             
## [169] Wrench_1.18.0               reshape2_1.4.4             
## [171] stats4_4.3.2                pkgload_1.3.2              
## [173] futile.options_1.0.1        evaluate_0.21              
## [175] latticeExtra_0.6-30         RcppParallel_5.1.7         
## [177] lambda.r_1.2.4              tzdb_0.4.0                 
## [179] foreach_1.5.2               httpuv_1.6.8               
## [181] tidyr_1.3.0                 purrr_1.0.1                
## [183] ggplot2_3.4.2               broom_1.0.4                
## [185] xtable_1.8-4                RSpectra_0.16-1            
## [187] rstatix_0.7.2               later_1.3.1                
## [189] viridisLite_0.4.2           ragg_1.2.5                 
## [191] rARPACK_0.11-0              tibble_3.2.1               
## [193] memoise_2.0.1               GenomicAlignments_1.36.0   
## [195] ellipse_0.4.3               IRanges_2.34.0             
## [197] cluster_2.1.4               here_1.0.1