### R code from vignette source 'SomatiCA.Rnw'

###################################################
### code chunk number 1: SomatiCA.Rnw:59-61
###################################################
library(SomatiCA)
SomatiCAUsersGuide()


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


###################################################
### code chunk number 3: SomatiCA.Rnw:71-86
###################################################
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)
data <- data.frame(chr, as.integer(position), as.character(zygo), as.integer(reads1), rawLAF, as.integer(reads2), germLAF) 



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


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


###################################################
### code chunk number 6: SomatiCA.Rnw:106-107
###################################################
seg


###################################################
### code chunk number 7: SomatiCA.Rnw:111-112 (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:128-130
###################################################
segmentwithRatio <- somaticRatio(seg$segment, input, method = "mle") 
segmentwithRatio


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


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


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


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


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


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


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


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


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


###################################################
### code chunk number 19: SomatiCA.Rnw:200-205 (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))
## }