①ROC曲线及根据cutoff值分组

输入文件为经多因素cox分析后得到的风险值文件

输入文件riskScore.txt.png

library(survivalROC)       #引用包
setwd("E:\\Master research")      #设置工作目录
rt=read.table("riskScore.txt", header=T, sep="\t", check.names=F, row.names=1)    #读取cox回归风险文件

#ROC曲线绘制
predictTime=1       #1年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker=rt$riskScore, predict.time =predictTime, method="KM")
pdf(file="ROC.pdf", width=5.5, height=5.5)
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col="black", 
    xlab="False positive rate", ylab="True positive rate",
    lwd = 2, cex.main=1.2, cex.lab=1.2, cex.axis=1.2, font=1.2)
polygon(x=c(0,roc$FP,1,0),y=c(0,roc$TP,1,0),col="#24B35D",border=NA)
text(0.85, 0.1, paste0("AUC=",sprintf("%.3f",roc$AUC)), cex=1.2)
segments(0,0,1,1,lty=2)
dev.off()

#获取最优cutoff
predictTime=1       #1年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker=rt$riskScore, predict.time =predictTime, method="KM")
sum=roc$TP-roc$FP
cutOp=roc$cut.values[which.max(sum)]
cutTP=roc$TP[which.max(sum)]
cutFP=roc$FP[which.max(sum)]
pdf(file="ROC.cutoff.pdf",width=5.5,height=5.5)
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col="black", 
    xlab="False positive rate", ylab="True positive rate",
    lwd = 2, cex.main=1.2, cex.lab=1.2, cex.axis=1.2, font=1.2)
polygon(x=c(0,roc$FP,1,0),y=c(0,roc$TP,1,0),col="#24B35D",border=NA)
segments(0,0,1,1,lty=2)
points(cutFP,cutTP, pch=20, col="red",cex=1.5)
text(cutFP+0.15,cutTP-0.05,paste0("Cutoff:",sprintf("%0.3f",cutOp)))
text(0.85, 0.1, paste0("AUC=",sprintf("%.3f",roc$AUC)), cex=1.2)
dev.off()

#输出风险文件
risk=as.vector(ifelse(rt$riskScore>cutOp,"high","low"))
outTab=cbind(rt, risk)
write.table(cbind(id=rownames(outTab),outTab),file="risk.txt",sep="\t",quote=F,row.names=F)

######多个时间的ROC曲线######
rocCol=c("red", "green", "blue")
aucText=c()
pdf(file="ROC.multiTime.pdf",width=6,height=6)
#绘制1年的ROC曲线
predictTime=1
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker=rt$riskScore, predict.time=predictTime, method="KM")
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[1], 
  xlab="False positive rate", ylab="True positive rate",
  lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
aucText=c(aucText,paste0("one year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
abline(0,1)
#绘制2年的ROC曲线
predictTime=2
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker=rt$riskScore, predict.time =predictTime, method="KM")
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[2],lwd = 2)
aucText=c(aucText,paste0("two year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
#绘制3年的ROC曲线
predictTime=3
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker=rt$riskScore, predict.time =predictTime, method="KM")
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[3],lwd = 2)
aucText=c(aucText,paste0("three year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
#绘制图例
legend("bottomright", aucText,lwd=2,bty="n",col=rocCol)
dev.off()

根据cutoff值分为高低风险组.png

②多指标ROC曲线

输入文件clinical.png
输入文件risk.png

library(survivalROC)         #引用包
riskFile="risk.txt"          #风险文件
cliFile="clinical.txt"       #临床数据文件
setwd("E:\\Master research\\T细胞和HNSC\\Tcell and HNSC-TCGA\\18.cliROC")     #修改工作目录

#读取风险文件
risk=read.table(riskFile, header=T, sep="\t", check.names=F, row.names=1)
risk=risk[,c("futime","fustat","riskScore")]

#读取临床数据文件
cli=read.table(cliFile,sep="\t",header=T,check.names=F,row.names=1)

#合并数据
samSample=intersect(row.names(risk), row.names(cli))
risk1=risk[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(risk1, cli)

#绘制risk打分的ROC曲线
rocCol=rainbow(ncol(rt)-2)
aucText=c()
pdf(file="cliROC.pdf", width=6, height=6)
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=risk$futime, status=risk$fustat, marker=risk$riskScore, predict.time=3, method="KM") #3年生存期预测
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[1], 
    xlab="False positive rate", ylab="True positive rate",
    lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
aucText=c(aucText,paste0("risk score"," (AUC=",sprintf("%.3f",roc$AUC),")"))
abline(0,1)

#绘制其他临床性状的ROC曲线
j=1
for(i in colnames(rt[,4:ncol(rt)])){
    roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,i], predict.time =1, method="KM")
    j=j+1
    lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[j],lwd = 2)
    aucText=c(aucText,paste0(i," (AUC=",sprintf("%.3f",roc$AUC),")"))
}
legend("bottomright", aucText, lwd=2, bty="n", col=rocCol)
dev.off()