rm(list = ls())
for (i in 1:1) {
###############################input data_1
dir_path <- "C:\\Users\\liyix\\OneDrive\\Desktop\\2021\\2021_CYP3A\\1.1_validate_fda\\substrate\\"
dir_path_name <- dir(dir_path,pattern = "*.",full.names = T)
dir_path_name
###############################merge data
data_train <- read.csv(grep("2019-11-15-X3A4.outcome-0.05_subs.csv",dir_path_name,value = T),header = T,stringsAsFactors = F)
data_test <- read.csv(grep("2021-01-07-data_for_model.csv",dir_path_name,value = T),header = T,stringsAsFactors = F)
dim(data_train) #[1] 2902 344
dim(data_test) #87 1026
#View(head(data_test))
unique(data_test$Col2)
unique(data_test$Col1)
data_test <- data_test[data_test$Col1 == "CYP3A Substrates", ]
dim(data_test) #[1] 41 1026
data_test$Col1 <- NULL
data_test <- data_test[,c(colnames(data_train)[-ncol(data_train)],"Col2")]
setdiff(colnames(data_test),colnames(data_train)) #[1] "Mapping.ID"
setdiff(colnames(data_train),colnames(data_test)) #[1] "X2020.02.24.ecfp4.txt"
dim(data_test) #[1] 2379 51
colnames(data_test)
###################################model
library(e1071)
library(ROCR)
library(pROC)
library(ROSE)
library(DMwR)
colnames(data_train)[ncol(data_train)] #"X2020.02.24.ecfp4.txt"
colnames(data_train)[ncol(data_train)] <- "num"
data_train <- data.frame(sapply(data_train, as.numeric))
data_train$num <- as.factor(data_train$num)
str(data_train$num)
############################################################################test_cross_validation
dim(data_train) #[1] 4981 51
#View(head(data_train))
colnames(data_train)
table(data_train$num) #4515 466
table(data_train[ncol(data_train)]) ##4515 466
testdata_ori <- data_train[sample(nrow(data_train),size = floor(nrow(data_train)*0.3),replace = F),]
traindata_ori <- data_train[-match(row.names(testdata_ori),row.names(data_train)),]
intersect(row.names(testdata_ori),row.names(traindata_ori)) #character(0)
table(testdata_ori$num) #1366 128
table(traindata_ori$num) #3149 338
#traindata_ori <- ROSE(num ~ ., data = traindata_ori, seed = 1)$data
dim(traindata_ori)
table(traindata_ori$num) #1795 1692
fit_svm_ori <- svm(num ~ .,data = traindata_ori, scale = F,probability=TRUE)
#### auc_svm
n_col_newdata <- ncol(data_train)
pre_auc_svm_ori <- predict(fit_svm_ori, testdata_ori[,-n_col_newdata], probability = TRUE, decision.values = TRUE)
#View(pre_auc_svm_ori)
pre_auc_svm_d_ori <- attr(pre_auc_svm_ori,"probabilities")
class(pre_auc_svm_d_ori) #[1] "matrix"
head(pre_auc_svm_d_ori)
#View(pre_auc_svm_d_ori)
pred_auc_svm_ori <- prediction(pre_auc_svm_d_ori[,grep("1",colnames(pre_auc_svm_d_ori))],testdata_ori[,n_col_newdata])
head(pre_auc_svm_d_ori)
auc_svm_ori <- performance(pred_auc_svm_ori,"auc")@y.values[[1]]
auc_svm_ori #[1] 0.8321537
pre_probability_ori <- data.frame(pre_auc_svm_ori)
pre_classification_ori <- predict(fit_svm_ori, testdata_ori[,-ncol(data_test)], type="class")
class(pre_classification_ori)
pre_probability_ori$num_1 <- as.character(pre_classification_ori)
pre_probability_ori$num <- testdata_ori$num
#head(pre_probability_ori)
table(pre_probability_ori$num,pre_probability_ori$num_1)
modelroc_svm_ori <- pROC::roc(testdata_ori[,n_col_newdata],as.numeric(pre_auc_svm_d_ori[,grep("1",colnames(pre_auc_svm_d_ori))]))
cutoff_svm_ori <- coords(modelroc_svm_ori, "best", ret = "threshold") #[1] 0.359122
predict.results_svm_ori <- ifelse(pre_auc_svm_d_ori[,grep("1",colnames(pre_auc_svm_d_ori))] > as.numeric(cutoff_svm_ori[1]),"1","0")
freq_default_svm_ori <- table(predict.results_svm_ori,testdata_ori[,n_col_newdata])
Sensitivity_svm_ori <- freq_default_svm_ori[1,1]/sum(freq_default_svm_ori[,1])
Specificity_svm_ori <- freq_default_svm_ori[2,2]/sum(freq_default_svm_ori[,2])
accuracy_svm_ori <- mean(testdata_ori[,n_col_newdata] == predict.results_svm_ori)
Balanced_accuracy_svm_ori <- (Sensitivity_svm_ori+Specificity_svm_ori)/2
Balanced_accuracy_svm_ori #[1] 0.7810957
library(caret)
confusionMatrix(factor(predict.results_svm_ori),factor(testdata_ori[,n_col_newdata]))
confusionMatrix(factor(pre_probability_ori$num_1),factor(testdata_ori[,n_col_newdata]))
#################################################################################################
##svm
#data_train <- ROSE(num ~ ., data = data_train, seed = 1)$data
#################################################probability
fit_probability <- svm(num ~ .,data = data_train, scale = F, probability=TRUE)
#### AUC
names(data_test)[ncol(data_test)]
data_test[,-ncol(data_test)] <- sapply(data_test[,-ncol(data_test)], as.numeric)
ncol(data_test) #[1] 81
pre_auc_probability <- predict(fit_probability, data_test[,-ncol(data_test)], probability=TRUE, decision.values = TRUE)
pre_probability <- attr(pre_auc_probability,"probabilities")
pre_probability <- data.frame(pre_probability)
#class(pre_auc_probability)
#pre_probability$num <- as.character(pre_auc_probability)
head(pre_probability)
colnames(pre_probability) <- c("probability_0","probability_1")
dim(pre_probability) #[1] 2379 3
pre_probability$name <- data_test$Col2
print(head(pre_probability))
class(pre_probability)
print(table(cut(pre_probability$probability_1,breaks = c(-0.1,.5,1))))
write.csv(pre_probability,paste0(dir_path,Sys.Date(),"-substrate.csv"),row.names = F)
table(pre_probability$classification)
}
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
## Loaded ROSE 0.0-3
## Loading required package: lattice
## Loading required package: grid
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Factor w/ 2 levels "0","1": 1 1 1 1 1 2 1 2 1 1 ...
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Warning in coords.roc(modelroc_svm_ori, "best", ret = "threshold"): The
## 'transpose' argument to FALSE by default since pROC 1.16. Set transpose = TRUE
## explicitly to revert to the previous behavior, or transpose = TRUE to silence
## this warning. Type help(coords_transpose) for additional information.
## Loading required package: ggplot2
## probability_0 probability_1 name
## 1 0.07707328 0.9229267 avanafil
## 2 0.14296841 0.8570316 Buspirone
## 3 0.16429482 0.8357052 conivaptan
## 4 0.26352409 0.7364759 darifenacin
## 5 0.49667293 0.5033271 darunavir
## 6 0.15967632 0.8403237 ebastine
##
## (-0.1,0.5] (0.5,1]
## 4 37
#0 1
#2003 376
#head(pre_probability)