exomeCopy.R
In exomeCopy: Copy number variant detection from exome sequencing read depth

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

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


###################################################
### code chunk number 2: exomeCopy.Rnw:55-77 (eval = FALSE)
###################################################
## library(exomeCopy)
## target.file <- "targets.bed"
## bam.files <- c("/path/to/file1.bam", "/path/to/file2.bam", "/path/to/file3.bam")
## sample.names <- c("sample1","sample2","sample3")
## reference.file <- "/path/to/reference_genome.fa"
## target.df <- read.delim(target.file, header = FALSE)
## target <- GRanges(seqname = target.df[,1], IRanges(start = target.df[,2] + 1, end = target.df[,3]))
## counts <- target
## for (i in 1:length(bam.files)) {
##   mcols(counts)[[sample.names[i]]] <- countBamInGRanges(bam.files[i], target)
## }
## counts$GC <- getGCcontent(target, reference.file)
## counts$GC.sq <- counts$GC^2
## counts$bg <- generateBackground(sample.names, counts, median)
## counts$log.bg <- log(counts$bg + .1)
## counts$width <- width(counts)
## fit.list <- lapply(sample.names, function(sample.name) {
##   lapply(seqlevels(target), function(seq.name) { 
##     exomeCopy(counts[seqnames(counts) == seq.name], sample.name, X.names = c("log.bg", "GC", "GC.sq","width"), S = 0:4, d = 2)
##   })
## })
## compiled.segments <- compileCopyCountSegments(fit.list)


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


###################################################
### code chunk number 4: exomeCopy.Rnw:99-106
###################################################
  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 5: exomeCopy.Rnw:112-119
###################################################
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 <- reduce(target, min.gapwidth=0)
target.sub <- subdivideGRanges(target)
target.sub


###################################################
### code chunk number 6: exomeCopy.Rnw:128-137
###################################################
bam.file <- system.file("extdata","mapping.bam",package="exomeCopy")
scanBamHeader(bam.file)[[1]]$targets
seqlevels(target.sub)
toy.counts <- target.sub
sample.df <- data.frame(samples="sample1",bam.files=bam.file,stringsAsFactors=FALSE)
for (i in 1:nrow(sample.df)) {
  mcols(toy.counts)[[sample.df$samples[i]]] <- countBamInGRanges(sample.df$bam.files[i],target.sub)
}
toy.counts


###################################################
### code chunk number 7: exomeCopy.Rnw:143-146
###################################################
reference.file <- system.file("extdata","reference.fa",package="exomeCopy")
toy.counts$GC <- getGCcontent(target.sub, reference.file)
toy.counts


###################################################
### code chunk number 8: exomeCopy.Rnw:155-158
###################################################
data(exomecounts)
length(exomecounts)
head(exomecounts)


###################################################
### code chunk number 9: exomeCopy.Rnw:164-165
###################################################
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 10: exomeCopy.Rnw:173-176
###################################################
chr1a <- exomecounts
seqlevels(chr1a) <- "chr1a"
example.counts <- c(exomecounts, chr1a)


###################################################
### code chunk number 11: exomeCopy.Rnw:183-187
###################################################
exome.samples <- grep("HG.+",colnames(mcols(example.counts)),value=TRUE)
example.counts$bg <- generateBackground(exome.samples, example.counts, median)
example.counts$log.bg <- log(example.counts$bg + .1)
example.counts$bg.var <- generateBackground(exome.samples, example.counts, var)


###################################################
### code chunk number 12: exomeCopy.Rnw:192-194
###################################################
summary(example.counts$bg)
example.counts <- example.counts[example.counts$bg > 0,]


###################################################
### code chunk number 13: exomeCopy.Rnw:199-201
###################################################
example.counts$GC.sq <- example.counts$GC^2
example.counts$width <- width(example.counts)


###################################################
### code chunk number 14: exomeCopy.Rnw:245-257
###################################################
simulateCNV <- function(x,indices,multiply,prob) {
  x[indices] <- x[indices] + multiply * rbinom(length(indices),prob=prob,size=x[indices])
  return(x)
}
set.seed(2)
cnv.probs <- rep(c(.99,.5,.5,.95),each=2)
cnv.mult <- rep(c(-1,1),each=4)
cnv.starts <- rep(c(1,301,601,901),each=2)
for (i in 1:8) {
  mcols(example.counts)[[exome.samples[i]]] <- simulateCNV(mcols(example.counts)[[exome.samples[i]]],cnv.starts[i]:(cnv.starts[i] + 99),multiply=cnv.mult[i],prob=cnv.probs[i])
  mcols(example.counts)[[exome.samples[i]]] <- simulateCNV(mcols(example.counts)[[exome.samples[i]]],1000 + cnv.starts[i]:(cnv.starts[i] + 99),multiply=cnv.mult[i],prob=cnv.probs[i])
}


###################################################
### code chunk number 15: exomeCopy.Rnw:266-268
###################################################
  fit <- exomeCopy(example.counts[seqnames(example.counts) == "chr1"],sample.name=exome.samples[3],X.names=c("log.bg","GC","GC.sq","width"),S=0:6,d=2)
  show(fit)


###################################################
### code chunk number 16: exomeCopy.Rnw:273-274
###################################################
  copyCountSegments(fit)


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


###################################################
### code chunk number 18: exomeCopy.Rnw:290-293
###################################################
runExomeCopy <- function(sample.name,seqs) {
  lapply(seqs,function(seq.name) exomeCopy(example.counts[seqnames(example.counts) == seq.name],sample.name,X.names=c("log.bg","GC","GC.sq","width"),S=0:4,d=2))
}


###################################################
### code chunk number 19: exomeCopy.Rnw:298-308
###################################################
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(example.counts))/1e3),"kb of target")
paste("~",round(time.elapsed / 60 / (8 * sum(width(example.counts))) * 32e6,1)," minutes for 1 sample, 32 Mb of target",sep="")


###################################################
### code chunk number 20: exomeCopy.Rnw:315-317
###################################################
compiled.segments <- compileCopyCountSegments(fit.list)
CNV.segments <- compiled.segments[compiled.segments$copy.count != 2]


###################################################
### code chunk number 21: exomeCopy.Rnw:325-328
###################################################
CNV.segments[1:6]
table(CNV.segments$nranges)
CNV.segments <- CNV.segments[CNV.segments$nranges > 5]


###################################################
### code chunk number 22: exomeCopy.Rnw:334-336
###################################################
CNV.overlaps.matrix <- as.matrix(findOverlaps(CNV.segments,drop.self=TRUE))
head(CNV.overlaps.matrix)


###################################################
### code chunk number 23: exomeCopy.Rnw:342-346
###################################################
par(mfrow=c(2,1),mar=c(4,3,2,1))
cnv.cols <- c("red","orange","black","deepskyblue","blue")
plotCompiledCNV(CNV.segments=CNV.segments, seq.name="chr1", col=cnv.cols)
plotCompiledCNV(CNV.segments=CNV.segments, seq.name="chr1a", col=cnv.cols)


###################################################
### code chunk number 24: exomeCopy.Rnw:354-356
###################################################
fit.list[[1]][[1]]@init.par$beta.hat
fit.list[[1]][[1]]@final.par$beta


###################################################
### code chunk number 25: exomeCopy.Rnw:361-362
###################################################
sessionInfo()