tumorVsNormal.R

tumorVsNormal.R

usage <- function() 
{
  writeLines("Usage:\n\tRscript tumorVsNormal.R tumor.bam normal.bam chromosome")
}

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]
chr = args[3]
tumor_gr <- getReadCountsFromBAM(tumor_bam, refSeqName=chr, WL=500, mode="paired", sampleNames="tumor")
normal_gr <- getReadCountsFromBAM(normal_bam, refSeqName=chr, WL=500, mode="paired", sampleNames="normal")
X <- tumor_gr
values(X) <- cbind(values(tumor_gr), values(normal_gr)) 
X <- normalizeGenome(X)
ref_analysis_norm <- referencecn.mops(X[,1], X[,2], norm=0, I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8,16,32,64), classes = c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6", "CN7","CN8","CN16","CN32","CN64","CN128"), segAlgorithm="DNAcopy")
ref_analysis_norm_file = paste("res_analysis.object_chr", chr, sep = "")
save(ref_analysis_norm, file=ref_analysis_norm_file)
if(length(cnvr(ref_analysis_norm))) {     
  # resCNMOPS <- calcIntegerCopyNumbers(ref_analysis_norm)
  resCNMOPS <-cn.mops:::.replaceNames(ref_analysis_norm, colnames(ref_analysis_norm@normalizedData),"tumor")
  segment_file=paste("cn.mops-segplot.chr", chr, ".png", sep = "") 
  png(segment_file)
  segplot(resCNMOPS)
  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()
}