---
title: "Dada2_Euka_surface"
author: "Sarah �rberg"
date: "Jan 2021"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r, results='hide'}
#install.packages('ape')
#install.packages('taxize')
#install.packages('ShortRead')

## try http:// if https:// URLs are not supported
#source("https://bioconductor.org/biocLite.R")
#biocLite("dada2")

#library(ShortRead)
#library(ape)
#library(taxize)
#library(dada2)
#library(tidyr)
#library(utils)
#library(dplyr)

```

CHANGE ME to the directory containing the fastq files after unzipping.

```{r}
path <- "U:/temp/Manuskripter/Arctic_surface_sediments/Data/AS_euka02_BR345/AS_euka02_cut" 
list.files(path)
```

Forward and reverse fastq filenames have format: SAMPLENAME_R1_001.fastq and SAMPLENAME_R2_001.fastq

```{r}

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

```

Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq

```{r}

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

```

# Visualize the quality profile of the forward and reverse reads
## 
```{r}
plotQualityProfile(fnFs[1]) #visualize the quality profile of the forward reads
plotQualityProfile(fnRs[1]) #visualize the quality profile of the reverse reads

```

## 
```{r}
plotQualityProfile(fnFs[2]) #visualize the quality profile of the forward reads
plotQualityProfile(fnRs[2]) #visualize the quality profile of the reverse reads

```

## 

```{r}
plotQualityProfile(fnFs[3]) #visualize the quality profile of the forward reads
plotQualityProfile(fnRs[3]) #visualize the quality profile of the reverse reads

```

```{r}
plotQualityProfile(fnFs[51]) #visualize the quality profile of the forward reads
plotQualityProfile(fnRs[51]) #visualize the quality profile of the reverse reads

```


# FILTER SEQUENCES

```{r}
filt_path <- file.path(path, "filtered") # Place filtered files in filtered/ subdirectory
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(90,90),
                    maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=FALSE) # On Windows set multithread=FALSE
head(out)
```

# LEARN THE ERROR RATE

```{r}
errF <- learnErrors(filtFs, multithread=FALSE)
errR <- learnErrors(filtRs, multithread=FALSE)

```

Plot errors from Forward

```{r}
plotErrors(errF, nominalQ=TRUE)

```

Plot errors from Reverse

```{r}
plotErrors(errR, nominalQ=TRUE)
```

# DEREPLICATE

```{r}

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

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaFs[[1]]

dadaRs <- dada(derepRs, err=errR, multithread=TRUE)
dadaRs[[1]]
```

#MERGE Reverse AND Forward

```{r}
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[19]])

#Construct sequence table
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

```

# REMOVE CHIMERAS

```{r}
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)

sum(seqtab.nochim)/sum(seqtab)
```

#TRACK reads through the pipeline

```{r}
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab), rowSums(seqtab.nochim))
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
rownames(track) <- sample.names
head(track)
```

# ASSIGN  - minBoot = 70

## Assign with 

```{r, eval = FALSE}
taxa <- assignTaxonomy(seqtab.nochim, minBoot = 70, "U:/temp/phD/Reference_databases_eDNA/SILVA_132_trimmed_euka02_dada2_names_with_ortega_sarah_ma_jan21.fasta", multithread=TRUE, taxLevels = c("Superkingdom","Kingdom","Phylum","Class","Order","Family","Genus","Species"))

```

```{r, eval = FALSE}
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
```

```{r, eval = FALSE}
head(taxa)
```

seqtab: is the sequence table with rownames=sample numbers and column names=unique sequences
taxonomy: rownames=unique sequences and columns are taxonomy 

transpose seq table and merge to tax table 
```{r, eval = FALSE}
tseq <- t(seqtab.nochim)
```

merge taxonomy and sequence tables because both now have sequences as rownames.

```{r, eval = FALSE}
newtable <- merge(taxa, tseq, by = "row.names")

#Isolate Algal taxa
Macroalgae <- subset(newtable, Phylum == "Chlorophyta" | Phylum == "Rhodophyta" | Class == "Florideophyceae" | Class == "Phaeophyceae" | Class == "Ulvophycea")

#Print csv to file
write.csv2(newtable, file = "U:/temp/Manuskripter/Arctic_surface_sediments/Output/All_samples/assign_tax_mixlib_updatedjan21_mb70.csv")

write.csv2(Macroalgae, file = "U:/temp/Manuskripter/Arctic_surface_sediments/Output/All_samples/assign_tax_mixlib_updatedjan21_macroalgae_mb70.csv")

                       
```
