long read Tutorial.Rmd

---
title: "Long Read Tutorial"
author: "Jorden Rabasco"
date: "5/19/2022"
output: html_document
---

```{r, include = FALSE}
library(knitr)
hook_output <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
  lines <- options$output.lines
  if (is.null(lines)) {
    return(hook_output(x, options))  # pass to default hook
  }
  x <- unlist(strsplit(x, "\n"))
  more <- "..."
  if (length(lines)==1) {        # first n lines
    if (length(x) > lines) {
      # truncate the output, but add ....
      x <- c(head(x, lines), more)
    }
  } else {
    x <- c(more, x[lines], more)
  }
  # paste these lines together
  x <- paste(c(x, ""), collapse = "\n")
  hook_output(x, options)
})
```

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

This tutorial is intended for those that wish to analyze their long read sequences in either the form of Pacbio reads or Loop seq reads. The output from this workflow will consist of assigned taxonomies of your ASVs generated by dada2. This tutorial also assumes a working version of R 4.2.0
Lets begin. If you first don’t have your own data you will want to get some sample data. This workflow will work off the Zymo mock community data set for both Pacbio reads and Loop seq reads which can be found here. It is worth noting that Pacbio reads and Loopseq reads are separate and you should not mix and match those discrete commands intended to analyze specifically Pacbio or Loopseq reads:
<https://www.zymoresearch.com/zymobiomics-community-standard>


##  Initial Set up

In the initial set up we will have to prep the reads for filtering and primer trimming. To do so, we will need to load libraries in R, but prior to this if Pacbio reads are being used you will need to run this code at the command line:

ccs --pbi --force --logLevel=DEBUG --numThreads=16 --minSnr=3.75 --minReadScore=0.65 --maxLength=7000 --minLength=10 --minPasses=3 --minZScore=-5 --maxDropFraction=0.34 --minPredictedAccuracy=0.999 subreads.bam ccs.bam


Prior to any R processing please make sure to install dada2 appropriately as well as the dependencies. This information can be found here: <https://benjjneb.github.io/dada2/dada-installation.html>

```{r, message = FALSE, warning = FALSE, eval=FALSE, include = FALSE}
install.packages("devtools", repos='http://cran.us.r-project.org', dependencies = TRUE)
library("devtools")
devtools::install_github("benjjneb/dada2", ref="v1.20")
```

Install all packages as described previously then load all appropriate R libraries. 
```{r, message = FALSE, warning = FALSE}
library(dada2);packageVersion("dada2")
library(Biostrings);packageVersion("Biostrings")
library(ShortRead);packageVersion("ShortRead")
library(ggplot2);packageVersion("ggplot2")
library(reshape2);packageVersion("reshape2")
library(RColorBrewer);packageVersion("RColorBrewer")

library(BiocGenerics);packageVersion("BiocGenerics")
library(S4Vectors);packageVersion("S4Vectors")
library(Biostrings);packageVersion("Biostrings")
library(Biobase);packageVersion("Biobase")
library(MatrixGenerics);packageVersion("MatrixGenerics")
```


## Pathing setup

This section will set up the paths and generate the appropriate folders in the desired working directory indicated by the variable "base_path" You will need to change this to where you want the output to be generated. You will also need to change the variable "fn" to the location of your input data.
```{r Pathing Set-up, message = FALSE, warning = FALSE}

existance_check <- function(pathway) {
    if(!file.exists(pathway)){
      dir.create(pathway)
    }
}
base_path <- "C:/Users/Jorden/Desktop/longread_tutorial" #output generation location
pac_path <- file.path(base_path, "Pacbio")
existance_check(pac_path)
pac_path_out <- file.path(pac_path, "Figures")
existance_check(pac_path_out )
pac_path_rds <- file.path(pac_path, "RDS")
existance_check(pac_path_rds)

loop_path <- file.path(base_path, "Loopseq")
existance_check(loop_path)
loop_path_out <- file.path(loop_path, "Figures")
existance_check(loop_path_out)
loop_path_rds <- file.path(loop_path, "RDS")
existance_check(loop_path_rds)

fn<-"C:/Users/Jorden/Desktop/longread_tutorial/zymo_CCS_99_9.fastq.gz" #sets locaiton of input data
```


## Identification of Primers

