机器学习筛选基因——随机生存森林

1.前言

1.1使用背景

对于主流的单因素cox筛选预后相关基因,以及目前最常见的最小绝对收缩和选择算法——Lasso(Least absolute shrinkage and selection operator)该方法是一种压缩估计。它通过构造一个惩罚函数得到一个较为精炼的模型,使得它压缩一些系数,同时设定一些系数为零。因此保留了子集收缩的优点,是一种处理具有复共线性数据的有偏估计。

Lasso算法最大的特点在于能识别高维变量中的过拟合(Overfitting)现象,通过一系列参数控制模型的复杂度,从而给出一个看似更好、更简易的模型。


1.2算法原理

今天介绍一下随机生存森林模型,看看根据随机生存森林在面对高维基因变量时,如何进行关键基因的筛选,以及森林生存森林建模

随机生存森林(Random survival forest,RSF)的原理为采用自助法(bootstrap)随机选取N个自助样本,对每个自助样本建立二元递归生存树,且每个自助样本保留37%的原始数据作为袋外数据(Out of bag data,OOB)用于模型的测试。最后对每棵树的预测结果进行综合。

RSF模型的变量筛选通常有两种方法,分别为VIMP法(Variable importance,VIMP)和最小深度法。RSF原理:在纳入一个新变量和不纳入该变量模型错误率的差值。VIMP为正值则说明该变量能增加模型的准确度,如果为负值则说明改变量会降低模型的准确度。最小深度法根据变量是否能区分更多的数据,越靠近根节点的变量越能区分更多变量。因此,最小深度法定义变量重要性根据变量所处节点与根节点之间的距离。换言之,值越小变量越重要。

2.方法

2.1.数据准备

此处的数据为一个转录组的生存数据,其特点是有生存时间,以及生存结局,其次拥有高维变量。

2.2.构建随机生存森林模型RSF

library(randomForestSRC)
rt=read.table(file = "tcga.expTime.txt", header=T, sep="\t", check.names=F)
rownames(rt)=rt[,1]
rt=rt[,-1]

#建立随机生存森林模型
rfsrc_pbcmy <- rfsrc(Surv(futime, fustat) ~ ., 
                     data = rt, 
                     nsplit = 10, 
                     na.action = "na.impute", 
                     tree.err = TRUE,
                     splitrule='logrank',
                     proximity = T,
                     forest = T,
                     #ntime = 60,#时间节点
                     #mtry = 10,#变量数
                     ntree=1000,#树的数量
                     importance = TRUE#变量重要性VIMP
                     #block.size=100,#第几颗树的错误率
)

 2.33.ggRandomForests变量重要性

library(ggRandomForests)
#首先查看VIMP法,排名前15的变量重要性
gg_dta <- gg_vimp(rfsrc_pbcmy, nvar=15)

#最小深度法查看变量重要性gg_dta <- gg_minimal_depth(rfsrc_pbcmy)plot(gg_dta)


#两种方法的结合  VIMP+min_depth
gg_dta <- gg_minimal_vimp(rfsrc_pbcmy)
plot(gg_dta)

 

这里我们要根据depth threshold 这个值,来确定筛选出的变量,当小于阈值时就得到了,在VIMP法和最小深度法结合的情况下筛选出的变量。

2.4.筛选变量,构建随机生存森林

filterExp_time是提取重要基因后的文件

rt=read.table(file = "filterExp_time.txt", header=T, sep="\t", row.names=1,check.names=F)

#模型构建
##构建随机生存森林
rsf_t <- rfsrc(Surv(futime,fustat)~.,data = rt,
               ntree = 1000,nodesize = 15,##该值建议多调整
               splitrule = 'logrank',
               importance = T,
               proximity = T,
               forest = T,
               seed = 123)
rsf_t
#树的棵数与错误率的曲线图
pdf("error.pdf",10,5)
plot(rsf_t)
dev.off()
library(survminer)
#KM+ROC 分析
## 获得每个样本的 riskscore,进一步进行下游分析
score_t <- data.frame(rt[,c(1,2)],Score=rsf_t$predicted)
cut <- surv_cutpoint(score_t,'futime','fustat','Score')
cut

#绘制高低风险组的阈值图
plot(cut)

## 生存分析
cat <- surv_categorize(cut)

fit <- survfit(Surv(futime,fustat)~Score,cat)
#mytheme <- theme_survminer(font.legend = c(14,'plain', 'black'),
#                           font.x = c(14,'plain', 'black'),
#                          font.y = c(14,'plain', 'black')) ## 自定义主题
diff=survdiff(Surv(futime, fustat) ~Score,data = cat)
pValue=1-pchisq(diff$chisq,df=1)
if(pValue<0.001){
        pValue="p<0.001"
}else{
        pValue=paste0("p=",sprintf("%.03f",pValue))
}

