probiotics_amplicons_tank.Rmd

---
title: "Probiotics_Trials-main"
author: "JM"
date: "12/2/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Load libraries

```{r}
library(dada2)
library(ShortRead)
library(ggplot2)
library(phyloseq)
library(vegan)
library(knitr)
library(ALDEx2)
library(CoDaSeq)
library(zCompositions)
library(igraph)
library(car)
library(grDevices)
library(propr)
library(cowplot)
library(randomcoloR)
library(dplyr)
library(reshape2)
library(tibble)
library(exactRankTests)
library(nlme)
library(data.table)
library(Rmisc)
writeLines(capture.output(sessionInfo()), "sessionInfo.txt")
```

## Quality-filter the sequencing reads and create Amplicon Sequence Variant (ASV) tables with DADA2

Put unjoined R1 and R2 fastq files, with adaptors and primers previously removed with cutadapt into a directory for DADA2. Here, our forward and reverse fastq filenames have format: SAMPLENAME_R1_cut.fastq.gz and SAMPLENAME_R2_cut.fastq.gz

*****If you have samples from multiple sequencing runs, you need to determine the sequence variants for each run separately, then merge the ASV tables.
Here is the dada2 page on merging runs: https://benjjneb.github.io/dada2/bigdata_paired.html


Make sure the full path is updated at the beginning and ending of this code chunk
```{r}

#### 1st sequencing run
path <- "~/Desktop/Probiotics_Trials-main/cutadapt_NS1968R_Tank" ##update path as needed
list.files(path)
fnFs <- sort(list.files(path, pattern="_R1_cut.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_cut.fastq.gz", full.names = TRUE))
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

# Perform filtering and trimming
filt_path <- file.path(path, "filtered") 
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(150,150),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE
head(out)
# Learn the Error Rates, it TAKES TIME!
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)
# Dereplicate the filtered fastq files
derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepFs) <- sample.names
names(derepRs) <- sample.names
# Infer the sequence variants in each sample
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
# Inspecting the dada-class object returned by dada:
dadaFs[[1]]
# Merge the denoised forward and reverse reads:
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])
# Construct sequence table
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled")
rownames(track) <- sample.names
head(track)
write.table(track, "dada_read_stats1.txt",sep="\t",col.names=NA)

saveRDS(seqtab, "~/Desktop/Probiotics_Trials-main/cutadapt_NS1968R_Tank/seqtab.rds") ##update path as needed
```



```{r}

#### 2nd sequencing run
path <- "~/Desktop/Probiotics_Trials-main/cutadapt_NS1995_Tank" ##update path as needed
list.files(path)
fnFs <- sort(list.files(path, pattern="_R1_cut.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_cut.fastq.gz", full.names = TRUE))
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

# Perform filtering and trimming
filt_path <- file.path(path, "filtered") 
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(150,150),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE
head(out)
# Learn the Error Rates, it TAKES TIME!
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)
# Dereplicate the filtered fastq files
derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepFs) <- sample.names
names(derepRs) <- sample.names
# Infer the sequence variants in each sample
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
# Inspecting the dada-class object returned by dada:
dadaFs[[1]]
# Merge the denoised forward and reverse reads:
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])
# Construct sequence table
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled")
rownames(track) <- sample.names
head(track)
write.table(track, "dada_read_stats2.txt",sep="\t",col.names=NA)

saveRDS(seqtab, "~/Desktop/Probiotics_Trials-main/cutadapt_NS1995_Tank/seqtab.rds") ##update path as needed
```


```{r}

#### 3rd sequencing run
path <- "~/Desktop/Probiotics_Trials-main/cutadapt_NS2051_Tank" ##update path as needed
list.files(path)
fnFs <- sort(list.files(path, pattern="_R1_cut.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_cut.fastq.gz", full.names = TRUE))
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

# Perform filtering and trimming
filt_path <- file.path(path, "filtered") 
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(150,150),
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE) # On Windows set multithread=FALSE
head(out)
# Learn the Error Rates, it TAKES TIME!
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)
# Dereplicate the filtered fastq files
derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepFs) <- sample.names
names(derepRs) <- sample.names
# Infer the sequence variants in each sample
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
# Inspecting the dada-class object returned by dada:
dadaFs[[1]]
# Merge the denoised forward and reverse reads:
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])
# Construct sequence table
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled")
rownames(track) <- sample.names
head(track)
write.table(track, "dada_read_stats3.txt",sep="\t",col.names=NA)

saveRDS(seqtab, "~/Desktop/Probiotics_Trials-main/cutadapt_NS2051_Tank/seqtab.rds") ##update path as needed
```

Now that I have determined the ASV tables for all 3 sequencing runs, I can merge the ASV tables and remove chimera sequences.


