#######################Use of cellMiner database#######################
#1. Preparation of Drug Data
##1.1 Read drug related data
library(readxl)
rt1 <- read_excel(path = "DTP_NCI60_ZSCORE.xls", skip = 7)
colnames(rt1) <- rt1[1,]
rt1 <- rt1[-1,-c(67,68)]
##1.2 Screening drug standards, #check drug FDA status
table(rt1$`FDA status`)  
rt1 <- rt1[rt1$`FDA status` %in% c("FDA approved", "Clinical trial"),]
rt1 <- rt1[,-c(1, 3:6)]
write.table(rt1, file = "drug.txt",sep = "\t",row.names = F,quote = F)

#2. Preparation of gene expression data
rt2 <- read_excel(path = "RNA__RNA_seq_composite_expression.xls", skip = 9)
colnames(rt2) <- rt2[1,]
rt2 <- rt2[-1,-c(2:6)]
write.table(rt2, file = "geneExp.txt",sep = "\t",row.names = F,quote = F)

#3. Drug Sensitivity Analysis
rm(list = ls())
##3.1 Reference packages
library(impute)
library(limma)

##3.2 Read drug input file
rt <- read.table("drug.txt",sep="\t",header=T,check.names=F)
rt <- as.matrix(rt)
rownames(rt) <- rt[,1]
drug <- rt[,2:ncol(rt)]
dimnames <- list(rownames(drug),colnames(drug))
data <- matrix(as.numeric(as.matrix(drug)),nrow=nrow(drug),dimnames=dimnames)
### Fill in gaps in drug data
mat <- impute.knn(data)
drug <- mat$data
drug <- avereps(drug)

##3.3 Read expression input file
exp <- read.table("geneExp.txt", sep="\t", header=T, row.names = 1, check.names=F)
dim(exp)
exp[1:4, 1:4]

#3.4 Extracting specific gene expression
gene <- read.table("gene.txt",sep="\t",header=F,check.names=F)
genelist <- as.vector(gene[,1])
genelist

genelist <- gsub(" ","",genelist)
genelist <- intersect(genelist,row.names(exp))
exp <- exp[genelist,]

##3.5 Calculation of drug susceptibility
outTab <- data.frame()
for(Gene in row.names(exp)){
  x <- as.numeric(exp[Gene,])
  for(Drug in row.names(drug)){
    y <- as.numeric(drug[Drug,])
    corT <- cor.test(x,y,method="pearson")
    cor <- corT$estimate
    pvalue <- corT$p.value
    if(pvalue < 0.01){
      outVector <- cbind(Gene,Drug,cor,pvalue)
      outTab <- rbind(outTab,outVector)
    }
  }
}
# output correlation analysis results
outTab <- outTab[order(as.numeric(as.vector(outTab$pvalue))),]
write.table(outTab, file="drugCor.txt", sep="\t", row.names=F, quote=F)

##3.6 Visualization
library(ggplot2)
library(ggpubr)

plotList_1 <- list()
corPlotNum <- 16
if(nrow(outTab)<corPlotNum){
  corPlotNum=nrow(outTab)
}
for(i in 1:corPlotNum){
  Gene <- outTab[i,1]
  Drug <- outTab[i,2]
  x <- as.numeric(exp[Gene,])
  y <- as.numeric(drug[Drug,])
  cor <- sprintf("%.03f",as.numeric(outTab[i,3]))
  pvalue=0
  if(as.numeric(outTab[i,4])<0.001){
    pvalue="p<0.001"
  }else{
    pvalue=paste0("p=",sprintf("%.03f",as.numeric(outTab[i,4])))
  }
  df1 <- as.data.frame(cbind(x,y))
  p1=ggplot(data = df1, aes(x = x, y = y))+
    geom_point(size=1)+
    stat_smooth(method="lm",se=FALSE, formula=y~x)+
    labs(x="Expression",y="IC50",title = paste0(Gene,", ",Drug),subtitle = paste0("Cor=",cor,", ",pvalue))+
    theme(axis.ticks = element_blank(), axis.text.y = element_blank(),axis.text.x = element_blank())+
    theme_bw()
  plotList_1[[i]]=p1
}

nrow <- ceiling(sqrt(corPlotNum))
ncol <- ceiling(corPlotNum/nrow)
ggarrange(plotlist=plotList_1,nrow=nrow,ncol=ncol)