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

print("Reading in datasets...")

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

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

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

print("...finished!")

print("Generating the design matrix...")
#Generate the design matrix
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
   }
}
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!")