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
creatCNarmfigure <- function(feature.file,tumor,cut.arm=0.1,fold.arm.gain=0.4,fold.arm.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 for arm copy number
  
  arm.copy <- feature.table[,grep("CN",colnames(feature.table))]
  arm.copy <- arm.copy[,!grepl("CLUS",colnames(arm.copy))]
  arm.copy <- arm.copy[,!grepl("mRNA",colnames(arm.copy))]
  arm.copy <- arm.copy[,!grepl("SMG",colnames(arm.copy))]
  arm.qval <- feature.qval[,match(colnames(arm.copy),colnames(feature.qval),nomatch=0)]
  colnames(arm.copy) <- as.character(sapply(colnames(arm.copy),function(x) paste(strsplit(x,"_")[[1]][2],strsplit(x,"_")[[1]][3],sep="-")))
  
  ##### pick significant copy number arm events data 
  arm.copy <- arm.copy[,arm.qval <= cut.arm]
  arm.qval <- arm.qval[arm.qval <= cut.arm]
  
  raw.arm.copy <- t(arm.copy)
  ##### calculate the absolute value of  copy number threshold 
  cut.arm.gain <- 2*2^fold.arm.gain-2
  cut.arm.loss <- 2*2^(fold.arm.loss)-2
  
  ######## charactrized the copy number data
  ann.arm.copy <- apply(raw.arm.copy,1:2,function(x) {
    if (is.na(x)) {
      return(5)
    } else if (nuch(x) >= cut.arm.gain) {
      return(2)
    } else if (nuch(x) <= cut.arm.loss) {
      return(1)
    } else {
      return(0)
    }
  }
  )
  rownames(ann.arm.copy) <- rownames(raw.arm.copy)
  colnames(ann.arm.copy) <- colnames(raw.arm.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.arm.copy,1,function(x) {
    x <- x[!is.na(x)]
    return(sum(x==2))
  }
  )
  sum.del <- apply(ann.arm.copy,1,function(x) {
    x <- x[!is.na(x)]
    return(sum(x==1))
  }
  )
  
  sum.arm <- data.frame(sum.amp,sum.del)
  colnames(sum.arm) <- c("amp","del")
  sum.arm[,"total"] <- rowSums(sum.arm)
  sum.arm[,"max"] <- apply(sum.arm,1,function(x) ifelse(x[1] >= x[2],x[1],x[2]))
  ann.arm.copy <- ann.arm.copy[sum.arm$max >= maxevent,]
  
  ########### remove samples with all missing values############## 
  index2 <- colSums(ann.arm.copy == 5) != nrow(ann.arm.copy)
  copy.arm.level <- data.frame(ann.arm.copy,check.names=F)
  copy.arm.level <- copy.arm.level[,index2]
  
  ########### create the event matrix 1:alterations, 0: no copy number change for cooccurrence figure
  t <- grep("Amp",rownames(copy.arm.level))
  if(length(t) > 0){
    arm.amp <- copy.arm.level[t,]
    arm.del <- copy.arm.level[-t,]
    event.arm.amp <- apply(arm.amp,1:2,function(x) ifelse(x == 2,1,0))
    event.arm.del <- apply(arm.del,1:2,function(x) ifelse(x == 1,1,0))
    event.arm <- rbind(event.arm.amp,event.arm.del)
  }else{
    event.arm <- apply(copy.arm.level,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.arm))
  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 arm event###########
  png(file=paste(tumor,'fisher.arm.copy.correlation.png',sep="."),width=1000,height=1000)
  get.mutual.heatmap(fisher[[1]],fisher[[2]],colnames(fisher[[1]]),cut.heatmap,"Correlations in Arm-level SCNAs",F,F,50)
  dev.off()
  return(event.arm)
}

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