### R code from vignette source 'vignettes/SomatiCA/inst/doc/SomatiCA.Rnw'

###################################################
### code chunk number 1: SomatiCA.Rnw:62-64
###################################################
library(SomatiCA)
SomatiCAUsersGuide()


###################################################
### code chunk number 2: SomatiCA.Rnw:70-71
###################################################
options(width=60)


###################################################
### code chunk number 3: SomatiCA.Rnw:74-89
###################################################
set.seed(1)
rawLAF <- c(rnorm(300, 0.2, 0.05), rnorm(300, 0.4, 0.05), 
            rnorm(200, 0.3, 0.05), rnorm(200, 0.2, 0.05), 
            rnorm(200, 0.3, 0.05), rnorm(250, 0.4, 0.05)) 
rawLAF <- ifelse(rawLAF>0.5, 1-rawLAF, rawLAF) 
germLAF <- c(rnorm(800+650, 0.4, 0.05)) 
germLAF <- ifelse(germLAF>0.5, 1-germLAF, germLAF) 
reads1 <- c(rpois(300, 25), rpois(300, 50), rpois(200, 60),  
            rpois(200, 25), rpois(200, 40), rpois(250, 50))
reads2 <- rpois(800+650,50)
chr <- c(rep("chr1", 800), rep("chr2", 650))
position <- c(seq(1, 16000000, by=20000), seq(1, 13000000, by=20000))
zygo <- rep("het", 800+650)
x <- data.frame(chr, as.integer(position), as.character(zygo), as.integer(reads1), rawLAF, as.integer(reads2), germLAF) 



###################################################
### code chunk number 4: SomatiCA.Rnw:92-96
###################################################
colnames(x) <- c("seqnames", "start", "zygosity", "tCount",
                 "LAF", "tCountN", "germLAF")            
input <- SomatiCAFormat(x, missing = F, verbose = T)
input


###################################################
### code chunk number 5: SomatiCA.Rnw:104-105
###################################################
seg <- larsCBSsegment(input, collapse.k = 0, ncores = 1, verbose = T, rss = FALSE)


###################################################
### code chunk number 6: SomatiCA.Rnw:109-110
###################################################
seg


###################################################
### code chunk number 7: SomatiCA.Rnw:114-115 (eval = FALSE)
###################################################
## plotSegment(seg$segment, input, k = 1, smooth = F)


###################################################
### code chunk number 8: fig1
###################################################
plotSegment(seg$segment, input, k = 1, smooth = F, dev.new=FALSE)


###################################################
### code chunk number 9: SomatiCA.Rnw:131-133
###################################################
segmentwithRatio <- somaticRatio(seg$segment, input, method = "mle") 
segmentwithRatio


###################################################
### code chunk number 10: SomatiCA.Rnw:138-141
###################################################
refined <- refineSegment(segmentwithRatio, input, method = "mle", 
                        adjust = TRUE, threshold1 = 0 , threshold2 = 0.05)   	
refined


###################################################
### code chunk number 11: SomatiCA.Rnw:148-150
###################################################
x <- admixtureRate(refined, mcmc = 5000, burnin = 1000, p = 0.01) 
admix <- x$admix


###################################################
### code chunk number 12: SomatiCA.Rnw:154-157
###################################################
y <- copynumberCorrected(segmentwithRatio, admix)
y



###################################################
### code chunk number 13: SomatiCA.Rnw:161-164
###################################################
data(GCcontent)
segmentGCcorrected <- segmentGCbiasRemoval(y, input, GCcontent)
segmentClonality <- subclonality(segmentGCcorrected, admix)


###################################################
### code chunk number 14: SomatiCA.Rnw:168-169
###################################################
merged <- MergeSegment(segmentClonality)


###################################################
### code chunk number 15: SomatiCA.Rnw:173-174 (eval = FALSE)
###################################################
## plotSubclonality(segmentClonality)


###################################################
### code chunk number 16: fig2
###################################################
plotSubclonality(segmentClonality, dev.new=FALSE)


###################################################
### code chunk number 17: SomatiCA.Rnw:188-189
###################################################
merged[elementMetadata(merged)[, "clonality"]=="clonal", ]


###################################################
### code chunk number 18: SomatiCA.Rnw:194-198 (eval = FALSE)
###################################################
## data(GCcontent)
## x <- SomatiCApipe(input, ncores = 1, collapse.k = 0, method = "mle", 
##                   mcmc = 50000, burnin = 10000, p = 0.001, GC = GCcontent)
## merged <- MergeSegment(x$segment)


###################################################
### code chunk number 19: SomatiCA.Rnw:204-209 (eval = FALSE)
###################################################
## chr <- paste("chr", c(1:22, "X"), sep="")
## url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/"
## GC <- foreach(i=chr, .combine=rbind)%dopar%{
##       return(GCcount(i, 10000, url))
## }