TCGA数据库构建生存预测模型之lasso回归

昨天我的COX分析运行了接近20个小时后，出了结果，AUC可以达到0.79，比一开始有提高，但是还不够好。

尽管我还看到一大票0.6的也在发文章。

比cox分析更快，更好的是用lasso回归来做。

我们先来看看以前的文章是怎么做的，这篇文章去年发表在Oncotarget上面

第一步，介绍一下TCGA纳入人群的基本信息

mark

第二步，把患者分成training组和testing 组，并给出基本信息

mark

第三步，把纳入的标本按照正常和癌症进行差异分析

mark

第四步，差异基因进行lasso回归得到几个关键基因

mark

第五步，按照构建的模型，把患者分为高风险组和低风险组

mark

第六步，使用构建的模型，分别在traning组和testing组，以及总组测试，这叫做内部验证

mark

因为没有认真看，很有可能跟给出的图有出入反正就是那个意思 

第七步，告诉人家，这个预测模型可以独立于临床相关信息，比如淋巴结，年龄这些，这样才有意义啊

mark

第八步，如果有机会，要拿点别人新的数据来测试啊，这个叫做普适性验证。

mark

再一次，这里的图只是占位置用。 到了这里，基本上一篇文章就结束了，当然如果条件允许，可以把这几个分子的表达在自己的标本里面跑一跑 如果再往下走，还有： 

第九步，关键基因的下游研究。

这些基因能预测生死，应该有厉害的功能才对啊，这里请参考这个帖子

课题设计：收不完的病人查不完的房，临床医生如何快速地设计一个靠谱的课题？

---

其中lasso回归这一步，基本上网上也没有什么教程，我也测试了一下，我自己的数据，最终发现他找出20个基因的模型，预测的AUC是0.788，跟我cox出来的差不多， 但是我的模型只要5个啊，所以，各有利弊。 

1.极速入门

我不能公开我的数据，所以就用公共数据记录一下： 首先我们安装R包， 加载R包

1. install.packages('glmnet')

2. library(glmnet)

加载测试数据，环境变量中出现，x和y，他们都是矩阵

data(CoxExample)

下面就开始了

1. fit = glmnet(x, y, family = 'cox')

2. plot(fit)

这么一搞，图就出来了

mark

再搞一搞

cvfit = cv.glmnet(x, y, family = 'cox')
plot(cvfit)

另外一张图就出来了

mark

图中有两根线，第一根线比较重要，后面的分析暗自用了第一根线的意义

下面这是第三个操作,就是找出来，哪几个基因被选中了

1. coef.min = coef(cvfit, s = 'lambda.min')

这边就是把这几个数据调取出来，包括名称，位置，系数

active.min = which(coef.min != 0)
index.min = coef.min[active.min]
index.min
coef.min

照着运行不会出错的话，会看到很多数字 我们看看哪些金榜题名

1. > row.names(coef.min)[active.min]

2.  [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万，运行起来相当缓慢，这时候可以用并行运算， 这个包支持的比较好，这样做：

1. require(doMC)

2. registerDoMC(cores=4)

3. system.time(cv.glmnet(x,y,family = 'cox'))

4. system.time(cv.glmnet(x,y,family = 'cox',parallel=TRUE))

其中cores=4 表示的是用4个核来算，我测试了一下，发现低次数比如1000次的时候，并行还要慢一点，没有测试100000次的，应该会快很多。

因为每次我都感觉到帖子很值钱，所以我关闭了赞赏，努力跟自己的气息对应起来。