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

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

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

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

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

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

第三步，把納入的標本按照正常和癌癥進行差異分析

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

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

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

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

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

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

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

第九步，關(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

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

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

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，默認是10萬，運行起來相當緩慢，這時候可以用并行運算， 這個包支持的比較好，這樣做：

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

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