Update and add scripts, files and notebook following reviews

After receiving reviews from our submission to the Environmental
Microbiology journal, we needed to modify some of our analysis
and perform new ones. This included generating new scripts and
files and updating scripts and notebook.

.env/R package dependencies.R

@@ -3,16 +3,18 @@ options("repos" = c(CRAN = "https://cran.revolutionanalytics.com"))
 library("devtools")
 library("BiocManager")
 
-install_github("jalvesaq/colorout", ref="7ea9440", upgrade="never")
+install_github("jalvesaq/colorout", ref = "7ea9440", upgrade = "never")
 
 # Microbiome analysis
-install_bioc("3.10/microbiome", upgrade="never")
-install("phyloseq",   update=FALSE)
-install_version("phylogram",  version="2.1.0",  upgrade="never")      # prune function
-install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis", ref="6e09713", upgrade="never")
-install_github("ropensci/taxa",   ref="49969dc", upgrade="never")
-install_github("jbisanz/qiime2R", ref="08e3b30", upgrade="never")
+install_bioc("3.10/microbiome", upgrade = "never")
+install("phyloseq", update = FALSE)
+install_bioc("3.8/DESeq2", upgrade = "never")
+install_version("phylogram", version = "2.1.0", upgrade = "never")      # prune function
+install_version("usedist",   version = "0.4.0", upgrade = "never")
+install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis", ref = "6e09713", upgrade = "never")
+install_github("ropensci/taxa",   ref = "49969dc", upgrade = "never")
+install_github("jbisanz/qiime2R", ref = "08e3b30", upgrade = "never")
 
 # Pathway analysis
-install_bioc("3.9/ALDEx2", upgrade="never")
-install_github("ggloor/CoDaSeq/CoDaSeq", ref="73f77ce", upgrade="never")
+install_bioc("3.9/ALDEx2", upgrade = "never")
+install_github("ggloor/CoDaSeq/CoDaSeq", ref = "73f77ce", upgrade = "never")

README.md

@@ -1,17 +1,15 @@
 # The hemolymph of *Biomphalaria* snail vectors of schistosomiasis supports a diverse microbiome
 
 ## Abstract
