诊断模型的构建与分析（随机森林与人工神经网络） by 生信技能树

StatQuest生物统计学视频是一个很优秀的生物统计学教程，教程作者是Josh Starmer （个人博客https://statquest.org/），生信菜鸟图很早之前就推过相关的学习资源。而且还组建过学习小分队，给视频写配套笔记：

但是这样的机器学习知识并没有很好的应用到生物信息学领域，尤其是数据挖掘方向。

关于机器学习，在生物信息学领域可能仍然是R语言应用更为广泛，简单的一个R包或者函数即可。在肿瘤数据库挖掘里面，主要是生存模型构建最广泛。我的4个小时TCGA肿瘤数据库知识图谱视频教程，中共使用了四种算法构建模型的R语言实践版本：

不管用了那种算法，核心都只是几句代码而已。不管是R语言还是python语言，都不会有太大的区别。

恰好看到了2020我们数据挖掘学习班优秀学员在自己的公众号分享了一个机器学习在数据挖掘的实例，所以申请转载给大家欣赏！

文初友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

下面是bio小本本公众号原文

今天，和大家分享一篇生信文章的解读与复现——随机森林与人工神经网络联合诊断心衰模型的构建与分析。

Part1文献解读

1摘要

心力衰竭是一个全球性的健康问题，影响到全世界约 2600 万人。由于传统的心力衰竭诊断技术在实践中存在诸多局限性，有必要开发新的诊断模型来补充现有的诊断方法。

近年来，随着基因测序技术的进步和进步，人们发现了更多与心力衰竭相关的基因。利用基因表达综合(GEO)数据库中现有的基因表达数据，筛选出心力衰竭患者差异表达基因(DEGS)，并通过随机森林分类器鉴定出 6 个关键基因（HMOX 2, SERPINA 3, LCN6, CSDC 2, FREM 1，和 ZMAT1)。

这些基因中，CSDC 2, FREM 1，和 ZMAT1 从未与心力衰竭有关联。利用人工神经网络成功地建立了心力衰竭的诊断模型，并在公共数据集上验证了该模型的有效性。

2前言

心力衰竭(HF)是所有类型心脏病常见的一种慢性疾病。心衰实质上是由心脏功能异常引起的病理生理状态，心脏不能满足正常心脏压力下正常代谢所需的抽吸速度。HF 可分为两类：一种是伴有射血分数降低的 HF，另一种是保留射血分数(HFpEF)的 HF，这两类 HF 的发生和发展机制明显不同。

HFpEF 常发生于压力超负荷肥大疾病。与 HFrEF 相比，HFpEF 更有可能降低心脏储备。在 HFpEF 的发病机制中，心肌细胞本身的凋亡程度较小，而其特征性变化是异常成纤维细胞的增殖和细胞基质蛋白的积累。这是 HFpEF 与 HFrEF 最显著的区别。

对于临床常用的 HF 诊断技术有几个限制。脑钠尿肽/N 端脯氨酸型钠尿肽在各种非心衰疾病如肺动脉高压、肝硬化腹水、急性或慢性肾功能衰竭、感染和炎症中也可能升高，但 HFpEF 患者正常。超声心动图是另一种常用的心功能评估技术，它更多地依赖于专家的个人操作能力和诊断经验，使检查的可重复性差。此外，单纯测量 EF 值，难以识别 HFpEF 患者。因此，有必要开发新的诊断模型来补充现有的诊断方法。

近年来，第二代测序技术的迅速发展为识别与多种疾病相关的标记基因提供了基础，为建立新的与基因相关的 HF 诊断模型奠定了坚实的基础。在本研究中，我们在 GEO 数据库中筛选了 HF 与正常心肌标本之间的差异表达基因(DEGS)。在这些 DEG 数据的基础上，我们采用随机森林算法对 HF 中的关键基因进行了识别。然后，我们将这些关键基因输入到人工神经网络中，以构建 HF 的遗传诊断模型.

3方法&结果