```{r}
st1 <- readRDS("~/Desktop/Probiotics_Trials-main/cutadapt_NS1968R_Tank/seqtab.rds") ##update path as needed
st2 <- readRDS("~/Desktop/Probiotics_Trials-main/cutadapt_NS2051_Tank/seqtab.rds") ##update path as needed
st3 <- readRDS("~/Desktop/Probiotics_Trials-main/cutadapt_NS2051_Tank/seqtab.rds") ##update path as needed
st.all <- mergeSequenceTables(st1, st2, st3, repeats="sum") # You will get the message "Duplicated sample names detected in the sequence table row names." to let you know that there are duplicate names across samples - it is not an error, just a message.

#Remove chimeric sequences:
seqtab.nochim <- removeBimeraDenovo(st.all, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(st.all)

# Combine read stats from 3 runs and add chimera summary - THIS ISN'T WORKING RIGHT, COME BACK TO LATER
stat1 <- read.table("dada_read_stats1.txt",sep="\t",header=TRUE, row.names=1)
stat2 <- read.table("dada_read_stats2.txt",sep="\t",header=TRUE, row.names=1)
stat3 <- read.table("dada_read_stats3.txt",sep="\t",header=TRUE, row.names=1)
stats.all<-bind_rows(stat1, stat2,stat3) 
#write.table(stats.all, "dada_read_stats_all.txt",sep="\t",col.names=NA)


# Track reads through the pipeline
# As a final check of our progress, we’ll look at the number of reads that made it through each step in the pipeline
rowSums(seqtab.nochim)
# need to write this out to add to dada read stats

# SAVE the non-chimeric sequence variant table SO YOU DON'T HAVE TO REPEAT ALL OF THE ABOVE STEPS
saveRDS(seqtab.nochim, file="~/Desktop/Probiotics_Trials-main/tank.rds") ##update path as needed
# RELOAD THE SAVED INFO FROM HERE (if you have closed the project):
#seqtab.nochim <- readRDS("~/Desktop/Probiotics_Trials-main/tank.rds") ##update path as needed

```

## Assign taxonomy in DADA2

Make sure the taxonomy reference database is in your working directory. Keep the database file gzipped. Adjust path name below. This step is very time consuming.

When taxonomy assignment is complete, we will use base R and phyloseq to clean up the taxonomy table. First, we will replace NAs and empty cells with the lowest taxonomy classification available. Second, we will use phyloseq to remove reads that are classified as Eukaryotes or unclassified at the domain level (ie, we are keeping only Bacteria and Archaea because that is what our primers target).

```{r}
taxa <- assignTaxonomy(seqtab.nochim, "~/Desktop/Probiotics_Trials-main/silva_nr99_v138.1_train_set.fa.gz", multithread=TRUE) ##update path as needed

# FIX the NAs in the taxa table
taxon <- as.data.frame(taxa,stringsAsFactors=FALSE)
taxon$Phylum[is.na(taxon$Phylum)] <- taxon$Kingdom[is.na(taxon$Phylum)]
taxon$Class[is.na(taxon$Class)] <- taxon$Phylum[is.na(taxon$Class)]
taxon$Order[is.na(taxon$Order)] <- taxon$Class[is.na(taxon$Order)]
taxon$Family[is.na(taxon$Family)] <- taxon$Order[is.na(taxon$Family)]
taxon$Genus[is.na(taxon$Genus)] <- taxon$Family[is.na(taxon$Genus)]
write.table(taxon,"silva_taxa_table.txt",sep="\t",col.names=NA)
write.table(seqtab.nochim, "silva_otu_table.txt",sep="\t",col.names=NA)

# Create phyloseq object from otu and taxonomy tables from dada2, along with the sample metadata.
otu <- read.table("silva_otu_table.txt",sep="\t",header=TRUE, row.names=1)
taxon <- read.table("silva_taxa_table.txt",sep="\t",header=TRUE,row.names=1)
samples<-read.table("metadata_tank.txt",sep="\t",header=T,row.names=1)
OTU = otu_table(otu, taxa_are_rows=FALSE)
taxon<-as.matrix(taxon)
TAX = tax_table(taxon)
sampledata = sample_data(samples)
ps <- phyloseq(otu_table(otu, taxa_are_rows=FALSE), 
               sample_data(samples), 
               tax_table(taxon))
ps #447 taxa and 24 samples

# remove chloroplasts and mitochondria and Eukaryota
get_taxa_unique(ps, "Family") #175
get_taxa_unique(ps, "Order") #114
get_taxa_unique(ps, "Kingdom") #3
ps <- subset_taxa(ps, Family !="Mitochondria")
ps <- subset_taxa(ps, Order !="Chloroplast")
ps <- subset_taxa(ps, Kingdom !="Eukaryota")
ps <- subset_taxa(ps, Kingdom !="NA")
get_taxa_unique(ps, "Family") #172
get_taxa_unique(ps, "Order") #112
get_taxa_unique(ps, "Kingdom") #2
ps #429 taxa and 24 samples

# Now export cleaned otu and taxa tables from phyloseq for future reference
otu = as(otu_table(ps), "matrix")
taxon = as(tax_table(ps), "matrix")
metadata = as(sample_data(ps), "matrix")
write.table(otu,"silva_nochloronomito_otu_table.txt",sep="\t",col.names=NA)
write.table(taxon,"silva_nochloronomito_taxa_table.txt",sep="\t",col.names=NA)
```


