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

###################################################
### code chunk number 1: exomeCopy.Rnw:18-19
###################################################
options(width = 70)


###################################################
### code chunk number 2: exomeCopy.Rnw:59-62
###################################################
library(exomeCopy)
gr <- GRanges(seqname="seq1",IRanges(start=1,end=345))
subdivideGRanges(gr)


###################################################
### code chunk number 3: exomeCopy.Rnw:68-75
###################################################
  plot(0,0,xlim=c(0,500),ylim=c(0,25),type="n",yaxt="n",ylab="",xlab="width of input GRanges object",main="Effect of subdivideGRanges")
  abline(v=1:5*100,col="grey")
  for (i in 1:24) {
    gr <- GRanges(seqname="chr1",IRanges(start=1,width=(i*20)))
    sbd.gr <- subdivideGRanges(gr)
    arrows(start(sbd.gr),rep(i,length(sbd.gr)),end(sbd.gr),rep(i,length(sbd.gr)),length=.04,angle=90,code=3)
  }


###################################################
### code chunk number 4: exomeCopy.Rnw:81-87
###################################################
target.file <- system.file("extdata","targets.bed",package="exomeCopy")
target.df <- read.delim(target.file,header=FALSE,col.names=c("seqname","start","end")) 
target <- GRanges(seqname=target.df$seqname,IRanges(start=target.df$start+1,end=target.df$end))
target
target.sub <- subdivideGRanges(target)
target.sub


###################################################
### code chunk number 5: exomeCopy.Rnw:96-105
###################################################
bam.file <- system.file("extdata","mapping.bam",package="exomeCopy")
scanBamHeader(bam.file)[[1]]$targets
seqlevels(target.sub)
rdata <- RangedData(space=seqnames(target.sub),ranges=ranges(target.sub))
sample.df <- data.frame(samples="sample1",bam.files=bam.file,stringsAsFactors=FALSE)
for (i in 1:nrow(sample.df)) {
  rdata[[sample.df$samples[i]]] <- countBamInGRanges(sample.df$bam.files[i],target.sub)
}
rdata


###################################################
### code chunk number 6: exomeCopy.Rnw:111-114
###################################################
reference.file <- system.file("extdata","reference.fa",package="exomeCopy")
rdata[["GC"]] <- getGCcontent(target.sub, reference.file)
rdata


###################################################
### code chunk number 7: exomeCopy.Rnw:123-126
###################################################
data(exomecounts)
dim(exomecounts)
exomecounts[1:5,1:3]


###################################################
### code chunk number 8: exomeCopy.Rnw:132-133
###################################################
plot(start(exomecounts),exomecounts$HG00551,xlim=c(0.8e6,1.8e6),xlab="genomic position",ylab="counts",main="HG00551 read counts in exonic ranges")


###################################################
### code chunk number 9: exomeCopy.Rnw:141-146
###################################################
chr1a <- exomecounts
chr1a.ranges <- ranges(chr1a)
names(chr1a.ranges) <- "chr1a"
chr1a <- RangedData(chr1a.ranges,as.data.frame(exomecounts)[,-c(1:4)])
exomecounts.plus <- c(exomecounts,chr1a)


###################################################
### code chunk number 10: exomeCopy.Rnw:153-156
###################################################
exome.samples <- grep("HG.+",colnames(exomecounts.plus),value=TRUE)
exomecounts.plus[["bg"]] <- generateBackground(exome.samples, exomecounts.plus, median)
exomecounts.plus[["bg.sd"]] <- generateBackground(exome.samples, exomecounts.plus, sd)


###################################################
### code chunk number 11: exomeCopy.Rnw:161-163
###################################################
exomecounts.plus[["GC.sq"]] <- exomecounts.plus$GC^2
exomecounts.plus[["width"]] <- width(exomecounts.plus)


