前面我们已经进行了单因素Cox回归，LASSO回归及多因素Cox回归建模，而我们想对建好的模型进行验证，本例通过内部验证组和外部验证组对Cox模型进行内部验证和外部验证，并对模型进行可视化。
7.内部验证组和外部验证组进行模型验证
#test1组的基因表达矩阵
test1GeneExp=test1[,CoxGenes]
#求test1组的风险值
Test1RiskScore=apply(test1GeneExp,1,riskFun)
#以中位值分为高风险组和低风险组
riskTest1=as.vector(ifelse(Test1RiskScore>TrainRiskMedian,"high","low"))
#构建生存分析模型
diffTest1=survdiff(Surv(time, status) ~riskTest1,data = test1)
#计算p值
Test1PValue=1-pchisq(diffTest1$chisq,df=1)
#做ROC分析
Test1ROC = survivalROC(Stime=test1$time, status=test1$status, marker = Test1RiskScore, predict.time =1, method="KM")
#test2组的基因表达矩阵
test2GeneExp=test2[,CoxGenes]
#求test2组的风险值
Test2RiskScore=apply(test2GeneExp,1,riskFun)
#以中位值分为高风险组和低风险组
riskTest2=as.vector(ifelse(Test2RiskScore>TrainRiskMedian,"high","low"))
#构建生存分析模型
diffTest2=survdiff(Surv(time, status) ~riskTest2,data = test2)
#计算p值
Test2PValue=1-pchisq(diffTest2$chisq,df=1)
#做ROC分析
Test2ROC = survivalROC(Stime=test2$time, status=test1$status, marker = Test2RiskScore, predict.time =1, method="KM")
8.输出结果文件
#判断p值和AUC值是否达到最低标准，一般p小于0.05且AUC大于0.65都是比较好的标准
if((TrainPValue<0.05) & (TrainROC$AUC>0.65) & (Test1PValue<0.05) & (Test1ROC$AUC>0.65) & (Test2PValue<0.05) & (Test2ROC$AUC>0.65)){
#输出分组的基因表达矩阵
trainExp=cbind(ID=row.names(trainGeneExp),trainGeneExp)
test1Exp=cbind(ID=row.names(test1GeneExp),test1GeneExp)
write.table(trainExp1,file="train_exp_time.txt",sep="\t",quote=F,row.names=F)
write.table(testExp1,file="test1_exp_time.txt",sep="\t",quote=F,row.names=F)
#这里讲下输出sep="\t"表示以空格分隔，quote=F表示输出的列名不用加双引号，row.name=F表示不输出行名，这里我们说一下，为什么都有了行名，cbind(ID=row.names(trainGeneExp),trainGeneExp)在这里我们还用cbind个行名到表格里，是因为如果用R给的方法输出行名，你就会发现，在用Excel打开txt文件的时候，列名会对不上列，因此我们自造行名，不用write.table函数自带的输出行名的功能
#输出单因素Cox结果
write.table(SingleCoxOut,file="SingleCoxOut.txt",sep="\t",quote=F,row.names=F)
#输出lasso回归结果
#绘制lambda图
pdf(file="lambda.pdf")
plot(r1, xvar = "lambda", label = TRUE)
dev.off()
#绘制交叉验证图
pdf(file="cvfit.pdf")
plot(cvr2)
abline(v=log(c(cvr2$lambda.min,cvr2$lambda.1se)),lty="dashed") #画出均方误差最小时的lambda值和距离均方误差最小时一个标准误的lambda值（非常拗口）
dev.off()
#输出多因素结果(这里也需要cbind，单独定一列行名，前面有讲过原因)
MultiCoxOut=cbind(ID=rownames(MultiCoxOut),MultiCoxOut)
write.table(MultiCoxOut,file="MultiCox.txt",sep="\t",quote=F,row.names=F)
#输出风险结果，需要用到前面的基因表达矩阵，高低风险分组信息
TrainRiskOut =cbind(trainGeneExp,riskTrain))
TrainRiskOut=cbind(ID=rownames(TrainRiskOut),TrainRiskOut)
write.table(TrainRiskOut,file="TrainRisk.txt",sep="\t",quote=F,row.names=F)
Test1rainRiskOut =cbind(test1GeneExp,riskTest1))
Test1rainRiskOut=cbind(ID=rownames(Test1RiskOut),Test1RiskOut)
write.table(Test1RiskOut,file="Test1Risk.txt",sep="\t",quote=F,row.names=F)
Test2rainRiskOut =cbind(test2GeneExp,riskTest2))
Test2rainRiskOut=cbind(ID=rownames(Test2RiskOut),Test2RiskOut)
write.table(Test2RiskOut,file="Test2Risk.txt",sep="\t",quote=F,row.names=F)
#如果结果已经满足了，就退出循环了
break
}
可视化部分之后更新，本教程到这里就结束啦，欢迎大家关注支持～误入BioInfor的大黄鸭，回复“TCGACox”获取完整版的代码
网友评论