零基础入门转录组数据分析——机器学习算法之boruta（训练模型）

---

您首先需要了解本贴是完全免费按实际案例分享基础知识和全部代码，希望能帮助到初学的各位更快入门，但是 尊重创作和知识才会有不断高质量的内容输出 ，如果阅读到最后觉得本贴确实对自己有帮助，希望广大学习者能够花点自己的小钱支持一下作者创作（条件允许的话一杯奶茶钱即可），感谢大家的支持~~~~~~ ^_^ !!!

注：当然这个并不是强制的哦，大家也可以白嫖~~，只是一点点小的期盼！！！

祝大家能够开心学习，轻松学习，在学习的路上少一些坎坷~~~

---

---

1. boruta基础知识

1.1 训练模型是什么？
训练模型指的是利用特定的算法和大量的数据，让计算机生成一个可以解决特定问题的模型，这个过程的目的是使模型能够对新的、未见过的数据进行准确的预测或分类

1.2 训练模型和上一章筛选特征基因有什么联系？
上一章筛选完特征基因后就可以用筛选后的那些特征基因作为输入来训练模型

1.3 训练模型能干什么？
其可以在筛选完特征基因后验证特征筛选的有效性，即通过在一个新的数据集（通常是测试集或验证集）上训练模型，并评估其性能，可以验证特征筛选是否真正提高了模型的泛化能力。同时还可以看看用这些特征基因构建模型对于预测和分类的准确性。

举个栗子： 在上一章中筛选出来了7个特征基因，想进一步了解下这7个特征基因在什么样的参数条件下能够构建出最优模型，以及构建出来的模型它的预测准确性是多少，就可以在训练模型中找到答案。

注意：关于特征基因的验证通常不是用模型去验证的，通常可以直接用实验的手段去检测基因的表达水平，因此在这一章中仅介绍如何训练模型，对于如何验证不做过多介绍。

这一章可以作为拓展知识大家了解即可（前面的处理方式都与上一章一致）

---

---

2. boruta（Rstudio）——代码实操

本项目以TCGA——肺腺癌为例展开分析
物种：人类（Homo sapiens）
R版本：4.2.2
R包：tidyverse，Boruta，caret，cowplot

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./11_boruta_train')){
  dir.create('./11_boruta_train')
} 
setwd('./11_boruta_train/') 

加载包：

library(tidyverse)
library(Boruta)
library(caret)
library(cowplot)

导入要分析的表达矩阵train_data ，并对train_data 的列名进行处理（这是因为在读入的时候系统会默认把样本id中的“-”替换成“.”，所以要给替换回去）

train_data <- read.csv("./data_fpkm.csv", row.names = 1, check.names = F)  # 行名为全部基因名，每列为样本名
colnames(train_data) <- gsub('.', '-', colnames(train_data), fixed = T)

导入分组信息表group

group <- read.csv("./data_group.csv", row.names = 1) # 为每个样本的分组信息（tumor和normal）
colnames(group) <- c('sample', 'group')

导入要筛选的基因hub_gene (8个基因)

hub_gene <- data.frame(symbol = gene <- c('ADAMTS2', 'ADAMTS4', 'AGRN', 'COL5A1', 'CTSB', 'FMOD', 'LAMB3', 'LAMB4'))
colnames(hub_gene) <- "symbol"

从全部的基因表达矩阵中取出这8个基因对应的表达矩阵，并且与之前准备的分组信息表进行合并

dat <- train_data[rownames(train_data) %in% hub_gene$symbol, ] %>%
  t() %>%
  as.data.frame() # 整理后行为样本名，列为基因名
dat$sample <- rownames(dat)
dat <- merge(dat, group, var = "sample")
dat <- column_to_rownames(dat, var = "sample") %>% as.data.frame()

table(dat$group)
dat$group <- factor(dat$group, levels = c('disease', 'control'))

2. 2 构建boruta模型

设置随机种子并运行Boruta特征选择

set.seed(123)
boruta.train <- Boruta(group~., data = dat, doTrace = 2, maxRuns = 500)

boruta.train结果如下图所示，运行了499次迭代，耗时20.38秒，其中有7个特征是确定的，有一个特征—LAMB4是不确定的。

