### R code from vignette source 'vignettes/TitanCNA/inst/doc/TitanCNA.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: TitanCNA.Rnw:75-80
###################################################
library(TitanCNA)
infile <- system.file("extdata", "test_alleleCounts_chr2.txt", 
                      package = "TitanCNA")
data <- loadAlleleCounts(infile)
names(data)


###################################################
### code chunk number 2: TitanCNA.Rnw:88-92
###################################################
numClusters <- 2
params <- loadDefaultParameters(copyNumber = 5, 
                                numberClonalClusters = numClusters)
params


###################################################
### code chunk number 3: TitanCNA.Rnw:105-107 (eval = FALSE)
###################################################
## params$ploidyParams$phi_0 <- 2 # for diploid or
## params$ploidyParams$phi_0 <- 4 # for tetraploid/ployploid


###################################################
### code chunk number 4: TitanCNA.Rnw:112-115
###################################################
K <- length(params$genotypeParams$alphaKHyper)
params$genotypeParams$alphaKHyper <- rep(500, K)
params$ploidyParams$phi_0 <- 1.5


###################################################
### code chunk number 5: TitanCNA.Rnw:123-129
###################################################
tumWig <- system.file("extdata", "test_tum_chr2.wig", package = "TitanCNA")
normWig <- system.file("extdata", "test_norm_chr2.wig", package = "TitanCNA")
gc <- system.file("extdata", "gc_chr2.wig", package = "TitanCNA")
map <- system.file("extdata", "map_chr2.wig", package = "TitanCNA")
cnData <- correctReadDepth(tumWig, normWig, gc, map)
head(cnData)


###################################################
### code chunk number 6: TitanCNA.Rnw:136-138
###################################################
logR <- getPositionOverlap(data$chr, data$posn, cnData)
data$logR <- log(2^logR)  #transform the log ratio to natural logs


###################################################
### code chunk number 7: TitanCNA.Rnw:144-146
###################################################
data <- filterData(data, 1:24, minDepth = 10, maxDepth = 200, 
                   positionList = NULL)


###################################################
### code chunk number 8: TitanCNA.Rnw:153-155 (eval = FALSE)
###################################################
## library(doMC)
## registerDoMC(cores = 10) #use 10 cores on a single machine


###################################################
### code chunk number 9: TitanCNA.Rnw:170-180
###################################################
convergeParams <- runEMclonalCN(data, gParams = params$genotypeParams, 
                                nParams = params$normalParams, 
                                pParams = params$ploidyParams, 
                                sParams = params$cellPrevParams, 
                                maxiter = 3, maxiterUpdate = 50, 
                                useOutlierState = FALSE, txnExpLen = 1e15, 
                                txnZstrength = 5e5, 
                                normalEstimateMethod = "map", 
                                estimateS = TRUE, estimatePloidy = TRUE)
names(convergeParams)


###################################################
### code chunk number 10: TitanCNA.Rnw:200-202
###################################################
optimalPath <- viterbiClonalCN(data, convergeParams)
head(optimalPath)


###################################################
### code chunk number 11: TitanCNA.Rnw:213-217
###################################################
results <- outputTitanResults(data, convergeParams, optimalPath,
                              filename = NULL, posteriorProbs = FALSE,
                              subcloneProfiles = TRUE)
head(results)


###################################################
### code chunk number 12: TitanCNA.Rnw:244-246 (eval = FALSE)
###################################################
## outparam <- paste("test_cluster02_params.txt", sep = "")
## outputModelParameters(convergeParams, results, outparam)


###################################################
### code chunk number 13: TitanCNA.Rnw:271-275
###################################################
ploidy <- tail(convergeParams$phi, 1)
ploidy
plotCNlogRByChr(results, chr = 2, ploidy = ploidy, ylim = c(-2, 2), cex = 0.25, 
                xlab = "", main = "Chr 2")


###################################################
### code chunk number 14: TitanCNA.Rnw:288-290
###################################################
plotAllelicRatio(results, chr = 2, ylim = c(0, 1), cex = 0.25, 
                 xlab = "", main = "Chr 2")


###################################################
### code chunk number 15: TitanCNA.Rnw:304-309
###################################################
norm <- tail(convergeParams$n, 1) 
norm # estimated normal contamination
1 - convergeParams$s[, ncol(convergeParams$s)] # estimated cellular prevalence 
plotClonalFrequency(results, chr = 2, normal = norm, ylim = c(0, 1), 
                    cex = 0.25, xlab = "", main = "Chr 2")


###################################################
### code chunk number 16: TitanCNA.Rnw:318-319
###################################################
plotSubcloneProfiles(results, chr = 2, cex = 2, spacing = 6, main = "Chr 2")


###################################################
### code chunk number 17: TitanCNA.Rnw:328-329
###################################################
toLatex(sessionInfo())