GEOQuery 软件包用于下载数据以获得芯片数据集的表达谱和临床表型数据:GSE57345, GSE42955，和 GSE84796 和 rna-seq 数据集:GSE141910 和 GSE116250。从 GEO 数据库中获取相应平台芯片探针的相应标注信息。在芯片探针 ID 和基因符号转换过程中，发现多个探针对应于 1 个基因符号，在这种情况下，平均探针表达作为基因表达水平。使用 org.Hs.eg.db 包 (3.7.0 版)对 rna-seq 表达式配置文件进行基因 ID 转换。

R 软件包 limma 对 GSE57345 数据集 136 个正常和 177 个心衰样本进行了差异分析。limma 软件包使用经典的贝叶斯数据分析来筛选 DEGS。DEGS 的显着性标准设置在 P 值小于 0.05， logFoldChang (LogFC) 大于 1.5。用 phatmap 软件包绘制 DEGS 的热图。

用 R 软件包绘制聚类剖面图对相关基因进行 GO 功能富集分析和 KEGG 富集分析，使用了 Benjamini–Hochberg 校正方法，阈值设置为 P 值<0.01，q 值<0.01。为了避免 GO 富集结果中的冗余，我们对 GO 富集项进行了去重复，消除了基因重叠>0.75（详见附录 2）。确定三种显著富集 GO 术语(P < 0.05) 和显著富集通路 (P < 0.05)。

在这些结果中，与 HF 相关的生物学过程包括细胞外基质的组织、心脏收缩、巨噬细胞活化和细胞-基质粘附。所涉及的细胞成分包括含有胶原的细胞外基质。分子功能包括整体结合和其他重要功能。图 3B 显示了部分 GO 富集的术语和显著表达的差异基因。我们还对 DEGS 进行了 kegg 通路富集分析，展示了所涉及的重要的富集生物途径的结果和相应的 DEG。

将 281 个 DEG 输入到随机森林分类器中。为了找到最优参数 mtry(即指定节点中二叉树的最优变量数)，我们对 1-281 个变量中的所有可能数进行了递归随机森林分类，并计算了模型的平均错误率。图 4A 显示所有变量被选中时的平均错误率。最后，选取 6 作为变量数的参数。变量数目越少越好，带外误差(out-of-band error)越小越好。引用模型误差与决策树数之间的关系图(图 4B)选取 2000 棵树作为最终模型的参数，模型误差稳定。

在建立随机森林模型的过程中，从降低精度和减小均方误差的角度出发，对输出结果的变量重要性(基尼系数法)进行了测量。然后，我们确定了 6 个重要性大于 2 的 DEGS 作为后续分析的候选基因。图 4C 表明在这六个变量中，HMOX 2 和 CSDC 2 是最重要的，其次是 ZMAT1, SERPINA3, FREM1，和 LCN6。基于这六个重要变量，我们对 GSE57345 数据集进行了 k-均值无监督聚类。图 4D 结果表明，在 GSE57345 数据集 313 例样本中，这 6 个基因可用于患者和正常人的鉴别。其中，ZMAT1 和 FREM 1 基因在正常组织中呈低表达，在疾病组织中呈高表达。另一方面，SERPINA3, LCN6, HMOX2，和 CSDC2 属于另一组，在正常标本中高表达，在疾病标本中低表达。

我们使用了另一个数据集 GSE141910 建立了基于 neuralnet 包的人工神经网络模型。第一步是对数据进行预处理，以实现数据的规范化。其次，选择最小-最大值方法[0，1]，在训练神经网络之前按下分离缩放数据。在开始计算之前，对最大和最小数据值进行标准化处理，并将隐藏层数设置为 5 层。在参数的选择上，不存在使用多少层和多个神经元的固定规则。神经元的数量应该介于输入层大小和输出层大小之间，通常是输入层大小的三分之二。为了更有效地评价神经网络模型的结果，我们选择了一种 5 折交叉验证方法。数据集随机分为训练集和验证集。训练集的目的是确定候选指标的权重。利用验证集验证了基于基因表达和基因权重构建的模型评分的分类效率。得到疾病神经网络模型分类评分。

