#TenFoldCVTimesTen_VirtualToxicityPanel_MMGBSA.R
#Requires the R package 'glmnet'

print("Reading in datasets...")

#Specify filepath in quotes
filepath <- "<your_filepath>"

#Read in GB re-scored CSV file as a data frame
dockScores <- read.csv(paste(filepath,"N560_p33_MMGBSAScore_Matrix.csv",sep=""),header=T)
dockScores <- as.matrix(dockScores[,5:ncol(dockScores)])

#Read in ADR group response matrix
Y <- read.csv(paste(filepath,"N560_p10_DrugADR_Matrix.csv",sep=""),header=T)

#Set binding threshold for MM/GBSA scores
#dockThresh <- -25.0

print("...finished!")
print("Generating the design matrix...")
#Generate the design matrix
dataMatrix <- matrix(nrow=nrow(dockScores),ncol=16)
i <- 0
while(i<nrow(dockScores)) {
   i <- i+1
   dataMatrix[i,1] <- dockScores[i,1]
   dataMatrix[i,2] <- (dockScores[i,2]+dockScores[i,3]+dockScores[i,4]+dockScores[i,5]+dockScores[i,6]+dockScores[i,7]+dockScores[i,8]+
                       dockScores[i,9]+dockScores[i,10]+dockScores[i,11]+dockScores[i,12])/11
   dataMatrix[i,3] <- (dockScores[i,13]+dockScores[i,14]+dockScores[i,15]+dockScores[i,16]+dockScores[i,17]+dockScores[i,18]+
                       dockScores[i,19]+dockScores[i,20])/8
   dataMatrix[i,4] <- dockScores[i,21]
   dataMatrix[i,5] <- dockScores[i,22]
   dataMatrix[i,6] <- dockScores[i,23]
   dataMatrix[i,7] <- dockScores[i,24]
   dataMatrix[i,8] <- dockScores[i,25]
   dataMatrix[i,9] <- dockScores[i,26]
   dataMatrix[i,10] <- dockScores[i,27]
   dataMatrix[i,11] <- dockScores[i,28]
   dataMatrix[i,12] <- dockScores[i,29]
   dataMatrix[i,13] <- dockScores[i,30]
   dataMatrix[i,14] <- dockScores[i,31]
   dataMatrix[i,15] <- dockScores[i,32]
   dataMatrix[i,16] <- dockScores[i,33]
#   j <- 0
#   while(j<16) {
#      j <- j+1
#      if (dataMatrix[i,j]<dockThresh) {
#         dataMatrix[i,j] <- 1
#      } else {
#         dataMatrix[i,j] <- 0
#      } 
#   }
}
dataMatrixNew <- matrix(nrow=nrow(dataMatrix),ncol=ncol(dataMatrix))
i <- 0
while(i<nrow(dataMatrix)) {
   i <- i+1
   dataMatrixNew[i,] <- dataMatrix[i,]
}
rm(dataMatrix)
dataMatrix <- dataMatrixNew
rm(dataMatrixNew)
dataMatrixNew <- matrix(nrow=nrow(dataMatrix),ncol=ncol(dataMatrix))
j <- 0
kcnt2 <- 0
while(j<ncol(dataMatrix)) {
   j <- j+1
   if (sum(dataMatrix[,j])!=0) {
#      print(colnames(dockScores)[[j]])
       kcnt2 <- kcnt2+1
#       colMean <- mean(dataMatrix[,j])
#       colSD <- sd(dataMatrix[,j])
#       i <- 0
#       while(i<length(dataMatrix[,j])) {
#          i <- i+1
#          zScoreNow <- (dataMatrix[i,j]-colMean)/colSD
#          if (zScoreNow<(-1)) {dataMatrixNew[i,kcnt2] <- 1} else {dataMatrixNew[i,kcnt2] <- 0}
#       }
       dockThresh <- as.double(quantile(dataMatrix[,j],0.1))
#       dataMatrixNew[,kcnt2] <- dataMatrix[,j]
      i <- 0
      while(i<length(dataMatrix[,j])) {
         i <- i+1
         if (dataMatrix[i,j]<=dockThresh) dataMatrixNew[i,kcnt2] <- 1
         if (dataMatrix[i,j]>dockThresh) dataMatrixNew[i,kcnt2] <- 0
      }
   }
}
rm(dataMatrix)
dataMatrix <- dataMatrixNew
rm(dataMatrixNew)
dataMatrixNew <- matrix(nrow=nrow(dataMatrix),ncol=kcnt2)
j <- 0
while(j<kcnt2) {
   j <- j+1
   dataMatrixNew[,j] <- dataMatrix[,j]
}
rm(dataMatrix)
dataMatrix <- dataMatrixNew
rm(dataMatrixNew)
print("...finished!")

