介绍

方法学参考文献

Liu H, Zhang W, Zhang Y, et al. Mime: A flexible machine-learning framework to construct and visualize models for clinical characteristics prediction and feature selection. Comput Struct Biotechnol J. 2024.

原文地址

https://pubmed.ncbi.nlm.nih.gov/39055398/

教程正文

Mime包提供了一个用户友好的解决方案，用于从转录组数据构建基于机器学习的集成模型。

随着高通量测序技术的广泛使用，对生物学和癌症异质性的理解已经发生了革命性的变化。Mime简化了开发高精度预测模型的过程，利用复杂的数据集识别与疾病进展、患者结果和治疗反应相关的关键基因。

它提供四个主要应用：

使用10种机器学习算法建立预测模型。
使用7种机器学习算法构建二元响应模型。
使用8种机器学习方法进行与预后相关的核心特征选择。
可视化每个模型的性能。

小胖体会：批量机器学习的基准测试用这个包实现起来比较方便，其不仅可以针对转录组数据，任何的结构化数据都可以方便的实现。

安装

可以像这样安装Mime的开发版：

# options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

depens'GSEABase', 'GSVA', 'cancerclass', 'mixOmics', 'sparrow', 'sva' , 'ComplexHeatmap' )
for(i in 1:length(depens)){
  depen  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)

Quick Start

这是一个向您展示如何使用Mime的基本示例：

Mime需要包含多个转录测序数据的队列，其中包含生存或治疗临床反应的信息以及基因集作为输入。您可以下载外部数据文件中的示例数据以检查输入的格式。

测试数据下载地址

https://github.com/l-magnificence/Mime/tree/main/data

load("Example.cohort.Rdata")
list_train_vali_Data[["Dataset1"]][1:5,1:5]

load("Example.ici.Rdata")
list_train_vali_Data[["training"]][1:5,1:5]

names(list_train_vali_Data)

df1 "Dataset1"]]
df2 "Dataset2"]]

小胖体会：注意两个数据集的前列的表头是一样的，这也是函数要求的形式。如果是自己的数据，需要整理成一样的形式。

第一列ID是样本ID，第二至第三列OS.time和OS是患者的生存时间和状态，其他列是用log2（x+1）缩放的基因表达水平。Dataset1是训练集，而其他数据集是验证。

load("./Example.ici.Rdata")
list_train_vali_Data[["training"]][1:5,1:5]

第一列ID是样本ID，第二列Var是患者的治疗反应（N：无反应；Y：反应），其他列是用log2（x+1）缩放的基因表达水平。

load("./genelist.Rdata")

该基因集与来自MSigDB的Wnt/β-catenin信号转导相关。
我们建议训练集超过100个样本，基因集超过50个基因。

1.构建预后预测模型

1.1选择最优模型

library(Mime)
load("./Example.cohort.Rdata")
load("./genelist.Rdata")


# 预后模型的基准测试
res_prog   ML.Dev.Prog.Sig(
    train_data = list_train_vali_Data$Dataset1,#训练数据
    list_train_vali_Data = list_train_vali_Data,#整体数据集列表
    unicox.filter.for.candi = T,#利用单因素cox对特征进行过滤
    unicox_p_cutoff = 0.05,# 利用单因素过滤cox的P值截断值
    candidate_genes = genelist,# 目标基因集（应用的特征集合）
    mode = "all",#机器学习模式，目前应用10种机器学习模型及其组合
    nodesize =10,# 随机森林的参数
    seed = 5201314,# 种子数，保证可重复
    cores_for_parallel=5# 并行计算的核心数
    )

AI知识补充

在随机森林中，nodesize（节点大小）具有以下重要意义：

一、对决策树生长的影响

在随机森林由多棵决策树组成，而每棵决策树在生长过程中，nodesize 参数起着关键作用。当决策树进行分支时，如果一个节点中的样本数量小于 nodesize，那么这个节点将不再进行分裂。

例如，假设 nodesize 设置为 50。在构建决策树时，如果某个节点包含的样本数量小于 50，该节点就会成为叶节点，不再进一步分裂。这样可以防止决策树过度生长，避免过拟合问题。

二、对模型性能的影响

准确性

