【原】RNAseq｜Mime代碼版-終極101 種機器學習算法組合構建最優(yōu)預后模型

生信補給站 2024-07-18 發(fā)布于北京

展開全文

Mime1為構建基于機器學習的集成模型提供了一個用戶友好的解決方案，利用復雜的數(shù)據(jù)集來識別與預后相關的關鍵基因。

前面單獨介紹了Lasso ，randomForestSRC，Enet（Elastic Net），CoxBoost 和 SuperPC 構建生存模型的方法和參數(shù)，本文介紹如何使用Mime1包一體式完成文獻中的101種機器學習組合的分析，并輸出文獻級別的圖。

除此以外額外介紹一下（1）替換自己數(shù)據(jù)時注意的點 （2）如何提取指定模型下的riskscore結果進行后續(xù)分析和（3）如何對目標癌種進行模型比較（膠質瘤可以使用數(shù)據(jù)內(nèi)置的）（個人認為更重要）

一載入R包，數(shù)據(jù)

該包集合了10種機器學習的包，所以安裝上會稍微繁瑣一下，給點耐心缺什么下載什么。

# options("repos"= c(CRAN="https://mirrors.tuna./CRAN/"))# options(BioC_mirror="http://mirrors.tuna./bioconductor/")if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
depens<-c('GSEABase', 'GSVA', 'cancerclass', 'mixOmics', 'sparrow', 'sva' , 'ComplexHeatmap' )for(i in 1:length(depens)){  depen<-depens[i]  if (!requireNamespace(depen, quietly = TRUE))  BiocManager::install(depen,update = FALSE)}
if (!requireNamespace("CoxBoost", quietly = TRUE))  devtools::install_github("binderh/CoxBoost")
if (!requireNamespace("fastAdaboost", quietly = TRUE))  devtools::install_github("souravc83/fastAdaboost")
if (!requireNamespace("Mime", quietly = TRUE))  devtools::install_github("l-magnificence/Mime")
library(Mime1)

首先查看示例數(shù)據(jù)

1. 基因集數(shù)據(jù)

可以是差異基因，熱點通路（MsigDB），WGCNA或者PPI找到的hub gene ，或者某種單細胞亞型的marker gene 等等任何方式得到的目標基因集

load("./External data/genelist.Rdata")#> [1] "MYC"    "CTNNB1" "JAG2"   "NOTCH1" "DLL1"   "AXIN2"  "PSEN2"  "FZD1"   "NOTCH4" "LEF1"   "AXIN1"  "NKD1"   "WNT5B" #>[14] "CUL1"   "JAG1"   "MAML1"  "KAT2A"  "GNAI1"  "WNT6"   "PTCH1"  "NCOR2"  "DKK4"   "HDAC2"  "DKK1"   "TCF7"   "WNT1"  #>[27] "NUMB"   "ADAM17" "DVL2"   "PPARD"  "NCSTN"  "HDAC5"  "CCND2"  "FRAT1"  "CSNK1E" "RBPJ"   "FZD8"   "TP53"   "SKP2"  #>[40] "HEY2"   "HEY1"   "HDAC11"

2. 生存數(shù)據(jù)與基因表達信息

load("./External data/Example.cohort.Rdata") # 生存數(shù)據(jù)與基因表達信息list_train_vali_Data[["Dataset1"]][1:5,1:5]#                 ID    OS.time OS   MT-CO1   MT-CO3#60  TCGA.DH.A66B.01 1281.65322  0 13.77340 13.67931#234 TCGA.HT.7607.01   96.19915  1 14.96535 14.31857#42  TCGA.DB.A64Q.01  182.37755  0 13.90659 13.65321#126 TCGA.DU.8167.01  471.97707  0 14.90695 14.59776#237 TCGA.HT.7610.01 1709.53901  0 15.22784 14.62756

其中l(wèi)ist_train_vali_Data是含有2個數(shù)據(jù)集的列表，每個數(shù)據(jù)集的第一列為ID ，2-3列為生存信息（OS.time ，OS），后面為基因表達量。

二構建預后模型

1. 構建101機器學習模型組合

該包大大降低了學習成本，可以通過ML.Dev.Prog.Sig函數(shù)直接構建

res <- ML.Dev.Prog.Sig(train_data = list_train_vali_Data$Dataset1,                     list_train_vali_Data = list_train_vali_Data,                     unicox.filter.for.candi = T,                     unicox_p_cutoff = 0.05,                     candidate_genes = genelist,                     mode = 'all',nodesize =5,seed = 5201314 )

ML.Dev.Prog.Sig() 可選 all, single 和 double三種模式. all 為所有10種算法以及組合 . single 為用10種算法中的一種. double 為兩種算法的組合，一般情況下使用 all 模式.
默認情況下 unicox.filter.for.candi 為 T , 會先對訓練集進行單因素cox分析，unicox_p_cutoff 顯著的基因會用于構建預后模型.

如果使用自己數(shù)據(jù)的時候，需要注意：

（1）替換自己數(shù)據(jù)注意前三列的要求，且將多個數(shù)據(jù)集以列表形式存儲。

（2）分析之前最好先確認 所有數(shù)據(jù)集中是否有基因集列表中的所有基因，減少報錯。

（3）種子數(shù)確定好，會有一些小的影響。

好奇查看后發(fā)現(xiàn)示例數(shù)據(jù)Dataset2中缺少基因集中的幾個基因，但是為什么沒有報錯呢？

