inst/doc/BaalChIP.R

In BaalChIP: BaalChIP: Bayesian analysis of allele-specific transcription factor binding in cancer genomes

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