较小的 nodesize 可能会使决策树更加细致地对数据进行划分，从而更好地拟合训练数据。但是，这也可能导致过拟合，使得模型在训练集上表现良好，而在测试集上的性能下降。
较大的 nodesize 可以限制决策树的生长，使模型更加简单和鲁棒。然而，如果 nodesize 过大，可能会导致模型欠拟合，无法充分捕捉数据中的复杂关系。

泛化能力

合适的 nodesize 有助于提高模型的泛化能力。通过限制节点的最小样本数量，可以避免决策树对噪声和异常值过于敏感，从而在新的、未见过的数据上表现更好。

计算效率

较小的 nodesize 通常会导致决策树的深度增加，从而增加模型的训练时间和计算复杂度。
较大的 nodesize 可以减少决策树的深度和复杂度，提高模型的训练速度和预测效率。

三、参数调整策略

在实际应用中，可以通过调整 nodesize 来优化随机森林模型的性能。以下是一些调整策略：

交叉验证：使用交叉验证方法，尝试不同的 nodesize 值，并选择在验证集上表现最佳的参数。可以将数据集划分为多个子集，分别使用不同的 nodesize 进行训练和验证，然后比较它们的性能指标，如准确率、召回率、F1 值等。
网格搜索：结合其他超参数一起进行网格搜索。随机森林中有多个超参数，如树的数量、最大深度、特征子集大小等。可以将 nodesize 与其他超参数组合在一起，进行网格搜索，以找到最佳的参数组合。
可视化分析：通过可视化决策树的结构和性能指标，观察不同 nodesize 值对决策树的影响。可以使用工具如 Graphviz 来可视化决策树，分析节点大小与树的复杂度、准确性之间的关系。

总之，nodesize 是随机森林中的一个重要参数，它对决策树的生长、模型性能和计算效率都有显著影响。在实际应用中，需要根据具体问题和数据集的特点，合理调整 nodesize 值，以获得最佳的模型性能。

class(res_prog)
names(res_prog )

res_prog$Cindex.res %>% head()

res_prog$Sig.genes

res_prog$ml.res

res_prog$ml.res$RSF

res_prog$riskscore




    
res_prog$riskscore$RSF$Dataset1 %>% head()

ML.Dev. Prog.Sig（）提供了三种模式，包括all, single和double。all表示使用所有十个算法和组合。single表示仅使用十个算法中的一个。double表示使用两个算法的组合。在大多数案例中，我们通常会使用all模式来分析数据。
如果您将unicox. filter.for.candi设置为T（默认），Mime将首先在训练集中提供的基因中执行单变量cox回归，以筛选出预测变量，然后用于构建模型。

绘制每个模型的C指数：

cindex_dis_all(object = res_prog,# 基准测试结果
               validate_set = names(list_train_vali_Data)[-1],#测试数据集名称
               order =names(list_train_vali_Data),#显示数据集顺序
               width = 0.35)

绘制不同数据集之间特定模型的C指数：

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

根据特定模型计算的风险评分，在不同数据集中绘制患者的生存曲线：