### Pacbio
Set up the primer information for Pacbio reads
```{r Pacbio 1) Pathing Primer Identification, message = FALSE, warning = FALSE}
#path <- "C:/Users/jorde/OneDrive/Desktop/Long_read_tutorial/Pacbio" # CHANGE ME to location of the fastq file
#fn <- file.path(path, "zymo_CCS_99_9.fastq.gz")
pacbio_Fw <- "AGRGTTYGATYMTGGCTCAG"
pacbio_rev <- "RGYTACCTTGTTACGACTT"
rc <- dada2:::rc
theme_set(theme_bw())
genusPalette <- c(Bacillus="#e41a1c", Enterococcus="#377eb8", Escherichia="#4daf4a", Lactobacillus="#984ea3",
                  Listeria="#ff7f00", Pseudomonas="#ffff33", Salmonella="#a65628", Staphylococcus="#f781bf")
```
### Loop-seq
If Loop-seq reads are being analyzed this step will set up the appropriate primers, as  we would if we were analyzing pacbio data, but with the added step of primer orientation checking. This step is used as a sanity check to make sure that primers are in their expected orientation in the vast majority of reads. 
```{r Loop-seq 1) Pathing Identification and Orientation of Primers, message = FALSE, warning = FALSE}
loop_FWD <- "AGAGTTTGATCMTGGC" # Loop 16S forward primer
loop_REV <- "TACCTTGTTACGACTT" # Loop 16S reverse primer
allOrients <- function(primer) {
    # Create all orientations of the input sequence
    require(Biostrings)
    dna <- DNAString(primer)
    orients <- c(Forward = dna, Complement = complement(dna), Reverse = reverse(dna), 
        RevComp = reverseComplement(dna))
    return(sapply(orients, toString))  # Convert back to character vector
}
primerHits <- function(primer, fn) {
    # Counts number of reads in which the primer is found
    nhits <- vcountPattern(primer, sread(readFastq(fn)), fixed = FALSE)
    return(sum(nhits > 0))
}
rbind(loop_FWD.Primer = sapply(allOrients(loop_FWD), primerHits, fn = fn), 
      loop_REV.Primer = sapply(allOrients(loop_REV), primerHits, fn = fn))
```


## Filtering and Primer Removal

Primer removal and filtering are steps used to keep the quality of the reads high and to assist in downstream taxonomic assignments. 
### Pacbio
The primer removal and read filtering steps designed for Pacbio reads are as follows:
```{r Pacbio 2) Filtering and Primer Removal Pacbio, message = FALSE, warning = FALSE}
nop <- file.path(pac_path, "noprimers", basename(fn))
prim <- removePrimers(fn, nop, primer.fwd=pacbio_Fw, primer.rev=dada2:::rc(pacbio_rev), orient=TRUE,  verbose=TRUE)
filtered_storage<-file.path(pac_path,"filtered")
existance_check(filtered_storage)
filt <- file.path(pac_path,"filtered", basename(fn))
track <- fastqFilter(nop, filt, minQ=3, minLen=1000, maxLen=1600, maxN=0, rm.phix=FALSE, maxEE=2, verbose=TRUE)
```
### Loop-seq
The primer removal and read filtering steps designed for Loop-seq reads are as follows:
```{r Loop-seq 2) Filtering and Primer Removal Loop-seq, message = FALSE, warning = FALSE}
loop_nop <- file.path(loop_path, "nop", basename(fn))
loop_out <- removePrimers(fn, loop_nop, loop_FWD, rc(loop_REV), verbose=TRUE)
loop_filt <- file.path(loop_path, "filtered", basename(fn))
loop_track <- filterAndTrim(loop_nop, loop_filt, maxEE=2, minLen=1400, maxLen=1600, verbose=TRUE)
```


## Scan for contaminants in filtered reads

While the previous step was shared by both Pacbio and Loop-seq workflows, this next step is specific to Loop-seq reads. We’ll do a simple scan for potential contaminants in the filtered unique sequences by assigning them taxonomy, and attributing any sequences assigned to genera outside those expected as potential contaminants. To do this we will need a supporting file which can be found here: <https://zenodo.org/record/3986799#.YrnLqnbMIuU> or you can use the following link address to download the file directly via the command line; <https://zenodo.org/record/3986799/files/silva_nr99_v138_train_set.fa.gz?download=1> As an aside, if the file "tax_filt.rds" has not been generated on previous run throughs then this step will take a signifigant amount of time, and is nothing to worry about. 
### Loop-seq
```{r Loop-seq 3) Scan for contaminants in filtered reads(pt.1), message = FALSE, warning = FALSE}
train_data_loc_loop<-"C:/Users/jorde/OneDrive/Desktop/Long_read_tutorial/tax/silva_nr99_v138_train_set.fa.gz" #change this address to wherever the file exists on your system
drp.filt <- derepFastq(loop_filt)
tax.filt.rds <- file.path(loop_path_rds, "tax_filt.rds") # RDS save/load to speed up reruns of the code
if(!file.exists(tax.filt.rds)) {
  tax.filt <- assignTaxonomy(drp.filt, train_data_loc_loop, minBoot=80, multi=TRUE)
  saveRDS(tax.filt, tax.filt.rds)
}
tax.filt <- readRDS(tax.filt.rds)
if(!identical(getSequences(tax.filt), getSequences(drp.filt))) stop("Taxonomy mismatch.")
table(tax.filt[,"Genus"], useNA="ifany")

```
As a follow up We can then look up any unexpected NAs with the following code:
### Loop-seq
```{r Loop-seq 4) Scan for contaminants in filtered reads(pt.2), output.lines=6, message = FALSE, warning = FALSE}
dada2:::pfasta(getSequences(drp.filt)[is.na(tax.filt[,"Genus"])])
```


