ROC曲线

作者: 萍智医信 | 来源:发表于2021-08-12 00:02 被阅读0次

    ①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()
    

    相关文章

      网友评论

        本文标题:ROC曲线

        本文链接:https://www.haomeiwen.com/subject/wzoxbltx.html