survplot "list",2) 
for (i in c(1:2)) {
  print(survplot[[i]]                              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_sur函数绘制的生存曲线，以下说明rs_sur函数的具体用法

Mime1::rs_sur(object = res_ML.Dev.Prog.Sig, #基准测试的结果
             model_name = "StepCox[both] + plsRcox",#具体绘制机器算法名称
             dataset = names(list_train_vali_Data)[1],#需要绘制的数据
             #color=c("blue","green"),
             median.line = "hv",



    
             cutoff = 0.5,# 分组截断值
             conf.int = T,#是否绘制可信区间
             xlab="Day",#x坐标
             pval.coord=c(4000,0.9)# P值的位置坐标
             )

1.2 计算每个模型的 AUC 分数


all.auc.1y   cal_AUC_ml_res(
    res.by.ML.Dev.Prog.Sig = res_prog,# 基准测试结果
    train_data = list_train_vali_Data[["Dataset1"]],#训练数据集
    inputmatrix.list = list_train_vali_Data,#基准测试时数据集整体列表
    mode = 'all',#模式，对应基准测试时候的参数
    AUC_time = 1,#指定计算auc的时间
    auc_cal_method="KM")#计算auc的方法


all.auc.3y   cal_AUC_ml_res(
    res.by.ML.Dev.Prog.Sig = res_prog,# 基准测试结果
    train_data = list_train_vali_Data[["Dataset1"]],#训练数据集
    inputmatrix.list = list_train_vali_Data,#基准测试时数据集整体列表
    mode = 'all',#模式，对应基准测试时候的参数
    AUC_time = 3,#指定计算auc的时间
    auc_cal_method="KM")#计算auc的方法

all.auc.5y   cal_AUC_ml_res(
    res.by.ML.Dev.Prog.Sig = res_prog,# 基准测试结果
    train_data = list_train_vali_Data[["Dataset1"]],#训练数据集
    inputmatrix.list = list_train_vali_Data,#基准测试时数据集整体列表
    mode = 'all',#模式，对应基准测试时候的参数
    AUC_time = 5,#指定计算auc的时间
    auc_cal_method="KM")#计算auc的方法

cal_AUC_ml_res（）还提供了三种模式，它们应该与ML.Dev使用的模式一致。
AUC_time1年、2年、3年……，我们建议使用所有数据集中最短的生存时间。

auc_dis_all(object = all.auc.1y,#基准测试的结果
            dataset = names(list_train_vali_Data),#要展示AUC的数据
            validate_set=names(list_train_vali_Data)[-1],#验证集的名称
            order= names(list_train_vali_Data),#要展示AUC的数据顺序
            width = 0.35,
            year=1)

在不同数据集之间绘制特定模型的 ROC 曲线：

roc_vis(object = 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)

在不同数据集中绘制特定模型的1、3和5年AUC：




    
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))

1.3特定模型单变量cox回归的Meta分析

unicox.rs.res   cal_unicox_ml_res(
    res.by.ML.Dev.Prog.Sig = res_prog,
    optimal.model = "StepCox[forward] + plsRcox",
    type ='categorical')



metamodel   cal_unicox_meta_ml_res(input = unicox.rs.res)



meta_unicox_vis(metamodel,
                dataset = names(list_train_vali_Data))

AI知识补充

在 meta 分析中，随机效应模型和固定效应模型是两种常用的统计方法，用于合并多个研究的结果。它们有不同的假设和适用场景。

一、固定效应模型

基本假设：

假设所有纳入的研究都是来自同一总体，即研究间的差异仅由随机误差引起。
认为各个研究的效应量是固定的，只是由于抽样误差而有所不同。

计算方法：

通常采用加权平均的方法来合并效应量。权重通常是各个研究的方差的倒数，方差较小的研究给予较大的权重。
合并后的效应量是各个研究效应量的加权平均值，其方差是各个研究方差的加权和。

适用场景：

当研究间的异质性较小，即各个研究的结果比较一致时，固定效应模型较为合适。
如果有明确的理由认为所有研究都是在相同的条件下进行的，且研究间的差异可以忽略不计，也可以使用固定效应模型。

二、随机效应模型

基本假设：

假设纳入的研究来自不同的总体，研究间的差异不仅由随机误差引起，还存在真实的效应量差异。
认为各个研究的效应量是随机变量，服从某种分布。

计算方法：

除了考虑各个研究的内部方差外，还考虑了研究间的方差。
合并后的效应量是各个研究效应量的加权平均值，其权重是各个研究的内部方差和研究间方差的函数。

适用场景：

当研究间存在较大的异质性时，随机效应模型更为合适。
如果无法确定各个研究是否来自同一总体，或者研究间的差异可能较大，也应该使用随机效应模型。

三、选择策略

异质性检验：

在进行 meta 分析之前，通常需要进行异质性检验，以判断研究间的异质性程度。
常用的异质性检验方法有 Q 检验和 I²统计量。如果异质性检验结果显示研究间存在显著的异质性（如 I²>50%），则应考虑使用随机效应模型；如果异质性较小，则可以使用固定效应模型。

实际考虑：

除了异质性检验结果外，还应考虑实际情况。例如，如果纳入的研究来自不同的地区、人群或研究方法差异较大，即使异质性检验结果不显著，也可能需要使用随机效应模型。
另外，如果研究的目的是推广到更广泛的人群或情况，随机效应模型可能更合适，因为它考虑了研究间的差异。