-The microbiome – the microorganism community that is found on or within an organism’s body – is increasingly recognized to shape many aspects of its host biology and is a key determinant of health and disease. Microbiomes modulate the capacity of insect disease vectors (mosquitos, tsetse flies, sandflies) to transmit parasites and disease. We investigate the diversity and abundance of microorganisms within the hemolymph (i.e. blood) of *Biomphalaria* snails, the intermediate host for *Schistosoma mansoni*, using Illumina MiSeq sequencing of the V4 variable region of the bacterial 16S ribosomal DNA. We sampled hemolymph from 5 snails from 6 different laboratory populations of *B. glabrata* and one population of *B. alexandrina*. We observed 281.34 ± 0.81 amplicon sequence variants (ASVs) per snail, where ASVs are defined as rDNA sequences differing by >3%. There were significant differences in microbiome composition at the level of individual snails, snail populations and snail species. Snail microbiomes were dominated by bacteria from the phyla Proteobacteria or Bacteroidetes while water microbiomes from snail tanks were dominated by bacteria from the phylum Actinobacteria. We investigated the absolute bacterial load using qPCR: hemolymph samples contained between ~200 to ~60,000 16S rDNA.µL<sup>-1</sup> independently of the snail populations or species, equivalent to 2,620 ± 321 bacteria per µL. We speculate that the microbiome may represent a critical, but unexplored intermediary in the snail-schistosome interaction as hemolymph is in very close contact to the parasite at each step of its development.
-
-Full article is available [here](https://doi.org/XX)
+The microbiome – the microorganism community that is found on or within an organism’s body – is increasingly recognized to shape many aspects of its host biology and is a key determinant of health and disease. Microbiomes modulate the capacity of insect disease vectors (mosquitos, tsetse flies, sandflies) to transmit parasites and disease. We investigate the diversity and abundance of microorganisms within the hemolymph (i.e. blood) of *Biomphalaria* snails, the intermediate host for *Schistosoma mansoni*, using Illumina MiSeq sequencing of the bacterial 16S V4 rDNA. We sampled hemolymph from 5 snails from 6 different laboratory populations of *B. glabrata* and one population of *B. alexandrina*. We observed 279.84 ± 0.79 amplicon sequence variants (ASVs) per snail. There were significant differences in microbiome composition at the level of individual snails, snail populations and species. Snail microbiomes were dominated by Proteobacteria and Bacteroidetes while water microbiomes from snail tank were dominated by Actinobacteria. We investigated the absolute bacterial load using qPCR: hemolymph samples contained 2,784 ± 339 bacteria per µL. We speculate that the microbiome may represent a critical, but unexplored intermediary in the snail-schistosome interaction as hemolymph is in very close contact to the parasite at each step of its development.
 
 ## Code
 
 The code used to analyze the data presented in the mansucript is available in this repository.
 
 ## Prerequisites
 
-Two dependencies must be first installed before running the Jupyter notebook:
+Two dependencies must be installed before running the Jupyter notebook:
 * A [conda](https://docs.conda.io/en/latest/) distribution, like [miniconda](https://docs.conda.io/en/latest/miniconda.html),
 * [Jupyter Notebook](https://jupyter.readthedocs.io/en/latest/install.html) with a [bash kernel](https://github.com/takluyver/bash_kernel).
 

data/contaminants

@@ -0,0 +1,44 @@
+e50787fe1618dfaa8cc8e04234dc6791
+86684e6cdc85c7a5c58856b828bcb884
+e1577fcea6c47e98cea38de4370b5c1f
+074184aae53fd811aaca58c6c5265db0
+dd691275a5ea31f51da3087b72dd790a
+1bedf65dfacd928021f17a00c8612b66
+dbf58d31450edad4f8f9730a260506a3
+0689fc071a5c23924102830f03845214
+3b9e8bd7fdd289917b118dcc39b4e618
+23655f57cb92ac756a35c41de74a3a36
+dc81033655acba9e81a5979ca66c2757
+a12b17e567cf24472991424087af3475
+76c9efe521603bbd07685a394db8bf28
+2e353e92cf6952000207b2ab4205a520
+9e6ae76fc8a7239225f4599e9598568e
+a7dff99b3e4c718d62257f78a739b4dc
+514e94617ccbb4ba653bcd0350a006fd
+19552b6e6b0cb22ae961f6d9570e19d9
+ea58229144659c110a636abcf65c17e0
+c5f92d9214664e5be73f1ce66c9b5b3c
+e0b546c4f474e4c9bbf9af603d51d07b
+df9195bf04fb150775c756a2b590abe1
+0575272b42a9ef3a6df0d5e6ebd0bcce
+5db3904fa6a39ed6cb17eb991e00c302
+8e7335797ae974292389935004afc15b
+71bfdc2a09e86bc0370edff95f6acf42
+fb7ce3eb7f1994fc7478e4ec7230ad75
+2d12f175c5372ad1a8454991fcbb9787
+bdc46ef01e63f778ad51debd92b69e23
+56c12da5edd8287e6a2c95542f42d6fc
+70f21655b927e4a37014d975e05eaa88
+f8d02bc709958d9c1f4caf6ac443dc94
+b4b19e14617baaafa27936fb537a7193
+db9f192baa66535ba6886618a64f90e1
+15f6723f440b523f991f0d550de86f1f
+e2b34f0f52c81239fdfcc01376418972
+25f86eb314363cf858787ed9239821b1
+ef3da3c788b0e2b58350daa67923cb29
+15d42dc249618c7326229450b0f46b39
+e705f7c98eb78bdb260aaa71b40f9cf4
+320646b00aa077b31219f515a1a9f475
+24123f1ed7bffb358f47b0eacfb43e98
+af866d9d9c4276901df4c35646ad24ca
+e1b2b6ae4af6686429c21926cc0314a0

data/whole_snails_main_bacteria.tsv

@@ -0,0 +1,23 @@
+Name	Rank	Source	Reference
+Gemmatimonas aurantiaca	Genus	BgBS90	Allan et al. 2018
+Micavibrio aeruginosavorus	Genus	BgBS90	Allan et al. 2018
+Parabacteroides distasonis	Genus	BgBS90	Allan et al. 2018
+Halobacteriovorax marinus	Genus	BgBS90	Allan et al. 2018
+Pseudomonas	Genus	Caribbean and lab populations	Ducklow et al. 1979
+Acinetobacter	Genus	Caribbean and lab populations	Ducklow et al. 1979
+Aeromonas	Genus	Caribbean and lab populations	Ducklow et al. 1979
+Enterobacteriaceae	Family	Caribbean and lab populations	Ducklow et al. 1979
+Vibrio extorquens	Genus	Caribbean and lab populations	Ducklow et al. 1979
+Pasteurella	Genus	Caribbean and lab populations	Ducklow et al. 1979
+Pseudomonas	Genus	Caribbean populations	Ducklow et al. 1981
+Acinetobacter	Genus	Caribbean populations	Ducklow et al. 1981
+Aeromonas	Genus	Caribbean populations	Ducklow et al. 1981
+Moraxella	Genus	Caribbean populations	Ducklow et al. 1981
+Mycoplasma	Genus	lab snails	Huot et al. 2020
+Flavobacteriaceae	Family	Bg snails	Huot et al. 2020
+Pirellula	Genus	lab snails	Huot et al. 2020
+Planctomyces	Genus	lab snails	Huot et al. 2020
+Candidatus	Genus	lab snails	Huot et al. 2020
+Odyssella	Genus	lab snails	Huot et al. 2020
+Mesorhizobium	Genus	lab snails	Huot et al. 2020
+Pseudomonas	Genus	lab snails	Huot et al. 2020

scripts/library_stats.R

@@ -86,7 +86,7 @@ mylb.tb.s <- apply(mylb.tb, 1, function(x) c(mean(x), sd(x), max(x), min(x)) )
 rownames(mylb.tb.s) <- c("Mean", "SD", "Max", "Min")
 
 
-# Per populations and replicates (mean and proportion)
+# Per population and replicate (mean and proportion)
 mytb.s.mn <- split(mydata, paste(mydata[,"Population"], mydata[,"Replicate"])) %>% lapply(., function(x) sapply(x[,mycln], mean) ) %>% as.data.frame()
 mytb.s.sd <- split(mydata, paste(mydata[,"Population"], mydata[,"Replicate"])) %>% lapply(., function(x) sapply(x[,mycln], sd) )   %>% as.data.frame()
 mytb.s.p  <- apply(mytb.s.mn, 2, function(x) x/x[1])

scripts/microbiome_comparison.R

@@ -0,0 +1,115 @@
+#!/usr/bin/env Rscript
+# Title: microbiome_comparison.R
+# Version: 1.0
+# Author: Frédéric CHEVALIER <[email protected]>
+# Created in: 2020-09-02
+
+
+
+#==========#
+# Packages #
+#==========#
+
+cat("Loading packages...\n")
+
+suppressMessages({
+    # System
+    library("magrittr")
+    library("tidyr")
+
+    # Microbiome
+    library("qiime2R")
+    library("phyloseq")
+})
+
+
+
+#===========#
+# Variables #
+#===========#
+
+cat("Setting variables...\n")
+
+# Working directory
+setwd(file.path(getwd(), "scripts"))
+
+# Suppress warning messages
+options(warn=-1)
+
+# Qiime files
+asv.f  <- "../results/1-qiime/table.qza"
+taxa.f <- "../results/1-qiime/rep-seqs_taxa.qza"
+tree.f <- "../results/1-qiime/rooted-tree.qza"
+
+# Metadata file
+md.f   <- "../data/sample-metadata.tsv"
+
+# Comparison file
+comp.f <- "../data/whole_snails_main_bacteria.tsv"
+
+# Population order and color
+pop.data <- matrix(c("Ba",      "#ffac40",
+                     "Bg26",    "#bf4d4d",
+                     "Bg36",    "#a469d6",
+                     "Bg121",   "#ff61e8",
+                     "BgBRE",   "#6be3a7",
+                     "BgBS90",  "#87d687",
+                     "BgNMRI",  "#ab7a78",
+                     "Water",   "#7798e0"), byrow=TRUE, ncol=2)
+
+
+
+#=================#
+# Data processing #
+#=================#
+
+# Loading ASV data
+asv    <- read_qza(asv.f)$data
+asv.tb <- otu_table(asv, taxa_are_rows=TRUE)
+
+# Loading taxa data and reorganizing table
+taxa <- read_qza(taxa.f)$data
+
+## Split info
+cln.nm <- c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
+taxa   <- separate(taxa, 2, cln.nm, sep=";")
+
+## Cleaning names
+taxa[, (1:length(cln.nm))+1] <- apply(taxa[,(1:length(cln.nm))+1], 2, function(x) gsub(".*__","", x))
+
+## Polishing data
+rownames(taxa) <- taxa[,1]
+taxa <- taxa[,-c(1,ncol(taxa))]
+
+## Assign "Unassigned" to all level of unassigned ASV
+taxa[ taxa[,1] == "Unassigned", ] <- "Unassigned"
+
+## Convert table
+taxa.tb             <- tax_table(as.matrix(taxa))
+taxa_names(taxa.tb) <- rownames(taxa)
+
+
+# Loading tree
+tree <- read_qza(tree.f)$data
+
+# Loading metadata and reorder factors
+md <- read.delim(md.f, row.names=1)
+md[,"Population"] <- factor(md[,"Population"], pop.data[,1])
+
+# Building phyloseq object
+mybiom <- phyloseq(asv.tb, taxa.tb, tree, sample_data(md))
+
+
+# Load comparison file
+comp <- read.delim(comp.f, header=TRUE)
+
+# Generate names of genus or family from the table
+comp.ls        <- apply(comp, 1, function(x) strsplit(x[1], "\ ") %>% unlist() %>% .[1] %>% grep(., tax_table(mybiom)[,x[2]]))
+names(comp.ls) <- strsplit(as.character(comp[,1]), "\ ") %>% lapply(., function(x) x[1]) %>% unlist()
+
+# Identify populations carrying the genus or family
+comp.pop <- lapply(comp.ls, function(x) if (length(x) > 0) { (otu_table(mybiom)[x,] %>% colSums() > 0) %>% sample_data(mybiom)[., "Population"] %>% as.matrix() %>% as.vector() %>% unique()})
+comp.pop <- comp.pop[names(comp.pop) %>% unique()]
+
+# Print the results on screen
+print(comp.pop)

scripts/microbiome_density.R

@@ -1,15 +1,17 @@
 #!/usr/bin/env Rscript
 # Title: microbiome_density.R
-# Version: 1.0
+# Version: 1.1
 # Author: Frédéric CHEVALIER <[email protected]>
-# Created in: 2020-03-25
+# Created in: 2020-09-08
 
 
 
 #==========#
 # Packages #
 #==========#
 
+cat("Loading packages...\n")
+
 suppressMessages({
     library("qiime2R")
     library("magrittr")
@@ -41,6 +43,8 @@ q <- function(x) {
 # Variables #
 #===========#
 
+cat("Setting variables...\n")
+
 # Working directory
 setwd(file.path(getwd(), "scripts"))
 
@@ -87,6 +91,12 @@ mrkr.pred <- read.delim(picrust.f)
 # Note: no data cleaning because contaminant are also quantified with qPCR
 asv <- apply(taxa, 2, function (x) sum(x > 0) ) %>% data.frame(asv=.)
 
+# Update mrkr.pred with contaminants
+contam    <- ! rownames(taxa) %in% mrkr.pred[,1]
+contam.mt <- cbind(rownames(taxa)[contam], 1, 1)
+colnames(contam.mt) <- colnames(mrkr.pred)
+mrkr.pred <- rbind(mrkr.pred, contam.mt)
+
 # Reset 16S count prediction for ASV with nsis > 2 because of incertainty of annotation (see PiCRUST documentation)
 mrkr.pred[mrkr.pred[, 3] > 2, 2] <- 1
 
@@ -236,7 +246,7 @@ p.ls[[3]] <- ggplot(mydata, aes(x = V4, y = asv)) +
         stat_cor(method = "kendall", cor.coef.name = "tau", label.x.npc = "right", label.y.npc = "top", hjust = 1)
 
 p <- ggarrange(plotlist=p.ls, ncol=1, labels=LETTERS[1:length(p.ls)])
-ggsave(paste0(graph.d,"Fig. 6 - qPCR-diversity.pdf"), p, height=15, width=5)
+ggsave(paste0(graph.d,"Fig. 5 - qPCR-diversity.pdf"), p, height=15, width=5)
 
 
 # Clean tmp file