5 折交叉验证结果使用受试者工作曲线(ROC)显示模型分类性能(图 5A)。此外，还使用了一个混淆矩阵来评估预测的性能(表 1)。五折交叉验证结果的 ROC 曲线下面积接近 1(平均 AUC>0.99)，表明模型具有较强的鲁棒性(Robust)。因此，接下来，我们利用所有的数据来构建神经网络模型。

根据神经网络模型的输出结果(补充档案 4 和图 5B)，可以看出，整个培训是按照 1423 步骤进行的。终止条件为误差函数的绝对偏导数<0.01(达到阈值)。输出结果表明，该模型的权重范围为 −4.67~4.53。重量预测为 4.527373(HMOX 2)，−4.7670777(CSDC 2), 1.478590 (ZMAT1), 2.332519 (SERPINA 3)，−4.522891(FREM 1)，以及 1.940819(LCN6).

使用三个独立的验证数据集 GSE116250, GSE42955，和 GSE84796 在进行最大和最小标准化数据处理后，计算出三个评分，并对其分类效率进行评价，并对 AUC 进行比较。这三个评分结果如下：1)neuraHF，将本研究中所识别的 DGS 与神经网络中的权重相乘得到的分数；2)CD8K；3)TP 53，这是文献中报道的与 HF 疾病相关的特征基因。

