2024-06-16

生信大杀器之神奇的贝叶斯网络（6）

小师妹生信果 2024-03-07 19:00:15

转自公众号：R语言学徒
http://mp.weixin.qq.com/s?__biz=Mzg5MDk3Mzg4OA==&mid=2247487163&idx=1&sn=e474e0d41a8cd041f7e068a2f0a2c3c7

哈喽，小果又和大家见面啦，经过前面5章的学习，不知道大家有没有掌握贝叶斯网络这一机器学习方法呢，其在生物信息学上的应用是非常关键的，今天，小果就带大家以一个切实的例子来进行一遍完整的复习！

小果选取的数据集来自于开放数据集GEO，标题为厌氧-有氧转变后大肠杆菌的RNA-Seq，源于 von Wulffen 等人的研究，研究了大肠杆菌（E. coli）在从厌氧条件转变为好氧条件时的基因表达变化，研究使用大肠杆菌K12菌株W3110。将细胞在pH7和37°C的定义培养基中在搅拌的3升生物反应器中厌氧生长，直到培养物达到600的OD（3nm）。此时抽取第一个样品，随后在曝气开始后1.0、5、1、2和5分钟以10l/min开始曝气，抽取其他样品。在这个例子中，小果将带大家学习在变量全是连续数据的情况下该怎么通过贝叶斯网络建模。

简单来说，现在有名（简记）为a、b、c的三个RNA基因的表达状况，每一个都分别记录了在6个时刻下的不同表现，在今天的学习中，小果将带大家应用贝叶斯网络的知识，探索RNA基因的内在结构，也就是三个基因有无相互关联，并给出一个简单的模型，通过该模型，还可以应用其他基因来预测我们所需要的基因表达。

导入R包

首先，让我们载入需要的相应R包，其中edgeR用来创建所需要的基因数据框，bnlearn用来进行贝叶斯网络的学习和预测

1.# 加载所需的R包  2.library(edgeR)  3.library(bnlearn)    加载数据   现在我们开始加载需要的数据，通过read_excel（）函数读取相应的xlsx文件1.# 加载数据  2.# 使用readxl包中的read_excel函数读取Excel文件中的数据，文件路径需要根据你的实际情况修改。3.wulffenTable <- read_excel("C:/Users/17240/Desktop/GSE71562_E14R012_raw_counts.xlsx")  4.# 显示数据的前几行  5.head(wulffenTable)  6.  7.# 选择需要分析的数据，这里选择了第2到19列的数据  8.X = wulffenTable[2:19]  9.colnames(X)<-c("at0","at0.5","at1","at2","at5","at10","bt0","bt0.5","bt1","bt2","bt5","bt10","ct0","ct0.5","ct1","ct2","ct5","ct10")

根据先验知识预设贝叶斯网络结构

我们可以根据先验知识来预设贝叶斯网络的结构，例如我们知道在同一基因例如a的情况下，t0时刻应为t0.5时刻的predictor而不是相反，同理t0.5因为t1时刻的predictor…，在这种条件下，我们可以进一步简化模型，认为基因随着时间的变换具有一阶马尔可夫关系，也就是说:

.根据这个条件，我们可以预设贝叶斯网络的结构：

1.#我们已知基因表达的时间顺序，故可以预设马尔可夫条件关系  2.w=data.frame(from=c("at0","at0.5","at1","at2","at5","bt0","bt0.5","bt1","bt2","bt5","ct0","ct0.5","ct1","ct2","ct5"),      3.                    to=c("at0.5","at1","at2","at5","at10","bt0.5","bt1","bt2","bt5","bt10","ct0.5","ct1","ct2","ct5","ct10"))

创建贝叶斯网络并调整

现在，我们开始正式地通过数据学习贝叶斯网络的结构，常用地，这里我们选用bnlearn中的gs()函数进行学习，该函数采用的是贪婪算法，贪婪算法在贝叶斯网络的构建过程中被广泛应用，用于逐步选择合适的变量顺序，以便逐步构建网络结构。贪婪算法的核心思想是每次选择最佳的变量来构建网络，以最小化某种准则（如信息准则、似然函数等）的值。

