# Copyright (C) 2015 Patrick Praher and Klambauer Guenter 
# <patrick.praher@gmail.com> <klambauer@bioinf.jku.at> 


require("xlsx")
require("cn.mops")
source("normalizeTumor.R")

##### read target regions from .bed file as genomic ranges
segments <- read.table("nexterarapidcapture_exome_targetedregions_v1.2.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
seqlevels(gr) <- gsub("chr","",seqlevels(gr)) #remove chr from  chromosome names

#Get the read counts from bam file as genomic ranges
#bam_names contains the paths to the bam files and bai_files contains the paths to the bam index files
rc <- cn.mops::getSegmentReadCountsFromBAM(bam_files,GR=gr,basenames(bam_names),mode="paired",BAIFiles=bai_files)

#read the purity estimates from .xlsx file as a numeric vector
purity <- as.numeric(as.character(read.xlsx2("new/CCL/Purity.xlsx",1,header = F)$X2))

#default vector for the fold change in cn.mops
IDefault <- c(0.025,0.5,1,1.5,2,2.5,3,3.5,4)

#fold change in normal data (1 for every copy number)
ICN2 <- rep(1, 9)

#results vecotr
cnvs_cnmops <- list()

#Normalize the samples
start <- 1
temp <- normalizeTumor(rc[[i]],calcCN2Segments(gr, toGR(cnvs_CCL[start:(start+4)])))

#The data sets in this example consists of 5 tumor samples and 1 normal samples, 
for(j in 1:5){
  
  #Calc I
  I <- IDefault * purity[j] + ICN2 * (1-purity[j])
  
  classes <- c("CN0","CN1","CN2","CN3","CN4","CN5","CN6","CN7","CN8")
  
  #calc data from each tumor sample by index
  cases <- temp
  cases@elementMetadata <- DataFrame(cases@elementMetadata[,j])
  names(cases@elementMetadata) <- colnames(temp@elementMetadata)[j]
  
  #calc data fro normal at position 6 in the data set
  controls <- temp
  controls@elementMetadata <- DataFrame(controls@elementMetadata[,6])
  names(controls@elementMetadata) <- colnames(temp@elementMetadata)[6]
  
  #perform cn.mops
  res <- cn.mops::referencecn.mops(cases = cases, controls = controls,norm = 0,I=I,classes=classes,lowerThreshold = -0.25,upperThreshold = 0.25)
  
  #convert the cn.mops results
  res_cn <- cn.mops::calcIntegerCopyNumbers(res)
  res_cn_cnvs <- as.data.frame(cnvs(res_cn))
  res_cn_cnvs <- res_cn_cnvs[order(res_cn_cnvs$sampleName, res_cn_cnvs$seqnames),]
  
  #save the cn.mops results in a list
  cnvs_cnmops[[j]] <- res_cn_cnvs
}