---
title: "DADA2_workflow"
output: html_notebook
---
  
Processing of filteredQC files from JGI iTags using the DADA2 pipeline (primers pre-trimmed using cutadapt). (https://benjjneb.github.io/dada2/tutorial.html)
  
  
# Download and install dada2, phyloseq with Bioconductor
```{r}
# dada2
source("https://bioconductor.org/biocLite.R")
biocLite("dada2")
library(dada2); packageVersion("dada2") #Should be version 1.4 for this pipeline

# phyloseq
biocLite('phyloseq')
library(phyloseq); packageVersion("phyloseq")

# Load ggplot2 package - phyloseq figures can be edited with ggplot
library(ggplot2)
```
  
  
# 16S sequeunce processing
  
```{r}
path <- "16S_filtered/" # CHANGE ME to the directory containing the demultiplexed and pre-processed fastq files
filtpath <- file.path(path, "filtered") # Filtered files go into the filtered/ subdirectory
fns <- list.files(path, pattern=".fastq.gz") # CHANGE if different file extensions

fns <- sort(list.files(path, pattern=".fastq.gz", full.names = TRUE))
plotQualityProfile(fn2[1:2])

# Filtering
out <- filterAndTrim(file.path(path,fns), file.path(filtpath,fns), 
              trimLeft = 21, truncLen = 250, truncQ = 2, maxEE=2, rm.phix=TRUE, maxN=0,
              compress=TRUE, verbose=TRUE, multithread=TRUE)

# File parsing
filtpath <- "16S_filtered/filtered/"
filts <- list.files(filtpath, pattern="fastq", full.names = TRUE)
sample.names <- sapply(strsplit(basename(filts), ".filter-ITAG_trimmed"), '[', 1)
names(filts) <- sample.names

# Learn error rates
set.seed(100)
err <- learnErrors(filts, multithread = TRUE, randomize=TRUE)

# Infer seq variants
dds <- vector("list", length(sample.names))
names(dds) <- sample.names
for(sam in sample.names) {
  cat("Processing:", sam, "\n")
  derep <- derepFastq(filts[[sam]])
  dds[[sam]] <- dada(derep, err=err, multithread = TRUE)
}

# Construct sequeunce table and write to disk
seqtab_16S <- makeSequenceTable(dds)
saveRDS(seqtab_16S, "seqtab_16S.rds")
seqtab_16S <- readRDS("seqtab_16S.rds")
  dim(seqtab_16S)

# Assign taxonomy
set.seed(100)
Taxa_16S <- assignTaxonomy(seqtab_16S, "rdp_train_set_16.fa.gz", multithread=TRUE, minBoot = 50)

# Save final tax table
saveRDS(Taxa_16S, "Taxa_16S.rds")
Taxa_16S <- readRDS("Taxa_16S.rds")


###########   Hand-off to phyloseq    ###########   


Bacteria_ps <- phyloseq(otu_table(seqtab_16S, taxa_are_rows=FALSE), 
               tax_table(Taxa_16S))

# Define a data frame from the metadata
sampleDF = read.csv("metadata.csv", row.names=1)

# Check that the rownames match the sample names
all(rownames(sampleDF) %in% sample_names(Bacteria_ps))

# Convert to "sample_data" class
sampledata = sample_data(sampleDF)

# Merge
Bacteria_ps = merge_phyloseq(sampledata, Bacteria_ps)

# Pre-processing steps for phyloseq object
Bacteria_ps_assign <- subset_taxa(Bacteria_ps, Phylum != "-1") # keeping only the taxa that were assigned to a phylum
Bacteria_ps_filt <- filter_taxa(Bacteria_ps_assign, function(x) sum(x > 3) > (0.10*length(x)), TRUE) # remove taxa not seen more than 3 times in at least 5% of the samples
Bacteria_ps_filt_rel <- transform_sample_counts(Bacteria_ps_filt, function(OTU) OTU/sum(OTU)) # Transform to relative abundances

# Write ASV table with taxa assignments
SeqIDs_16S <- as.data.frame(tax_table(Bacteria_ps_filt_rel)) #store the sequences in a dataframe
taxa_names(Bacteria_ps_filt_rel) <- paste0("ASV", seq(ntaxa(Bacteria_ps_filt_rel))) #replaces actual sequeunce with SV# id
ASV_table = as(otu_table(Bacteria_ps_filt_rel), "matrix")
Taxa_table = as(tax_table(Bacteria_ps_filt_rel), "matrix")
if(taxa_are_rows(Bacteria_ps_filt_rel)){ASV_table <- t(ASV_table)}
ASVdf = as.data.frame(ASV_table)
ASVdf_16S_taxa <- cbind(t(ASVdf), SeqIDs_16S, Taxa_table)
write.csv(ASVdf_16S_taxa, file="ASV_Taxa_table_16S.csv")
```
  
  
  
## ITS2 sequeunce processing
    
```{r}
path <- "ITS_filtered/" # CHANGE ME to the directory containing the demultiplexed and pre-processed fastq files
filtpath <- file.path(path, "filtered") # Filtered files go into the filtered/ subdirectory
fns <- list.files(path, pattern=".fastq.gz") # CHANGE if different file extensions

fns <- sort(list.files(path, pattern=".fastq.gz", full.names = TRUE))
plotQualityProfile(fns[1:2])

# Filtering
out <- filterAndTrim(file.path(path,fns), file.path(filtpath,fns), 
              trimLeft = 10, truncQ = 2, maxEE=2, rm.phix=TRUE, maxN=0,
              compress=TRUE, verbose=TRUE, multithread=TRUE)

# File parsing
filtpath <- "ITS_filtered/filtered/"
filts <- list.files(filtpath, pattern="fastq", full.names = TRUE)
sample.names <- sapply(strsplit(basename(filts), ".filter-ITAG_trimmed"), '[', 1)
names(filts) <- sample.names

# Learn error rates
set.seed(100)
err <- learnErrors(filts, multithread = TRUE, randomize=TRUE)

# Infer seq variants
dds <- vector("list", length(sample.names))
names(dds) <- sample.names
for(sam in sample.names) {
  cat("Processing:", sam, "\n")
  derep <- derepFastq(filts[[sam]])
  dds[[sam]] <- dada(derep, err=err, multithread = TRUE)
}

# Construct sequeunce table and write to disk
seqtab_ITS <- makeSequenceTable(dds)
saveRDS(seqtab_ITS, "seqtab_ITS.rds")
seqtab_ITS <- readRDS("seqtab_ITS.rds")
  dim(seqtab_ITS)

# Assign taxonomy
set.seed(100)
Taxa_ITS <- assignTaxonomy(seqtab_ITS, "sh_general_release_dynamic_28.06.2017.fasta", multithread=TRUE, minBoot = 50, tryRC = TRUE)

# Save final tax table
saveRDS(Taxa_ITS, "Taxa_ITS.rds")
Taxa_ITS <- readRDS("Taxa_ITS.rds")


###########   Hand-off to phyloseq    ###########   


Fungi_ps <- phyloseq(otu_table(seqtab_ITS, taxa_are_rows=FALSE), 
               tax_table(Taxa_ITS))

# Define a data frame from the metadata
sampleDF = read.csv("~/Desktop/PANG_JGI_CSP/ITS/JGI_CSP_ITS_metadata.csv", row.names=1)

# Check that the rownames match the sample names
all(rownames(sampleDF) %in% sample_names(Fungi_ps))

# Convert to "sample_data" class
sampledata = sample_data(sampleDF)

# Merge
Fungi_ps = merge_phyloseq(sampledata, Fungi_ps)

# Pre-processing steps for phyloseq object
Fungi_ps_assign <- subset_taxa(Fungi_ps, Phylum != "-1")
Fungi_ps_filt = filter_taxa(Fungi_ps_assign, function(x) sum(x > 3) > (0.05*length(x)), TRUE) # remove taxa not seen more than 3 times in at least 5% of the samples
Fungi_ps_filt_rel <- transform_sample_counts(fung_ps_filt, function(OTU) OTU/sum(OTU)) # Transform to relative abundances

# Write ASV table with taxa assignments
SeqIDs_ITS <- as.data.frame(tax_table(Fungi_ps_filt_rel)) #store the sequences in a dataframe
taxa_names(Fungi_ps_filt_rel) <- paste0("ASV", seq(ntaxa(Fungi_ps_filt_rel))) #replaces actual sequeunce with SV# id
ASV_table = as(otu_table(Fungi_ps_filt_rel), "matrix")
Taxa_table = as(tax_table(Fungi_ps_filt_rel), "matrix")
if(taxa_are_rows(Fungi_ps_filt_rel)){ASV_table <- t(ASV_table)}
ASVdf = as.data.frame(ASV_table)
ASVdf_ITS_taxa <- cbind(t(ASVdf), SeqIDs_ITS, Taxa_table)
write.csv(ASVdf_ITS_taxa, file="ASV_Taxa_table_ITS.csv")
```