修正Boruta结果（这一步是可选的，主要用于修正Boruta算法可能产生的某些不稳定的特征选择结果）

final.boruta <- TentativeRoughFix(boruta.train)

修正后的final.boruta结果如下图所示，其余没什么变化，唯独那个不确定的LAMB4变成了不重要的

------------------------------注意：从这一步开始就和上一章不一致了----------------------------

定义一个函数用来生成用于随机森林模型调参的mtry值列表。

• mtry值——是随机森林（Random Forest）算法中的一个重要参数，它决定了每个决策树在进行特征选择时考虑的特征数量。具体来说，mtry参数定义了每个决策树节点可用于分割的特征的最大数量

那么选择什么样的mtry值是合适的？（mtry值的选择）

• （1）经验值： 通常，mtry的取值范围是0到p，其中p是特征的总数。取特征总数的平方根是一个常见的经验值，但这并不是绝对的，实际选择时需要根据数据和任务进行调整。
• （2）交叉验证： 为了找到最佳的mtry值，可以通过交叉验证等方法进行调优。例如，在R的randomForest包中，可以使用rfcv()函数来执行随机森林交叉验证，从而选择最佳的mtry值。
• （3）自适应选择： 有些随机森林实现支持自适应地选择mtry的值。在每次分裂时，算法会动态调整mtry的大小，以便更好地适应数据。

# 交叉验证选择参数并拟合模型,定义一个函数生成一些列用来测试的mtry (一系列不大于总变量数的数值)。
generateTestVariableSet <- function(num_toal_variable){
  max_power <- ceiling(log10(num_toal_variable))
  tmp_subset <- c(unlist(sapply(1:max_power, function(x) (1:10)^x, simplify = F)), ceiling(max_power/3))
  #return(tmp_subset)
  base::unique(sort(tmp_subset[tmp_subset<num_toal_variable]))
}

函数的具体作用分解如下：

• （1）输入num_toal_variable——这个表示特征的总数
• （2）计算最大幂次——通过 max_power <- ceiling(log10(num_toal_variable)) 计算特征总数的以10为底的对数的上整，这个值用于确定在生成候选 mtry 值时要考虑的幂次范围
• （3）生成幂次序列——使用 sapply 函数和匿名函数 (1:10)^x 生成一个包含从1的幂到max_power的幂（每个幂次都取1到10的幂次方）的列表。这个列表通过 unlist 展开，以便形成一个单一的向量
• （4）筛选和排序——从生成的幂次序列中筛选出小于或等于特征总数的值，以确保 mtry 值不会超过可用的特征数量。然后，使用 sort 函数对这些值进行排序，并使用 unique 函数去除任何重复的值，以生成一个有序的、无重复的 mtry 候选值列表。
• （5）输出结果——最终返回这个排序且去重的 mtry 候选值列表

定义好计算mtry值的函数后就要准备输入的数据（从boruta模型中提取出通过筛选的基因）

#提取重要的变量和可能重要的变量
boruta.finalVarsWithTentative <- data.frame(Item = getSelectedAttributes(final.boruta, withTentative = T), 
                                            Type = "Boruta_with_tentative")

boruta.finalVarsWithTentative如下图所示，第一列item就是那7个通过检验的特征基因，第二列type就是对应的类型

之后从全部的基因表达矩阵中提取出这7个基因对应的表达矩阵

boruta_train_data <- dat[, boruta.finalVarsWithTentative$Item]

boruta_train_data 如下图所示，行名为样本名，每一列对应一个基因。

之后将这个boruta_train_data 作为参数，输入到自定义函数中去计算相应的mtry值列表

boruta_mtry <- generateTestVariableSet(ncol(boruta_train_data))

可以看出输出mtry值列表有6个数值，分别是123456，不同的mtry值计算出来的结果不一样

接下来就要为训练模型做准备，首先是设置交叉验证方式

这里采用trainControl函数去设置：
（1）method=“repeatedcv”——指定了交叉验证的方法为“重复交叉验证”（repeatedcv）
（2）number=10——表示每次交叉验证将数据集分成10
（3）repeats=5——表示这个过程将重复5次

trControl <- trainControl(method = "repeatedcv", number = 10, repeats = 5)