总之，在 meta 分析中，选择固定效应模型还是随机效应模型需要综合考虑异质性检验结果和实际情况。如果研究间异质性较小，可以使用固定效应模型；如果异质性较大或存在不确定性，随机效应模型更为合适。

1.4 与先前发布模型的比较

rs.glioma.lgg.gbm   cal_RS_pre.prog.sig(use_your_own_collected_sig = F,#是否用自己收集的特征集
                      type.sig = c('LGG','GBM','Glioma'),
                      list_input_data = list_train_vali_Data)
                      
rs.glioma.lgg.gbm$Chen.33591634$Dataset1 %>% head()

cal_RS_pre. prog.sig（）将根据以前论文计算风险分数。
如果use_your_own_collected_sig设置为T，您应该提供一个包含签名信息的数据框。数据框的列名是模型、PMID、癌症、作者、Coef和符号。模型由第一作者的名字和论文的PMID组成。癌症使用缩写，如TCGA的格式。作者是第一作者的名字。Coef是每个变量的系数。符号是基因名称。否则，我们使用我们收集的胶质瘤模型。

小胖体会：实测按照上述的做法并不能使用自身的特征集合进行计算。

为了解决上述问题，我查看了作者原始使用的其构建的特征矩阵情况，下载地址如下：

https://github.com/l-magnificence/Mime/blob/main/data/pre.prog.sig.rda

pre.prog.sig$Glioma %>% head(10)

以下查看原始基因集文件的具体情况


pre.prog.sig$Glioma %>% 
  distinct(PMID,.keep_all = T) %>% 
  pull(PMID)

(symbol30810537   pre.prog.sig$Glioma %>% 
  filter(PMID == 30810537) %>% 
  pull(symbol))


(symbol31907037   pre.prog.sig$Glioma %>% 
  filter(PMID == 31907037) %>% 
  pull(symbol))

(symbol33879726     pre.prog.sig$Glioma %>% 
    filter(PMID == 33879726) %>% 
    pull(symbol))

# 利用测试数据的一个数据框作为自建数据集
sig_table   pre.prog.sig$LGG

rs.glioma.lgg.col   cal_RS_pre.prog.sig(use_your_own_collected_sig = T,#是否用自己收集的特征集
                      collected_sig_table = sig_table,
                      list_input_data = list_train_vali_Data)

小胖体会：实测代码可以运行，看来需要输入函数的是一个包含每个基因集基因coef系数的长数据

将特定模型的HR与以前发布的模型进行比较：

HR_com(object = rs.glioma.lgg.col,
       object2 = res_prog,
       model_name="StepCox[forward] + plsRcox",
       dataset=names(list_train_vali_Data),
       type = "categorical")

cal_cindex_pre.prog.sig() 函数将根据先前论文中的特征计算 C 指数，就像 cal_RS_pre.prog.sig() 函数一样。

将特定模型的C指数与以前发布的模型进行比较：

cindex_comp(object = cc.glioma.lgg.gbm,
            object2 = res_prog,
            model_name="StepCox[forward] + plsRcox",
            dataset=names(list_train_vali_Data))

auc.glioma.lgg.gbm.1                                             type.sig = c('Glioma','LGG','GBM'),
                                            list_input_data = list_train_vali_Data,AUC_time = 1,
                                            auc_cal_method = 'KM')

cal_auc_pre. prog.sig（）将根据以前论文中的特征集计算AUC，例如.prog.sig（）cal_RS_pre函数。 AUC_time参数与cal_AUC_ml_res（）的要求相同。

将特定模型的 AUC 与先前发布的模型进行比较：

auc_comp(auc.glioma.lgg.gbm.1,
         all.auc.1y,
         model_name="StepCox[forward] + plsRcox",
         dataset=names(list_train_vali_Data))

1.5 免疫浸润分析

完成风险分组后，用户可以对分组后的数据进行下游分析。在这里，我们将Mime与R包immunedeconv结合，以帮助用户快速预览免疫浸润情况。如果用户需要更精确的免疫浸润分析，他们可以自己微调参数。

devo

