library(gplots)
library(RColorBrewer)
### this function is used to calcuate the fisher exact test within a matrix
source("./script/fisher.exact.test.combined.R")
### this function is used to generate the heatmap figures 
source("./script/heatmap.3.modified.R")
#### this function is used to read the feature table and obtain the copy number and mutation info and then do preprocess the matrix, 
#####calculate the fisher exact test and then draw the covarance figure

#### function to make the copy number data as numerical
nuch <- function(x) {return(as.numeric(as.character(x)))}

#### function to pull out mutation data from the feature table and plot the co-occuring figure of somatic muation genes
creatCNfocalfigure <- function(feature.file,tumor,cut.focal=0.1,fold.focal.gain=0.4,fold.focal.loss=-0.4,maxevent=5){
  ## read the feature file 
  feature.table.original <- read.table(feature.file,header=T,comment='',sep='\t',as.is=T,na.strings='NA',stringsAsFactors=F,check.names=F)
  ## create the feature table only with samples and features
  sample.index <- nchar(feature.table.original[,1]) == 12
  feature.table <- feature.table.original[sample.index,]
  ## pull out q value row
  t <- grep("all_q",feature.table.original[,1])
  if(length(t) > 0) {
    feature.qval <- feature.table.original[t,]
  }else {
    return(NULL)
  }
  
  sample.name <- feature.table[,1]
  rownames(feature.table) <- sample.name
  rownames(feature.qval) <- feature.qval[,1]
  
  ##### preprocess focal copy number section
  amp.peak <- feature.table[,grep("Amp_",colnames(feature.table))]
  amp.peak <- amp.peak[,!grepl("mRNA",colnames(amp.peak))]
  del.peak <- feature.table[,grep("Del_",colnames(feature.table))]
  del.peak <- del.peak[,!grepl("mRNA",colnames(del.peak))]
  focal.copy <- cbind(amp.peak,del.peak)
  focal.qval <- feature.qval[,match(colnames(focal.copy),colnames(feature.qval),nomatch=0)]
  
  ##### pick significant copy number arm events data 
  focal.copy <- focal.copy[,focal.qval <= cut.focal]
  
  raw.focal.copy <- t(focal.copy)
  
  ##### calculate the absolute value of  copy number threshold 
  cut.focal.gain <- 2*2^fold.focal.gain-2
  cut.focal.loss <- 2*2^(fold.focal.loss)-2
  
  ######## charactrized the copy number data
  ann.focal.copy <- apply(raw.focal.copy,1:2,function(x) {
    if (is.na(x)) {
      return(5)
    } else if (nuch(x) >= cut.focal.gain) {
      return(2)
    } else if (nuch(x) <= cut.focal.loss) {
      return(1)
    } else {
      return(0)
    }
  }
  )
  
  rownames(ann.focal.copy) <- rownames(raw.focal.copy)
  colnames(ann.focal.copy) <- colnames(raw.focal.copy)
  
  ########## calculate the total number of amplification and deletion seperately for each event and then only pick the events with at least 5 alterations
  sum.amp <- apply(ann.focal.copy,1,function(x) {
    x <- x[!is.na(x)]
    return(sum(x %in% c(3,4)))
  }
  )
  sum.del <- apply(ann.focal.copy,1,function(x) {
    x <- x[!is.na(x)]
    return(sum(x %in% c(1,2)))
  }
  )
  sum.focal <- data.frame(sum.amp,sum.del)
  colnames(sum.focal) <- c("amp","del")
  sum.focal[,"total"] <- rowSums(sum.focal)
  sum.focal[,"max"] <- apply(sum.focal,1,function(x) ifelse(x[1] >= x[2],x[1],x[2]))
  ann.focal.copy <- ann.focal.copy[sum.focal$max >=maxevent,]
  
  ########### remove samples with all missing values############## 
  index3 <- colSums(ann.focal.copy == 5) != nrow(ann.focal.copy)
  copy.focal <- data.frame(ann.focal.copy,check.names=F)
  copy.focal <- copy.focal[,index3]
  
  ########### create the event matrix 1:alterations, 0: no copy number change for cooccurrence figure
  tt <- grep("Amp",rownames(copy.focal))
  if(length(tt) > 0){
    copy.focal.amp <- copy.focal[tt,]
    copy.focal.del <- copy.focal[-tt,]
    event.focal.amp <- apply(copy.focal.amp,1:2,function(x) ifelse(x == 2,1,0))
    event.focal.del <- apply(copy.focal.del,1:2,function(x) ifelse(x == 1,1,0))
    event.focal <- rbind(event.focal.amp,event.focal.del)
  }else{
    event.focal <- apply(copy.focal,1:2,function(x) ifelse(x == 1,1,0))
  }
  
  #+++++++ calculate correlation based on event matrix +++++++++++
  cut.sample <- 0.01
  cut.heatmap.default <- 4
  ## fisher exact test is used to calculate the association withon two pair of mutations: the pvalue with exclusiveness and cocurrence will be returned.
  fisher <- fisher.simple.test(t(event.focal))
  maxp <- max(-log10(unlist(fisher)))
  cut.heatmap <- ifelse(maxp >= cut.heatmap.default,cut.heatmap.default,round(maxp))
  
  ##### draw exclusiveness and cocurrence figure on copy number focal event###########
  png(file=paste(tumor,'fisher.focal.copy.correlation.png',sep="."),width=1000,height=1000)
  get.mutual.heatmap(fisher[[1]],fisher[[2]],colnames(fisher[[1]]),cut.heatmap,"Correlations in Focal SCNAs",F,F,50)
  dev.off()
  return(event.focal)
}

#creatCNfocalfigure(feature.file="./CESC-TP.samplefeatures.txt",tumor="CESC-TP",cut.focal=0.1,fold.focal.gain=0.4,fold.focal.loss=-0.4,maxevent=5)