tumorVsNormals.captureRegions.R

tumorVsNormals.captureRegions.R

usage <- function()
{
writeLines("Usage:\n\tRscript tumorVsNormal.R tumor.bam normal.bam capture.bed")
}
 
args <- commandArgs(trailingOnly = TRUE)
 
if(length(args) != 3) {
usage()
quit()
}
 
library(cn.mops)
library(snow)
library(RUnit)
 
tumor_bam = args[1]
normal_bam = args[2]
capture_bed = args[3]


segments <- (read.table(capture_bed,sep="\t",as.is=TRUE))
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))

X.nor <- getSegmentReadCountsFromBAM(normal_bam, GR=gr, mode="paired")
X.tum <- getSegmentReadCountsFromBAM(tumor_bam, GR=gr, mode="paired")

X <- X.tum
values(X) <- cbind(values(X.tum), values(X.nor))
X <- normalizeGenome(X)

ref_analysis_norm <- referencecn.mops(X[,1], X[,2], segAlgorithm="DNAcopy")
ref_analysis_norm <- cn.mops:::.replaceNames(ref_analysis_norm, colnames(ref_analysis_norm@normalizedData),"tumor"))
ref_analysis_norm_file = "ref_analysis_norm.Robject"
save(ref_analysis_norm, file=ref_analysis_norm_file)


segment_file="cn.mops-segplot.png"
png(segment_file)
segplot(ref_analysis_norm)
abline(h=log2(1/2), col = "green")
abline(h=log2(3/2), col = "blue")
abline(h=log2(4/2), col = "yellow")
legend("bottomright", inset=.05, title="Copy Number difference (tumor - normal)", c("1","2","-1"), horiz=TRUE, col = c("blue", "yellow", "green"), lty=1)
dev.off()