###################################################
### code chunk number 12: exomeCopy.Rnw:208-220
###################################################
simulateCNV <- function(x,indices,multiply,prob) {
  x[indices] <- x[indices] + multiply * rbinom(length(indices),prob=prob,size=x[indices])
  return(x)
}
rdata <- exomecounts.plus
set.seed(2)
cnv.probs <- rep(c(.99,.5,.5,.95),each=2)
cnv.mult <- rep(c(-1,1),each=4)
for (i in 1:8) {
  rdata[[exome.samples[i]]] <- simulateCNV(rdata[[exome.samples[i]]],(i*100 + 1):((i+1)*100),multiply=cnv.mult[i],prob=cnv.probs[i])
  rdata[[exome.samples[i]]] <- simulateCNV(rdata[[exome.samples[i]]],1000 + (i*100 + 1):((i+1)*100),multiply=cnv.mult[i],prob=cnv.probs[i])
}


###################################################
### code chunk number 13: exomeCopy.Rnw:229-231
###################################################
  fit <- exomeCopy(rdata["chr1"],sample.name=exome.samples[3],X.names=c("bg","GC","GC.sq","width"),S=0:6,d=2)
  show(fit)


###################################################
### code chunk number 14: exomeCopy.Rnw:236-237
###################################################
  copyCountSegments(fit)


###################################################
### code chunk number 15: exomeCopy.Rnw:243-245
###################################################
  cnv.cols <- c("red","orange","black","deepskyblue","blue","blue2","blue4")
  plot(fit,col=cnv.cols)


###################################################
### code chunk number 16: exomeCopy.Rnw:253-256
###################################################
runExomeCopy <- function(sample.name,seqs) {
  lapply(seqs,function(seq.name) exomeCopy(rdata[seq.name],sample.name,X.names=c("bg","GC","GC.sq","width"),S=0:4,d=2))
}


###################################################
### code chunk number 17: exomeCopy.Rnw:261-271
###################################################
seqs <- c("chr1","chr1a")
names(seqs) <- seqs
samples <- exome.samples[1:8]
names(samples) <- samples
t0 <- as.numeric(proc.time()[3])
fit.list <- lapply(samples,runExomeCopy,seqs)
t1 <- as.numeric(proc.time()[3])
time.elapsed <- as.numeric(t1 - t0)
paste(round(time.elapsed),"seconds for",length(samples),"samples,",round(sum(width(rdata))/1e3),"kb of target")
paste("~",round(time.elapsed / 60 / (8 * sum(width(rdata))) * 32e6,1)," minutes for 1 sample, 32 Mb of target",sep="")


###################################################
### code chunk number 18: exomeCopy.Rnw:278-280
###################################################
compiled.segments <- compileCopyCountSegments(fit.list)
CNV.segments <- compiled.segments[compiled.segments$copy.count != 2,]


###################################################
### code chunk number 19: exomeCopy.Rnw:288-290
###################################################
CNV.segments[1:6,]
CNV.segments <- CNV.segments[CNV.segments$nranges > 5,]


###################################################
### code chunk number 20: exomeCopy.Rnw:297-299
###################################################
cnv.cols <- c("red","orange","black","deepskyblue","blue")
plotCompiledCNV(CNV.segments=CNV.segments, seq.name="chr1", col=cnv.cols)


###################################################
### code chunk number 21: exomeCopy.Rnw:302-303
###################################################
plotCompiledCNV(CNV.segments=CNV.segments, seq.name="chr1a", col=cnv.cols)


###################################################
### code chunk number 22: exomeCopy.Rnw:311-313
###################################################
fit.list[[1]][[1]]@init.par$beta.hat
fit.list[[1]][[1]]@final.par$beta


###################################################
### code chunk number 23: exomeCopy.Rnw:319-324
###################################################
par(mfrow=c(3,3),mar=c(3,5,1,1))
for (i in 1:8) {
  barplot(rbind(fit.list[[i]][[1]]@final.par$beta,fit.list[[i]][[1]]@init.par$beta.hat),horiz=TRUE,las=1,beside=TRUE,col=c("white","darkgrey"),xlim=c(-150,350))
  legend("topright",legend=c("initial","final"),fill=c("darkgrey","white"))
}


###################################################
### code chunk number 24: exomeCopy.Rnw:330-331
###################################################
sessionInfo()