零基础入门转录组数据分析——机器学习算法之xgboost（筛选特征基因）

零基础入门转录组数据分析——机器学习算法之xgboost（筛选特征基因）

---

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

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

---

---

1. xgboost基础知识

1.1 xgboost是什么？
XGBoost（eXtreme Gradient Boosting）是一种基于梯度提升（Gradient Boosting）框架的高效、灵活且广泛使用的机器学习算法，通过多次迭代后旨在找到最优预测模型。

1.2 xgboost的算法原理是什么？
XGBoost的核心思想是通过组合多个弱学习器（通常是决策树）来形成一个强学习器。在每一次迭代中，XGBoost都会基于当前模型的预测结果和真实值之间的差异（残差）来训练一个新的弱学习器，然后将这个新的弱学习器加入到已有的模型集合中，以期望能够减少整体的预测误差，从而找到一个最优模型。

1.3 xgboost的优势？

• 高效性： XGBoost在算法层面进行了多项优化，包括使用高效的树学习算法、支持并行和分布式计算、利用CPU缓存进行加速等，这使得它虽然需要经过多次迭代但是速度相对会很快
• 准确性： XGBoost通过引入二阶泰勒展开来近似目标函数，并使用正则化项来控制模型的复杂度，从而能够更准确地拟合数据，提高模型的预测精度。
• 灵活性： XGBoost支持多种类型的损失函数和弱学习器，同时也支持自定义评估函数，这使得它能够适应不同的机器学习任务。
• 鲁棒性： XGBoost通过加入正则化项和使用贪心算法进行节点分裂，来防止过拟合，提高模型的泛化能力

1.4 xgboost的缺点？
在特征选择的时候只能输出每个特征的重要性，无法自动筛选出关键特征，需要人为选择关键特征，给最终结果添加部分主观性。

1.5 xgboost筛选基因的本质是什么？
XGBoost在训练过程中会自动评估每个特征的重要性（重要性通常基于该特征在构建决策树时对模型性能提升的贡献程度），之后人为的手动筛选关键特征。

举个栗子： 有8个基因构建xgboost模型来预测患病的概率，算法结果会输出每个基因的重要性评分，这样就知道哪些基因对预测结果影响最大，之后就人为选择一些关键基因。

综上所述： xgboost先通过多次迭代达到最优模型，然后输出每个特征的重要性，这时就可以根据重要性排序去获取那些比较重要的特征。

---

---

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

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

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(tidyverse)
library(Matrix)
library(xgboost) 
library(ggplot2)

导入要分析的表达矩阵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)

train_data 如下图所示，行为基因名（symbol），列为样本名

导入分组信息表group

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

group 如下图所示，第一列sample为样本名，第二列为样本对应的分组 （分组为二分类变量：disease和control）

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

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

hub_gene 如下图所示，只有一列：8个基因的基因名

从全部的基因表达矩阵中取出这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'))

dat 如下图所示，行为基因名，前8列为基因对应的表达矩阵，第9列为合并的分组信息表

2. 2 构建boruta模型

前置准备：将前8列转换成matrix矩阵

# 将trainset的1-8列（自变量）转换为矩阵
traindata1 <- data.matrix(dat[, -9])

traindata1 如下图所示，和dat类似，区别就在于traindata1没有第9列分组信息，并且是个矩阵，不是数据框。

之后将traindata1转换成稀疏矩阵traindata2，这一步是为了进一步减少内存使用。

# 利用Matrix函数，将sparse参数设置为TRUE，转化为稀疏矩阵
traindata2 <- Matrix(traindata1, sparse = T)

指定dat的最后一列分组信息作为因变量，并转为数值型

# 将因变量转换为numeric类型，-1是为了从0开始计数
train_y <- as.numeric(dat[, 9])-1

将准备好的稀疏矩阵traindata2和因变量train_y打包成一个列表（这么做是为了方便后续用数据构建xgboost模型）

# 将自变量和因变量拼接为list
traindata <- list(data = traindata2, label = train_y)