data2 <- data2 %>%    dplyr::select(ID , OS.time , OS, genelist )#Error in `dplyr::select()`:#! Can't select columns that don't exist.#? Columns `JAG1`, `DKK4`, and `WNT1` don't exist.#Run `rlang::last_trace()` to see where the error occurred.

通過View(ML.Dev.Prog.Sig) 檢查函數(shù)，設置unicox.filter.for.candi = T 后會先做單因素cox分析，單因素顯著的基因才會作為機器學習101模型組合的候選基因進行后續(xù)分析，而下圖紅框中單基因顯著的基因，恰好沒有dataset2中缺少的`JAG1`, `DKK4`, and WNT1 基因，因此沒有報錯。

但是在不確定哪些基因單因素預后顯著的前提下，分析自己數(shù)據(jù)時候還是先確保訓練集和驗證集均有基因列表的所有基因。

2. C-index 展示

示例數(shù)據(jù)list_train_vali_Data 為2個數(shù)據(jù)集的list，結果圖中隊列為2個，最后兩列為Cindex的均值，這也就是機器學習模型組合文獻中的主圖。

cindex_dis_all(res,               validate_set = names(list_train_vali_Data)[-1],               order = names(list_train_vali_Data),               width = 0.35)

3. 查看指定模型的結果

假設我們選擇第一個模型（StepCox[forward] + plsRcox），可以單獨查看該模型下各個數(shù)據(jù)集的cindex表現(xiàn)

cindex_dis_select(res,                  model="StepCox[forward] + plsRcox",                  order= names(list_train_vali_Data))

也可以查看該模型下各個數(shù)據(jù)集的KM曲線情況

survplot <- vector("list",2) for (i in c(1:2)) {  print(survplot[[i]]<-rs_sur(res, model_name = "StepCox[forward] + plsRcox",                              dataset = names(list_train_vali_Data)[i],                              #color=c("blue","green"),                              median.line = "hv",                              cutoff = 0.5,                              conf.int = T,                              xlab="Day",pval.coord=c(1000,0.9))  )}aplot::plot_list(gglist=survplot,ncol=2)

提取模型RS結果

這里有個很重要的點是要提取指定模型下的RS結果，然后就可以根據(jù)自己的需求重新繪制KM 以及獨立預后分析，森林圖，列線圖等其他分析了。

結果都在res中，根據(jù)str(res)知道對應的信息，提取即可

head(res$riskscore$`StepCox[forward] + plsRcox`[[1]])head(res$riskscore$`StepCox[forward] + plsRcox`[[2]])

4. AUC結果

計算每個模型的1年，3年，5年的 auc值，并可視化所有模型的1年auc結果

all.auc.1y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,train_data = list_train_vali_Data[["Dataset1"]],                             inputmatrix.list = list_train_vali_Data,mode = 'all',AUC_time = 1,                             auc_cal_method="KM")all.auc.3y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,train_data = list_train_vali_Data[["Dataset1"]],                             inputmatrix.list = list_train_vali_Data,mode = 'all',AUC_time = 3,                             auc_cal_method="KM")all.auc.5y <- cal_AUC_ml_res(res.by.ML.Dev.Prog.Sig = res,train_data = list_train_vali_Data[["Dataset1"]],                             inputmatrix.list = list_train_vali_Data,mode = 'all',AUC_time = 5,                             auc_cal_method="KM")
auc_dis_all(all.auc.1y,            dataset = names(list_train_vali_Data),            validate_set=names(list_train_vali_Data)[-1],            order= names(list_train_vali_Data),            width = 0.35,            year=1)

同樣可以繪制選定模型下的auc曲線

roc_vis(all.auc.1y,        model_name = "StepCox[forward] + plsRcox",        dataset = names(list_train_vali_Data),        order= names(list_train_vali_Data),        anno_position=c(0.65,0.55),        year=1)

auc_dis_select(list(all.auc.1y,all.auc.3y,all.auc.5y),               model_name="StepCox[forward] + plsRcox",               dataset = names(list_train_vali_Data),               order= names(list_train_vali_Data),               year=c(1,3,5))

5. 模型比較

該包還提供了和之前文獻報道的預后模型比較的函數(shù)，當然只提供了膠質瘤的。

那如果你做的是其他癌種呢？可以通過查看函數(shù)了解是怎樣的輸入形式，然后就做對應的替換后就可以分析了。（很重要）

cc.glioma.lgg.gbm <- cal_cindex_pre.prog.sig(use_your_own_collected_sig = F,                                             type.sig = c('Glioma','LGG','GBM'),                                             list_input_data = list_train_vali_Data)cindex_comp(cc.glioma.lgg.gbm,            res,            model_name="StepCox[forward] + plsRcox",            dataset=names(list_train_vali_Data))

查看函數(shù)，找到內(nèi)置模型的形式

type.sig = c('Glioma','LGG','GBM')pre.prog.sig <- Mime1::pre.prog.sig    if (all(type.sig %in% names(pre.prog.sig))) {      if (length(type.sig) == 1) {        sig.input <- pre.prog.sig[[type.sig[1]]]      }      else {        sig.input <- pre.prog.sig[[type.sig[1]]]        for (i in 2:length(type.sig)) {          sig.input <- rbind(sig.input, pre.prog.sig[[type.sig[i]]])        }      }    }    View(sig.input)