昨天我的COX分析运行了接近20个小时后,出了结果,AUC可以达到0.79,比一开始有提高,但是还不够好。
尽管我还看到一大票0.6的也在发文章。
比cox分析更快,更好的是用lasso回归来做。
我们先来看看以前的文章是怎么做的,这篇文章去年发表在Oncotarget上面
第一步,介绍一下TCGA纳入人群的基本信息
第二步,把患者分成training组和testing 组,并给出基本信息
第三步,把纳入的标本按照正常和癌症进行差异分析
第四步,差异基因进行lasso回归得到几个关键基因
第五步,按照构建的模型,把患者分为高风险组和低风险组
第六步,使用构建的模型,分别在traning组和testing组,以及总组测试,这叫做内部验证
因为没有认真看,很有可能跟给出的图有出入反正就是那个意思
第七步,告诉人家,这个预测模型可以独立于临床相关信息,比如淋巴结,年龄这些,这样才有意义啊
第八步,如果有机会,要拿点别人新的数据来测试啊,这个叫做普适性验证。
再一次,这里的图只是占位置用。 到了这里,基本上一篇文章就结束了,当然如果条件允许,可以把这几个分子的表达在自己的标本里面跑一跑 如果再往下走,还有:
第九步,关键基因的下游研究。
这些基因能预测生死,应该有厉害的功能才对啊,这里请参考这个帖子
课题设计:收不完的病人查不完的房,临床医生如何快速地设计一个靠谱的课题?
其中lasso回归这一步,基本上网上也没有什么教程,我也测试了一下,我自己的数据,最终发现他找出20个基因的模型,预测的AUC是0.788,跟我cox出来的差不多, 但是我的模型只要5个啊,所以,各有利弊。
1.极速入门
我不能公开我的数据,所以就用公共数据记录一下: 首先我们安装R包, 加载R包
install.packages('glmnet')
library(glmnet)
加载测试数据,环境变量中出现,x和y,他们都是矩阵
data(CoxExample)
下面就开始了
fit = glmnet(x, y, family = 'cox')
plot(fit)
这么一搞,图就出来了
再搞一搞
cvfit = cv.glmnet(x, y, family = 'cox')
plot(cvfit)
另外一张图就出来了
图中有两根线,第一根线比较重要,后面的分析暗自用了第一根线的意义
下面这是第三个操作,就是找出来,哪几个基因被选中了
coef.min = coef(cvfit, s = 'lambda.min')
这边就是把这几个数据调取出来,包括名称,位置,系数
active.min = which(coef.min != 0)
index.min = coef.min[active.min]
index.min
coef.min
照着运行不会出错的话,会看到很多数字 我们看看哪些金榜题名
> row.names(coef.min)[active.min]
[1] 'V1' 'V2' 'V3' 'V4' 'V5' 'V6' 'V7' 'V8' 'V9' 'V10' 'V13' 'V17' 'V21' 'V22' 'V25' 'V27' 'V30'
因为是测试数据,显示的是V1,V1,实际上如果是真实数据,显示的是基因名称 基本上模型就做好了,然后用predict就可以算出风险值,往下做就全部出来了。
2.练手材料
下面的数据用来练手,需要注意的点是两个,
第一,x,y最终都是矩阵,其中包含time和status的y,我用survival包的Surv功能让他们合在一起
第二,测试数据的'VignetteExample.rdata'需要以这个字样检索,自行下载,放在同一个工作目录才能使用
library('glmnet')
library('survival')
load('VignetteExample.rdata')
x <- patient.data$x
y <- data.matrix(Surv(patient.data$time,patient.data$status))
cv.fit <- cv.glmnet(x, y, family='cox', maxit = 1000)
plot(cv.fit)
fit <- glmnet(x, y, family = 'cox', maxit = 1000)
plot(fit)
Coefficients <- coef(fit, s = cv.fit$lambda.min)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index
Active.Coefficients
row.names(Coefficients)[Active.Index]
3.并行化处理
上面有一步是写的maxit=1000,默认是10万,运行起来相当缓慢,这时候可以用并行运算, 这个包支持的比较好,这样做:
require(doMC)
registerDoMC(cores=4)
system.time(cv.glmnet(x,y,family = 'cox'))
system.time(cv.glmnet(x,y,family = 'cox',parallel=TRUE))
其中cores=4 表示的是用4个核来算,我测试了一下,发现低次数比如1000次的时候,并行还要慢一点,没有测试100000次的,应该会快很多。
因为每次我都感觉到帖子很值钱,所以我关闭了赞赏,努力跟自己的气息对应起来。
联系客服