%\VignetteIndexEntry{VAR-Seq Workflow Template}
%\VignetteDepends{rjson, ggplot2, limma, edgeR, GOstats, GO.db, annotate, pheatmap, VariantAnnotation}
%\VignetteKeywords{compute cluster, pipeline, reports}
%\VignetteEngine{knitr::knitr}
%\VignettePackage{systemPipeR}

% Generate vignette with knitr
% R CMD Sweave --engine=knitr::knitr --pdf systemPipeVARseq.Rnw

\documentclass{article}
%<<style, eval=TRUE, echo=FALSE, results=tex>>=
%BiocStyle::latex(use.unsrturl=FALSE)
%@
<<style-knitr, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex(use.unsrturl=FALSE)
@

\usepackage[authoryear,round]{natbib}
\bibliographystyle{plainnat}
\def\bibsection{\section{References}}

\usepackage{graphicx}
\usepackage{color}
\usepackage{hyperref}
\usepackage{url}
\usepackage{float}

%\newcommand{\comment}[1]{}
%\newcommand{\Rfunction}[1]{{\texttt{#1}}}
%\newcommand{\Robject}[1]{{\texttt{#1}}}
%\newcommand{\Rpackage}[1]{{\textit{#1}}}
%\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
%\newcommand{\Rclass}[1]{{\textit{#1}}}
	
% Define header and footer area with fandyhdr package (see: http://www.ctan.org/tex-archive/macros/latex/contrib/fancyhdr/fancyhdr.pdf)
\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhead{}
\fancyfoot{}
\rhead{\nouppercase{\leftmark}}
\lhead{\textit{systemPipeR VAR-Seq Workflow}}
\rfoot{\thepage}

\begin{document}

<<setup, include=FALSE, cache=FALSE>>=
library(knitr)
# set global chunk options for knitr
opts_chunk$set(comment=NA, warning=FALSE, message=FALSE, fig.path='figure/systemPipeR-')
options(formatR.arrow=TRUE, width=95)
unlink("test.db")
@

\title{VAR-Seq workflow template: Some Descriptive Title}
\author{Project ID: VARseq\_PI\_Name\_Organism\_Jul2015 \\
	Project PI: First Last (first.last@inst.edu)\\
        Author of Report: First Last (first.last@inst.edu)}
\maketitle

\tableofcontents

\section{Introduction}
\subsection{Background and objectives}

This report describes the analysis of a VAR-Seq project studying the genetic differences among several strains ... from \textit{organism} .... 

\subsection{Experimental design}
Typically, users want to specify here all information relevant for the analysis of their NGS study. This includes detailed descriptions of FASTQ files, experimental design, reference genome, gene annotations, etc.  

\section{Load workflow environment}
\subsection{Load packages and sample data}
The \Rpackage{systemPipeR} package needs to be loaded to perform the analysis steps shown in this report \citep{Girke2014-oy}. 
<<eval=TRUE>>=
library(systemPipeR)
@

Load workflow environment with sample data into your current working directory. The sample data are described  \href{http://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#load-sample-data-and-workflow-templates}{\textcolor{blue}{here}}.
<<eval=FALSE>>=
library(systemPipeRdata)
genWorkenvir(workflow="varseq")
setwd("varseq")
@

In the workflow environments generated by \Rfunction{genWorkenvir} all data inputs are stored in a \Robject{data/} directory and all analysis results will be written to a separate \Robject{results/} directory, while the \Robject{systemPipeVARseq.Rnw} script and the \Robject{targets} file are expected to be located in the parent directory. The R session is expected to run from this parent directory. Additional parameter files are stored under \Robject{param/}.

To work with real data, users want to organize their own data similarly and substitute all test data for their own data. To rerun an established workflow on new data, the initial \Robject{targets} file along with the corresponding FASTQ files are usually the only inputs the user needs to provide.

If applicable users can load custom functions not provided by \Rpackage{systemPipeR}. Skip this step if this is not the case.
<<eval=FALSE>>=
source("systemPipeVARseq_Fct.R")
@

\subsection{Experiment definition provided by \Robject{targets} file}
The \href{run:targetsPE.txt}{\Robject{targets}} file defines all FASTQ files and sample comparisons of the analysis workflow. 
<<eval=TRUE>>=
targetspath <- system.file("extdata", "targetsPE.txt", package="systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")[,1:5]
targets
@

\section{Read preprocessing}
\subsection{Read quality filtering and trimming}
The following removes reads with low quality base calls (here Phred scores below 20) from all FASTQ files.
<<eval=FALSE>>=
args <- systemArgs(sysma="param/trimPE.param", mytargets="targetsPE.txt")[1:4] # Note: subsetting!
filterFct <- function(fq, cutoff=20, Nexceptions=0) {
    qcount <- rowSums(as(quality(fq), "matrix") <= cutoff)
    fq[qcount <= Nexceptions] # Retains reads where Phred scores are >= cutoff with N exceptions
}
preprocessReads(args=args, Fct="filterFct(fq, cutoff=20, Nexceptions=0)", batchsize=100000)
writeTargetsout(x=args, file="targets_PEtrim.txt", overwrite=TRUE)
@

\subsection{FASTQ quality report}
The following \Rfunction{seeFastq} and \Rfunction{seeFastqPlot} functions generate and plot a series of
useful quality statistics for a set of FASTQ files including per cycle quality
box plots, base proportions, base-level quality trends, relative k-mer
diversity, length and occurrence distribution of reads, number of reads above
quality cutoffs and mean quality distribution. The results are written to a PDF file named 
\href{run:./results/fastqReport.pdf}{\Robject{fastqReport.pdf}}.

<<eval=FALSE>>=
args <- systemArgs(sysma="bwa.param", mytargets="targets.txt")
fqlist <- seeFastq(fastq=infile1(args), batchsize=100000, klength=8)
pdf("./results/fastqReport.pdf", height=18, width=4*length(fqlist))
seeFastqPlot(fqlist)
dev.off()
@
\begin{figure}[H]
  \centering
   \includegraphics[width=18cm]{fastqReport.pdf}
   \caption{QC report for 18 FASTQ files.}
   \label{fig:fastqreport}
\end{figure}

\section{Alignments}
\subsection{Read mapping with \Rfunction{BWA}}
The NGS reads of this project are aligned against the reference genome sequence using the highly variant tolerant short read aligner \Robject{BWA} \citep{Li2013-oy, Li2009-oc}. The parameter settings of the aligner are defined in the \Robject{bwa.param} file.
<<eval=FALSE>>=
args <- systemArgs(sysma="bwa.param", mytargets="targets.txt")
sysargs(args)[1] # Command-line parameters for first FASTQ file
@

Runs the alignments sequentially (\textit{e.g.} on a single machine)
<<eval=FALSE>>=
bampaths <- runCommandline(args=args)
@

Alternatively, the alignment jobs can be submitted to a compute cluster, here using 72 CPU cores (18 \Robject{qsub} processes each with 4 CPU cores).
<<eval=FALSE>>=
moduleload(modules(args))
system("bwa index -a bwtsw ./data/tair10.fasta")
resources <- list(walltime="20:00:00", nodes=paste0("1:ppn=", cores(args)), memory="10gb")
reg <- clusterRun(args, conffile=".BatchJobs.R", template="torque.tmpl", Njobs=18, runid="01", 
                  resourceList=resources)
waitForJobs(reg)
@
Check whether all BAM files have been created
<<eval=FALSE>>=
file.exists(outpaths(args))
@

\subsection{Read mapping with \Rfunction{gsnap}}
An alternative variant tolerant aligner is \Rfunction{gsnap} from the \Rpackage{gmapR} package \citep{Wu2010-iq}. The following code shows how to run this aligner on multiple nodes of a comute cluster that uses Torque as scheduler.
<<eval=FALSE>>=
library(gmapR); library(BiocParallel); library(BatchJobs)
gmapGenome <- GmapGenome(reference(args), directory="data", name="gmap_tair10chr", create=TRUE)
args <- systemArgs(sysma="gsnap.param", mytargets="targetsPE.txt")
f <- function(x) {
    library(gmapR); library(systemPipeR)
    args <- systemArgs(sysma="gsnap.param", mytargets="targetsPE.txt")
    gmapGenome <- GmapGenome(reference(args), directory="data", name="gmap_tair10chr", create=FALSE)
    p <- GsnapParam(genome=gmapGenome, unique_only=TRUE, molecule="DNA", max_mismatches=3)
    o <- gsnap(input_a=infile1(args)[x], input_b=infile2(args)[x], params=p, output=outfile1(args)[x])
}
funs <- makeClusterFunctionsTorque("torque.tmpl")
param <- BatchJobsParam(length(args), resources=list(walltime="20:00:00", nodes="1:ppn=1", memory="6gb"), cluster.functions=funs)
register(param)
d <- bplapply(seq(along=args), f)
writeTargetsout(x=args, file="targets_gsnap_bam.txt")
@
\subsection{Read and alignment stats}
The following generates a summary table of the number of reads in each sample and how many of them aligned to the reference.
<<eval=FALSE>>=
read_statsDF <- alignStats(args=args) 
write.table(read_statsDF, "results/alignStats.xls", row.names=FALSE, quote=FALSE, sep="\t")
@

\subsection{Create symbolic links for viewing BAM files in IGV}
The \Rfunction{symLink2bam} function creates symbolic links to view the BAM alignment files in a genome browser such as IGV. The corresponding URLs are written to a file with a path specified under \Robject{urlfile}, here \href{run:./results/IGVurl.txt}{IGVurl.txt}.
<<eval=FALSE>>=
symLink2bam(sysargs=args, htmldir=c("~/.html/", "somedir/"), 
            urlbase="http://biocluster.ucr.edu/~tgirke/", 
	    urlfile="./results/IGVurl.txt")
@

\section{Variant calling}
The following performs variant calling with \Robject{GATK}, \Robject{BCFtools} and \Rpackage{VariantTools} in parallel mode on a compute cluster \citep{McKenna2010-ql, Li2011-ll}. If a cluster is not available, the \Rfunction{runCommandline()} function can be used to run the variant calling with \Robject{GATK} and \Robject{BCFtools} for each sample sequentially on a single machine, or \Rfunction{callVariants} in case of \Rpackage{VariantTools}. Typically, the user would choose here only one variant caller rather than running several ones.   
\subsection{Variant calling with GATK}
The following creates in the inital step a new targets file (\Robject{targets\_bam.txt}). The first column of this file gives the paths to the BAM files created in the alignment step. The new targets file and the parameter file \Robject{gatk.param} are used to create a new \Robject{SYSargs} instance for running GATK. Since GATK involves many processing steps, it is executed by a bash script \Robject{gatk\_run.sh} where the user can specify the detailed run parameters. All three files are expected to be located in the current working directory. Samples files for \Robject{gatk.param} and \Robject{gatk\_run.sh} are available in the subdirectory \Robject{./inst/extdata/} of the source file of the \Rpackage{systemPipeR} package. Alternatively, they can be downloaded directly from \href{https://github.com/tgirke/systemPipeR/tree/master/inst/extdata}{here}.
<<eval=FALSE>>=
writeTargetsout(x=args, file="targets_bam.txt")
system("java -jar CreateSequenceDictionary.jar R=./data/tair10.fasta O=./data/tair10.dict")
# system("java -jar /opt/picard/1.81/CreateSequenceDictionary.jar R=./data/tair10.fasta O=./data/tair10.dict")
args <- systemArgs(sysma="gatk.param", mytargets="targets_bam.txt")
resources <- list(walltime="20:00:00", nodes=paste0("1:ppn=", 1), memory="10gb")
reg <- clusterRun(args, conffile=".BatchJobs.R", template="torque.tmpl", Njobs=18, runid="01",
                  resourceList=resources)
waitForJobs(reg)
writeTargetsout(x=args, file="targets_gatk.txt")
@

\subsection{Variant calling with BCFtools}
The following runs the variant calling with BCFtools. This step requires in the current working directory the parameter file \Robject{sambcf.param} and the bash script \Robject{sambcf\_run.sh}.  
<<eval=FALSE>>=
args <- systemArgs(sysma="sambcf.param", mytargets="targets_bam.txt")
resources <- list(walltime="20:00:00", nodes=paste0("1:ppn=", 1), memory="10gb")
reg <- clusterRun(args, conffile=".BatchJobs.R", template="torque.tmpl", Njobs=18, runid="01",
                  resourceList=resources)
waitForJobs(reg)
writeTargetsout(x=args, file="targets_sambcf.txt")
@

\subsection{Variant calling with \Rpackage{VariantTools}}
<<eval=FALSE>>=
library(gmapR); library(BiocParallel); library(BatchJobs)
args <- systemArgs(sysma="vartools.param", mytargets="targets_gsnap_bam.txt")
f <- function(x) {
    library(VariantTools); library(gmapR); library(systemPipeR)
    args <- systemArgs(sysma="vartools.param", mytargets="targets_gsnap_bam.txt")
    gmapGenome <- GmapGenome(systemPipeR::reference(args), directory="data", name="gmap_tair10chr", create=FALSE)
    tally.param <- TallyVariantsParam(gmapGenome, high_base_quality = 23L, indels = TRUE)
    bfl <- BamFileList(infile1(args)[x], index=character())
    var <- callVariants(bfl[[1]], tally.param)
    sampleNames(var) <- names(bfl)
    writeVcf(asVCF(var), outfile1(args)[x], index = TRUE)
}
funs <- makeClusterFunctionsTorque("torque.tmpl")
param <- BatchJobsParam(length(args), resources=list(walltime="20:00:00", nodes="1:ppn=1", memory="6gb"), cluster.functions=funs)
register(param)
d <- bplapply(seq(along=args), f)
writeTargetsout(x=args, file="targets_vartools.txt")
@

\section{Filter variants}
The function \Rfunction{filterVars} filters VCF files based on user definable quality parameters. It sequentially imports each VCF file into R, applies the filtering on an internally generated \Robject{VRanges} object and then writes the results to a new subsetted VCF file. The filter parameters are passed on to the corresponding argument as a character string. The function applies this filter to the internally generated \Robject{VRanges} object using the standard subsetting syntax for two dimensional objects such as: \Robject{vr[filter, ]}. The parameter files (\Robject{filter\_gatk.param}, \Robject{filter\_sambcf.param} and \Robject{filter\_vartools.param}), used in the filtering steps, define the paths to the input and output VCF files which are stored in new \Robject{SYSargs} instances.  

\subsection{Filter variants called by \Rpackage{GATK}}
The below example filters for variants that are supported by $\ge$\texttt{x} reads and $\ge$80\% of them support the called variants. In addition, all variants need to pass $\ge$\texttt{x} of the soft filters recorded in the VCF files generated by GATK. Since the toy data used for this workflow is very small, the chosen settings are unreasonabley relaxed. A more reasonable filter setting is given in the line below (here commented out).  
<<eval=FALSE>>=
library(VariantAnnotation)
args <- systemArgs(sysma="filter_gatk.param", mytargets="targets_gatk.txt")
filter <- "totalDepth(vr) >= 2 & (altDepth(vr) / totalDepth(vr) >= 0.8) & rowSums(softFilterMatrix(vr))>=1"
# filter <- "totalDepth(vr) >= 20 & (altDepth(vr) / totalDepth(vr) >= 0.8) & rowSums(softFilterMatrix(vr))==6"
filterVars(args, filter, varcaller="gatk", organism="A. thaliana")
writeTargetsout(x=args, file="targets_gatk_filtered.txt")
@

\subsection{Filter variants called by \Rpackage{BCFtools}}
The following shows how to filter the VCF files generated by \Rpackage{BCFtools} using similar parameter settings as in the previous filtering of the GATK results.
<<eval=FALSE>>=
args <- systemArgs(sysma="filter_sambcf.param", mytargets="targets_sambcf.txt")
filter <- "rowSums(vr) >= 2 & (rowSums(vr[,3:4])/rowSums(vr[,1:4]) >= 0.8)"
# filter <- "rowSums(vr) >= 20 & (rowSums(vr[,3:4])/rowSums(vr[,1:4]) >= 0.8)"
filterVars(args, filter, varcaller="bcftools", organism="A. thaliana")
writeTargetsout(x=args, file="targets_sambcf_filtered.txt")
@

\subsection{Filter variants called by \Rpackage{VariantTools}}
The following shows how to filter the VCF files generated by \Rpackage{VariantTools} using similar parameter settings as in the previous filtering of the GATK results.
<<eval=FALSE>>=
args <- systemArgs(sysma="filter_vartools.param", mytargets="targets_vartools.txt")
filter <- "(values(vr)$n.read.pos.ref + values(vr)$n.read.pos) >= 2 & (values(vr)$n.read.pos / (values(vr)$n.read.pos.ref + values(vr)$n.read.pos) >= 0.8)"
# filter <- "(values(vr)$n.read.pos.ref + values(vr)$n.read.pos) >= 20 & (values(vr)$n.read.pos / (values(vr)$n.read.pos.ref + values(vr)$n.read.pos) >= 0.8)"
filterVars(args, filter, varcaller="vartools", organism="A. thaliana")
writeTargetsout(x=args, file="targets_vartools_filtered.txt")
@

\section{Annotate filtered variants}
The function \Rfunction{variantReport} generates a variant report using utilities provided by the \Rpackage{VariantAnnotation} package. The report for each sample is written to a tabular file containing genomic context annotations (\textit{e.g.} coding or non-coding SNPs, amino acid changes, IDs of affected genes, etc.) along with confidence statistics for each variant. The parameter file \Robject{annotate\_vars.param} defines the paths to the input and output files which are stored in a new \Robject{SYSargs} instance. 

\subsection{Annotate filtered variants called by GATK}

<<eval=FALSE>>=
library("GenomicFeatures")
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_gatk_filtered.txt")
txdb <- loadDb("./data/tair10.sqlite")
fa <- FaFile(systemPipeR::reference(args))
variantReport(args=args, txdb=txdb, fa=fa, organism="A. thaliana")
@

\subsection{Annotate filtered variants called by BCFtools}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_sambcf_filtered.txt")
txdb <- loadDb("./data/tair10.sqlite")
fa <- FaFile(systemPipeR::reference(args))
variantReport(args=args, txdb=txdb, fa=fa, organism="A. thaliana")
@

\subsection{Annotate filtered variants called by \Rpackage{VariantTools}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_vartools_filtered.txt")
txdb <- loadDb("./data/tair10.sqlite")
fa <- FaFile(systemPipeR::reference(args))
variantReport(args=args, txdb=txdb, fa=fa, organism="A. thaliana")
@

\section{Combine annotation results among samples}
To simplify comparisons among samples, the \Rfunction{combineVarReports} function combines all variant annotation reports referenced in a \Robject{SYSargs} instance (here \Robject{args}). At the same time the function allows to consider only certain feature types of interest. For instance, the below setting \Rfunarg{filtercol=c(Consequence="nonsynonymous")} will include only nonsysynonymous variances listed in the \Rfunarg{Consequence} column of the annotation reports. To omit filtering, one can use the setting \Rfunarg{filtercol="All"}

\subsection{Combine results from \Rpackage{GATK}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_gatk_filtered.txt")
combineDF <- combineVarReports(args, filtercol=c(Consequence="nonsynonymous"))
write.table(combineDF, "./results/combineDF_nonsyn_gatk.xls", quote=FALSE, row.names=FALSE, sep="\t")
@

\subsection{Combine results from \Rpackage{BCFtools}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_sambcf_filtered.txt")
combineDF <- combineVarReports(args, filtercol=c(Consequence="nonsynonymous"))
write.table(combineDF, "./results/combineDF_nonsyn_sambcf.xls", quote=FALSE, row.names=FALSE, sep="\t")
@

\subsection{Combine results from \Rpackage{VariantTools}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_vartools_filtered.txt")
combineDF <- combineVarReports(args, filtercol=c(Consequence="nonsynonymous"))
write.table(combineDF, "./results/combineDF_nonsyn_vartools.xls", quote=FALSE, row.names=FALSE, sep="\t")
@

\section{Summary statistics of variants}
The function \Rfunction{varSummar} counts the number of variants for each feature type included in the anntation reports. 
\subsection{Summary for \Rpackage{GATK}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_gatk_filtered.txt")
write.table(varSummary(args), "./results/variantStats_gatk.xls", quote=FALSE, col.names = NA, sep="\t")
@

\subsection{Summary for \Rpackage{BCFtools}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_sambcf_filtered.txt")
write.table(varSummary(args), "./results/variantStats_sambcf.xls", quote=FALSE, col.names = NA, sep="\t")
@

\subsection{Summary for \Rpackage{VariantTools}}
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_vartools_filtered.txt")
write.table(varSummary(args), "./results/variantStats_vartools.xls", quote=FALSE, col.names = NA, sep="\t")
@

\section{Venn diagram of variants}
The venn diagram utilities defined by the \Rpackage{systemPipeR} package can be used to identify common and unique variants reported for different samples and/or variant callers. The below generates a 4-way venn diagram comparing four sampes for each of the two variant callers.
<<eval=FALSE>>=
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_gatk_filtered.txt")
varlist <- sapply(names(outpaths(args))[1:4], function(x) as.character(read.delim(outpaths(args)[x])$VARID))
vennset_gatk <- overLapper(varlist, type="vennsets")
args <- systemArgs(sysma="annotate_vars.param", mytargets="targets_sambcf_filtered.txt")
varlist <- sapply(names(outpaths(args))[1:4], function(x) as.character(read.delim(outpaths(args)[x])$VARID))
vennset_bcf <- overLapper(varlist, type="vennsets")
pdf("./results/vennplot_var.pdf")
vennPlot(list(vennset_gatk, vennset_bcf), mymain="", mysub="GATK: red; BCFtools: blue", colmode=2, ccol=c("blue", "red"))
dev.off()
@

\begin{figure}[H]
  \centering
   \includegraphics[width=10cm]{vennplot_var.pdf}
   \caption{Venn Diagram for 4 samples from GATK and BCFtools.}
   \label{fig:vennplot}
\end{figure}

\section{Version Information}
<<sessionInfo, results='asis'>>=
toLatex(sessionInfo())
@

\section{Funding}
This project was supported by funds from the National Institutes of Health (NIH).

\bibliography{bibtex}

\end{document}