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

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

比cox分析更快，更好的是用lasso回歸來做,。

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

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

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

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

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

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

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

因為沒有認(rèn)真看,，很有可能跟給出的圖有出入反正就是那個意思 

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

第八步,，如果有機(jī)會,，要拿點別人新的數(shù)據(jù)來測試啊，這個叫做普適性驗證,。

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

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

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

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

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

1.極速入門

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

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

data(CoxExample)

下面就開始了

這么一搞,，圖就出來了

再搞一搞

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

另外一張圖就出來了

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

下面這是第三個操作,就是找出來,，哪幾個基因被選中了

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

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

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

因為是測試數(shù)據(jù),，顯示的是V1，V1,，實際上如果是真實數(shù)據(jù),，顯示的是基因名稱 基本上模型就做好了，然后用predict就可以算出風(fēng)險值,，往下做就全部出來了,。

2.練手材料

下面的數(shù)據(jù)用來練手，需要注意的點是兩個,， 

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

第二,，測試數(shù)據(jù)的'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,，默認(rèn)是10萬，運(yùn)行起來相當(dāng)緩慢,，這時候可以用并行運(yùn)算,， 這個包支持的比較好，這樣做：

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

因為每次我都感覺到帖子很值錢,，所以我關(guān)閉了贊賞，努力跟自己的氣息對應(yīng)起來,。