--- a
+++ b/Analysis_2.R
@@ -0,0 +1,127 @@
+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')
+
+