如果你想使用这个功能，你应该提前安装软件包immunedeconv。
TME_deconvolution_all（）包括10个反卷积方法（"quantseq"、"xcell"、"epic"、"abis"、"mcp_counter"、"估计"、"cibersor"、"cibersort_abs"、"estimate"、"consensus_tme"）来自::deconvolution_methods。默认情况下，反卷积方法设置为（"xcell"、"epic"、"abis"、estimate"、"cibersor"、"cibersort_abs"）。

展示结果

immuno_heatmap(res,
               devo,
               model_name="StepCox[backward] + plsRcox",
               dataset="Dataset1")

2. 构建二分类响应预测模型

load("Example.ici.Rdata")
load("genelist.Rdata")

list_train_vali_Data$training->df

res.ici   ML.Dev.Pred.Category.Sig(
    train_data = list_train_vali_Data$training,
    list_train_vali_Data = list_train_vali_Data,
    candidate_genes = genelist,
    methods = c('nb','svmRadialWeights','rf','kknn','adaboost','LogitBoost','cancerclass'),
    seed = 5201314,
    cores_for_parallel = 10
    )

ML.Dev.Pred.Category.Sig() 利用机器学习算法为二元变量开发预测模型。

在不同数据集中绘制不同方法的AUC：

auc_vis_category_all(res.ici,dataset = c("training","validation"),
                     order= c("training","validation"))

在不同数据集中绘制具体方法的ROC：

plot_listmethods 'nb','svmRadialWeights','rf','kknn','adaboost','LogitBoost','cancerclass')
for (i in methods) {
  plot_list[[i]]"training","validation"),
                                   order= c("training","validation"),
                                   anno_position=c(0.4,0.25))
}
aplot::plot_list(gglist=plot_list,ncol=3)

将 AUC 与其他已发表的与免疫治疗反应相关的模型进行比较：

auc.other.pre                                       train_data = list_train_vali_Data$training,
                                      cores_for_parallel = 32)

cal_auc_previous_sig（）将根据以前论文中免疫治疗反应的签名计算AUC。 •cores_for_parallel意味着您可以选择用于并行操作的内核。如果多核条件错误，请将cores_for_parallel设置为1。

auc.other.pre                                       train_data = list_train_vali_Data$training,
                                      cores_for_parallel = 32)

cal_auc_previous_sig（）将根据以前论文中免疫治疗反应的签名计算AUC。
cores_for_parallel意味着您可以选择用于并行操作的内核。如果多核条件错误，请将cores_for_parallel设置为1。

与之前免疫治疗反应的AUC做对比

auc_category_comp(res.ici,
                  auc.other.pre,
                  model_name="svmRadialWeights",
                  dataset=names(list_train_vali_Data))

3.核心特征选择

load("./Example.cohort.Rdata")
load("./genelist.Rdata")
res.feature.all $Dataset1,
                                              candidate_genes = genelist,
                                              mode = "all",nodesize =5,seed = 5201314 )

全部模式是指使用所有八种方法（RSF、Enet、Boruta、Xgboop、SVM-REF、Lasso、CoxBoost、StepCox）进行选择。单一模式是指只使用一种方法进行运行。如果使用单一模式，则应在八种方法中指定single_ml。由于SVM花费的时间太长，我们将用于选择的其他七种方法定义为all_without_SVM模式。
输出基因与患者预后密切相关，更高的筛查频率意味着更关键。

不同方法过滤的基因的upset图：

core_feature_select(res.feature.all)

绘制通过不同方法筛选的基因的等级

core_feature_rank(res.feature.all, top=20)

在这里，我们随机选择前两个基因来分析它们的相关性：

dataset_col"#3182BDFF","#E6550DFF")



    
corplot for (i in c(1:2)) {
  print(corplot[[i]]                               dataset=names(list_train_vali_Data)[i],
                               color = dataset_col[i],
                               feature1="PSEN2",
                               feature2="WNT5B",
                               method="pearson"))
}
aplot::plot_list(gglist=corplot,ncol=2)

根据不同数据集中特定基因的中位表达水平绘制患者生存曲线：

survplot "list",2) 
for (i in c(1:2)) {
  print(survplot[[i]]"PSEN2", 
                                        InputMatrix=list_train_vali_Data[[i]],
                                        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)