#绘制训练集的K-M曲线
pdf("KM_train.pdf",width = 6.5,height =5.5)
ggsurvplot(fit,cat,
           conf.int=T,
           pval=pValue,
           pval.size=6,
           legend.title="Risk",
           legend.labs=c("High risk", "Low risk"),
           xlab="Time(years)",
           break.time.by = 1,
           palette=c("#ED0000FF" ,"#0099B4FF" ),
           risk.table=TRUE,
           risk.table.title="",
           risk.table.height=.25)
dev.off()
## ROC
library(timeROC)
col <- c("#0099B4FF","#925E9FFF","#FDAF91FF") ## 自定义颜色
ROC_rt <- timeROC(score_t$futime,score_t$fustat,score_t$Score,
              cause = 1,weighting = 'marginal',
              times=c(1,3,5), ROC=TRUE)#,iid = T
ROC_rt

#绘制训练集的ROC曲线
pdf("ROC_train.pdf",5,5)
par(cex.axis=1,cex.lab=1,font.axis=2,font.lab=2)
plot(ROC_rt,time=1,col="#00468BFF",title=FALSE,lwd=3)
plot(ROC_rt,time=3,col="#ED0000FF",add=TRUE,title=FALSE,lwd=3)
plot(ROC_rt,time=5,col="#42B540FF",add=TRUE,title=FALSE,lwd=3)
legend('bottomright',
       c(paste0('AUC at 1 years: ',sprintf("%.03f",ROC_rt$AUC[1])),
         paste0('AUC at 3 years: ',sprintf("%.03f",ROC_rt$AUC[2])),
         paste0('AUC at 5 years: ',sprintf("%.03f",ROC_rt$AUC[3]))),
       col=c("#00468BFF", "#ED0000FF", "#42B540FF"),lwd=3,bty = 'n',text.font=2)
dev.off()



## 测试集
testdata=read.table(file = "testdata.txt", header=T, sep="\t", row.names=1,check.names=F)
testdata$futime=testdata$futime/12
rsf_v<- predict(rsf_t,newdata = testdata,proximity=T)
score_v <- data.frame(testdata[,c(1,2)],Score=rsf_v$predicted)
#将测试集根据cutoff值划分为高风险组和低风险组
score_v$Group <- ifelse(score_v$Score>cut$cutpoint$cutpoint,'high','low')
fit <- survfit(Surv(futime,fustat)~Group,score_v)
diff=survdiff(Surv(futime, fustat) ~Group,data = score_v)
pValue=1-pchisq(diff$chisq,df=1)
if(pValue<0.001){
        pValue="p<0.001"
}else{
        pValue=paste0("p=",sprintf("%.03f",pValue))
}
#绘制测试集K-M曲线
pdf("KM_tes.pdf",width = 6.5,height =5.5)
ggsurvplot(fit,score_v,
           conf.int=T,
           pval=pValue,
           pval.size=6,
           legend.title="Risk",
           legend.labs=c("High risk", "Low risk"),
           xlab="Time(years)",
           break.time.by = 1,
           palette=c("#ED0000FF" ,"#0099B4FF" ),
           risk.table=TRUE,
           risk.table.title="",
           risk.table.height=.25)
dev.off()

ROC_rt <- timeROC(score_v$futime,score_v$fustat,score_v$Score,
                  cause = 1,weighting = 'marginal',
                  times=c(1,3,5), ROC=TRUE)#,iid = T
ROC_rt
#绘制测试集的ROC曲线
pdf("ROC_test.pdf",5,5)
par(cex.axis=1,cex.lab=1,font.axis=2,font.lab=2)
plot(ROC_rt,time=1,col="#00468BFF",title=FALSE,lwd=3)
plot(ROC_rt,time=3,col="#ED0000FF",add=TRUE,title=FALSE,lwd=3)
plot(ROC_rt,time=5,col="#42B540FF",add=TRUE,title=FALSE,lwd=3)
legend('bottomright',
       c(paste0('AUC at 1 years: ',sprintf("%.03f",ROC_rt$AUC[1])),
         paste0('AUC at 3 years: ',sprintf("%.03f",ROC_rt$AUC[2])),
         paste0('AUC at 5 years: ',sprintf("%.03f",ROC_rt$AUC[3]))),
       col=c("#00468BFF" ,"#ED0000FF" ,"#42B540FF"),lwd=3,bty = 'n',text.font=2)
dev.off()

根据预测评分划分高低风险组

ROC曲线

K-M曲线

以上就是随机生存森林,筛选基因以及构建生存预测模型的整个过程。

值得注意的是,当我们进行第一次变量筛选的时候,可能很难达到数目很小的基因个数,对于这种情况,我们可以进行二次迭代,就是将得到的基因再次放进去进行筛选,两次迭代后,基因数目大概在5-15的时候,就可以用来建模。

  • 25
    点赞
  • 35
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

小庞统计

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值