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

## ----load-purecn, echo=FALSE, message=FALSE------------------------------
library(PureCN)

## ----example_files, message=FALSE, warning=FALSE, results='hide'---------

gatk.normal.file <- system.file("extdata", "example_normal.txt",
    package="PureCN") 
gatk.normal2.file <- system.file("extdata", "example_normal2.txt", 
    package="PureCN") 
gatk.normal.files <- c(gatk.normal.file, gatk.normal2.file) 
gatk.tumor.file <- system.file("extdata", "example_tumor.txt", 
    package="PureCN") 
vcf.file <- system.file("extdata", "example_vcf.vcf", package="PureCN")

## ----example_files2------------------------------------------------------
gc.gene.file <- system.file("extdata", "example_gc.gene.file.txt", 
    package="PureCN")

## ----gccorrect-----------------------------------------------------------
gatk.normal.file <- system.file("extdata", "example_normal.txt", 
    package = "PureCN")
gc.gene.file <- system.file("extdata", "example_gc.gene.file.txt", 
    package = "PureCN")
coverage <- correctCoverageBias(gatk.normal.file, gc.gene.file, "example_normal_loess.txt")

## ----normaldb------------------------------------------------------------
normalDB <- createNormalDatabase(gatk.normal.files)

## ----normaldbpca---------------------------------------------------------
# get the best normal
gatk.best.normal.file <- findBestNormal(gatk.tumor.file, normalDB)

## ----normaldbpcapool-----------------------------------------------------
# get the best 2 normals and average them
gatk.best.normal.files <- findBestNormal(gatk.tumor.file, normalDB, 
    num.normals=2)
pool <- poolCoverage(lapply(gatk.best.normal.files, readCoverageGatk),
     remove.chrs=c('chrX', 'chrY'))

## ----exonweightfile1-----------------------------------------------------
exon.weight.file <- "exon_weights.txt"
createExonWeightFile(gatk.tumor.file, gatk.normal.files, exon.weight.file)

## ----exonweightfile2-----------------------------------------------------
#createExonWeightFile(gatk.normal.files[1], gatk.normal.files[-1], 
#   "exon_weights.txt")

## ----snpblacklist--------------------------------------------------------
#mutect.normal.files <- dir("poolofnormals", pattern="vcf$", full.names=TRUE) 
#snp.blacklist <- createSNPBlacklist(mutect.normal.files)
#write.csv(snp.bl[[1]], file="SNP_blacklist.csv")
#write.csv(snp.bl[[2]], file="SNP_blacklist_segmented.csv", row.names=FALSE, 
#   quote=FALSE)

## ----runpurecn-----------------------------------------------------------
ret <-runAbsoluteCN(gatk.normal.file=pool, gatk.tumor.file=gatk.tumor.file, 
    vcf.file=vcf.file, sampleid='Sample1', gc.gene.file=gc.gene.file, 
#args.filterVcf=list(snp.blacklist=snp.blacklist, stats.file=mutect.stats.file), 
    args.segmentation=list(exon.weight.file=exon.weight.file), 
    post.optimize=FALSE, plot.cnv=FALSE)

## ----createoutput--------------------------------------------------------
file.rds <- 'Sample1_PureCN.rds'
saveRDS(ret, file=file.rds)
pdf('Sample1_PureCN.pdf', width=10, height=12)
plotAbs(ret, type='all')
dev.off()

## ----figureexample1, fig.show='hide'-------------------------------------
plotAbs(ret, type="overview")

## ----figureexample2, fig.show='hide'-------------------------------------
plotAbs(ret, 1, type="hist")

## ----figureexample3, fig.show='hide'-------------------------------------
plotAbs(ret, 1, type="BAF")

## ----figureexample4, fig.show='hide'-------------------------------------
plotAbs(ret, 1, type="AF")

## ----curationfile--------------------------------------------------------
createCurationFile(file.rds) 

## ----readcurationfile----------------------------------------------------
ret <- readCurationFile(file.rds)

## ----curationfileshow----------------------------------------------------
read.csv('Sample1_PureCN.csv')

## ----output1-------------------------------------------------------------
names(ret)
head(ret$results[[1]]$gene.calls, 3)

## ----output2-------------------------------------------------------------
head(ret$results[[1]]$SNV.posterior$beta.model$posteriors, 3)

## ----output3-------------------------------------------------------------
head(predictSomatic(ret), 3)

## ----sessioninfo, results='asis', echo=FALSE-----------------------------
toLatex(sessionInfo())