## Denoising

The next step is arguable the most important part of the tutorial; denoising via dada2 package. This was already loaded into the R environment in the beginning of the workflow so there is no need to worry. This will produce error plots as well as an ASV table which we can then use for taxonomic assingment as well as other helpful information. 
### Pacbio
The denoising step designed for Pacbio data is as follows:
```{r Pacbio 3) Denoising Pacbio, message = FALSE, warning = FALSE}
drp <- derepFastq(filt, verbose=TRUE)
err <- learnErrors(drp, BAND_SIZE=32, multithread=TRUE, errorEstimationFunction=dada2:::PacBioErrfun) # 10s of seconds
plotErrors(err)
dd <- dada(drp, err=err, BAND_SIZE=32, multithread=TRUE) # seconds
cbind(ccs=prim[,1], primers=prim[,2], filtered=track[[2]], denoised=sum(dd$denoised))
dd$clustering[,-1]
```
### Loop-seq
The denoising step designed for loop-seq data is as follows:
```{r Loop-seq 5) Denoising Loop-seq, message = FALSE, warning = FALSE}
loop_err <- learnErrors(loop_filt, multi=TRUE, verbose=0)
plotErrors(loop_err, nominalQ=TRUE)
loop_dd <- dada(loop_filt, loop_err, multi=TRUE, OMEGA_C=0)
loop_dd
```


## Assigning Taxonomy

### Pacbio
For Pacbio data; to assign the ASV output of the dada2 step a taxonomic assingment we will need to download a supporting file. This is used as the taxonomic informaiton that will be used in this step. The download for the file can be found here <https://zenodo.org/record/801832#.YrnEgnbMIuV> or you can use the following link address to download the file directly via the command line; <https://zenodo.org/record/801832/files/silva_nr_v128_train_set.fa.gz?download=1>
```{r Pacbio 4) Assigning Taxonomy Pacbio (pt.1), message = FALSE, warning = FALSE}
train_data <- "C:/Users/Jorden/Desktop/longread_tutorial/tax/silva_nr_v128_train_set.fa.gz" #change this address to wherever the file exists on your system
tax <- assignTaxonomy(dd, train_data, multithread=TRUE)
tax[,"Genus"] <- gsub("Escherichia/Shigella", "Escherichia", tax[,"Genus"]) # Reformat to be compatible with other data sources
unname(tax)
```
###Pacbio
we can also use the following code to see if any chimeras are present. 
```{r Pacbio 5) Assigning Taxonomy Pacbio (pt.2), message = FALSE, warning = FALSE}
bim <- isBimeraDenovo(dd, minFoldParentOverAbundance=3.5) 
# Higher MFPOA to avoid flagging intra-genomic variants
table(bim)
```
As you can see there are 29!
### Loop-seq
The taxonomic assinment stPaep for loop-seq is signifigantly smaller than for pacBio data as much of the processing was do in Loop-seq step 2.The last piece to needed for the assingment was the ASV output of the dada2 package which we have just accomplished. Therefore to finish up this tutorial the steps desinged for loop-seq taxonomic analysis are included below.
```{r Loop-seq 6) Assigning Taxonomy loop-seq, message = FALSE, warning = FALSE}
tax <- tax.filt[getSequences(loop_dd),]
if(!identical(getSequences(tax), getSequences(loop_dd))) stop("Taxonomy mismatch.")
table(tax[,"Genus"], useNA="ifany")
```
Note: This chimera check step for loop-seq reads was incorporated into the “Scan for contaminants in filtered reads” step when you blast the NAs in the dataframe to determine their origin.