print("Begin 10-fold CV...")
kClass <- 10
Y <- as.matrix(Y[,5:ncol(Y)])
adrNameString <- c("neoplasms","bloodAndLymph","immuneSystem","endocrineDisorders","psychDisorders","cardiacDisorders","vascularDisorders","gastroDisorders","hepatoDisorders","renalDisorders")
varNameString <- c("lambda_min","cvlo","cvm","cvup","nzero")
tenAUCMatrix <- matrix(nrow=10,ncol=kClass)
tenCVResultsMatrix <- matrix(nrow=kClass,ncol=5)
lambdaMinVect <- vector(length=10)
cvloVect <- vector(length=10)
cvmVect <- vector(length=10)
cvupVect <- vector(length=10)
nzeroVect <- vector(length=10)
betas <- matrix(nrow=10,ncol=(ncol(dataMatrix)+1))
i <- 0
while(i<kClass) {
   i <- i+1
   print("Beginning 10-fold CV cycle:")
   j <- 0
   while(j<10) {
      j <- j+1
      print(j)
      myglmnetcv <- cv.glmnet(x=dataMatrix,y=Y[,i],family="binomial",type.measure="auc",nfolds=10)
      minLambdaIndx <- which(myglmnetcv$lambda==myglmnetcv$lambda.min)
      lambdaMinVect[[j]] <- myglmnetcv$lambda.min
      cvloVect[[j]] <- myglmnetcv$cvlo[[minLambdaIndx]]
      cvmVect[[j]] <- myglmnetcv$cvm[[minLambdaIndx]]
      cvupVect[[j]] <- myglmnetcv$cvup[[minLambdaIndx]]
      nzeroVect[[j]] <- myglmnetcv$nzero[[minLambdaIndx]]
      betaCoeff <- as.vector(t(coef(myglmnetcv,s="lambda.min")))
      l <- 0;while(l<length(betaCoeff)) {l <- l+1;betas[j,l] <- betaCoeff[[l]]}
      tenAUCMatrix[j,i] <- myglmnetcv$cvm[[minLambdaIndx]]
   }
   medianDiffVect <- abs(median(cvmVect)-cvmVect)
   chosenOneIndx <- which(medianDiffVect==min(medianDiffVect))
   if (length(chosenOneIndx)>1) chosenOneIndx <- chosenOneIndx[[1]]
   tenCVResultsMatrix[i,1] <- lambdaMinVect[[chosenOneIndx]]
   tenCVResultsMatrix[i,2] <- cvloVect[[chosenOneIndx]]
   tenCVResultsMatrix[i,3] <- cvmVect[[chosenOneIndx]]
   tenCVResultsMatrix[i,4] <- cvupVect[[chosenOneIndx]]
   tenCVResultsMatrix[i,5] <- nzeroVect[[chosenOneIndx]]
   betaVect <- vector(length=length(betas[chosenOneIndx,]))
   nn <- 0;while(nn<length(betas[chosenOneIndx,])) {nn <- nn+1;betaVect[[nn]] <- betas[chosenOneIndx,nn]}
   betaVect <- as.character(betaVect)
   dumString <- adrNameString[[i]]
   dumString <- paste(dumString,"_Betas.txt",sep="")
   fileString <- paste(filepath,dumString,sep="")
   write(x=betaVect,file=fileString)
}
colnames(tenAUCMatrix) <- adrNameString
xdat1 <- data.frame(tenAUCMatrix)
write.csv(x=xdat1,file=paste(filepath,"tenAUCMatrix.csv",sep=""))
rownames(tenCVResultsMatrix) <- adrNameString
colnames(tenCVResultsMatrix) <- varNameString
xdat2 <- data.frame(tenCVResultsMatrix)
write.csv(x=xdat2,file=paste(filepath,"tenCVResultsMatrix.csv",sep=""))
print("...finished!")