将traindata列表里的表达矩阵traindata2和因变量train_y传入到xgb.DMatrix函数中转换成一个DMatrix对象（这个对象是xgboost算法内部用于存储数据和标签的优化格式，可以提高计算效率）

dtrain <- xgb.DMatrix(data = traindata$data, label = traindata$label)

之后就是构建xgboost模型

• data = dtrain —— 这个参数指定了要训练的数据，传入的应该是个DMatrix对象
• max_depth = 5 ——这个参数定义了树的最大深度，增加这个值会使模型更加复杂，可以关注到更多细节，但同时也会导致过拟合。
• eta = 0.3 ——这个参数通常被称为学习率/收缩率。它控制了在每一步中如何更新权重，从而减少每次迭代对模型的影响。较小的eta值就意味着更小的权重更新，可以使模型训练变慢，从而得到更好的结果和更好的泛化能力。
• objective = ‘binary:logistic’ ——这个参数指定了学习任务的类型和相应的损失函数。binary:logistic表明是一个二分类问题，并且使用的是逻辑回归损失函数。xgboost算法支持很多不同的objective，适用于不同的情况，具体在这里不展开介绍，感兴趣的小伙伴可以自行查询函数帮助。
• nround = 25 —— 这个参数指定了迭代次数，数值越大迭代次数越多，会使模型变得更复杂，同时也会增加过拟合风险和训练时间。

res.xgb <- xgboost(data = dtrain, max_depth = 5, eta=0.3, objective='binary:logistic', nround = 25)

构建完模型之后此时就可以从模型中提取变量的贡献度了

xgb_importance <- xgb.importance(train_matrix@Dimnames[[2]], model = res.xgb)    ##特征重要度

xgb_importance 如下图所示

• Feature——这一列就是基因的名称
• Gain——这一列表示每个特征在模型构建过程中平均增益的贡献，较高的增益值意味着该特征对于提升模型预测性能有更大的作用（重要） 。
• Cover——这一列可以理解成平均覆盖度（不需要关注）
• Frequency——这一列表示每个特征在树的构建过程中被用作分裂节点的频率，简单来说，值越大表明频率越高，也说明该特征在模型中被频繁使用，同样的对模型的预测能力有重要的贡献。（不需要关注）

接下来从构建好的xgboost模型中提取出最重要的5个基因并输出

## 结果
hub_gene <- xgb_importance[c(1 : 5), ]
hub_gene$Feature <- gsub('.','-',hub_gene$Feature,fixed = T )
write.csv(hub_gene, '01.hub_gene.csv')
write.csv(xgb_importance, '02.xgb_importance.csv')

2. 3 xgboost结果简单可视化

接下来一步就是要对xgboost结果进行简单可视化，毕竟文章里是要放图的，并且图片展现的效果会更好!!!

# xgboost结果简单可视化（ggplot2函数）
ggplot(xgb_importance, aes(x= reorder( Feature,Gain), y=Gain,fill=Feature)) +geom_bar(stat="identity") +theme_classic() +guides(fill=FALSE)+#theme(legend.position = )+scale_fill_manual(values=c("#8DD3C7","#FFFFB3","#BEBADA","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F","#377EB8","#100EB2","#FDB999"))+coord_flip()+theme_bw()+ggtitle('XGBoost')+theme(plot.title = element_text(size=24,color='black', face = "bold",family='Times'),axis.title.x =element_text(size=18,color='black', face = "bold",family='Times'),axis.text.x =element_text(size=16, color='black', face = "bold",family='Times'),axis.title.y =element_blank(),axis.text.y=element_text(size=16,   color='black',face = "bold",family='Times'),legend.title=element_text(size=20, color='black', face = "bold",family='Times'),legend.text=element_text(size=18, color='black', face = "bold",family='Times'),title=element_text(size=20, color='black', face = "bold",family='Times'),strip.text = element_text(size = 14,family = "Times", face = "bold"))+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+labs(x="gene",y="Gain",fill="")

结果如下图所示横坐标为贡献度（也可以叫成重要性），纵坐标为基因，柱子越高表明该基因越重要。

---

---

结语：

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

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

---

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