library(openxlsx)

library(plyr)

library(dplyr)

pt = read.xlsx("MEGA PT Data for Lu%2c August 21%2c 2017.xlsx", sheet = 1)

pt$mutID = paste(pt$Chrom, pt$Position,pt$BaseFrom,pt$BaseTo)

pt$Matched.Plasma = as.character(pt$Matched.Plasma)

cancer = c("CRC","Lung","Breast","Pancreas","Ovarian","Esophagus","Liver","Stomach","Small Intestine","Gastric","Ovary","pancreas")

uidThr=200

## threshold for omega value

thres = 1.9

final_cv = data.frame()

fp_cv = data.frame()

for(m in 1:10){

  allResults = c()

  cv = matrix(0, nrow = 10, ncol = 3)

  colnames(cv) = c("FP","Sens", "Con")

  for(c in 1:10){

    load(paste0("PvalueRatio1209cv_5perc_bl_",m,"_", c,"_1208data.rda"))

    result$fail = (result$UID1<uidThr)+(result$UID2<uidThr)+(result$UID3<uidThr)+

      (result$UID4<uidThr)+(result$UID5<uidThr)+(result$UID6<uidThr)

    result = result[result$fail<=2,]

    ## difference between average MAF in the test and the max MAF in the normal controls

    result$diff = result$aveMAF-result$max

    result$diff_r = (result$aveMAF-result$max)/(result$max+10e-6)

    result$PT.Avg.MAF[is.na(result$PT.Avg.MAF)] = -1

    result$PT.Avg.MAF = as.numeric(result$PT.Avg.MAF)

    ratio = matrix(c(result$r1, result$r2,result$r3,result$r4, result$r5,result$r6), nrow = nrow(result), ncol = 6)

    uid = matrix(c(result$UID1,result$UID2, result$UID3,result$UID4,result$UID5,result$UID6), nrow = nrow(result), ncol = 6)

    order.ratio = t(apply(ratio, 1, order))

    for(i in 1:nrow(ratio)){

      ratio[i,] = ratio[i,][order.ratio[i,]]

      uid[i,] = uid[i,][order.ratio[i,]]

    }

    uid[uid<uidThr] = 0

    ## eliminate the min and max wells in the testSet

    result$omega = rowSums(log(ratio[,-c(1,6)])*uid[,-c(1,6)])/rowSums(uid[,-c(1,6)])

    result$omega[is.na(result$omega)] = Inf

    result$class = FALSE

    result$class = (result$omega>=thres)

    result$iteration = m

    result$fold = c

    ## normal plasmas in the testSet

    nltemp = setdiff(nlpls,nlt)

    ## summarizing results

    summ = ddply(result, .(Template,Sample.Category), summarise, PT = max(PT.Avg.MAF*class), class = max(class))

    cv[c,1] = sum(summ$class[summ$Template %in% nltemp])/length(nltemp)

    cv[c,2] = sum(summ$class[summ$Sample.Category %in% cancer])/nrow(summ[summ$Sample.Category %in% cancer,])

    check = summ[summ$class==TRUE & (summ$Sample.Category %in% cancer) & (summ$Template %in% pt$Matched.Plasma),]

    cv[c,3] = sum(check$PT>=1)/sum(check$PT>=0)

    allResults = rbind(allResults,result)

  }

  save(allResults, file = paste0("allResults_fromLuMethod_",m,"_20171209.rda"))

  print(m)

  cv = data.frame(cv)

  final_cv = rbind.data.frame(final_cv,cv)

}

sampleTable_ALL = data.frame()

for(m in 1:10){

  load(file = paste0("allResults_fromLuMethod_",m,"_20171209.rda"))

  # subset by sample omega values and mutations

  sampleOmegaList = tapply(allResults[,'omega'], INDEX = allResults[,'Template'], FUN = function(x)return(x))

  sampleMutList = tapply(allResults[,'mutID'], INDEX = allResults[,'Template'], FUN = function(x)return(x))

  sampleCosmicList = tapply(allResults[,'CosmicCount'], INDEX = allResults[,'Template'], FUN = function(x)return(x))

  sampleGeneList = tapply(allResults[,'Gene'], INDEX = allResults[,'Template'], FUN = function(x)return(x))

  sampleAmpList = tapply(allResults[,"ampMatchName"], INDEX = allResults[,'Template'], FUN = function(x)return(x))

  # find max value

  sampleMaxS = unlist(tapply(allResults[,'omega'], INDEX = allResults[,'Template'], FUN = which.max, simplify = T))

  maxS = sapply(names(sampleMaxS), function(i) sampleOmegaList[[i]][sampleMaxS[i]])

  maxSmut = sapply(names(sampleMaxS), function(i) sampleMutList[[i]][sampleMaxS[i]])

  cosm = sapply(names(sampleMaxS), function(i) sampleCosmicList[[i]][sampleMaxS[i]])

  gene = sapply(names(sampleMaxS), function(i) sampleGeneList[[i]][sampleMaxS[i]])

  amp = sapply(names(sampleMaxS), function(i) sampleAmpList[[i]][sampleMaxS[i]])

  sampleTable = allResults[,c('Template','Sample.Category',"iteration","fold")]

  sampleTable = sampleTable[!duplicated(sampleTable[,1]),]

  rownames(sampleTable) = sampleTable[,1]

  sampleTable = cbind(sampleTable,maxOmega = maxS[rownames(sampleTable)], maxO_mut = maxSmut[rownames(sampleTable)],

                      CosmicCount = cosm[rownames(sampleTable)], gene = gene[rownames(sampleTable)], ampMatchName = amp[rownames(sampleTable)])

  sampleTable$remove = paste(sampleTable$Template,sampleTable$maxO_mut)

  sampleTable = sampleTable[,-1]

  uid_matrix = allResults[,c(1:51)]

  uid_matrix_subset = uid_matrix[match(sampleTable$remove,uid_matrix$remove),]

  sampleTable = cbind(sampleTable,uid_matrix_subset)

  #save(sampleTable, file = 'omegaPerSample_20171019.rda')

  # add maximum diff and diff_r per sample

  diff = unlist(tapply(allResults[,'diff'], INDEX = allResults[,'Template'], FUN = max, na.rm = T, simplify = T))

  diff_r = unlist(tapply(allResults[,'diff_r'], INDEX = allResults[,'Template'], FUN = max, na.rm = T, simplify = T))

  sampleTable_diff = cbind(sampleTable,maxDiff = diff[rownames(sampleTable)],maxDiff_r = diff_r[rownames(sampleTable)])

  sampleTable_ALL = rbind.data.frame(sampleTable_ALL, sampleTable_diff)

}

save(sampleTable_ALL, file = 'maxValuesPerSample_20171209_FORJosh.rda')