Bacteria

Prepare environment

library(tidyverse); packageVersion('tidyverse')
library(reshape2); packageVersion('reshape2')
library(ampvis2); packageVersion('ampvis2')
library(phyloseq); packageVersion('phyloseq')
library(microbiomeMarker); packageVersion('microbiomeMarker')
library(microbiome); packageVersion('microbiome')
library(RColorBrewer); packageVersion('RColorBrewer')
library(vegan); packageVersion('vegan')
library(ggpubr); packageVersion('ggpubr')
library(MicEco); packageVersion('MicEco')
library(RVAideMemoire); packageVersion('RVAideMemoire')

path <- "...set path to your read files..."
list.files(path)

fnFs <- sort(list.files(path, pattern="_R1_001_trimmed.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001_trimmed.fastq", full.names = TRUE))

sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

Quality profiles of raw reads

plotQualityProfile(fnFs[1:2])

plotQualityProfile(fnRs[1:2])

Filter raw reads

filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(225,220),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=FALSE)

out

Quality profiles of filtered reads

plotQualityProfile(filtFs[1:2])

plotQualityProfile(filtRs[1:2])

Assess and plot errors

errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

plotErrors(errF, nominalQ=TRUE)

plotErrors(errR, nominalQ=TRUE)

Dereplicate

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

# Name the derep-class objects by the sample names
names(derepFs) <- sample.names
names(derepRs) <- sample.names

Sample inference

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)

Generate sequence table

seqtab <- makeSequenceTable(mergers)

# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

seqtab_subset <- seqtab[,nchar(colnames(seqtab)) %in% seq(426,428)]

# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab_subset)))

seqtab_subset.nochim <- removeBimeraDenovo(seqtab_subset, 
                                           method="consensus", 
                                           multithread=TRUE, 
                                           verbose=TRUE)

dim(seqtab_subset.nochim)

sum(seqtab_subset.nochim)/sum(seqtab_subset)

Track reads through the pipeline

getN <- function(x) sum(getUniques(x))

track <- cbind(out, 
               sapply(dadaFs, getN), 
               sapply(dadaRs, getN), 
               sapply(mergers, getN), 
               rowSums(seqtab_subset.nochim))

colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names

head(track)

Assign taxonomy

taxa <- assignTaxonomy(seqtab_subset.nochim, "../silva_nr_v132_train_set.fa.gz", multithread=TRUE)

taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)

write.csv2(seqtab.nochim, "../output/16S-asvmat.csv", row.names = T)
write.csv2(taxa, "../output/16S-taxmat.csv", row.names = T)

Ordering and colors

# Default order for ID and Stage variables
order_id <- c('GS', 'GU', 'LH', 'CP', 'FA', 'WS', 'EA', 'EO', 'EC', 'ES')
order_stage <- c('Larva', 'Pupa', 'Adult', 'Eggs')

# Default color palettes for ID, Stage, and Label variables
cols_stage <- c('Larva' = '#C6AC8F', 'Pupa' = '#5E503F', 'Adult' = '#0A0908', 'Eggs' = '#EAE0D5')

cols_label <- c('Larval haemolymph' = '#C7928F', 'Larval gut unsterile' = '#A35752', 
                'Larval gut sterile' = '#DDBDBB', 'Pupal cell pulp' = '#5E413F', 
                'Female abdomen' = '#4F4740', 'Wash solution' = '#70655C',
                'Eggs cage' = '#F1E0D0', 'Eggs ovarium' = '#E3C1A1', 
                'Eggs oviposition apparatus' = '#D09762', 'Eggs sterile' = '#C78243')

cols_id <- c('LH' = '#C7928F', 'GU' = '#A35752', 'GS' = '#DDBDBB', 'CP' = '#5E413F', 
             'FA' = '#4F4740', 'WS' = '#70655C', 'EC' = '#F1E0D0', 'EO' = '#E3C1A1', 
             'EA' = '#D09762', 'ES' = '#C78243')

Load the data

asvmat <- read.csv("../data/dada2-output/16S-asvmat.csv", sep = ";", row.names = 1)
taxmat <- read.csv("../data/dada2-output/16S-taxmat.csv", sep = ";", row.names = 1) 
metmat <- read.csv("../data/metadata.csv", sep = ";")

# Compare and check sequences, assign ASV labels to replace sequences in tables
seq_check <- data.frame(ASV_SEQS = colnames(otumat),
                        TAX_SEQS = rownames(taxmat),
                        ASV = paste0('OTU', seq(1:nrow(taxmat)))) %>% 
  mutate(TEST = ifelse(.$ASV_SEQS == .$TAX_SEQS, 1, 0))

# Sanity check to see if sequences in ASV and taxonomy table match
any(seq_check$TEST == 0) # If any values are 0, seqs do not match

Prepare tables

asvmat <- data.frame(t(asvmat)) %>% 
  rownames_to_column('ASV_SEQS') %>% 
  left_join(seq_check %>% select(ASV_SEQS, ASV)) %>% 
  column_to_rownames('ASV') %>% 
  select(-ASV_SEQS)

