## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"------------------
BiocStyle::latex(use.unsrturl=FALSE)

## ----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")

## ----eval=TRUE-----------------------------------------------------------
library(systemPipeR)

## ----eval=FALSE----------------------------------------------------------
## library(systemPipeRdata)
## genWorkenvir(workflow="varseq")
## setwd("varseq")

## ----eval=FALSE----------------------------------------------------------
## source("systemPipeVARseq_Fct.R")

## ----eval=TRUE-----------------------------------------------------------
targetspath <- system.file("extdata", "targetsPE.txt", package="systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")[,1:5]
targets

## ----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)

## ----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()

## ----eval=FALSE----------------------------------------------------------
## args <- systemArgs(sysma="bwa.param", mytargets="targets.txt")
## sysargs(args)[1] # Command-line parameters for first FASTQ file

## ----eval=FALSE----------------------------------------------------------
## bampaths <- runCommandline(args=args)

## ----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)

## ----eval=FALSE----------------------------------------------------------
## file.exists(outpaths(args))

## ----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")

## ----eval=FALSE----------------------------------------------------------
## read_statsDF <- alignStats(args=args)
## write.table(read_statsDF, "results/alignStats.xls", row.names=FALSE, quote=FALSE, sep="\t")

## ----eval=FALSE----------------------------------------------------------
## symLink2bam(sysargs=args, htmldir=c("~/.html/", "somedir/"),
##             urlbase="http://biocluster.ucr.edu/~tgirke/",
## 	    urlfile="./results/IGVurl.txt")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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")

## ----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()

## ----sessionInfo, results='asis'-----------------------------------------
toLatex(sessionInfo())