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

print("Reading in datasets...")

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

#Read in VinaLC CSV file as a data frame
dockScores <- read.csv(paste(filepath,"N560_p409_DockScore_Matrix.csv",sep=""),header=T)

#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 <- -8.0

print("...finished!")

print("Generating the design matrix...")
#Generate a data matrix 'X' that consists of only the docking scores for drugs that are found in the SIDER dataset
dataMatrix <- matrix(nrow=nrow(dockScores),ncol=(ncol(dockScores)-4))
i <- 0
while(i<nrow(dockScores)) {
   i <- i+1
   j <- 4
   while(j<ncol(dockScores)) {
      j <- j+1
      scoreNow <- as.double(dockScores[i,j])  
      dataMatrix[i,j-4] <- scoreNow
#      if (scoreNow<=dockThresh) dataMatrix[i,j-4] <- 1
#      if (scoreNow>dockThresh) dataMatrix[i,j-4] <- 0
   }
}
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(dockScoreNames[[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<(-2)) {dataMatrixNew[i,kcnt2] <- 1} else {dataMatrixNew[i,kcnt2] <- 0}
       #}
       #      dockThresh <- as.double(quantile(dataMatrix[,j],0.05))
       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!")