taxmat <- taxmat %>% 
  rownames_to_column('TAX_SEQS') %>% 
  left_join(seq_check %>% select(TAX_SEQS, ASV)) %>% 
  column_to_rownames('ASV') %>% 
  select(-TAX_SEQS)

Create phyloseq object

ps <- phyloseq(
  otu_table(as.matrix(asvmat), taxa_are_rows = T),
  tax_table(as.matrix(taxmat)),
  sample_data(metmat %>% column_to_rownames('SampleID')))

sample_variables(ps) # check available metadata variables
sample_names(ps) # check sample names
sample_data(ps) # check metadata table

# Make sure the variables are plotted in the right order
sample_data(ps)$ID <- factor(sample_data(ps)$ID, levels = order_id)
sample_data(ps)$Stage <- factor(sample_data(ps)$Stage, levels = order_stage)

Rarefaction

rarecurve(t(otu_table(ps)), step=50, cex=0.5)

Remove outlier and rarefy

# Check smallest sample size
min(sample_sums(ps))

# Remove smallest sample
ps.sub <- subset_samples(ps, sample_names(ps) != 'EA1')

# Check again for smallest sample size and remaining number of samples
min(sample_sums(ps.sub))
nsamples(ps.sub)

# Rarefy abundance data to even depth and plot rarefaction curves
ps.rarefied <- rarefy_even_depth(ps.sub, 
                                 rngseed = 1000, 
                                 sample.size = min(sample_sums(ps.sub)), 
                                 replace = F)

rarecurve(t(otu_table(ps.rarefied)), step=50, cex=0.5)

amp.rarefied <- phyloseq_to_ampvis2(ps.rarefied)

Alpha diversity

alpha <- alpha(ps.rarefied, index = 'all')

# Extract metadata
ps.rarefied.meta <- meta(ps.rarefied)

# Add selected alpha diversity measures to metadata table
ps.rarefied.meta$Shannon <- alpha$diversity_shannon
ps.rarefied.meta$Evenness <- alpha$evenness_pielou

# Set order and grouping for statistics
levels_stage <- levels(as.factor(ps.rarefied.meta$Stage))
pairs_stage <- combn(seq_along(levels_stage), 2, simplify = F, FUN = function(i)levels_stage[i])

Plot alpha diversity measures and add statics

• Shannon
• Pielou’s evenness
• Rank abundance

# Set shapes for each group of samples
shapes_id <- c('LH' = 16, 'GU' = 17, 'GS' = 18, 'CP' = 16, 'FA' = 16, 'WS' = 17, 'EC' = 16, 'EO' = 17, 'EA' = 18, 'ES' = 15)

# Plot alpha diversity values
shannon_p <- ggplot(ps.rarefied.meta, aes(x = Stage, y = Shannon)) +
  geom_violin(aes(fill = Stage), 
              colour = NA, trim = F, alpha = 0.4) +
  geom_point(aes(colour = ID, shape = ID), 
             size = 4, position = position_jitter(width = 0.15, seed = 12)) +
  geom_boxplot(aes(fill = Stage), 
               colour = 'white', width = 0.2, alpha = 0.25) +
  scale_fill_manual(values = cols_stage) +
  scale_color_manual(values = cols_id) +
  scale_shape_manual(values = shapes_id) +
  labs(x = 'Developmental stage',
       y = 'Shannon index',
       colour = '',
       shape = '') +
  guides(fill = 'none') +
  theme_classic()

# Add statistics
shannon_p <- shannon_p + stat_compare_means(aes(label = ..p.signif..), comparisons = pairs_stage, ref.group = '0.5', method = 'wilcox.test')

# Assign sample groups to developmental stages to subset legend
larva <- c('LH', 'GU', 'GS')
pupa <- c('CP')
adult <- c('FA', 'WS')
eggs <- c('EC', 'EO', 'EA', 'ES')

# Create separate legends for each developmental stage
legend_larva <- get_legend(
  ggplot(ps.rarefied.meta %>% filter(Stage == 'Larva'), 
         aes(x = Stage, y = Shannon)) +
    geom_point(aes(colour = ID, shape = ID), 
               size = 4, position = position_jitter(width = 0.15, seed = 12)) +
    labs(fill = 'Larva', colour = 'Larva', shape = 'Larva') +
    scale_fill_manual(values = cols_stage[larva]) + 
    scale_color_manual(values = cols_id[larva]) + 
    scale_shape_manual(values = shapes_id[larva]))

legend_pupa <- get_legend(
  ggplot(ps.rarefied.meta %>% filter(Stage == 'pupa'), 
         aes(x = Stage, y = Shannon)) +
    geom_point(aes(colour = ID, shape = ID), 
               size = 4, position = position_jitter(width = 0.15, seed = 12)) +
    labs(fill = 'Pupa', colour = 'Pupa', shape = 'Pupa') +
    scale_fill_manual(values = cols_stage[pupa]) + 
    scale_color_manual(values = cols_id[pupa]) + 
    scale_shape_manual(values = shapes_id[pupa]))