图 6 显示了三个独立验证数据集的三个分数的比较。在 GSE116250 数据集(图 6A neuraHF、 CD8K 和 TP 53 的 AUC 分别为 0.991、0.683 和 0.597。neuraHF 的敏感性为 100%，特异性为 96%。在 GSE 42955 数据集(图 6B 神经高频、CD8K 和 TP 53 的 AUC 分别为 0.858、0.517 和 0.65。neuraHF 的敏感性为 80%，特异性为 95.8%。在核查结果中 GSE 84796 (图 6C neuraHF、 CD8K 和 TP 53 的 AUC 分别为 0.871、0.586 和 0.486。neuraHF 的敏感性和特异性分别为 85.7%和 80%。

4讨论

在本研究中，我们首次计算了与 HF 相关的 DEGS，并通过随机森林分类器获得了 6 个重要候选 DEGS。利用神经网络模型确定了相关基因的预测权重，构建了与心衰相关的 neuraHF 分类模型，并在三个独立样本数据集上评价了模型评分的分类效率。与其他两种 HF 相关的生物标志物相比，AUC 的分类效率较高，且具有较好的分类效率。然而，由于 rna-seq 在构建神经网络模型时所预测的权重在理论上更适用于 rna-seq 数据的疾病分类，GSE116250 数据集在验证结果中显示出最佳性能。同时，由于 GEO 数据库中缺乏 HFpEF 的基因数据，HFpEF 的遗传特征没有被纳入诊断模型的构建中，从而影响了 HFpEF 模型的诊断效果。

更有趣的是，在构建 HF 诊断模型的分析过程中，我们首次发现了三个关键基因(CSDC 2, FREM 1，和 ZMAT1)可能在 HF 的发病机制中起一定作用。

心脏标本的获取困难可能会减少 HF 的潜在应用。然而，我们目前的研究并不打算完全取代现有的诊断和治疗方法，而是为了补充这些方法。通常，目前的诊断标准和程序是基于来自 HFrEF 患者的数据。然而，仍不清楚这些是否完全适用于 HFpEF 患者。例如，使用这些非侵入性方法很难诊断 HFpEF 的轻微症状。然而，从我们的研究中得出的诊断模型可以通过及时的心脏活检来确定心力衰竭的可能性。因此，我们的方法具有一定的临床应用价值。显然，根据我们目前的研究结果，该模型的精度还有待进一步研究。

5总结

亮点：

1.使用了机器学习来筛选 hub gene

2.使用了深度学习来构建诊断模型

3.使用了多个数据集进行数据分析

4.文章思路清晰，逻辑明了

不足：

1.文章里包含明显错误如文中实际 logFC 是>0.5 而不是 1.5,GSE141910 在 GEO 中仅包含 366 个样本而非 399 个样本

2.流程图未描述 kegg,带外误差的英文为 out-of-bag error

3.本文临床意义不大

Part2文章复现

step1 start GEO

#数据下载
rm(list = ls())
library(GEOquery)
gse = "GSE57345"
Sys.setenv("VROOM_CONNECTION_SIZE" = 2621440*2)
eSet <- getGEO(gse,
               destdir = '.',
               getGPL = F)
#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
#boxplot(exp)
#exp = log2(exp+1)
#(2)提取临床信息
pd <- pData(eSet[[1]])
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl <- eSet[[1]]@annotation
save(eSet,gse,pd,exp,gpl,file = "step1output.Rdata")

因为 GSE57345 是常规的微阵列芯片，所以我们按照标准流程代码来就好，得到了 GSE57345 的表达矩阵。

step2 group id

#1.group_list(实验分组)和ids(芯片注释)
rm(list = ls())
load("step1output.Rdata")
library(stringr)
library(dplyr)
gr = pd$characteristics_ch1.1
k1 = str_detect(gr,"yes");table(k1)
k2 = str_detect(gr,"no");table(k2)

group = ifelse(k1,"disease",'normal');table(group)
group = factor(group,levels = c("normal","disease"))

#2.探针注释----
#library(devtools)
#install_github("jmzeng1314/idmap1")
library(idmap1)
temp=getIDs("GPL11532")

rowname <- data.frame(probe_id=rownames(exp))
gene <- left_join(rowname,temp,by='probe_id')
gene <- na.omit(gene)
exp <- exp[match(gene$probe_id,rownames(exp)),]
# 多个探针对应一个基因的话，取平均表达量最大的那个探针
result <- data.frame(ex=rowMeans(exp),
                        symbol=gene$symbol)

result <- arrange(result,symbol,desc(ex))
de <- result[!duplicated(result$symbol),]
exp <- exp[match(rownames(de),rownames(exp)),]
exp_sy <- data.frame(symbol=de$symbol,
                     exp)

save(eSet,pd,exp_sy,group,gpl,file = "step2output.Rdata")

将样本分组，用曾老师的 idmap1 包给探针注释，这里开始和文章不一样了。我们经常会遇到多个探针对应一个基因的情况，但是处理方法有很多种，我们选择取平均表达量最大的那个探针。得到 18837 个基因。

step3 pca

rm(list = ls())
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
exp = exp_sy[,-1]
group_list = group

dat=as.data.frame(t(exp))
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra)
# pca的统一操作走起
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group_list, # color by groups
                         #palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot

save(pca_plot,file = "pca_plot.Rdata")

做主成分分析图，看看该芯片中对照组和实验组的差异性。（它们两个分组并不是泾渭分明，所以才需要构建诊断模型）

step4 DEG

rm(list = ls())
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和group_list，不需要改
exp =exp_sy[,-1]
group_list = group
library(limma)
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
#2.加symbol列，火山图要用
deg <- cbind(deg,symbol=exp_sy$symbol)
head(deg)

#3.加change列,标记上下调基因
logFC_t=0.5
P.Value_t = 0.05
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
change = ifelse(k1,
                "down",
                ifelse(k2,
                       "up",
                       "stable"))
table(change)
deg <- mutate(deg,change)
ex_deg <- deg[k1|k2,]
save(group_list,deg,logFC_t,P.Value_t,ex_deg,
     file = "step4output.Rdata")

做差异分析，logfc 取 0.5，结果如下，得到 426 个差异基因。

step5 random forest

rm(list = ls())
library(tidyverse)
load(file = "step2output.Rdata")
load(file = "step4_output.Rdata")

df <- exp_sy[match(ex_deg$symbol,exp_sy$symbol),]
rownames(df) <- df$symbol
df <- data.frame(t(df[,-1])) %>%
  cbind(group,.)

library(randomForest)
set.seed(1423)
fit.forest <- randomForest(group~., data=df,
                           ntree=3000)
fit.forest
plot(fit.forest)
#找出使模型准确率达到最优所需要的树的数量
which.min(fit.forest$err.rate[,1])

set.seed(1423)
fit.forest_2 <- randomForest(group~., data=df,
                           ntree=49)
fit.forest_2
(127+173)/313
# 模型正确率为95.85%

varImpPlot(fit.forest_2)
result_df <- data.frame(importance(fit.forest_2, type=2)) %>%
  arrange(desc(.))

impo_df <- filter(result_df,MeanDecreaseGini>4)
impo_df
# 选取了比较重要的几个基因放入后续分析，
# 但是我觉得更好的方法应该是使用Boruta包进行特征选择。
# Kursa M., Rudnicki W. (2010), Feature Selection with the Boruta Package,
# Journal of Statistical Software, 36(11), 1 - 13

impo_exp <- exp_sy[match(rownames(impo_df),exp_sy$symbol),]
rownames(impo_exp) <- impo_exp[,1]

n=impo_exp[,-1]
dim(n)
library(pheatmap)
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(ggplotify)
heatmap_plot <- as.ggplot(pheatmap(n,show_colnames =F,
                                   show_rownames = T,
                                   cluster_cols = T,
                                   scale = 'row',
                                   annotation_col=annotation_col))
heatmap_plot

rownames(impo_exp)
# "ALG3"     "AGFG1"    "ADAMTS15" "ACE"      "ADGRD1"

save(list = ls(),file = "step5output.Rdata")

我们设随机种子为1423，然后做随机森林，确定最佳随机树数目为49

看看变量重要性统计图

选取了基尼指数平均减少量大于4的变量作为后续分析的重要变量。

并画了个热图。

step6 neuralnet

rm(list = ls())

filelist <- dir('gse141910')
file_path <- paste('gse141910',filelist,sep = '/')
data_ne <-  read.csv(file_path[1])

for (n in 2:length(file_path)) {
  a <- read.csv(file_path[n])
  data_ne <- merge(data_ne,a)
}

# 转换基因id
library(clusterProfiler)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)## 查看org.Hs.eg.db包含的数据
# browseVignettes("clusterProfiler") #查看教程
eg = bitr(data_ne$X, fromType="ENSEMBL", 
          toType="SYMBOL", OrgDb="org.Hs.eg.db")
impo_sy <- c("ALG3","AGFG1","ADAMTS15","ACE","ADGRD1")
nt <- eg[match(impo_sy,eg$SYMBOL),]
nt_exp <- data_ne[match(nt$ENSEMBL,data_ne$X),] %>% 
  .[,-1] %>% 
  scale()

rownames(nt_exp) <- impo_sy
library(tidyverse)
k <- str_detect(colnames(nt_exp),'P')
table(k)
# 240名病人，126名正常人
nt_gr <- rbind(nt_exp,group=k) %>% 
  t(.) %>% 
  data.frame()

nt_gr$group <- factor(nt_gr$group)

library(neuralnet)

fit <- neuralnet(group~.,nt_gr,
                 hidden = 5)
fit$result.matrix

plot(fit)
# 我也不知道做出来的结果对不对，反正看得不是很懂，
# 不知道那篇文章中的权重是怎么得出来的。
# 但是如果真要应用深入学习的方法来研究问题时，
# 应该多尝试一些不同的算法或r包来比较一下结果。
# 比如TensorFlow,MXNet，H2O包等。

save(list = ls(),file = "step6output.Rdata")

GSE141910是rna-seq芯片，处理原始数据可能有点麻烦，所以我就直接用了芯片上传者处理好的表达矩阵来做后续分析。

用Y叔的clusterProfiler包，和org.Hs.eg.db数据库来对探针进行注释，并提取出之前的5个重要变量。

然后构建神经网络模型，隐藏神经元数量设为5，和文章里的图比起来就丑好多了。

step7 auc

rm(list = ls())
load('step6output.Rdata')
library(ggplot2)
library(pROC)
source("ROC_Plot.R")
source("Theme_Publication.R")

HUB <- nt_gr
colnames(HUB)
# "ALG3","AGFG1","ADAMTS15","ACE","ADGRD1","group"
x1 <- HUB$ALG3
x2 <- HUB$AGFG1
x3 <- HUB$ADAMTS15
x4 <- HUB$ACE
x5 <- HUB$ADGRD1
y <- HUB$group

R1 <- roc(y,x1,ci=T); R2 <- roc(y,x2,ci=T); R3 <- roc(y,x3,ci=T)
R4 <- roc(y,x4,ci=T);R5 <- roc(y,x5,ci=T)

dat1 <- data.frame(Sensitivities = R1$sensitivities, 
                   FalsePositiveRate = (1-R1$specificities), 
                   Model = rep("ARG1",length(R1$sensitivities)))
dat2 <- data.frame(Sensitivities = R2$sensitivities, 
                   FalsePositiveRate = (1-R2$specificities), 
                   Model = rep("CXCL1",length(R2$sensitivities)))
dat3 <- data.frame(Sensitivities = R3$sensitivities, 
                   FalsePositiveRate = (1-R3$specificities), 
                   Model = rep("LRG1",length(R3$sensitivities)))
dat4 <- data.frame(Sensitivities = R4$sensitivities, 
                   FalsePositiveRate = (1-R4$specificities), 
                   Model = rep("MMP9",length(R4$sensitivities)))
dat5 <- data.frame(Sensitivities = R5$sensitivities, 
                   FalsePositiveRate = (1-R5$specificities), 
                   Model = rep("OSCAR",length(R5$sensitivities)))

dat <- rbind(dat1,dat2,dat3,dat4,dat5)
colnames(dat)[2] <- ("1-Specificities")

ggplot(dat,aes(`1-Specificities`,Sensitivities,colour = Model)) + geom_roc_plot()+
  scale_colour_Publication() + theme_Publication()+
  theme(panel.border = element_rect(colour = 'black'),legend.title = element_blank(),
        legend.position = c(0.9,0.15),legend.direction = 'vertical')

#做表
sy <- (c("ALG3","AGFG1","ADAMTS15","ACE","ADGRD1"))
auc <- data.frame(rn= c('low','auc','up'))
roc <- list(R1,R2,R3,R4,R5)

for (n in 1:length(sy)) {
  auc1 <- data.frame(roc[[n]]$ci)
  auc <- cbind(auc,auc1)
}
rownames(auc) <- auc[,1]
auc <- auc[,-1]
colnames(auc) <- sy
write.csv(auc,'auc.csv')

save(list = ls(),file = "step7output.Rdata")

# 这里只是举了个画auc的例子，并没有像文章一样去独立验证了，
# 因为我不知道神经网络模型怎么在其他芯片里验证并做roc曲线。

我并没有像原文一样去做auc，原文的思想是用GSE57345和GSE141910两个大样本芯片作为训练集，再用三个小样本GSE42995，GSE84796，GSE116250作为验证集。最近有点其他事耽误了，就没有做这一步，但是还是找了一下之前做auc的代码放上来了。

还有富集分析忘了做了，大家可以参考一下其他推文：批量GO-KEGG富集分析注释何必自己写脚本，看Y叔的神器

Part3最后

本文最关键的方法是随机森林和神经网络，这两种方法主要参考参考资料里的那两本书来的，其中的原理本文并未细讲，如要用这两种方法来分析自己的数据，最好还是看一下那两本书所对应的章节，这样才可以灵活运用。

另外，本文所讲文献及文献附件，还有代码、数据都放原创者的公众号（bio小本本）上了，关注后台回复：20211024 即可获取。

参考资料

[1]

Robert I. Kabacoff著，《R语言实战》第2版，人民邮电出版社。

[2]

Cory Lesmeister著，《精通机器学习：基于R》，人民邮电出版社。

[3]

生信技能树微信公众号。

文末友情推荐