下面是贝叶斯网络中的贪婪算法的一般步骤：

1.选择起始节点：算法从一个空网络开始，选择一个变量作为起始节点，可以是根据领域知识或随机选择。

2.逐步添加变量：算法根据一定的评价准则，逐步添加其他变量到网络中。在每一步中，选择一个未加入网络的变量，将其添加到网络中，并与已存在的变量建立概率依赖关系。

3.选择依赖关系：在每次添加新变量时，算法需要确定该变量与哪些已存在的变量之间建立概率依赖关系。这通常涉及评估变量之间的相关性，以及使用统计方法估计概率分布。

4.评价网络结构：在每次添加变量后，使用某种评价准则来度量当前网络结构的好坏。常用的评价准则包括贝叶斯信息准则（BIC）、最大似然估计（MLE）等。算法会尝试不同的变量顺序和依赖关系来优化这些准则。

5.迭代优化：重复步骤2到步骤4，直到满足一定的终止条件，例如网络结构不再发生改变或达到预定的变量数量。

1.pdag = gs(X,whitelist = w)  2.pdag    3.pp = graphviz.plot(pdag)

上图就是我们学习到的贝叶斯网络结构，但是小伙伴们有没有注意到，bt0.5在ct0上面，根据我们的常识，一个在0.5时刻的变量是不会出现在0时刻之上的，所以相应的，我们需要做出一定的修改，在gs方法中加入黑名单，屏蔽bt0.5到ct0的一段有向弧：

1.pdag = gs(X,whitelist = w,blacklist = c("bt0.5","ct0"))  2.pdag4.pp = graphviz.plot(pdag)

上图是经我们修改过的贝叶斯网络，已经差不多了！但是还是存在一个小问题，那就是at5和ct5、bt10和ct10之间分别存在一条无向弧，无向弧的存在使我们无法估计贝叶斯网络的具体参数，由于这两对变量之间不存在明显的因果关系，故而可以先预设其中弧的方向：

1.pdag = set.arc(pdag, from = "at5", to = "ct5")  2.pdag = set.arc(pdag, from = "bt10", to = "ct10")  3.pp = graphviz.plot(pdag)

这样，我们需要的一个完整的贝叶斯网络就成功创建了，通过该网络，可以直观地反应基因之间的关系，其中仅包含a、b、c的三条子图各自构成一条马尔科夫链，而a、b、c三者之间也存在着一定的联系，这就是贝叶斯网络对数据的结构的挖掘作用，不止于此，我们还可以通过学习贝叶斯网络的参数进而对所需要的节点进行预测！

参数学习

1.#拆分数据为训练数据集和测试数据集  2.wul.train=X[1:4000,]  3.wul.test=X[4001:4319,]  4.set.seed(123)5.fitted = bn.fit(pdag, wul.train)  6.fitted

通过bn.fit函数学习贝叶斯网络的参数，其中第一个输入是我们已有的贝叶斯网络对象，第二个是需要放入的训练数据

预测

在这里，小果运用ct5节点像大家演示，即通过除了ct5外的其他变量的取值和贝叶斯网络结构来预测ct5变量的取值，这里采用的是蒙特卡洛后验推理的方法进行预测：

1.pred = predict(fitted, node = "ct5", data = wul.test, method = "bayes-lw")  2.cbind(wul.test$ct5,pred)  3.plot(wul.test$ct5-pred)

根据如上所示的表格和预测残差图可以看出，预测值大部分浮动在真实值两侧，残差基本稳定在0附近，这说明我们采用的简单贝叶斯网络结构所作的预测是合理的，预测的效果显著。

那么贝叶斯网络的学习就到此结束啦，在小果的陪伴下，相信小伙伴们已经基本掌握了贝叶斯网络的相关知识，也会通过R语言编写相应的程序，小果要和大家说再见咯，一定要自己练习一下哦，同时如果大家想要继续了解更多有关R语言内容可以持续关注小果哦～～或者也可以关注我们的官网也会持续更新的哦～ http://www.biocloudservice.com/home.html

小果还提供思路设计、定制生信分析、文献思路复现；有需要的小伙伴欢迎直接扫码咨询小果，竭诚为您的科研助力！