legend_adult <- get_legend(
  ggplot(ps.rarefied.meta %>% filter(Stage == 'adult'), 
         aes(x = Stage, y = Shannon)) +
    geom_point(aes(colour = ID, shape = ID), 
               size = 4, position = position_jitter(width = 0.15, seed = 12)) +
    labs(fill = 'Adult', colour = 'Adult', shape = 'Adult') +
    scale_fill_manual(values = cols_stage[adult]) + 
    scale_color_manual(values = cols_id[adult]) + 
    scale_shape_manual(values = shapes_id[adult])) 

legend_eggs <- get_legend(
  ggplot(ps.rarefied.meta %>% filter(Stage == 'eggs'), 
         aes(x = Stage, y = Shannon)) +
    geom_point(aes(colour = ID, shape = ID), 
               size = 4, position = position_jitter(width = 0.15, seed = 12)) +
    labs(fill = 'Eggs', colour = 'Eggs', shape = 'Eggs') +
    scale_fill_manual(values = cols_stage[eggs]) + 
    scale_color_manual(values = cols_id[eggs]) + 
    scale_shape_manual(values = shapes_id[eggs]))

# Arrange legends
legends <- ggarrange(legend_larva, legend_pupa, legend_adult, legend_eggs, nrow = 4, heights = c(0.4, 0.2, 0.3, 0.5))

# Combine the alpha diversity plot with its legends
shannon_p <- ggarrange(shannon_p + theme(legend.position = 'none'), ggarrange(legends,nrow = 2, heights = c(0.6, 0.4)), ncol = 2, widths = c(1.8, 0.3))

shannon_p

click figure to see larger version

pielou_p <- ggplot(ps.rarefied.meta, aes(x = Stage, y = Evenness)) +
  geom_violin(aes(fill = Stage), colour = NA, trim = F, alpha = 0.4) +
  geom_point(aes(colour = ID, shape = Tissue), size = 4, position = position_jitter(width = 0.15, seed = 12)) +
  geom_boxplot(aes(fill = Stage), colour = 'white', width = 0.2, alpha = 0.25) +
  scale_fill_manual(values = cols_stage) +
  scale_color_manual(values = cols_id) +
  labs(x = 'Developmental stage',
       y = 'Shannon index',
       colour = '',
       shape = '') +
  guides(fill = 'none') +
  theme_classic()

pielou_p <- pielou_p + stat_compare_means(aes(label = ..p.signif..), comparisons = pairs_stage, ref.group = '0.5', method = 'wilcox.test')

pielou_p

click figure to see larger version

amp_rankabundance(amp.rarefied, group_by = 'Stage', showSD = TRUE, log10_x = TRUE) +
  scale_colour_manual(values = cols_stage) +
  scale_fill_manual(values = cols_stage)

click figure to see larger version

Heatmap

Plot heatmap

amp_heatmap(
  amp.rarefied,
  group_by = 'SampleID',
  facet_by = 'Stage', 
  tax_aggregate = 'Genus',
  tax_add = 'Phylum',
  tax_show = 25, 
  normalise = T, 
  showRemainingTaxa = T, 
  plot_values_size = 3, 
  color_vector = c('white', '#A2708A'),
  plot_colorscale = 'sqrt',
  plot_values = F) +
  theme(axis.text.x = element_text(angle = 45, size = 8, vjust = 1),
        axis.text.y = element_text(size = 8),
        legend.position='right')

click figure to see larger version

Ordination

Canonical Correspondence Analysis

ordinationresult <- amp_ordinate(
  amp.rarefied, 
  type = 'CCA',
  constrain = 'ID',
  transform = 'Hellinger',
  sample_color_by = 'Stage',
  sample_colorframe = TRUE, 
  sample_colorframe_label_size = 4, 
  sample_colorframe_label = 'Stage', 
  sample_label_by = 'ID', 
  sample_label_size = 3.5,
  repel_labels = T,
  detailed_output = T)

• CCA
• Screeplot

ordinationresult$plot +
  scale_fill_manual(values = cols_stage) +
  scale_colour_manual(values = cols_stage) +
  labs(fill = '', colour = '') +
  theme_classic() +
  theme(legend.position = 'none')

click figure to see larger version

ordinationresult$screeplot

click figure to see larger version

LefSe

lefse <- run_lefse(ps.rarefied,
                   group = 'Stage',
                   taxa_rank = 'Genus')

Plot

plot_ef_bar(x) +
  scale_fill_manual(values = cols_stage) +
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.position = c(0.8, 0.1))

click figure to see larger version

Network

library(ggrepel)

set.seed(100)

ig <- make_network(ps.rarefied, max.dist = 0.5, distance = 'bray')

plot_network(ig, ps.rarefied, color = 'Stage', point_size = 4, type = 'samples', label = NA) +
  geom_text_repel(aes(label = value), size = 3.5, box.padding = 0.5) +
  scale_colour_manual(values = cols_stage) +
  theme(legend.position = c(0.8, 0.2))

click figure to see larger version