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

% Generate vignette with knitr
% R CMD Sweave --engine=knitr::knitr --pdf systemPipeChIPseq.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 ChIP-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{ChIP-Seq workflow template: Some Descriptive Title}
\author{Project ID: ChIPseq\_PI\_Name\_Organism\_Jun2015 \\
	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 several ChIP-Seq experiments studying the DNA binding patterns of the transcriptions factors ... 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="chipseq")
setwd("chipseq")
@

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{systemPipeChIPseq.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("systemPipeChIPseq_Fct.R")
@

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

\section{Read preprocessing}
\subsection{Read quality filtering and trimming}
The following example shows how one can design a custom read preprocessing function using utilities provided by the \Rpackage{ShortRead} package, and then apply it with \Rfunction{preprocessReads} in batch mode to all FASTQ samples referenced in the corresponding \Robject{SYSargs} instance (\Robject{args} object below). More detailed information on read preprocessing is provided in \Rpackage{systemPipeR}'s main vignette.
<<eval=FALSE, messages=FALSE, warning=FALSE, cache=TRUE>>=
args <- systemArgs(sysma="param/trim.param", mytargets="targets_chip.txt")
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_chip_trim.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="param/bowtieSE.param", mytargets="targets_chip_trim.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{Bowtie2}}
The NGS reads of this project will be aligned with \Robject{Bowtie2} against the reference genome sequence \citep{Langmead2012-bs}. The parameter settings of the aligner are defined in the \Robject{bowtieSE.param} file. In ChIP-Seq experiments it is usually more appropriate to eliminate reads mapping to multiple locations. To achieve this, users wants to remove the argument setting '\Rfunarg{-k 50 --non-deterministic}' in the \Robject{bowtieSE.param} file.
<<eval=FALSE>>=
args <- systemArgs(sysma="param/bowtieSE.param", mytargets="targets_chip_trim.txt")
sysargs(args)[1] # Command-line parameters for first FASTQ file
moduleload(modules(args)) # Skip if a module system is not used
system("bowtie2-build ./data/tair10.fasta ./data/tair10.fasta") # Indexes reference genome
runCommandline(args)
writeTargetsout(x=args, file="targets_bam.txt", overwrite=TRUE)
@
Check whether all BAM files have been created
<<eval=FALSE>>=
file.exists(outpaths(args))
@

\subsection{Read and alignment stats}
The following provides an overview 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")
read.delim("results/alignStats.xls")
@

\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 without moving these large files to a local system. 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{Peak calling with MACS2} 
\subsection{Merge BAM files of replicates prior to peak calling}
Merging BAM files of technical and/or biological replicates can improve the sensitivity of the peak calling by increasing the depth of read coverage. The \Rfunction{mergeBamByFactor} function merges BAM files based on grouping information specified by a \Rfunction{factor}, here the \Rfunarg{Factor} column of the imported targets file. It also returns an updated \Rfunarg{SYSargs} object containing the paths to the merged BAM files as well as to any unmerged files without replicates. This step can be skipped if merging of BAM files is not desired. 
<<eval=FALSE>>=
args <- systemArgs(sysma=NULL, mytargets="targets_bam.txt")
args_merge <- mergeBamByFactor(args, overwrite=TRUE)
writeTargetsout(x=args_merge, file="targets_mergeBamByFactor.txt", overwrite=TRUE)
@

\subsection{Peak calling without input/reference sample}
MACS2 can perform peak calling on ChIP-Seq data with and without input samples \citep{Zhang2008-pc}. The following performs peak calling without input on all samples specified in the corresponding \Robject{args} object. Note, due to the small size of the sample data, MACS2 needs to be run here with the '\Rfunarg{--nomodel}' setting. For real data sets, users want to remove this parameter in the corresponding \Robject{*.param} file(s).
<<eval=FALSE>>=
args <- systemArgs(sysma="param/macs2_noinput.param", mytargets="targets_mergeBamByFactor.txt")
sysargs(args)[1] # Command-line parameters for first FASTQ file
runCommandline(args)
file.exists(outpaths(args))
writeTargetsout(x=args, file="targets_macs.txt", overwrite=TRUE)
@

\subsection{Peak calling with input/reference sample}
To perform peak calling with input samples, they can be most conveniently specified in the \Rfunarg{SampleReference} column of the initial targets file. The function \Rfunction{writeTargetsRef} uses this information to create a \Robject{targets} file intermediate for running MACS2 with the corresponding input samples.  
<<eval=FALSE>>=
writeTargetsRef(infile="targets_mergeBamByFactor.txt", outfile="targets_bam_ref.txt", silent=FALSE, overwrite=TRUE)
args <- systemArgs(sysma="param/macs2.param", mytargets="targets_bam_ref.txt")
sysargs(args)[1] # Command-line parameters for first FASTQ file
runCommandline(args)
file.exists(outpaths(args))
writeTargetsout(x=args, file="targets_macs.txt", overwrite=TRUE)
@

The peak calling results from MACS2 are written for each sample to separate files in the \Robject{results} directory. They are named after the corresponding \Robject{*.bam} files with extensions used by MACS2.

\section{Annotate peaks with genomic context}
\subsection{Annotation with \Rpackage{ChIPpeakAnno} package}
The following annotates the identified peaks with genomic context information using the \Rpackage{ChIPpeakAnno} and \Rpackage{ChIPseeker} packages, respectively \citep{Zhu2010-zo, Yu2015-xu}.
<<eval=FALSE>>=
library(ChIPpeakAnno); library(GenomicFeatures)
args <- systemArgs(sysma="param/annotate_peaks.param", mytargets="targets_macs.txt")
txdb <- loadDb("./data/tair10.sqlite")
ge <- genes(txdb, columns=c("tx_name", "gene_id", "tx_type")) 
for(i in seq(along=args)) {
    peaksGR <- as(read.delim(infile1(args)[i], comment="#"), "GRanges")
    annotatedPeak <- annotatePeakInBatch(peaksGR, AnnotationData=genes(txdb))
    df <- data.frame(as.data.frame(annotatedPeak), as.data.frame(values(ge[values(annotatedPeak)$feature,])))
    write.table(df, outpaths(args[i]), quote=FALSE, row.names=FALSE, sep="\t")
}
writeTargetsout(x=args, file="targets_peakanno.txt", overwrite=TRUE)
@
<<eval=FALSE, include=FALSE>>=
## Perform previous step with full genome annotation from Biomart
# txdb <- makeTxDbFromBiomart(biomart="ENSEMBL_MART_PLANT", dataset="athaliana_eg_gene")
# tx <- transcripts(txdb, columns=c("tx_name", "gene_id", "tx_type")) 
# ge <- genes(txdb, columns=c("tx_name", "gene_id", "tx_type")) # works as well 
# seqlevels(ge) <- c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "ChrC", "ChrM")
# table(mcols(tx)$tx_type)
# tx <- tx[!duplicated(unstrsplit(values(tx)$gene_id, sep=",")) # Keeps only first transcript model for each gene]
# annotatedPeak <- annotatePeakInBatch(macsOutput, AnnotationData = tx)
@

The peak annotation results are written for each peak set to separate files in the \Robject{results} directory. They are named after the corresponding peak files with extensions specified in the \Robject{annotate\_peaks.param} file, here '\Robject{*.peaks.annotated.xls}'.

\subsection{Annotation with \Rpackage{ChIPseeker} package}
Same as in previous step but using the \Rpackage{ChIPseeker} package for annotating the peaks.
<<eval=FALSE>>=
library(ChIPseeker)
txdb <- loadDb("./data/tair10.sqlite")
for(i in seq(along=args)) {
    peakAnno <- annotatePeak(infile1(args)[i], TxDb=txdb, verbose=FALSE)
    df <- as.data.frame(peakAnno)
    write.table(df, outpaths(args[i]), quote=FALSE, row.names=FALSE, sep="\t")
}
writeTargetsout(x=args, file="targets_peakanno.txt", overwrite=TRUE)
@

Summary plots provided by the \Rpackage{ChIPseeker} package. Here applied only to one sample for demonstration purposes.
<<eval=FALSE>>=
peak <- readPeakFile(infile1(args)[1])
covplot(peak, weightCol="X.log10.pvalue.")
peakHeatmap(outpaths(args)[1], TxDb=txdb, upstream=1000, downstream=1000, color="red")
plotAvgProf2(outpaths(args)[1], TxDb=txdb, upstream=1000, downstream=1000, xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
@

\section{Count reads overlapping peak regions}
The \Rfunction{countRangeset} function is a convenience wrapper to perform read counting iteratively over serveral range sets, here peak range sets. Internally, the read counting is performed with the \Rfunction{summarizeOverlaps} function from the \Rpackage{GenomicAlignments} package. The resulting count tables are directly saved to files, one for each peak set.

<<eval=FALSE>>=
library(GenomicRanges)
args <- systemArgs(sysma="param/count_rangesets.param", mytargets="targets_macs.txt")
args_bam <- systemArgs(sysma=NULL, mytargets="targets_bam.txt")
bfl <- BamFileList(outpaths(args_bam), yieldSize=50000, index=character())
countDFnames <- countRangeset(bfl, args, mode="Union", ignore.strand=TRUE)
writeTargetsout(x=args, file="targets_countDF.txt", overwrite=TRUE)
@

\section{Differential binding analysis of peaks}
The function \Rfunction{runDiff} performs differential binding analysis in batch mode for several count tables using \Rpackage{edgeR} or \Rpackage{DESeq2} \citep{Robinson2010-uk, Love2014-sh}. Internally, it calls the functions \Rfunction{run\_edgeR} and \Rfunction{run\_DESeq2}. It also returns the filtering results and plots from the downstream \Rfunction{filterDEGs} function using the fold change and FDR cutoffs provided under the \Rfunarg{dbrfilter} argument.
 
<<eval=FALSE>>=
args_diff <- systemArgs(sysma="param/rundiff.param", mytargets="targets_countDF.txt")
cmp <- readComp(file=args_bam, format="matrix") 
dbrlist <- runDiff(args=args_diff, diffFct=run_edgeR, targets=targetsin(args_bam), 
                    cmp=cmp[[1]], independent=TRUE, dbrfilter=c(Fold=2, FDR=1))
writeTargetsout(x=args_diff, file="targets_rundiff.txt", overwrite=TRUE)
@

\section{GO term enrichment analysis}
The following performs GO term enrichment analysis for each annotated peak set. 
<<eval=FALSE>>=
args <- systemArgs(sysma="param/macs2.param", mytargets="targets_bam_ref.txt")
args_anno <- systemArgs(sysma="param/annotate_peaks.param", mytargets="targets_macs.txt")
annofiles <- outpaths(args_anno)
gene_ids <- sapply(names(annofiles), function(x) unique(as.character(read.delim(annofiles[x])[,"gene_id"])))
load("data/GO/catdb.RData")
BatchResult <- GOCluster_Report(catdb=catdb, setlist=gene_ids, method="all", id_type="gene", CLSZ=2, cutoff=0.9, gocats=c("MF", "BP", "CC"), recordSpecGO=NULL)
@

\section{Motif analysis}
\subsection{Parse DNA sequences of peak regions from genome}
Enrichment analysis of known DNA binding motifs or \textit{de novo} discovery of novel motifs requires the DNA sequences of the identified peak regions. To parse the corresponding sequences from the reference genome, the \Rfunction{getSeq} function from the \Rpackage{Biostrings} package can be used. The following example parses the sequences for each peak set and saves the results to separate FASTA files, one for each peak set. In addition, the sequences in the FASTA files are ranked (sorted) by increasing p-values as expected by some motif discovery tools, such as \Rpackage{BCRANK}.

<<eval=FALSE>>=
library(Biostrings); library(seqLogo); library(BCRANK)
args <- systemArgs(sysma="param/annotate_peaks.param", mytargets="targets_macs.txt")
rangefiles <- infile1(args)
for(i in seq(along=rangefiles)) {
    df <- read.delim(rangefiles[i], comment="#")
    peaks <- as(df, "GRanges")
    names(peaks) <- paste0(as.character(seqnames(peaks)), "_", start(peaks), "-", end(peaks))
    peaks <- peaks[order(values(peaks)$X.log10.pvalue, decreasing=TRUE)]
    pseq <- getSeq(FaFile("./data/tair10.fasta"), peaks)
    names(pseq) <- names(peaks)
    writeXStringSet(pseq, paste0(rangefiles[i], ".fasta")) 
}
@

\subsection{Motif discovery with \Rpackage{BCRANK}}
The Bioconductor package \Rpackage{BCRANK} is one of the many tools available for \textit{de novo} discovery of DNA binding motifs in peak regions of ChIP-Seq experiments. The given example applies this method on the first peak sample set and plots the sequence logo of the highest ranking motif.

<<eval=FALSE>>=
set.seed(0)
BCRANKout <- bcrank(paste0(rangefiles[1], ".fasta"), restarts=25, use.P1=TRUE, use.P2=TRUE)
toptable(BCRANKout)
topMotif <- toptable(BCRANKout, 1)
weightMatrix <- pwm(topMotif, normalize = FALSE)
weightMatrixNormalized <- pwm(topMotif, normalize = TRUE)
pdf("results/seqlogo.pdf")
seqLogo(weightMatrixNormalized)
dev.off()
@
\begin{figure}[H]
  \centering
   \includegraphics[width=9cm]{seqlogo.pdf}
   \caption{One of the motifs identified by \Rpackage{BCRANK}}
   \label{fig:seqlogo}
\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) and the National Science Foundation (NSF).

\bibliography{bibtex}

\end{document}