Now, time to explore the data.

```{r}
# Read in data and create phyloseq object
otu <- read.table("silva_nochloronomito_otu_table.txt",sep="\t",header=TRUE, row.names=1)
taxon <- read.table("silva_nochloronomito_taxa_table.txt",sep="\t",header=TRUE,row.names=1)
samples<-read.table("metadata_tank.txt",sep="\t",header=T,row.names=1)
OTU = otu_table(otu, taxa_are_rows=FALSE)
taxon<-as.matrix(taxon)
TAX = tax_table(taxon)
sampledata = sample_data(samples)
ps <- phyloseq(otu_table(otu, taxa_are_rows=FALSE), 
               sample_data(samples), 
               tax_table(taxon))
ps #429 taxa and 24 samples
ntaxa(ps) #429
ps1<-filter_taxa(ps, function(x) mean(x) >1, TRUE)
ntaxa(ps1) #148
ps2<-filter_taxa(ps, function(x) mean(x) >2, TRUE)
ntaxa(ps2) #72
ps5<-filter_taxa(ps, function(x) mean(x) >5, TRUE)
ntaxa(ps5) #21
ps10<-filter_taxa(ps, function(x) mean(x) >10, TRUE)
ntaxa(ps10) #9
get_taxa_unique(ps, "Genus") #263
get_taxa_unique(ps1, "Genus") #104
get_taxa_unique(ps2, "Genus") #57
get_taxa_unique(ps5, "Genus") #16
get_taxa_unique(ps10, "Genus") #6

##### bar charts
ps_ra<-transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
otu = as(otu_table(ps_ra), "matrix")
write.table(otu,"silva_nochloronomito_otu_table_ra.txt",sep="\t",col.names=NA)
get_taxa_unique(ps_ra, "Order") #115
get_taxa_unique(ps_ra, "Class") #60
get_taxa_unique(ps_ra, "Genus") #263

n <- 60
palette <- distinctColorPalette(n)

sample_data(ps_ra)$coral<-factor(sample_data(ps_ra)$coral,levels=c("coral_25","coral_26","coral_27","coral_28"))
sample_data(ps_ra)$date<-factor(sample_data(ps_ra)$date,levels=c("pre","1d","3d","7d","21d","28d"))

pdf("barchart_Class.pdf",width=13)
p1=plot_bar(ps_ra, "date", fill="Class", facet_grid=.~coral)+
  geom_bar(aes(fill=Class), stat="identity",position="stack")+
  theme_bw()+
  theme(strip.text=element_text(face="bold", size=12))+
  theme(axis.text.x=element_text(size=12))+
  theme(axis.text.y=element_text(size=12))+
  scale_fill_manual(values=palette)+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "bottom")
p1
dev.off()

### Plot genus level with filtered dataset
ps2_ra<-transform_sample_counts(ps2, function(OTU) OTU/sum(OTU))
sample_data(ps2_ra)$coral<-factor(sample_data(ps2_ra)$coral,levels=c("coral_25","coral_26","coral_27","coral_28"))
sample_data(ps2_ra)$date<-factor(sample_data(ps2_ra)$date,levels=c("pre","1d","3d","7d","21d","28d"))
get_taxa_unique(ps2_ra, "Genus") #57

n2 <- 57
palette2 <- distinctColorPalette(n2)

pdf("barchart_Genus.pdf",width=13)
p2=plot_bar(ps2_ra, "date", fill="Genus", facet_grid=.~coral)+
  geom_bar(aes(fill=Genus), stat="identity",position="stack")+
  theme_bw()+
  theme(strip.text=element_text(face="bold"))+
  theme(axis.text.x=element_text(angle = 90))+
  scale_fill_manual(values=palette2)+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "bottom")
p2
dev.off()


```

I really want to know about the Pseudoalteromonas and if I can detect McH1-7


