Простой способ построения средней / средней кривой ROC в R - PullRequest
Новая
       40

Простой способ построения средней / средней кривой ROC в R

1 голос
/ 01 июля 2019

Я занимаюсь молекулярной стыковкой, проверяя 4 различных метода стыковки.У меня есть 4 таблицы данных: в каждой есть столбцы Имя лиганда, Оценка связывания и Активность соединения (активно = 1, не активно = 0).Я построил кривые ROC для каждого метода на графике, но я не могу понять, как построить для них «среднее» или среднее значение ROC.

Вот изображение моей кривой ROC plot.

Я не знаю, возможно ли это с имеющимися у меня диаграммами данных или мне нужно создать новый файл с ROCданные прогноза.Если да, то как я могу это сделать?

Я новичок в R, поэтому было бы очень полезно получить более простое объяснение! 

#load ROCR
library(ROCR)
library(pROC)

#load ligands and decoys
lig <- unique(read.table("C:/ROC-Curve/Input-Files/actives")[,1]);
dec <- unique(read.table("C:/ROC-Curve/Input-Files/decoys")[,1]);

#load data file from docking
m10n <- read.table("C:/ROC-Curve/Input-Files/med10-nmr",header=T);
m10x <- read.table("C:/ROC-Curve/Input-Files/med10-xray",header=T);
m30n <- read.table("C:/ROC-Curve/Input-Files/med30-nmr-results",header=T);
m30x <- read.table("C:/ROC-Curve/Input-Files/med30-xray-results",header=T);

#change colnames
colnames(m10n)[1]="LigandName";
colnames(m10n)[2]="BindingScore"

colnames(m10x)[1]="LigandName";
colnames(m10x)[2]="BindingScore"

colnames(m30n)[1]="LigandName";
colnames(m30n)[2]="BindingScore"

colnames(m30x)[1]="LigandName";
colnames(m30x)[2]="BindingScore"

#add column with ligand/decoy info
m10n$IsActive <- as.numeric(m10n$LigandName %in% lig)
m10n$IsActive = ifelse(m10n$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), 
m10n$IsActive+1,m10n$IsActive*0) 

m10x$IsActive <- as.numeric(m10x$LigandName %in% lig)
m10x$IsActive = ifelse(m10x$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m10x$IsActive+1,m10n$IsActive*0)

m30n$IsActive <- as.numeric(m30n$LigandName %in% lig)
m30n$IsActive = ifelse(m30n$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m30n$IsActive+1,m10n$IsActive*0)

m30x$IsActive <- as.numeric(m30x$LigandName %in% lig)
m30x$IsActive = ifelse(m30x$LigandName %in% c("active-kc01", 'active-440-11', "active-ak32515", "active-kf7", 'active-kf-11', "active-kh48", "active-kh49a", "active-kh57a", "active-tmi6", "active-zt1826"), m30x$IsActive+1,m10n$IsActive*0)

#define ROC parameters 
#here pred is selected to compare between ligands using AutoDock
PREDm10n <- prediction(m10n$BindingScore*-1, m10n$IsActive {class(x) = "Prediction"; return(x)})
PERFm10n <- performance(PREDm10n, 'tpr','fpr')

PREDm10x <- prediction(m10x$BindingScore*-1, m10x$IsActive)
PERFm10x <- performance(PREDm10x, 'tpr','fpr')

PREDm30n <- prediction(m30n$BindingScore*-1, m30n$IsActive)
PERFm30n <- performance(PREDm30n, 'tpr','fpr')

PREDm30x <- prediction(m30x$BindingScore*-1, m30x$IsActive)
PERFm30x <- performance(PREDm30x, 'tpr','fpr')

jpeg("Medium Runs.jpg")
#plot m10n
m10nsemilog=PERFm10n@x.values[[1]]
m10nsemilog[m10nsemilog < 0.0005]=0.0005
plot(m10nsemilog,PERFm10n@y.values[[1]],type="l",lwd=2,xlab="False Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="blue",main="Medium Runs")
axis(1, c(0,0.001,0.01,0.1,1))
x<-seq(0,1,0.001)
points(x,x,col="gray",type="l")

#plotm10x
m10xsemilog=PERFm10x@x.values[[1]]
m10xsemilog[m10xsemilog < 0.0005]=0.0005
lines(m10xsemilog,PERFm10x@y.values[[1]],lwd=2,type="l",xlab="False 
Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", 
col="red",main="Short Runs")
axis(1, c(0,0.001,0.01,0.1,1))
x<-seq(0,1,0.001)
points(x,x,col="gray",type="l")

#plotm30n
m30nsemilog=PERFm30n@x.values[[1]]
m30nsemilog[m30nsemilog < 0.0005]=0.0005
lines(m30nsemilog,PERFm30n@y.values[[1]],lwd=2,type="l",xlab="False Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="green",main="Short Runs")
axis(1, c(0,0.001,0.01,0.1,1))
x<-seq(0,1,0.001)
points(x,x,col="gray",type="l")

#plotm30x
m30xsemilog=PERFm30x@x.values[[1]]
m30xsemilog[m30xsemilog < 0.0005]=0.0005
lines(m30xsemilog,PERFm30x@y.values[[1]],lwd=2,type="l",xlab="False 
Positive Rate", ylab="True Positive Rate",xaxt="n", log="x", col="purple",main="Short Runs")
axis(1, c(0,0.001,0.01,0.1,1))
x<-seq(0,1,0.001)
points(x,x,col="gray",type="l")

#make legend
legend("topleft", c("medium10-nmr (AUC = 0.586)", "medium10-xray (AUC = 0.564)", "medium30-nmr (AUC = 0.577)", "medium30-xray (AUC= 0.543)", "random"), lty=1, 
   col = c("red", "blue", "green", "purple", "gray"), bty="n")

dev.off()
...