## ----style, echo=FALSE, results="asis", message=FALSE-------------------------
BiocStyle::markdown()
knitr::opts_chunk$set(tidy = FALSE,
                      warning = FALSE,
                      message = FALSE)

## ----echo=FALSE, results='hide', message=FALSE--------------------------------
library(BaalChIP)

## ----eval=TRUE----------------------------------------------------------------
wd <- system.file("test",package="BaalChIP") #system directory

samplesheet <- file.path(wd, "exampleChIP1.tsv")
hets <- c("MCF7"="MCF7_hetSNP.txt", "GM12891"="GM12891_hetSNP.txt")
res <- BaalChIP(samplesheet=samplesheet, hets=hets)

## ----quick1, eval=FALSE-------------------------------------------------------
# res <- BaalChIP(samplesheet=samplesheet, hets=hets)
# res <- BaalChIP.run(res, cores=4, verbose=TRUE) #cores for parallel computing

## ----quick2, eval=FALSE-------------------------------------------------------
# 
# #create BaalChIP object
# samplesheet <- file.path(wd, "exampleChIP1.tsv")
# hets <- c("MCF7"="MCF7_hetSNP.txt", "GM12891"="GM12891_hetSNP.txt")
# res <- BaalChIP(samplesheet=samplesheet, hets=hets)
# 
# #Now, load some data
# data(blacklist_hg19)
# data(pickrell2011cov1_hg19)
# data(UniqueMappability50bp_hg19)
# 
# #run one at the time (instead of BaalChIP.run)
# res <- alleleCounts(res, min_base_quality=10, min_mapq=15, verbose=FALSE)
# res <- QCfilter(res,
#                 RegionsToFilter=c("blacklist_hg19", "pickrell2011cov1_hg19"),
#                 RegionsToKeep="UniqueMappability50bp_hg19",
#                 verbose=FALSE)
# res <- mergePerGroup(res)
# res <- filter1allele(res)
# res <- getASB(res,
#               Iter=5000,
#               conf_level=0.95,
#               cores=4, RMcorrection = TRUE,
#               RAFcorrection=TRUE)

## -----------------------------------------------------------------------------
gDNA <- list("TaT1"=
               c(file.path(wd, "bamFiles/TaT1_1_gDNA.test.bam"),
                 file.path(wd, "bamFiles/TaT1_2_gDNA.test.bam")),
              "AMOVC"=
               c(file.path(wd, "bamFiles/AMOVC_1_gDNA.test.bam"),
                 file.path(wd, "bamFiles/AMOVC_2_gDNA.test.bam")))

## ----quick3, eval=FALSE-------------------------------------------------------
# wd <- system.file("test",package="BaalChIP") #system directory
# 
# samplesheet <- file.path(wd, "exampleChIP2.tsv")
# hets <- c("TaT1"="TaT1_hetSNP.txt", "AMOVC"="AMOVC_hetSNP.txt")
# res2 <- BaalChIP(samplesheet=samplesheet, hets=hets, CorrectWithgDNA=gDNA)
# res2 <- BaalChIP.run(res2, cores=4, verbose=TRUE) #cores for parallel computing

## -----------------------------------------------------------------------------
samplesheet <- read.delim(file.path(wd,"exampleChIP1.tsv"))
samplesheet

## -----------------------------------------------------------------------------
samplesheet <- read.delim(file.path(wd,"exampleChIP2.tsv"))
samplesheet

## -----------------------------------------------------------------------------
head(read.delim(file.path(system.file("test",package="BaalChIP"),"MCF7_hetSNP.txt")))

## ----eval=TRUE----------------------------------------------------------------
samplesheet <- file.path(wd,"exampleChIP1.tsv")
hets <- c("MCF7"="MCF7_hetSNP.txt", "GM12891"="GM12891_hetSNP.txt")
res <- BaalChIP(samplesheet=samplesheet, hets=hets)
res

## -----------------------------------------------------------------------------
BaalChIP.get(res, what="samples")

## ----eval=TRUE----------------------------------------------------------------
#run alleleCounts
res <- alleleCounts(res, min_base_quality=10, min_mapq=15, verbose=FALSE)

## -----------------------------------------------------------------------------
data(blacklist_hg19)
data(pickrell2011cov1_hg19)
data(UniqueMappability50bp_hg19)