```{r}

Pseudoalt<-subset_taxa(ps_ra, Genus=="Pseudoalteromonas")
Pseudoalt
otu.Pseudoalt = as(otu_table(Pseudoalt), "matrix")
taxon.Pseudoalt = as(tax_table(Pseudoalt), "matrix")
meta.Pseudoalt = as(sample_data(Pseudoalt), "matrix")
otu.Pseudoalt<-as.data.frame(otu.Pseudoalt)
otu.Pseudoalt<-rownames_to_column(otu.Pseudoalt,var="Sample")
#export ASV sequences for supplemental table
write.table(otu.Pseudoalt,"Pseudoalteromonas_ASVs.txt",sep="\t",col.names=NA)

```


Did communities shift with treatment?


```{r}
###################### Perform center-log-ratio transformation on ASVs and calculate Aitchison Distance and principal components
                  
otu <- read.table("silva_nochloronomito_otu_table.txt",sep="\t",header=TRUE, row.names=1)
taxon <- read.table("silva_nochloronomito_taxa_table.txt",sep="\t",header=TRUE,row.names=1)
samples<-read.table("metadata_tank.txt",sep="\t",header=T,row.names=1)

# First, replace 0 values with an estimate (because normalization is taking log, can't have 0)
# Also transposing here, need samples as rows
d.czm <- cmultRepl(t(otu), method="CZM", label=0)
# Perform the center-log-ratio (CLR) transformation 
d.clr <- codaSeq.clr(d.czm)
# transpose matrix of CLR transformed data for ordination and dendrogram
E.clr <- t(d.clr)
# plot compositional PCA biplot (perform a singular value decomposition)
d.pcx <- prcomp(E.clr)
# calculate percent variance explained for the axis labels
pc1 <- round(d.pcx$sdev[1]^2/sum(d.pcx$sdev^2),2)
pc2 <- round(d.pcx$sdev[2]^2/sum(d.pcx$sdev^2),2)
xlab <- paste("PC1: ", pc1, sep="")
ylab <- paste("PC2: ", pc2, sep="")
biplot(d.pcx, cex=c(0.6,0.4), var.axes=F,scale=1, xlab=xlab, ylab=ylab)
summary(d.pcx)
str(d.pcx)
screeplot(d.pcx)

# replot PCA with ggplot2 (showing samples only)
df_out <- as.data.frame(d.pcx$x)
theme_set(theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()))
cols<-c("pre"="#000000","1d"="#999999","3d"="#D55E00","7d"="#E69F00","21d"="#0072B2","28d"="#56B4E9")
samples$coral<-factor(samples$coral,levels=c("coral_25","coral_26","coral_27","coral_28"))
samples$date<-factor(samples$date,levels=c("pre","1d","3d","7d","21d","28d"))

pdf("PCA.pdf",width=8.5)
p<-ggplot(df_out,aes(x=PC1,y=PC2,fill=samples$date,shape=samples$coral))
p<-p+geom_point(size=5)+theme(axis.title = element_text(size=14))+theme(axis.text=element_text(size=12))+
  theme(legend.title = element_text(size=14))+theme(legend.text = element_text(size=12))+
  scale_fill_manual(values=cols)+
  scale_shape_manual(values=c(21,22,23,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))
p + labs(x=xlab, y=ylab, fill="Date", shape="Coral") + coord_fixed()
dev.off()

###################### Use phyloseq/vegan to perform ANOSIM/PERMANOVA
                  
# set metadata as factors for anosim
coral<-as.character(samples$coral)
date<-as.character(samples$date)

# permanova between groups using Aitchison distance
dist.clr <- dist(E.clr)
perm<-adonis(dist.clr~coral*date,as(sample_data(ps),"data.frame"))
print(perm)


##### Beta Dispersion

# Calculate multivariate dispersions based on date
mod<-betadisper(dist.clr,date)
anova(mod)
plot(mod)
boxplot(mod)

# merge beta dispersion data and metadata to make a prettier boxplot
tapply(mod$distances, date, mean)
dis <- mod$distances
dis.melt <- melt(dis)
dis.melt$Sample <- rownames(dis.melt)
samples$Sample <- rownames(samples)
dis.treat <- merge(samples, dis.melt)
colnames(dis.treat)[4] <- "distance"

#run linear model to test significance
distlm <-lm(distance~date*coral, data=dis.treat)
summary(distlm)
anova(distlm)

dis.treat$date<-factor(dis.treat$date,levels=c("pre","1d","3d","7d","21d","28d"))
pdf("DistanceToCentroid.pdf",width=6, height=4)
p3<-ggplot(dis.treat,aes(x=date,y=distance))+
  geom_boxplot()+
  theme_bw()+
  geom_point(aes(color=coral),size=3)+
  scale_color_manual(values=c("#F0E442","#009E73","#CC79A7","#666666"),name="Coral")+
  theme(axis.title.x=element_blank())+
  #theme(legend.position="none")+
  theme(text=element_text(size=14))+
  theme(strip.text.y=element_text(face="italic",size=14))+
  ylab("Distance to Centroid")
p3
dev.off()

```