昨天我的COX分析運(yùn)行了接近20個(gè)小時(shí)后，出了結(jié)果，AUC可以達(dá)到0.79，比一開(kāi)始有提高，但是還不夠好。

盡管我還看到一大票0.6的也在發(fā)文章。

比cox分析更快，更好的是用lasso回歸來(lái)做。

我們先來(lái)看看以前的文章是怎么做的，這篇文章去年發(fā)表在Oncotarget上面

第一步，介紹一下TCGA納入人群的基本信息

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

第三步，把納入的標(biāo)本按照正常和癌癥進(jìn)行差異分析

第四步，差異基因進(jìn)行l(wèi)asso回歸得到幾個(gè)關(guān)鍵基因

第五步，按照構(gòu)建的模型，把患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組

第六步，使用構(gòu)建的模型，分別在traning組和testing組，以及總組測(cè)試，這叫做內(nèi)部驗(yàn)證

因?yàn)闆](méi)有認(rèn)真看，很有可能跟給出的圖有出入反正就是那個(gè)意思 

第七步，告訴人家，這個(gè)預(yù)測(cè)模型可以獨(dú)立于臨床相關(guān)信息，比如淋巴結(jié)，年齡這些，這樣才有意義啊

第八步，如果有機(jī)會(huì)，要拿點(diǎn)別人新的數(shù)據(jù)來(lái)測(cè)試啊，這個(gè)叫做普適性驗(yàn)證。

再一次，這里的圖只是占位置用。 到了這里，基本上一篇文章就結(jié)束了，當(dāng)然如果條件允許，可以把這幾個(gè)分子的表達(dá)在自己的標(biāo)本里面跑一跑 如果再往下走，還有： 

第九步，關(guān)鍵基因的下游研究。

這些基因能預(yù)測(cè)生死，應(yīng)該有厲害的功能才對(duì)啊，這里請(qǐng)參考這個(gè)帖子

課題設(shè)計(jì)：收不完的病人查不完的房，臨床醫(yī)生如何快速地設(shè)計(jì)一個(gè)靠譜的課題？

其中l(wèi)asso回歸這一步，基本上網(wǎng)上也沒(méi)有什么教程，我也測(cè)試了一下，我自己的數(shù)據(jù)，最終發(fā)現(xiàn)他找出20個(gè)基因的模型，預(yù)測(cè)的AUC是0.788，跟我cox出來(lái)的差不多， 但是我的模型只要5個(gè)啊，所以，各有利弊。 

1.極速入門

我不能公開(kāi)我的數(shù)據(jù)，所以就用公共數(shù)據(jù)記錄一下： 首先我們安裝R包， 加載R包

加載測(cè)試數(shù)據(jù)，環(huán)境變量中出現(xiàn)，x和y，他們都是矩陣

data(CoxExample)

下面就開(kāi)始了

這么一搞，圖就出來(lái)了

再搞一搞

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

另外一張圖就出來(lái)了

圖中有兩根線，第一根線比較重要，后面的分析暗自用了第一根線的意義

下面這是第三個(gè)操作,就是找出來(lái)，哪幾個(gè)基因被選中了

這邊就是把這幾個(gè)數(shù)據(jù)調(diào)取出來(lái)，包括名稱，位置，系數(shù)

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

照著運(yùn)行不會(huì)出錯(cuò)的話，會(huì)看到很多數(shù)字 我們看看哪些金榜題名

因?yàn)槭菧y(cè)試數(shù)據(jù)，顯示的是V1，V1，實(shí)際上如果是真實(shí)數(shù)據(jù)，顯示的是基因名稱 基本上模型就做好了，然后用predict就可以算出風(fēng)險(xiǎn)值，往下做就全部出來(lái)了。

2.練手材料

下面的數(shù)據(jù)用來(lái)練手，需要注意的點(diǎn)是兩個(gè)， 

第一，x,y最終都是矩陣，其中包含time和status的y，我用survival包的Surv功能讓他們合在一起 

第二，測(cè)試數(shù)據(jù)的'VignetteExample.rdata'需要以這個(gè)字樣檢索，自行下載，放在同一個(gè)工作目錄才能使用

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，默認(rèn)是10萬(wàn)，運(yùn)行起來(lái)相當(dāng)緩慢，這時(shí)候可以用并行運(yùn)算， 這個(gè)包支持的比較好，這樣做：

其中cores=4 表示的是用4個(gè)核來(lái)算，我測(cè)試了一下，發(fā)現(xiàn)低次數(shù)比如1000次的時(shí)候，并行還要慢一點(diǎn)，沒(méi)有測(cè)試100000次的，應(yīng)該會(huì)快很多。

因?yàn)槊看挝叶几杏X(jué)到帖子很值錢，所以我關(guān)閉了贊賞，努力跟自己的氣息對(duì)應(yīng)起來(lái)。