2021.11.12 初版
2021.11.16 更新單因素COX輸出結果，多了半邊括號...
本系列主要用于通過使用訓練集（train）、內部驗證集（test1組）、外部驗證集（test2）聯(lián)合對感興趣的基因進行Cox回歸模型的構建，并用于預測患者的臨床預后情況，通過生存曲線、ROC曲線、生存狀態(tài)圖對模型效能進行可視化
胡亂調侃，與正題無關，趕時間的友友直接刷掉哈~
誤入BioInfor的大黃鴨 --一個喜歡把教程寫著寫著寫成科普的本科臨床醫(yī)學生
大黃鴨自從誤入了生信的大坑之后，也走了很多彎路（畢竟自己也只是個臨床醫(yī)學的苦逼大三生），花了很多錢不說，電腦游戲軟件圖標逐漸被Rstudio等軟件給取代。。。甚至被舍友給逼問為啥不直接報計算機專業(yè)？？？誠然，，自學生信花的時間也太多了吧。。從R說起，卑微的我隨后在淘寶上開了個小店（雜貨鋪？），以別人R語言報錯改代碼為生，就以這點微薄收入度日，主要不是為了掙錢，還是為了認識更多師兄師姐，而且?guī)蛣e人改代碼真的學到了很多很多，R語言水平也提升了很多很多（離譜的是還有人為了8塊錢改個代碼跟我講過價。。。）之后就在想為啥不把這些經驗分享給友友們，好讓友友們學生信少走點彎路（嘻嘻有一說一，其實還是想讓大家關注一下我，以后我修改代碼的話也會實時在里面更新教程，有學員親測看了我改代碼過程一個月之后，他R語言啥報錯再也沒怕過，還說想把我微信給刪了）
（本文第一次試水，多有不足，歡迎大家批評指正）
首先我們先來講包和輸入文件的準備
目錄
- 0.腦子的準備
- 1.包的準備
- 2.輸入文件的準備
- 3.分組
- 4.單因素Cox分析
- 5.LASSO回歸
- 6.多因素Cox回歸（構建Cox模型）
- 7.模型驗證
- 8.輸出結果文件
- 9.可視化
Cox回歸模型介紹：參見百度百科
算了還是簡單講講吧。。又是廢話一通，請各位知道Cox回歸模型是啥東西的大神自行略過。
在我們日常觀看的電視電影里，經常看到此景：“醫(yī)生，我父親還能活幾年?。俊?，“10-15年”
為啥是10-15年呢？是醫(yī)生隨便說的，沒有任何依據(jù)嗎？我亂報個數(shù)，說20年行不行？
不是的。
這時Cox回歸模型就開始刷它的存在感了
Cox模型是由英國統(tǒng)計學家D.R.Cox(1972)年提出的一種半?yún)?shù)回歸模型，它以生存結局和生存時間為因變量，可同時分析眾多因素對生存期的影響。在這里因素可以是臨床因素，比如說腫瘤各種分期，年齡，患者體征，甚至可以是患者的生活習慣等，也可以是各種基因的表達量，可以納入你所感興趣的基因進行模型的構建，它是一個系數(shù)×表達量的公式，通過計算得出的風險值來預測患者的預后情況。本例以基因表達作為例子。
本例基本思路：
機器學習一般可將數(shù)據(jù)分為訓練集（train）、內部驗證集（test1組）、外部驗證集（test2），在本教程中我們在TCGA獲取的數(shù)據(jù)進行隨機分組，隨機分為訓練集和內部驗證集，在ICGC數(shù)據(jù)庫中我們獲取數(shù)據(jù)對模型進行外部驗證。訓練集負責訓練模型，而內部驗證集我們用于對Cox模型的效能進行檢驗，隨機分組則體現(xiàn)出數(shù)據(jù)的隨機性及普遍性，而為了進一步說明模型的可用性，我們使用其他數(shù)據(jù)庫中的數(shù)據(jù)，使用模型檢驗效果，數(shù)據(jù)的來源可以是其他數(shù)據(jù)庫中的基因表達和臨床數(shù)據(jù)，也可以使用我們自己實驗中獲得的患者基因表達和臨床數(shù)據(jù)。
單因素Cox分析主要是分析其因素單獨有沒有對生存有影響，跟Kaplan-Meier生存曲線差不多，算法不一樣
1.包的準備
install.packages("glmnet")
install.packages("caret")
install.packages("survival")
install.packages("survivalROC")
install.packages("survminer")
直接跑代碼就完事了，后續(xù)我更新一個方便批量安裝包的教程哦~（瘋狂暗示-關注一波大黃鴨嘻嘻）
2.輸入文件的準備
我們需要的數(shù)據(jù)文件是TCGA數(shù)據(jù)庫中FPKM的表達數(shù)據(jù)文件，以及臨床數(shù)據(jù)文件（包含患者生存狀態(tài)及生存時間的），將其整理到一個文件中
表格樣式為：
| ID | time | status | Gene1 | Gene2 | ... | GeneN |
|---|---|---|---|---|---|---|
| TCGA樣本名 | 生存時間（天） | 生存狀態(tài)（0/1 活著/死亡） | Gene1exp | Gene2exp | ... | GeneNexp |
樣本名用TCGA-XX-XXXX的格式，可以用Excel的（mid函數(shù)提取前12位）免得到時候再代碼里刪，因為后面的部分為樣本的信息，前面的才是患者信息，如果遇到重復行名（rowname），需要在Excel上尋找兩個一樣的，刪掉其中一個，生存時間一定要用年??！以天為單位的要除以365，這些工作直接用Excel就完事了?？梢杂肰LOOKUP函數(shù)，也可以手動整理，工作量一般般吧，或者等我后續(xù)更新的教程也行，關注大黃鴨我就完事了哈哈哈~
當然，train集和test1集都是在同一個表達矩陣中進行隨機分組的，而test2集需要通過外部數(shù)據(jù)庫（ICGC或GEO等）或者自己實驗的數(shù)據(jù)（由于本例納入計算的數(shù)據(jù)為TCGA和ICGC數(shù)據(jù)庫的Fpkm值，其他的標準可能與模型不適用，建議還是以同一標準處理的RNA-Seq量作為外部驗證組的驗證（比如說用edgeR標準化的TCGA Counts矩陣就不能用Fpkm標準化矩陣進行外部驗證））整理出另外一個矩陣
接下來就是讀取數(shù)據(jù)啦，要記得更改工作路徑哦~（setwd()走一波）
我們這文件保存為exp_time1.txt（train集和test1集），exp_time2.txt（test2集）（timeexp不好看。。兩個e給整重疊了），使用代碼讀?。?/p>
exptime=read.table("expTime1.txt" , header=T , row.names=1 , sep="\t" , check.names=F)
test2=read.table("expTime2.txt" , header=T , row.names=1 , sep="\t" , check.names=F)
header，row.names表示設置表頭和行名，sep的這個\t表示識別空格為分隔符。。畢竟俺們保存為txt的時候，每個詞之間是變成空格的嘛，check.names這個如果不設置成F，到時候列名前面多個“X.”的話別來找我
3.分組
咱前面說了，這個TCGA得隨機分組，也就是分成訓練集和內部驗證集，這一塊有人55分的，也有人64分的，73分的都有，咱這用的是73分，就是訓練集7份，內部驗證集3份，當然怎么分還是得看您的心情哈哈哈
隨機分組，結果不一定好啊，所以這里咱得用循環(huán)，結果不好的話就重新再隨機分組，一直到出現(xiàn)好的結果就輸出。所以，這里得弄個循環(huán)，憑咱電腦配置，咱這設置為最多循環(huán)800次，也夠跑半天了，當然您電腦差點的話，稍微調低一點也行，咱也不敢反對，否則怕您到時候電腦出問題啥的，或者跑了半天也沒出結果，來找我麻煩也不太好。如果您電腦是幾顆銳龍3990x+128G多少頻率的內存，調太高也不好，畢竟分組幾萬次才出個好結果，這也太僥幸了吧，太過違背隨機原則了。
for(i in 1:800){
...
}
構建一個這樣的循環(huán)，本文下面的所有代碼都放到里面去就好了。
咱這用了createDataPartition這個函數(shù)，當然您用隨機數(shù)生成的seed也行，只是這個函數(shù)專門為機器學習分組的，真的很香啊，他可以從我們樣本中隨機抽出train和test，死亡和存活的樣本都有，而且比例不是很過分，而用隨機數(shù)的話，假如倒霉，train組大部分抽到死亡的，test組大部分抽到生存的，就不好了。所以建議你們還是用這個函數(shù)，真的很方便機器學習這塊。。當然隨機數(shù)生成seed也行（畢竟那么倒霉的幾率很?。?，這塊還得友友們親自去百度，后續(xù)如果我覺得價值大可能我也會更新。
grouplist<-createDataPartition(y=exptime[,2] , p=0.7 , list=F)
train<-exptime[grouplist,]
test1<-exptime[-grouplist,]
代碼解析：y=exptime[,2]是為了以第二排作為判斷分類的標準，生存0，死亡1，以免出現(xiàn)train組大部分抽到死亡的，test組大部分抽到生存的
p=0.7當然就是抽70%了
list=F是不輸出為list這個數(shù)據(jù)類型（有點難用），直接輸出為矩陣
下面就是生成train組，test1組矩陣了
本教程就先講到這啦，后續(xù)隨后更新，歡迎大家關注支持～大家關注一下我公眾號：誤入BioInfor的大黃鴨，回復“TCGACox”獲取完整版的代碼