接下来就要根据前面计算出来的mtry值列表生成一个网格，这个网格将用于在训练过程中测试不同的mtry值，以找到最优的mtry

tuneGrid <- expand.grid(mtry = boruta_mtry)

做好上述准备后就可以开始通过train函数尝试训练模型

• x——表示输入的特征数据（boruta_train_data）
• y——表示指定的目标变量（可以理解成分组）
• method——使用的模型类型（这里是"rf"，即随机森林）。
• tuneGrid——用于调优的网格（这里是前面计算出来的mtry不同值）
• metric——用于评估模型性能的度量标准
• trControl——交叉验证的控制参数

train_data_group <- as.factor(dat$group)
borutaConfirmed_rf_default <- train(x = boruta_train_data, 
                                    y = train_data_group, 
                                    method = "rf", 
                                    tuneGrid = tuneGrid, 
                                    metric = "Accuracy", 
                                    trControl = trControl)

注意：在metric参数中，如果分组是一个二元分类问题的目标变量，需要选择"Accuracy"（准确率）或"ROC"（对于二分类问题的AUC值）。如果分组是一个具有多个类别的分类问题的目标变量，应该选择"Kappa"或"Accuracy"。

borutaConfirmed_rf_default如下图所示，介绍训练模型相对应的信息，比如：7个特征（7 predictor），二分类变量（2 classes），以及不同mtry对应的模型精确性，最后发现当mtry = 2时，模型为最优模型，其预测准确性为0.9613（最高）。

明确了最优模型之后，就要从模型中进一步提取出筛选后的基因及其对应的重要性并保存输出

gene_importance <- varImp(borutaConfirmed_rf_default)
gene_importance <- gene_importance$importance
gene_importance$symbol <- rownames(gene_importance)
colnames(gene_importance) <- c('importance', 'symbol')
gene_importance <- dplyr::select(gene_importance, symbol, importance)
gene_importance <- gene_importance[order(gene_importance$importance, decreasing = T), ]
write.csv(gene_importance, file = 'gene_importance_all.csv')

2. 3 训练模型结果简单可视化

接下来一步就是要对训练模型结果进行简单可视化，毕竟文章里是要放图的!!!

{
  p1 <- plot(borutaConfirmed_rf_default, xlab = "", xaxt = "n")
  
  p2 <- ggplot(gene_importance, 
               aes(x = reorder(symbol, importance), y = importance , fill=importance)) + #x、y轴定义；根据ONTOLOGY填充颜色
    geom_bar(stat="identity", width=0.8) +  #柱状图宽度
    scale_fill_continuous(low = '#8ECFC9', high = '#72B063')+  #柱状图填充颜色
    coord_flip() +  #让柱状图变为纵向
    xlab("") +  #x轴标签
    ylab("Importance") +  #y轴标签
    labs(title = "Boruta_analyse")+  #设置标题
    theme_bw() +
    theme(panel.grid = element_blank(),
          axis.text.x = element_text(family="Times", face = "bold", color="gray50", size = 14),
          axis.text.y = element_text(family="Times", face = "bold", color="black", size = 14),
          legend.title = element_text(family="Times", face = "bold", size = 15), # 图例标题大小
          legend.text = element_text(family="Times", face = "bold", size = 15), # 图例标签文字大小
          axis.title = element_text(family="Times", face = "bold", size = 15), # 坐标轴标题大小
          legend.position = "right",
          plot.title = element_text(family="Times", face = "bold", size = 15, hjust = 0.45), # 标题设置
          panel.background = element_rect(fill = "#f2f2f2") # 面板背景颜色
    )
  
  combined_plot <- plot_grid(p1, p2, ncol = 2, rel_heights = c(1, 4), align = 'v')
  
}

结果如下图所示，左边的图就是在训练模型的时候不同mtry对用的预测准确性，可以直观的看出当mtry = 2的时候此时模型预测准确性最高，右边的图就是输出了用于构建模型那7个基因的重要性。

---

---

结语：

以上就是boruta算法筛选完关键基因后训练模型的所有过程，如果有什么需要补充或不懂的地方，大家可以私聊我或者在下方评论。

如果觉得本教程对你有所帮助，点赞关注不迷路！！！

---

• 目录部分跳转链接：零基础入门生信数据分析——导读