```{r lib2, echo=T, results='show', message=FALSE, warning=FALSE}
library(dada2) # Bioconductor
library(phyloseq) # Bioconductor
library(Biostrings) # Bioconductor
library(tidyverse)
library(ggthemes)
theme_set(theme_bw(base_size = 12))
```

```{r load_dada2_data, echo=T, results='show', message=FALSE, warning=FALSE}

path2dada2 <- "MiSeq_SOP" # Name directory with .fastq files
# Forward and reverse fastq filenames have format: SAMPLENAME_R1_001.fastq and SAMPLENAME_R2_001.fastq
fnFs <- sort(list.files(path2dada2, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path2dada2, pattern="_R2_001.fastq", full.names = TRUE))
# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

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

# Filter and trim input .fastq files
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
              maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
              compress=TRUE, multithread=TRUE) 

```

```{r error_rates, echo=T, results='show', message=FALSE, warning=FALSE}
# These function "learn" errror rates. By using the distribution of nucleotides
# paired with underlying quality scores, dada2 can make a best guess at PCR or sequencing errors
# using machine learning.
errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)
plotErrors(errF, nominalQ=TRUE)
```

```{r dereplicate, echo=T, results='show', message=FALSE, warning=FALSE}
# Dereplication cuts down later computation steps by combining identical sequencing 
# reads into into “unique sequences” with a corresponding “abundance” to cut down computation
# time. This step is no longer necessary in the newest version of dada2. 
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
```

```{r sample_inference, echo=T, results='show', message=FALSE, warning=FALSE}
# Sample inference (group reads into 16s amplicons)
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
saveRDS(dadaFs,
        "dadaFs.RDS")
save(dadaFs,
     dadaRs,
     # derepFs,
     # derepRs,
     file = "dadaFs.RDS")
dadaFs[["F3D0"]] # Summary for first read pair of first read (F3D0)
dadaFs[[1]] # Same as previous line
```

```{r merge, echo=T, results='show', message=FALSE, warning=FALSE}
# Merge paired end sequences
mergers <- 
  mergePairs(dadaFs, 
             derepFs, 
             dadaRs, 
             derepRs, 
             verbose=TRUE) 

head(mergers[[1]]) # Inspect the merger date from first sample
```

```{r seqtab, echo=T, results='show', message=FALSE, warning=FALSE}
seqtab <- makeSequenceTable(mergers) # Makes table with counts for amplicon sequence variants
dim(seqtab) # How many species (amplicon sequence variants) do we see?
table(nchar(getSequences(seqtab))) # Inspect distribution of sequence lengths
```

```{r chimeras, echo=T, results='show', message=FALSE, warning=FALSE}
# Remove chimeras (i.e. r1 from one sequence, r2 from another)
seqtab.nochim <- 
  removeBimeraDenovo(seqtab, 
                     method="consensus", 
                     multithread=TRUE, 
                     verbose=TRUE)
sum(seqtab.nochim)/sum(seqtab)
dim(seqtab.nochim)
```

```{r get_counts, echo=T, results='show', message=FALSE, warning=FALSE}
# Track read counts through pipeline
getN <- function(x) sum(getUniques(x)) # Function to get unique reads
track <- 
  cbind(out, 
        sapply(dadaFs, getN), 
        sapply(dadaRs, getN), 
        sapply(mergers, getN), 
        rowSums(seqtab.nochim)) # Get unique reads for each setp

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

```{r taxons, echo=T, results='show', message=FALSE, warning=FALSE}
# Assign taxonomy
taxa <- assignTaxonomy(seqtab.nochim, 
                       refFasta =
                         "taxonomy/silva_nr_v138_train_set.fa.gz", 
                       multithread=TRUE)
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
```

```{r phyloseq_format, echo=T, results='hide', message=FALSE, warning=FALSE}
# Format data for phyloseq
samples.out <- rownames(seqtab.nochim)
subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
gender <- substr(subject,1,1)
subject <- substr(subject,2,999)
day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
samdf$When <- "Early"
samdf$When[samdf$Day>100] <- "Late"
rownames(samdf) <- samples.out
```

```{r phyloseq_object, echo=T, results='show', message=FALSE, warning=FALSE}
# Construct phyloseq object
ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE), 
               sample_data(samdf), 
               tax_table(taxa))
ps <- prune_samples(sample_names(ps) != "Mock", ps) # Remove mock sample


# Store full DNA sequences and give ASVs short names to more easily visualize
dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- taxa_names(ps)
ps <- merge_phyloseq(ps, dna)
taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))
ps
```

```{r alpha_diversity, echo=T, results='show', message=FALSE, warning=FALSE}
# Alpha diversity
plot_richness(ps, x="Day", measures=c("Shannon", "Simpson"), color="When")
```

```{r ordinate, echo=T, results='hide', message=FALSE, warning=FALSE}
# Transform data to proportions as appropriate for Bray-Curtis distances
# (conceptually similar to PCA)
ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")
plot_ordination(ps.prop, ord.nmds.bray, color="When", title="Bray NMDS")
```

```{r barplot, echo=T, results='show', message=FALSE, warning=FALSE}
top20 <- names(sort(taxa_sums(ps), decreasing=TRUE))[1:20]
ps.top20 <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
plot_bar(ps.top20, x="Day", fill="Family") + facet_wrap(~When, scales="free_x")
```


```{r extra, echo=T, results='show', message=FALSE, warning=FALSE}



results_df2 <- results_df %>%
  filter(!is.na(padj)) %>%
  mutate(Direction = case_when(padj >= 0.05 | abs(log2FoldChange)<=1 ~ "No Change",
                               log2FoldChange < -1 ~ "Down",
                               log2FoldChange > 1 ~ "Up"),
         trans_p = if_else(is.infinite(-log10(padj)),
                                  300,-log10(padj)))

cbbPalette <- c("No Change"="#999999", 
                "A"="#E69F00", 
                "B"="#56B4E9", 
                "Down"="#009E73", 
                "Up"="#F0E442", 
                "Q"="#0072B2", 
                "W"="#D55E00", 
                "E"="#CC79A7")
p <- ggplot(results_df2) + 
  geom_point(aes(x = log2FoldChange,
                 y = trans_p,
                 colour = Direction)) +
  scale_y_continuous(name = "-log10( padj )") +
  scale_colour_manual(values = cbbPalette) +
  theme(legend.position = c(0.2,0.8),
        legend.box.background = element_rect(fill="white",
                                         colour="black"))
# library(cowplot)
# save_plot("logFC.png",
#           p,
#           base_asp = 1.3)
```