PCM_Aromatase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Mtitle: Unraveling the origin of aromatase inhibitory activity via proteochemometrics\n",
    "  modeling\n",
    "author: \"Saw Simeon, Sunanta Nabu, Naravut Suvannang, Ola Spjuth, Maris Lapins, Virapong Prachayasittikul, Jar E.S. Wikberg, Chanin Nantasenamat\"\n",
    "date: \"September 18, 2015\"\n",
    "output: html_document\n",
    "---\n",
    "\n",
    "Initialize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "library(readxl); packageVersion(\"readxl\")\n",
    "library(prospectr); packageVersion(\"prospectr\")\n",
    "library(Rcpi); packageVersion(\"Rcpi\")\n",
    "library(caret); packageVersion(\"caret\")\n",
    "library(pls); packageVersion(\"pls\")\n",
    "library(genefilter); packageVersion(\"genefilter\")\n",
    "library(ggplot2); packageVersion(\"ggplot2\")\n",
    "library(grid); packageVersion(\"grid\")\n",
    "library(cowplot); packageVersion(\"cowplot\")\n",
    "library(ridge); packageVersion(\"ridge\")\n",
    "library(randomForest); packageVersion(\"randomForest\")\n",
    "library(XLConnect); packageVersion(\"XLConnect\")\n",
    "library(XLConnectJars); packageVersion(\"XLConnectJars\")\n",
    "library(RCurl); packageVersion(\"RCurl\")\n",
    "library(paran); packageVersion(\"paran\")\n",
    "library(dplyr); packageVersion(\"dplyr\")\n",
    "library(doSNOW); packageVersion(\"doSNOW\")\n",
    "library(doParallel); packageVersion(\"doParallel\")\n",
    "library(parallel); packageVersion(\"parallel\")\n",
    "library(gtools); packageVersion(\"gtools\")\n",
    "library(reshape); packageVersion(\"reshape\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "url <- \"https://raw.githubusercontent.com/Rnewbie/PCM_Aromatase/master/aromatase_data.xlsx\"\n",
    "f = CFILE(\"data.xlsx\", mode = \"wb\")\n",
    "curlPerform(url = url, writedata = f@ref, ssl.verifypeer = FALSE)\n",
    "close(f)\n",
    "out <- readWorksheetFromFile(file = \"data.xlsx\", sheet = \"Sheet1\", header = TRUE,\n",
    "                             startRow = 11, startCol = 2, endCol = 15, endRow = 35)\n",
    "head(out)\n",
    "data <- read_excel(\"data.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preapring the data set for 13 PCM models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "### preparing data for 13 models \n",
    "compound <- data[, 5:311]\n",
    "protein <- data[, 313:351]\n",
    "pIC50 <- data[2]\n",
    "set.seed(23)\n",
    "compound = compound[, -nearZeroVar(compound)]\n",
    "compound[compound == \"0\"] <- -1\n",
    "CxP <- getCPI(compound, protein, type = \"tensorprod\")\n",
    "CxP <- as.data.frame(CxP)\n",
    "dfcompound <- names(data.frame(compound[,1:39]))\n",
    "dfprotein <- names(data.frame(protein[,1:39]))\n",
    "compoundNamecross <- rep(dfcompound, each = 39)\n",
    "proteinNamecross <- rep(dfprotein, times = 39)\n",
    "label <- paste(compoundNamecross, proteinNamecross, sep=\"_\")\n",
    "colnames(CxP) <- label\n",
    "PxP <- getCPI(protein, protein, type = \"tensorprod\")\n",
    "proteinName2 <- rep(dfprotein, times = 39)\n",
    "proteinName1 <- rep(dfprotein, each = 39)\n",
    "label_protein <- paste(proteinName1, proteinName2, sep = \"_\")\n",
    "colnames(PxP) <- label_protein\n",
    "index <- seq(1, 1521, by = 40)\n",
    "protein_selfcross <- PxP[, -index]\n",
    "transposedIndexed_protein <- t(protein_selfcross)\n",
    "index1 <- which(duplicated(transposedIndexed_protein))\n",
    "removed_duplicated_protein <- transposedIndexed_protein[-index1, ]\n",
    "PxP <- t(removed_duplicated_protein)\n",
    "\n",
    "CxC <- getCPI(compound, compound, type = \"tensorprod\")\n",
    "compoundName2 <- rep(dfcompound, times = 39)\n",
    "compoundName1 <- rep(dfcompound, each = 39)\n",
    "label <- paste(compoundName1, compoundName2, sep = \"_\")\n",
    "colnames(CxC) <- label\n",
    "index3 <- seq(1, 1369, by = 38)\n",
    "compound_selfcross <- CxC[, -index3]\n",
    "transposedIndexed_compound <- t(compound_selfcross)\n",
    "index4 <- which(duplicated(transposedIndexed_compound))\n",
    "removed_compound <- transposedIndexed_compound[-index4, ]\n",
    "compound_finalcrossterms <- t(removed_compound)\n",
    "CxC <- compound_finalcrossterms\n",
    "\n",
    "C <- compound\n",
    "P <- protein\n",
    "C_P <- cbind(C, P)\n",
    "C_P_CxP_data_block_scale <- cbind(C, P, CxP) * (1/sqrt(length(C)+length(P)+length(CxP)))\n",
    "C_P_CxC_data_block_scale <- cbind(C, P, \n",
    "                                  CxC) * (1/sqrt(length(C)+length(P)+length(CxC)))\n",
    "C_P_PxP_data_block_scale <- cbind(C, P,\n",
    "                                  PxP) * (1/sqrt(length(C)+length(P)+length(PxP)))\n",
    "C_P_CxP_CxC_data_block_scale <- cbind(C, P, CxP,\n",
    "                                      CxC) * (1/sqrt(length(C)+length(P)+length(CxP)+length(CxC)))\n",
    "C_P_CxP_PxP_data_block_scale <- cbind(C, P, CxP,\n",
    "                                      PxP) * (1/sqrt(length(C)+length(P)+length(CxP)+length(PxP)))\n",
    "C_P_CxC_PxP_data_block_scale <- cbind(C, P, CxC,\n",
    "                                      PxP) * (1/sqrt(length(C)+length(P)+length(CxC)+length(PxP)))\n",
    "C_P_CxP_CxC_PxP_data_block_scale <- cbind(C, P, CxP, CxC, PxP) * (1/sqrt(length(C)+length(P)+\n",
    "                                                                           length(CxC)+length(PxP)))\n",
    "\n",
    "C <- cbind(pIC50, C)\n",
    "P <- cbind(pIC50, P)\n",
    "CxP <- cbind(pIC50, CxP)\n",
    "CxC <- cbind(pIC50, CxC)\n",
    "PxP <- cbind(pIC50, PxP)\n",
    "C_P <- cbind(pIC50, C_P)\n",
    "C_P_CxP <- cbind(pIC50, C_P_CxP_data_block_scale)\n",
    "C_P_CxC <- cbind(pIC50, C_P_CxC_data_block_scale)\n",
    "C_P_PxP <- cbind(pIC50, C_P_PxP_data_block_scale)\n",
    "C_P_CxP_CxC <- cbind(pIC50, C_P_CxP_CxC_data_block_scale)\n",
    "C_P_CxP_PxP <- cbind(pIC50, C_P_CxP_PxP_data_block_scale)\n",
    "C_P_CxC_PxP <- cbind(pIC50, C_P_CxC_PxP_data_block_scale)\n",
    "C_P_CxP_CxC_PxP <- cbind(pIC50, C_P_CxP_CxC_PxP_data_block_scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for trianing PLS Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "pls_training <- function(x){\n",
    "  ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50 ~., data = x,  method = \"pls\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    model <- plsr(pIC50~., data = Train, ncomp = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, Train, ncomp = tune$bestTune[[1]])\n",
    "    value <- data.frame(obs = Train$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "mean_and_sd <- function(x) {\n",
    "  c(round(rowMeans(x, na.rm = TRUE), digits = 4),\n",
    "    round(genefilter::rowSds(x, na.rm = TRUE), digits = 4))\n",
    "}\n",
    "\n",
    "pls_train <- function(x) {\n",
    "  ok <- pls_training(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for trianing PLS Model 10-fold cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "pls_cross_validation <- function(x){\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    myData <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    k = 10\n",
    "    index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "    folds <- 1:k\n",
    "    myRes <- data.frame()\n",
    "    ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    for (j in 1:k)\n",
    "      training <- subset(myData, index %in% folds[-j])\n",
    "    testing <- subset(myData, index %in% c(j))\n",
    "    tune <- train(pIC50 ~., data = training,  method = \"pls\", tuneLength = 10,\n",
    "                  trControl = ctrl)\n",
    "    model <- plsr(pIC50~., data = training, ncomp = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, testing, ncomp = tune$bestTune[[1]])\n",
    "    ok <- data.frame(obs = testing$pIC50, pred = prediction)\n",
    "    value <- rbind(myRes, ok)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "pls_10_CV <- function(x) {\n",
    "  ok <- pls_cross_validation(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for PLS Model testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "pls_cross_validation <- function(x){\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    myData <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    k = 10\n",
    "    index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "    folds <- 1:k\n",
    "    myRes <- data.frame()\n",
    "    ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    for (j in 1:k)\n",
    "      training <- subset(myData, index %in% folds[-j])\n",
    "    testing <- subset(myData, index %in% c(j))\n",
    "    tune <- train(pIC50 ~., data = training,  method = \"pls\", tuneLength = 10,\n",
    "                  trControl = ctrl)\n",
    "    model <- plsr(pIC50~., data = training, ncomp = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, testing, ncomp = tune$bestTune[[1]])\n",
    "    ok <- data.frame(obs = testing$pIC50, pred = prediction)\n",
    "    value <- rbind(myRes, ok)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "pls_10_CV <- function(x) {\n",
    "  ok <- pls_cross_validation(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for ridge Model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "ridge_training <- function(x){\n",
    "  library(ridge)\n",
    "  ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50 ~., data = x,  method = \"foba\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    model <- linearRidge(pIC50~., \n",
    "                         data = Train, scaling = \"none\",\n",
    "                         lambda = tune$bestTune$lambda[[1]])\n",
    "    prediction <- predict(model, Train)\n",
    "    value <- data.frame(obs = Train$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "mean_and_sd <- function(x) {\n",
    "  c(round(rowMeans(x, na.rm = TRUE), digits = 4),\n",
    "    round(genefilter::rowSds(x, na.rm = TRUE), digits = 4))\n",
    "}\n",
    "\n",
    "ridge_train <- function(x) {\n",
    "  ok <- ridge_training(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for ridge model 10-fold CV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "ridge_cross_validation <- function(x){\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    myData <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    k = 10\n",
    "    index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "    folds <- 1:k\n",
    "    myRes <- data.frame()\n",
    "    ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    for (j in 1:k)\n",
    "      training <- subset(myData, index %in% folds[-j])\n",
    "    testing <- subset(myData, index %in% c(j))\n",
    "    tune <- train(pIC50 ~., data = training,  method = \"foba\", tuneLength = 10,\n",
    "                  trControl = ctrl)\n",
    "    model <- linearRidge(pIC50~., scaling = \"none\", \n",
    "                         data = training, lambda = tune$bestTune$lambda[[1]])\n",
    "    prediction <- predict(model, testing)\n",
    "    ok <- data.frame(obs = testing$pIC50, pred = prediction)\n",
    "    value <- rbind(myRes, ok)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "ridge_10_CV <- function(x) {\n",
    "  ok <- ridge_cross_validation(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for ridge model testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "ridge_testing <- function(x){\n",
    "  ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50 ~., data = x,  method = \"foba\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    model <- linearRidge(pIC50~., data = Train, scaling = \"none\",\n",
    "                         lambda = tune$bestTune$lambda[[1]])\n",
    "    prediction <- predict(model, Test)\n",
    "    value <- data.frame(obs = Test$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "ridge_test <- function(x) {\n",
    "  ok <- ridge_testing(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for Random Forest model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "randomForest_training <- function(x){\n",
    "  library(randomForest)\n",
    "  ctrl <- trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50 ~., data = x,  method = \"rf\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "  results <- list(100)\n",
    "  for (i in 1:100) {\n",
    "    sel <- naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    model <- randomForest(pIC50~., \n",
    "                         data = Train,\n",
    "                         mtry = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, Train)\n",
    "    value <- data.frame(obs = Train$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "}\n",
    "\n",
    "mean_and_sd <- function(x) {\n",
    "  c(round(rowMeans(x, na.rm = TRUE), digits = 4),\n",
    "    round(genefilter::rowSds(x, na.rm = TRUE), digits = 4))\n",
    "}\n",
    "\n",
    "randomForest_train <- function(x) {\n",
    "  ok <- randomForest_training(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for Random Forest 10-fold Cross Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "randomForest_cross_validation <- function(x){\n",
    "  library(doSNOW)\n",
    "  library(foreach)\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(8)\n",
    "  registerDoSNOW(cl)\n",
    "  \n",
    "  results <- list(100)\n",
    "  results <- foreach (i=1:100) %dopar% {\n",
    "    sel <- prospectr::naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    myData <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    k = 10\n",
    "    index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "    folds <- 1:k\n",
    "    myRes <- data.frame()\n",
    "    ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    for (j in 1:k)\n",
    "      training <- subset(myData, index %in% folds[-j])\n",
    "    testing <- subset(myData, index %in% c(j))\n",
    "    tune <- caret::train(pIC50 ~., data = training,  method = \"rf\", tuneLength = 10,\n",
    "                  trControl = ctrl)\n",
    "    model <- randomForest::randomForest(pIC50~., scaling = \"none\", \n",
    "                         data = training, mtry = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, testing)\n",
    "    ok <- data.frame(obs = testing$pIC50, pred = prediction)\n",
    "    value <- rbind(myRes, ok)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- caret::defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "  stopCluster(cl)\n",
    "}\n",
    "\n",
    "randomForest_10_CV <- function(x) {\n",
    "  ok <- randomForest_cross_validation(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function for Random Forest external testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "randomForest_testing <- function(x){\n",
    "  ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- caret::train(pIC50 ~., data = x,  method = \"rf\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "  library(doSNOW)\n",
    "  library(foreach)\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(8)\n",
    "  registerDoSNOW(cl)\n",
    "  \n",
    "  results <- list(100)\n",
    "  results <- foreach (i=1:100) %dopar% {\n",
    "    sel <- prospectr::naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    model <- randomForest::randomForest(pIC50~., data = Train,\n",
    "                         mtry = tune$bestTune[[1]])\n",
    "    prediction <- predict(model, Test)\n",
    "    value <- data.frame(obs = Test$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- caret::defaultSummary(value)\n",
    "  }\n",
    "  return(results)\n",
    "  stopCluster(cl)\n",
    "}\n",
    "\n",
    "randomForest_test <- function(x) {\n",
    "  ok <- randomForest_testing(x)\n",
    "  data <- data.frame(ok)\n",
    "  result <- mean_and_sd(data)\n",
    "  df <- data.frame(result)\n",
    "  R2_and_RMSE <- t(df)\n",
    "  label <- c(\"RMSE_Mean\", \"Rsquared_Mean\", \"RMSE_SD\", \"Rsquared_SD\")\n",
    "  colnames(R2_and_RMSE) <- label\n",
    "  return(R2_and_RMSE)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Component Analysis Score Plot (Compound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "data <- read_excel(\"data.xlsx\")\n",
    "compound <- data[, 5:311]\n",
    "set.seed(122)\n",
    "compound = compound[, -nearZeroVar(compound)]\n",
    "compoundname <- as.character(unique(data$Compound))\n",
    "dataframe <- cbind(Compound = data$Compound, compound)\n",
    "sample_compound <- function(x) {\n",
    "  myRef  <- data.frame()\n",
    "  for (i in compoundname) {\n",
    "    index_data <- which(x$Compound == i)\n",
    "    ok <- x[index_data, ]\n",
    "    sample <- sample_n(ok, size = 1, replace = TRUE)\n",
    "    myRef <- rbind(myRef, sample)\n",
    "  }\n",
    "  return(myRef)\n",
    "}\n",
    "df <- sample_compound(dataframe)[, 2:length(dataframe)]\n",
    "set.seed(2000)\n",
    "pca <- prcomp(df, retx=TRUE,scale.=TRUE)\n",
    "scores <- pca$x[,1:5]\n",
    "loadings <- pca$rotation[,1:5]\n",
    "km <- kmeans(scores, center=2, nstart=5)\n",
    "ggdata <- data.frame(scores, Cluster=km$cluster)\n",
    "set.seed(23)\n",
    "x <- ggplot(ggdata, aes(x = PC1, y = PC2, colour = Cluster)) +\n",
    "  geom_point(aes(fill=factor(Cluster)), size=5, shape=20, pch = 21, alpha = 0.8) +\n",
    "  ggtitle(\" \") +\n",
    "  stat_ellipse(aes(fill=factor(Cluster)), colour = \"black\", \n",
    "               geom=\"polygon\", level=0.95, alpha=0.2) +\n",
    "  guides(color=guide_legend(\"Cluster\"),fill=guide_legend(\"Cluster\")) +\n",
    "  #geom_text(aes(label=compoundnumber), size=7, hjust=0.5, vjust= 1.5, alpha=0.45) +\n",
    "  theme(\n",
    "    legend.position=(\"none\"),\n",
    "    #plot.title = element_text(size=20, face=\"bold\", colour=\"black\", vjust = 2, hjust=-0.07),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1),\n",
    "    axis.text.y = element_text(size = 15),\n",
    "    axis.ticks.length = unit(0.3, \"cm\"),\n",
    "    axis.text.x = element_text(size = 15),\n",
    "    legend.title=element_blank(),\n",
    "    axis.title.x = element_text(color=\"black\", size=20),\n",
    "    axis.title.y = element_text(color=\"black\", size=20)) +\n",
    "  coord_cartesian(ylim = c(-15, 15), xlim = c(-15, 15))\n",
    "print(x)\n",
    "# MUST MANUALLY SAVE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Component Analysis Loading Plot (Compound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "data <- read_excel(\"data.xlsx\")\n",
    "compound <- data[, 5:311]\n",
    "set.seed(122)\n",
    "compound = compound[, -nearZeroVar(compound)]\n",
    "\n",
    "compoundname <- as.character(unique(data$Compound))\n",
    "dataframe <- cbind(Compound = data$Compound, compound)\n",
    "sample_compound <- function(x) {\n",
    "  myRef  <- data.frame()\n",
    "  for (i in compoundname) {\n",
    "    index_data <- which(x$Compound == i)\n",
    "    ok <- x[index_data, ]\n",
    "    sample <- sample_n(ok, size = 1, replace = TRUE)\n",
    "    myRef <- rbind(myRef, sample)\n",
    "  }\n",
    "  return(myRef)\n",
    "}\n",
    "df <- sample_compound(dataframe)[, 2:length(dataframe)]\n",
    "set.seed(2000)\n",
    "pca <- prcomp(df, retx=TRUE,scale.=TRUE)\n",
    "scores <- pca$x[,1:5]\n",
    "set.seed(300)\n",
    "loadings <- pca$rotation[,1:5]\n",
    "km <- kmeans(loadings, center=3, nstart=5)\n",
    "ggdata <- data.frame(loadings, Cluster=km$cluster)\n",
    "labelcompound <- rownames(loadings)\n",
    "ggdata <- cbind(labelcompound, ggdata)\n",
    "### paper numbering\n",
    "library(grid)\n",
    "set.seed(23)\n",
    "a <- ggplot(ggdata, aes(x = PC1, y = PC2, colour = Cluster)) +\n",
    "  geom_point(aes(fill=factor(Cluster)), size=5, shape=20, pch = 21, alpha = 0.8) +\n",
    "  ggtitle(\" \") +\n",
    "  stat_ellipse(aes(fill=factor(Cluster)), colour = \"black\", \n",
    "               geom=\"polygon\", level=0.95, alpha=0.2) +\n",
    "  guides(color=guide_legend(\"Cluster\"),fill=guide_legend(\"Cluster\")) +\n",
    "  #geom_text(aes(label=labelcompound), size=7, hjust=0.5, vjust= 1.5, alpha=0.45) +\n",
    "  theme(\n",
    "    legend.position=(\"none\"),\n",
    "    #plot.title = element_text(size=20, face=\"bold\", colour=\"black\", vjust = 2, hjust=-0.07),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1),\n",
    "    axis.text.y = element_text(size = 15),\n",
    "    axis.ticks.length = unit(0.3, \"cm\"),\n",
    "    axis.text.x = element_text(size = 15),\n",
    "    legend.title=element_blank(),\n",
    "    axis.title.x = element_text(color=\"black\", size=20),\n",
    "    axis.title.y = element_text(color=\"black\", size=20)) +\n",
    "  coord_cartesian(ylim = c(-0.8, 0.8), xlim = c(-0.8, 0.8))\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Component Analysis Score Plot (Protein)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "data <- read_excel(\"data.xlsx\")\n",
    "protein <- data[, 313:351]\n",
    "set.seed(200)\n",
    "proteinname <- as.character(unique(data$Protein))\n",
    "dataframe <- cbind(Protein = data$Protein, protein)\n",
    "sample_protein <- function(x) {\n",
    "  myRef  <- data.frame()\n",
    "  for (i in proteinname) {\n",
    "    index_data <- which(x$Protein == i)\n",
    "    ok <- x[index_data, ]\n",
    "    sample <- sample_n(ok, size = 1, replace = TRUE)\n",
    "    myRef <- rbind(myRef, sample)\n",
    "  }\n",
    "  return(myRef)\n",
    "}\n",
    "df <- sample_protein(dataframe)[, 2:length(dataframe)]\n",
    "set.seed(100)\n",
    "df = df[, -nearZeroVar(df)]\n",
    "pca <- prcomp(df, retx=TRUE,scale.=TRUE)\n",
    "scores <- pca$x[,1:5]\n",
    "loadings <- pca$rotation[,1:5]\n",
    "km <- kmeans(scores, center=1, nstart=5)\n",
    "ggdata <- data.frame(scores, Cluster=km$cluster)\n",
    "ggdata <- cbind(proteinname, ggdata)\n",
    "### paper numbering\n",
    "library(grid)\n",
    "set.seed(2300)\n",
    "y <- ggplot(ggdata, aes(x = PC1, y = PC2, colour = Cluster)) +\n",
    "  geom_point(aes(fill=factor(Cluster)), size=5, shape=20, pch = 21, alpha = 0.8) +\n",
    "  ggtitle(\" \") +\n",
    "  stat_ellipse(aes(fill=factor(Cluster)), colour = \"black\", \n",
    "               geom=\"polygon\", level=0.95, alpha=0.2) +\n",
    "  guides(color=guide_legend(\"Cluster\"),fill=guide_legend(\"Cluster\")) +\n",
    "  #geom_text(aes(label=proteinname), size=7, hjust=0.5, vjust= 1.5, alpha=0.45) +\n",
    "  theme(\n",
    "    legend.position=(\"none\"),\n",
    "   # plot.title = element_text(size=20, face=\"bold\", colour=\"black\", vjust = 2, hjust=-0.07),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1),\n",
    "    axis.text.y = element_text(size = 15),\n",
    "    axis.ticks.length = unit(0.3, \"cm\"),\n",
    "    axis.text.x = element_text(size = 15),\n",
    "    legend.title=element_blank(),\n",
    "    axis.title.x = element_text(color=\"black\", size=20),\n",
    "    axis.title.y = element_text(color=\"black\", size=20)) +\n",
    "  coord_cartesian(ylim = c(-5, 5), xlim = c(-5, 5))\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Component Analysis Loading Plot (Protein)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "### loading plot\n",
    "protein <- data[, 313:351]\n",
    "set.seed(200)\n",
    "proteinname <- as.character(unique(data$Protein))\n",
    "dataframe <- cbind(Protein = data$Protein, protein)\n",
    "sample_protein <- function(x) {\n",
    "  myRef  <- data.frame()\n",
    "  for (i in proteinname) {\n",
    "    index_data <- which(x$Protein == i)\n",
    "    ok <- x[index_data, ]\n",
    "    sample <- sample_n(ok, size = 1, replace = TRUE)\n",
    "    myRef <- rbind(myRef, sample)\n",
    "  }\n",
    "  return(myRef)\n",
    "}\n",
    "df <- sample_protein(dataframe)[, 2:length(dataframe)]\n",
    "set.seed(100)\n",
    "df = df[, -nearZeroVar(df)]\n",
    "column_names <- colnames(df)\n",
    "label_names2 <- gsub(\"zscl\", \"z\", column_names)\n",
    "colnames(df) <- label_names2\n",
    "pca <- prcomp(df, retx=TRUE,scale.=TRUE)\n",
    "scores <- pca$x[,1:5]\n",
    "loadings <- pca$rotation[,1:5]\n",
    "km <- kmeans(loadings, center=2, nstart=5)\n",
    "ggdata <- data.frame(loadings, Cluster=km$cluster)\n",
    "labelprotein <- rownames(loadings)\n",
    "ggdata <- cbind(labelprotein, ggdata)\n",
    "### paper numbering\n",
    "library(grid)\n",
    "set.seed(23)\n",
    "b <- ggplot(ggdata, aes(x = PC1, y = PC2, colour = Cluster)) +\n",
    "  geom_point(aes(fill=factor(Cluster)), size=5, shape=20, pch = 21, alpha = 0.8) +\n",
    "  ggtitle(\" \") +\n",
    "  stat_ellipse(aes(fill=factor(Cluster)), colour = \"black\", \n",
    "               geom=\"polygon\", level=0.95, alpha=0.2) +\n",
    "  guides(color=guide_legend(\"Cluster\"),fill=guide_legend(\"Cluster\")) +\n",
    "  #geom_text(aes(label=labelprotein), size=7, hjust=0.5, vjust= 1.5, alpha=0.45) +\n",
    "  theme(\n",
    "    legend.position=(\"none\"),\n",
    "    #plot.title = element_text(size=20, face=\"bold\", colour=\"black\", vjust = 2, hjust=-0.07),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1),\n",
    "    axis.text.y = element_text(size = 15),\n",
    "    axis.ticks.length = unit(0.3, \"cm\"),\n",
    "    axis.text.x = element_text(size = 15),\n",
    "    legend.title=element_blank(),\n",
    "    axis.title.x = element_text(color=\"black\", size=20),\n",
    "    axis.title.y = element_text(color=\"black\", size=20)) +\n",
    "  coord_cartesian(ylim = c(-1.5, 1.5), xlim = c(-1.5, 1.5))\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot with interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "data <- read_excel(\"data.xlsx\")\n",
    "### preparing data for 13 models \n",
    "compound <- data[, 5:311]\n",
    "protein <- data[, 313:351]\n",
    "pIC50 <- data[2]\n",
    "set.seed(23)\n",
    "compound = compound[, -nearZeroVar(compound)]\n",
    "compound[compound == \"0\"] <- -1\n",
    "CxP <- getCPI(compound, protein, type = \"tensorprod\")\n",
    "CxP <- as.data.frame(CxP)\n",
    "dfcompound <- names(data.frame(compound[,1:39]))\n",
    "dfprotein <- names(data.frame(protein[,1:39]))\n",
    "compoundNamecross <- rep(dfcompound, each = 39)\n",
    "proteinNamecross <- rep(dfprotein, times = 39)\n",
    "label <- paste(compoundNamecross, proteinNamecross, sep=\"_\")\n",
    "colnames(CxP) <- label\n",
    "PxP <- getCPI(protein, protein, type = \"tensorprod\")\n",
    "proteinName2 <- rep(dfprotein, times = 39)\n",
    "proteinName1 <- rep(dfprotein, each = 39)\n",
    "label_protein <- paste(proteinName1, proteinName2, sep = \"_\")\n",
    "colnames(PxP) <- label_protein\n",
    "index <- seq(1, 1521, by = 40)\n",
    "protein_selfcross <- PxP[, -index]\n",
    "transposedIndexed_protein <- t(protein_selfcross)\n",
    "index1 <- which(duplicated(transposedIndexed_protein))\n",
    "removed_duplicated_protein <- transposedIndexed_protein[-index1, ]\n",
    "PxP <- t(removed_duplicated_protein)\n",
    "\n",
    "CxC <- getCPI(compound, compound, type = \"tensorprod\")\n",
    "compoundName2 <- rep(dfcompound, times = 39)\n",
    "compoundName1 <- rep(dfcompound, each = 39)\n",
    "label <- paste(compoundName1, compoundName2, sep = \"_\")\n",
    "colnames(CxC) <- label\n",
    "index3 <- seq(1, 1369, by = 38)\n",
    "compound_selfcross <- CxC[, -index3]\n",
    "transposedIndexed_compound <- t(compound_selfcross)\n",
    "index4 <- which(duplicated(transposedIndexed_compound))\n",
    "removed_compound <- transposedIndexed_compound[-index4, ]\n",
    "compound_finalcrossterms <- t(removed_compound)\n",
    "CxC <- compound_finalcrossterms\n",
    "\n",
    "\n",
    "C <- compound\n",
    "P <- protein\n",
    "C_P <- cbind(C, P)\n",
    "C_P_CxP_data_block_scale <- cbind(C, P, CxP) * (1/sqrt(length(C)+length(P)+length(CxP)))\n",
    "\n",
    "C_P_CxC_data_block_scale <- cbind(C, P, \n",
    "                                  CxC) * (1/sqrt(length(C)+length(P)+length(CxC)))\n",
    "\n",
    "C_P_PxP_data_block_scale <- cbind(C, P,\n",
    "                                  PxP) * (1/sqrt(length(C)+length(P)+length(PxP)))\n",
    "\n",
    "C_P_CxP_CxC_data_block_scale <- cbind(C, P, CxP,\n",
    "                                      CxC) * (1/sqrt(length(C)+length(P)+length(CxP)+length(CxC)))\n",
    "\n",
    "C_P_CxP_PxP_data_block_scale <- cbind(C, P, CxP,\n",
    "                                      PxP) * (1/sqrt(length(C)+length(P)+length(CxP)+length(PxP)))\n",
    "\n",
    "C_P_CxC_PxP_data_block_scale <- cbind(C, P, CxC,\n",
    "                                      PxP) * (1/sqrt(length(C)+length(P)+length(CxC)+length(PxP)))\n",
    "\n",
    "C_P_CxP_CxC_PxP_data_block_scale <- cbind(C, P, CxP, CxC, PxP) * (1/sqrt(length(C)+length(P)+\n",
    "                                                                           length(CxC)+length(PxP)))\n",
    "\n",
    "\n",
    "cor_removed <- function(x) {\n",
    "  x <- x[, !apply(x, 2, function(x) length(unique(x)) ==1)]\n",
    "  raw <- cor(x)\n",
    "  raw_2 <- raw[1: ncol(raw), 1: ncol(raw)]\n",
    "  high <- findCorrelation(raw_2, cutoff = 0.7)\n",
    "  remove <- x[, -high]\n",
    "  input <- cbind(pIC50, remove)\n",
    "  return(input)\n",
    "}\n",
    "### Applying function to all the 13 descriptors block\n",
    "\n",
    "C <- cor_removed(C)\n",
    "P <- cor_removed(P)\n",
    "CxP <- cor_removed(CxP)\n",
    "CxC <- cor_removed(CxC)\n",
    "PxP <- cor_removed(PxP)\n",
    "C_P <- cor_removed(C_P)\n",
    "C_P_CxP <- cor_removed(C_P_CxP_data_block_scale)\n",
    "C_P_CxC <- cor_removed(C_P_CxC_data_block_scale)\n",
    "C_P_PxP <- cor_removed(C_P_PxP_data_block_scale)\n",
    "C_P_CxP_CxC <- cor_removed(C_P_CxP_CxC_data_block_scale)\n",
    "C_P_CxP_PxP <- cor_removed(C_P_CxP_PxP_data_block_scale)\n",
    "C_P_CxC_PxP <- cor_removed(C_P_CxC_PxP_data_block_scale)\n",
    "C_P_CxP_CxC_PxP <- cor_removed(C_P_CxC_PxP_data_block_scale)\n",
    "\n",
    "file <- function(x) {\n",
    "  sel <- prospectr::naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    Train <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- caret::train(pIC50 ~., data = Train,  method = \"rf\", tuneLength = 10,\n",
    "                trControl = ctrl)\n",
    "    model <- randomForest::randomForest(pIC50~., data = Train,\n",
    "                         mtry = tune$bestTune[[1]])\n",
    "    prediction_Internal <- predict(model, Train)\n",
    "    value <- data.frame(obs = Train$pIC50, pred = prediction_Internal)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    value$Label <- c(\"Internal\")\n",
    "    prediction_External <- predict(model, Test)\n",
    "    value_external <- data.frame(obs = Test$pIC50, pred = prediction_External)\n",
    "    colnames(value_external) <- labeling\n",
    "    value_external$Label <- c(\"External\")\n",
    "    results <- rbind(value, value_external)\n",
    "    return(results)\n",
    "}\n",
    "\n",
    "get_interval <- function(x) {\n",
    "  file <- file(x)\n",
    "  x <- file[,1]\n",
    "  y <- file[,2]\n",
    "  label <- file[3]\n",
    "  fit <- lm(y~x)\n",
    "  pred.int <- predict(fit, interval = \"prediction\")\n",
    "  pred.lower = pred.int[,2]\n",
    "  pred.upper = pred.int[,3]\n",
    "  df <- cbind(x, y, label, pred.lower, pred.upper)\n",
    "  return(df)\n",
    "}\n",
    "\n",
    "plot_graph <- function(x) {\n",
    "  ok <- get_interval(x)\n",
    "  good <- ggplot(ok, aes(x = x)) +\n",
    "    geom_point(size = 7, colour = \"black\", pch = 21, alpha = 0.8,\n",
    "               aes(y= y, fill = factor(Label))) +\n",
    "    geom_abline(intercept = ok$pred.lower[[1]], linetype = 2, size = 1.5, colour = \"grey\") +\n",
    "  geom_abline(intercept = ok$pred.upper[[1]], linetype = 2, size = 1.5, colour = \"grey\") +\n",
    "    theme(\n",
    "      panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill =NA, size = 1),\n",
    "      axis.text.y = element_text(colour = \"black\", face = \"bold\"),\n",
    "      axis.text.x = element_text(colour = \"black\", face = \"bold\"),\n",
    "      axis.title= element_text(colour = \"black\", face = \"bold\"),\n",
    "     # axis.title.x = element_text(colour = \"black\", face = \"bold\"),\n",
    "      legend.position = (\"none\")) +\n",
    "    labs(list(y = expression(Predicted~pIC[50])), face = \"bold\") +\n",
    "    labs(list(x = expression(Experimental~pIC[50])), face = \"bold\") +\n",
    "    coord_cartesian(ylim = c(-2, 4), xlim = c(-2, 4))\n",
    "  return(good)\n",
    "}\n",
    "\n",
    "a <- plot_graph(C)\n",
    "print(a)\n",
    "b <- plot_graph(P)\n",
    "print(b)\n",
    "c <- plot_graph(CxP)\n",
    "print(c)\n",
    "d <- plot_graph(CxC)\n",
    "print(d)\n",
    "e <- plot_graph(PxP)\n",
    "print(e)\n",
    "f <- plot_graph(C_P)\n",
    "print(f)\n",
    "g <- plot_graph(C_P_CxP)\n",
    "print(g)\n",
    "h <- plot_graph(C_P_CxC)\n",
    "print(h)\n",
    "i <- plot_graph(C_P_PxP)\n",
    "print(i)\n",
    "j <- plot_graph(C_P_CxP_CxC)\n",
    "print(j)\n",
    "k <- plot_graph(C_P_CxP_PxP)\n",
    "print(k)\n",
    "l <- plot_graph(C_P_CxC_PxP)\n",
    "print(l)\n",
    "m <- plot_graph(C_P_CxP_CxC_PxP)\n",
    "print(m)\n",
    "\n",
    "\n",
    "#plot_grid(a, b, c, d, e, f, g, h, i, j, k, l, m,\n",
    "#          labels = c(\"1\", \"2\", \"3\", \"4\", \"5\", \"6\", \"7\", \"8\", \"9\", \n",
    "#                     \"10\", \"11\", \"12\", \"13\"),\n",
    "#          ncol = 4, align = \"h\", label_size = 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feature Importance Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "randomForest_feature_importance <- function(x) {\n",
    "  library(doSNOW)\n",
    "  library(foreach)\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(8)\n",
    "  registerDoSNOW(cl)\n",
    "  \n",
    "  results <- list(100)\n",
    "  results <- foreach (i = 1:100) %dopar% {\n",
    "    sel <- prospectr::naes(x, k = 90, pc = 5, iter.max = 100)\n",
    "    myData <- x[sel$model, ]\n",
    "    Test <- x[sel$test, ]\n",
    "    ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    tune <- caret::train(pIC50 ~., data = myData,  method = \"rf\", tuneLength = 10,\n",
    "                         trControl = ctrl)\n",
    "    model <- randomForest::randomForest(pIC50~., data = myData, importance = TRUE,\n",
    "                                        mtry = tune$bestTune[[1]])\n",
    "    importance <- model$importance\n",
    "    results[[i]] <- importance\n",
    "  }\n",
    "  return(results)\n",
    "  stopCluster(cl)\n",
    "}\n",
    "\n",
    "results_feature_importance_RF <- randomForest_feature_importance(C_P_CxP_CxC_PxP)\n",
    "\n",
    "df <- as.data.frame(results_feature_importance_RF)\n",
    "df <- t(apply(df, 1, round, digits = 4))\n",
    "index <- seq(2, 200, by = 2)\n",
    "df <- df[, index]\n",
    "geni <- mean_and_sd(df)\n",
    "geni <- as.data.frame(geni)\n",
    "mean <- data.frame(apply(df, 1, mean))\n",
    "sd <- data.frame(apply(df, 1, sd))\n",
    "geni_index <- cbind(mean, sd)\n",
    "colnames(geni_index) <- c(\"mean\", \"sd\")\n",
    "geni_index <- geni_index[order(mean, decreasing = TRUE),]\n",
    "geni_index <- head(geni_index, 20)\n",
    "set.seed(200)\n",
    "a <- cbind(descriptors = rownames(geni_index), geni_index)\n",
    "a$descriptors <- factor(a$descriptors, levels = a[order(a$mean), \"descriptors\"])\n",
    "### plot\n",
    "z <- ggplot(a, aes(x = mean, y = descriptors)) +\n",
    "  geom_point(size = 4, colour = \"black\", fill = \"red\", pch = 21) +\n",
    "  geom_errorbarh(aes(xmin = mean-sd, xmax = mean+sd), colour = \"black\") +\n",
    "  ggtitle(\"\") + xlab(\"Gini index\") + ylab(\"\") +\n",
    "  theme(\n",
    "    plot.margin = unit(c(1, 1, 1, 1), \"cm\"),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1),\n",
    "    axis.title.x = element_text(size = 20, face = \"bold\", colour = \"black\"))\n",
    "z + background_grid(major = \"xy\", minor = \"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scrambling Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "scrambling_R2 <- function(x) {\n",
    "  library(doSNOW)\n",
    "  library(foreach)\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(8)\n",
    "  registerDoSNOW(cl)\n",
    "  \n",
    "  results <- list(50)\n",
    "  results <- foreach (i = 1:50) %dopar% {\n",
    "    pIC50 <- gtools::permute(x$pIC50)\n",
    "    pIC50 <- data.frame(pIC50)\n",
    "    fake_data <- cbind(pIC50, x[, 2:ncol(x)])\n",
    "    ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    tune <- caret::train(pIC50~., data = fake_data, method = \"rf\",\n",
    "                  trControl  = ctrl, tuneLength = 10)\n",
    "    \n",
    "    fit <- randomForest::randomForest(pIC50~., data = fake_data, mtry = tune$bestTune[[1]])\n",
    "    prediction <- predict(fit, x)\n",
    "    value <- data.frame(obs = x$pIC50, pred = prediction)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- caret::defaultSummary(value)\n",
    "   }\n",
    "  R2 <- data.frame(results)\n",
    "  R2 <- t(R2)\n",
    "  R2 <- as.numeric(R2[,2])\n",
    "  R2 <- round(R2, digits = 2)\n",
    "  return(R2)\n",
    "  stopCluster(cl)\n",
    "}\n",
    "real_R2 <- function(x) {\n",
    "  ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50~., data = x, method = \"rf\",\n",
    "                  trControl  = ctrl, tuneLength = 10)\n",
    "  fit <- randomForest(pIC50~., data = x, mtry = tune$bestTune[[1]])\n",
    "  prediction <- predict(fit, x)\n",
    "  value <- data.frame(obs = x$pIC50, pred = prediction)\n",
    "  labeling <- c(\"obs\", \"pred\")\n",
    "  colnames(value) <- labeling\n",
    "  result <- defaultSummary(value)\n",
    "  R2 <- as.data.frame(result)\n",
    "  R2 <- R2[2, ]\n",
    "  R2 <- round(R2, digits = 2)\n",
    "  return(R2)\n",
    "}\n",
    "\n",
    "scrambling_Q2 <- function(x) {\n",
    "  library(doSNOW)\n",
    "  library(foreach)\n",
    "  library(parallel)\n",
    "  cl <- makeCluster(8)\n",
    "  registerDoSNOW(cl)\n",
    "  \n",
    "  results <- list(50)\n",
    "  results <- foreach (i = 1:50) %dopar% {\n",
    "    pIC50 <- gtools::permute(x$pIC50)\n",
    "    pIC50 <- data.frame(pIC50)\n",
    "    myData <- cbind(pIC50, x[, 2:ncol(x)])\n",
    "  k = 10\n",
    "  index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "  folds <- 1:k\n",
    "  myRes = data.frame()\n",
    "  for (j in 1:k)\n",
    "    training <- subset(myData, index %in% folds[-j])\n",
    "  testing <- subset(myData, index %in% c(j))\n",
    "  #pIC50 <- gtools::permute(training$pIC50)\n",
    "  #pIC50 <- data.frame(pIC50)\n",
    "  #fake_data <- cbind(pIC50, training[2:ncol(training)])\n",
    "  ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "    tune <- caret::train(pIC50~., data = training, method = \"rf\",\n",
    "                  trControl  = ctrl, tuneLength = 10)\n",
    "    \n",
    "    fit <- randomForest::randomForest(pIC50~., data = training, mtry = tune$bestTune[[1]])\n",
    "    prediction <- predict(fit, testing)\n",
    "    value <- data.frame( obs = testing$pIC50, pred = prediction)\n",
    "    value <- rbind(myRes, value)\n",
    "    labeling <- c(\"obs\", \"pred\")\n",
    "    colnames(value) <- labeling\n",
    "    results[[i]] <- caret::defaultSummary(value)\n",
    "  }\n",
    "  Q2 <- data.frame(results)\n",
    "  Q2 <- t(Q2)\n",
    "  Q2 <- as.numeric(Q2[,2])\n",
    "  Q2 <- round(Q2, digits = 2)\n",
    "  return(Q2)\n",
    "  stopCluster(cl)\n",
    "}\n",
    "\n",
    "real_Q2 <- function(x) {\n",
    "  myData <- x\n",
    "  k = 10\n",
    "  index <- sample(1:k, nrow(myData), replace = TRUE)\n",
    "  folds <- 1:k\n",
    "  myRes <- data.frame()\n",
    "  for (j in 1:k)\n",
    "    training <- subset(myData, index %in% folds[-j])\n",
    "  testing <- subset(myData, index %in% c(j))\n",
    "  ctrl <- caret::trainControl(method = \"repeatedcv\", number = 10, repeats = 1)\n",
    "  tune <- train(pIC50~., data = training, method = \"rf\",\n",
    "                  trControl  = ctrl, tuneLength = 10)\n",
    "  fit <- randomForest(pIC50~., data = training, mtry = tune$bestTune[[1]])\n",
    "  prediction <- predict(fit, testing)\n",
    "  value <- data.frame(obs = testing$pIC50, pred = prediction)\n",
    "  value <- rbind(myRes, value)\n",
    "  labeling <- c(\"obs\", \"pred\")\n",
    "  colnames(value) <- labeling\n",
    "  result <- defaultSummary(value)\n",
    "  Q2 <- as.data.frame(result)\n",
    "  Q2 <- Q2[2, ]\n",
    "  Q2 <- round(Q2, digits = 2)\n",
    "  return(Q2)\n",
    "}\n",
    "\n",
    "data_pre <- function(x) {\n",
    "  fake_R2 <- scrambling_R2(x)\n",
    "  real_R2 <- real_R2(x)\n",
    "  fake_Q2 <- scrambling_Q2(x)\n",
    "  real_Q2 <- real_Q2(x)\n",
    "  fake_R2 <- as.data.frame(fake_R2)\n",
    "  fake_R2$Label <- \"Fake\"\n",
    "  real_R2 <- as.data.frame(real_R2)\n",
    "  real_R2$Label <- \"Real\"\n",
    "  fake_Q2 <- as.data.frame(fake_Q2)\n",
    "  fake_Q2$Label <- \"Fake\"\n",
    "  real_Q2 <- as.data.frame(real_Q2)\n",
    "  real_Q2$Label <- \"Real\"\n",
    "  fake <- cbind(fake_R2, fake_Q2)\n",
    "  real <- cbind(real_R2, real_Q2)\n",
    "  combine <- data.frame(fake, real)\n",
    "  return(combine)\n",
    "}\n",
    "\n",
    "plot_scrambling <- function(x){\n",
    "  ok <- data_pre(x)\n",
    "colnames(ok) <- c(\"R2\", \"Label\", \"Q2\", \"Label\", \"R2\", \"Label\", \"Q2\", \"Label\")\n",
    "R2 <- ok[c(1, 5)]\n",
    "R2 <- reshape2::melt(R2)\n",
    "R2 <- R2$value\n",
    "R2 <- data.frame(R2)\n",
    "Q2 <- ok[c(3, 7)]\n",
    "Q2 <- reshape2::melt(Q2)\n",
    "Q2 <- Q2$value\n",
    "Q2 <- data.frame(Q2)\n",
    "Label <- c(\"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\", \"Fake\",\n",
    "           \"Fake\", \"Fake\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\", \"Real\",\n",
    "           \"Real\", \"Real\", \"Real\", \"Real\")\n",
    "data <- cbind(R2, Q2, Label)\n",
    "plot <- ggplot(data, aes(x = R2, y = Q2, colour = Label)) +\n",
    "  geom_point(size = 7, colour = \"black\", aes(fill = factor(Label)), pch = 21, alpha = 0.8) +\n",
    "  theme(\n",
    "    legend.position = (\"none\"),\n",
    "    panel.border = element_rect(linetype = \"solid\", colour = \"black\", fill = NA, size = 1)) +\n",
    "  labs(y = expression(paste(italic(Q^2)))) +\n",
    "  labs(x = expression(paste(italic(R^2)))) +\n",
    "  coord_cartesian(ylim = c(-0.2, 1.2), xlim = c(-0.2, 1.2))\n",
    "return(plot)\n",
    "}\n",
    "\n",
    "a <- plot_scrambling(C)\n",
    "print(a)\n",
    "b <- plot_scrambling(P)\n",
    "print(b)\n",
    "c <- plot_scrambling(CxP)\n",
    "print(c)\n",
    "d <- plot_scrambling(CxC)\n",
    "print(d)\n",
    "e <- plot_scrambling(PxP)\n",
    "print(e)\n",
    "f <- plot_scrambling(C_P)\n",
    "print(f)\n",
    "g <- plot_scrambling(C_P_CxP)\n",
    "print(g)\n",
    "h <- plot_scrambling(C_P_CxC)\n",
    "print(h)\n",
    "i <- plot_scrambling(C_P_PxP)\n",
    "print(i)\n",
    "j <- plot_scrambling(C_P_CxP_CxC)\n",
    "print(j)\n",
    "k <- plot_scrambling(C_P_CxP_PxP)\n",
    "print(k)\n",
    "l <- plot_scrambling(C_P_CxC_PxP)\n",
    "print(l)\n",
    "m <- plot_scrambling(C_P_CxP_CxC_PxP)\n",
    "print(m)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 0
}