## ----eval=TRUE----------------------------------------------------------------
#run QC filter
res <- QCfilter(res, 
                RegionsToFilter=list("blacklist"=blacklist_hg19, "highcoverage"=pickrell2011cov1_hg19), 
                RegionsToKeep=list("UniqueMappability"=UniqueMappability50bp_hg19), 
                verbose=FALSE)
res <- mergePerGroup(res)
res <- filter1allele(res)

## ----eval=FALSE, message=FALSE, error=FALSE, warning=FALSE--------------------
# res <- filterIntbias(res,
#                      simul_output="directory_name",
#                      tmpfile_prefix="prefix_name",
#                      simulation_script = "local",
#                      alignmentSimulArgs=c("picard-tools-1.119",
#                                           "bowtie-1.1.1",
#                                           "genomes_test/male.hg19",
#                                           "genomes_test/maleByChrom")
#                      verbose=FALSE)

## -----------------------------------------------------------------------------
preComputed_output <- system.file("test/simuloutput",package="BaalChIP")
list.files(preComputed_output)

## ----eval=TRUE, message=FALSE-------------------------------------------------
res <- filterIntbias(res, simul_output=preComputed_output, 
                     tmpfile_prefix="c67c6ec6c433", 
                     skipScriptRun=TRUE,
                     verbose=FALSE)

## -----------------------------------------------------------------------------
res <- mergePerGroup(res)

## -----------------------------------------------------------------------------
res <- filter1allele(res)

## ----eval=FALSE, message=FALSE------------------------------------------------
# res <- getASB(res, Iter=5000, conf_level=0.95, cores = 4,
#               RMcorrection = TRUE,
#               RAFcorrection=TRUE)

## ----eval=TRUE, echo=FALSE, include=FALSE-------------------------------------
data(baalObject)
res <- BaalObject

## -----------------------------------------------------------------------------
res

## -----------------------------------------------------------------------------
result <- BaalChIP.report(res)
head(result[["MCF7"]])

## -----------------------------------------------------------------------------
data(ENCODEexample)
ENCODEexample

## -----------------------------------------------------------------------------
a <- summaryQC(ENCODEexample)

## -----------------------------------------------------------------------------
summaryQC(ENCODEexample)[["filtering_stats"]]

## -----------------------------------------------------------------------------
summaryQC(ENCODEexample)[["average_stats"]]

## ----plotENCODE1--------------------------------------------------------------
plotQC(ENCODEexample, what="barplot_per_group")

## ----plotENCODE2--------------------------------------------------------------
plotQC(ENCODEexample, what="boxplot_per_filter")

## ----plotENCODE3--------------------------------------------------------------
plotQC(ENCODEexample, what="overall_pie")

## ----plotSimul----------------------------------------------------------------
plotSimul(ENCODEexample)

## -----------------------------------------------------------------------------
#ENCODE example
a <- BaalChIP.get(ENCODEexample, "assayedVar")[["MCF7"]]
a[1:5,1:5]

#FAIRE exmaple
a <- BaalChIP.get(FAIREexample, "assayedVar")[["MDA134"]]
a[1:5,]

## -----------------------------------------------------------------------------
summaryASB(ENCODEexample)

## -----------------------------------------------------------------------------
summaryASB(FAIREexample)

## ----adjENCODE, fig.width=12--------------------------------------------------
adjustmentBaalPlot(ENCODEexample)

## ----adjFAIRE1, fig.width=7, fig.height=2.5-----------------------------------
adjustmentBaalPlot(FAIREexample)

## ----adjFAIRE2, fig.width=7, fig.height=2.5-----------------------------------
adjustmentBaalPlot(FAIREexample, col=c("cyan4","chocolate3"))

## -----------------------------------------------------------------------------
biastable <- BaalChIP.get(ENCODEexample, "biasTable")
head(biastable[["K562"]])

## -----------------------------------------------------------------------------
biastable <- BaalChIP.get(FAIREexample, "biasTable")
head(biastable[["T47D"]])

## -----------------------------------------------------------------------------
result <- BaalChIP.report(FAIREexample)[["T47D"]]

#show ASB SNPs
result[result$isASB==TRUE,]

## ----echo=FALSE---------------------------------------------------------------
sessionInfo()