从基因到表型（3）-实战之一文搞定时空生长模拟【含安装与一键式分析脚本】

阅前须知

在本期中将会使用到gapseq与bacarena，这两款软件均用到了sybilSBML库，但由于R的CRAN库管理要求，sybilSBML的作者并不想在测试服务器上安装最新版本的libsbml，导致sybil从CRAN上下架，无法通过install.package的方式安装。如果不能成功安装sybilSBML库，则会导致gapseq无法正常输出xml格式的模型，BacArena也无法正确的读取xml的模型。不过，如果你只会使用gapseq与bacarena做上下游的连贯分析，两者可以完全通过读取RDS文件来解决。所以在全文的最后，我提供了付费解锁的Gapseq与BacArena快速安装与部署的自制脚本（包含sybilSBML）。当然，如果有兴趣的朋友也可以自行参考两款软件的官方指南进行安装，丰俭由人。

前言与背景

在前两期内容中，我介绍了目前基因组规模代谢模型的一些基本知识和常见软件的选择策略。本期将基于Gapseq和BacArena进行交叉喂养分析，这两款软件目前在google scholar上的引用量分别是119和237。能够实现从基因组建模到单菌或多菌互作的时空模拟，并且对菌株（群）的生长过程和交叉喂养状态进行可视化。在文章中，我特别标明了Gapseq与BacArena的常用分析命令，帮助大家快速定位和使用。并且提供了一个基于真实数据的抑菌预测，为大家提供一些分析应用上的思路。

创建一个基因组规模代谢模型草图

gapseq的分析流程

根据官方文档介绍，gapseq具有两种用法，一是用来预测代谢途径，二是可以重建代谢网络，并且两种功能实际上是互相关联的。

下图则总结了该流程的主要用法：

gapseq的输入通常是以基因组DNA序列，文件可以是直接的fa、fna或者压缩后的基因组文件；（请注意同样是代谢重建软件的Carveme，要求的输入文件是蛋白质序列）

gapseq在接收到文件输入后，首先运行的是代谢途径和转运蛋白的预测gapseq find 与gapseq find-transport ；

预测完成代谢产物后，gapseq便可进行代谢模型草图的搭建gapseq draft和对草图模型的间隙填充gapseq fill，以解决草图模型无法在常见生长培养基上生长的问题。

需要说明的是，目前填充算法实现的模型离替代手动管理模型还有很长的距离，如果你的目标是精细化管理而不是探索分析，那么更建议在获取到模型之后，再使用基于matlab的DEMETER或者基于python的refinegems对手头的模型进行更进一步的调教。

gapseq的用法示例

---

1.预测代谢途径

在代谢途径预测中，可以使用参数-p来指定全部代谢途径的预测或是特定的产物预测

toy/ecoli.fna.gz 文件可在gapseq安装目录下找到

# 预测所有代谢途径
gapseq find -p all toy/ecoli.fna.gz

# 预测特定的化合物(几丁质)
gapseq find -p chitin toy/ecoli.fna.gz

# 预测酶的可用性(细胞色素C氧化酶)
gapseq find -e 1.9.3.1 toy/ecoli.fna.gz

# 按名称搜索酶
gapseq find -r "dye-decolorizing peroxidase" toy/ecoli.fna.gz

# 搜索转运蛋白
gapseq find-transport toy/myb71.fna.gz

使用find命令输出结果文件为以Pathways.tbl或Reactions.tbl结尾，使用find-transport命令输出结果文件为以Transporter.tbl结尾，其表头与含义可查看官方文档获取详细的说明：https://gapseq.readthedocs.io/en/latest/usage/output.html

👉 基于这些表格，我们可以迅速比对不同菌株之间的代谢特征的异同

2.代谢网络草图重建和间隙填充

myb71.fna.gz文件同样是gapseq安装目录下提供的测试文件

# 创建草图模型
gapseq draft -r toy/myb71-all-Reactions.tbl -t toy/myb71-Transporter.tbl -p toy/myb71-all-Pathways.tbl -c toy/myb71.fna.gz

# 间隙填充
gapseq fill -m toy/myb71-draft.RDS -c toy/myb71-rxnWeights.RDS -g toy/myb71-rxnXgenes.RDS -n dat/media/TSBmed.csv

draft命令中，则会输出以rxnWeights.RDS、rxnXgenes.RDS、draft.RDS和draft.xml为尾缀的文件，其中draft.xml是只有在sybilSBML库安装时才会创建。

fill命令中，则会输出以RDS和xml（同样，只有在sybilSBML库安装时才会创建）为尾缀的文件，这两个文件是进行下游分析的重要文件。

3.全自动化分析

gapseq doall toy/myb71.fna.gz # 单样本运行完成需要4小时左右

👉 使用doall命令，可以一键自动化生成所有文件

4.反应追踪

在全自动化分析中，gapseq会记录每个反应、为什么以及在哪个阶段将其添加到网络中，这可以有效的帮助用户了解其输出网络以及手动矫正模型的工作。我们可以通过以下简单的代码来获取反应的跟踪情况：

library(sybil)

mod <- readRDS("model.RDS") # 需要调试的模型

# 获取反应的属性表（这个表很大）
mod@react_attr

# 仅检索反应 ID、名称、ec编号和反应的来源
mod@react_attr[,c("seed", "name", "ec", "gs.origin")]

# 将结果输出成一个文件
write.csv(mod@react_attr[,c("seed", "name", "ec", "gs.origin")],
          file = "model_reaction_tracing.csv")

其中，gs.origin包含了较多的反应信息，具体可查看官方文档

---

对单菌/多菌在特定条件下进行时空模拟

BacArena的分析流程

BacArena旨在模拟细菌在群落中的行为；下图是我自制的一个BacArena基本分析的简易流程图，简单来说，我们在使用BacArena进行时空模拟时，需要将候选菌株的模型放置在一个设定好的环境中，通过给予不同的底物进行模拟培养，观察菌株在不同底物下的生长情况，多菌株之间的喂养关系或者底物生产情况等等。

BacArena的用法示例

1.模拟单菌生长

# 使用sybil包自带的Escherichia coli（大肠杆菌）核心代谢模型进行简单模拟
# 使用大肠杆菌核心代谢模型进行简单模拟
library("BacArena")

data("Ec_core")

# 将代谢模型数据转化为Bac类的对象
# （参数 limit_growth=FALSE 设置为禁用生长限制。当空间有限时，这一点很重要，因为细菌可以停止生长或积累大量生物量。）现在我们必须创建一个有机体可以相互作用的环境
bac <- Bac(Ec_core, limit_growth=FALSE)

# 这里，默认情况下我们选择20乘以20作为环境的网格大小。如果未指定，则会自动设置环境的物理区域，在本例中为 0.005cm 乘以 0.005cm。详细信息可以通过以下方式检查
arena <- Arena(n=20, m=20)

arena

# 接下来，我们想通过以下方式将我们创造的有机体置于其环境中
arena <- addOrg(arena,bac,amount=5)
# 通过这条命令，我们添加了 5 个细菌个体。现在我们可以将这些物质添加到环境中

# 在这里，我们添加了含有 0.5 mM 葡萄糖的最小培养基。通过再次调用arena对象，可以检查哪些物质的添加量是多少
arena <- addSubs(arena, smax=0.5, mediac="EX_glc(e)", unit="mM")
arena <- addSubs(arena, smax=1, mediac=c("EX_pi(e)", "EX_h2o(e)", "EX_o2(e)", "EX_nh4(e)"), unit="mM")

# 最后，我们可以开始模拟 12 小时大肠杆菌生长的计算机实验
eval <- simEnv(arena,time=12)

# 在模拟过程中，将打印一些关于增长的基本统计数据。对象评估存储所有 12 个模拟步骤。
plotSpecActivity(eval)
# 检索 eval 对象后，我们可以开始分析数据。为了简单了解发生了什幺，我们可以看看具有高变异性的物质
getVarSubs(eval)

# 选择一种物质并检查其随时间的变化
getSubHist(eval, "EX_glc(e)")

# 从图形分析开始，让我们结合物质变化来研究生长曲线
pdf('curve.pdf')
par(mfrow=c(1,2))
print(plotCurves2(eval, legendpos = "right"))
dev.off()

# 空间和时间变化可以通过以下方式显示
pdf('spatial_time_curve.pdf')
# 由于是12个小时生长，我们以2个小时为间隔，总共就得到6张生长分布图，因此需要将我们的画布设置的可以容纳6张图片
par(mfrow=c(3,3))
evalArena(eval, show_legend = FALSE, time=seq(1,12,2))
dev.off()

结果参考：

2.模拟突变体生长

library("BacArena")
data("Ec_core")

# 现在我们想在环境中多种生物或生物类型。
# 为此，我们创建了两种不同类型的大肠杆菌核心代谢模型：野生型大肠杆菌和无法使用有氧呼吸的营养不良突变体。

bac1 <- Bac(Ec_core,type="ecoli_wt")

# 现在，我们使用 sybil 包的基本命令创建营养不良突变体。
ecore_aux <- changeBounds(Ec_core, "EX_o2(e)",lb=0)
bac2 <- Bac(ecore_aux,type="ecoli_aux", setExInf=FALSE)

# 我们再次设置一个环境并插入生物体和物质
arena <- Arena(n=20, m=20)
arena <- addOrg(arena,bac1,amount=5)
arena <- addOrg(arena,bac2,amount=5)
arena <- addSubs(arena, smax=0.5, mediac="EX_glc(e)", unit="mM")
arena <- addSubs(arena, smax=1, mediac=c("EX_pi(e)", "EX_h2o(e)", "EX_o2(e)", "EX_nh4(e)"), unit="mM")
eval <- simEnv(arena,time=10)

# 在这里，我们将我们创建的两种生物类型（由它们的 x 位置给出）放在环境中，然后开始模拟 10 个时间步长。
# 接下来，我们再次执行所有评估步骤

pdf('multiple_growth_curves.pdf')
par(mfrow=c(1,2))
plotCurves2(eval)
dev.off()

# 以及以不同方式绘制的含有葡萄糖、氧气和乙醇物质的社区空间模式：
pdf('multiple_spatial_time_curves.pdf')
# par(mfrow=c(1,4))
# evalArena(eval,c("Population","EX_glc(e)","EX_o2(e)","EX_etoh(e)"), time=seq(1,10,2))
par(mfrow=c(2,2))
evalArena(eval, show_legend = FALSE,time=seq(1,10,2))
plotSubDist2(eval, sub = c("EX_etoh(e)"), times = c(1,5,10))
dev.off()

# 在这里，第一张图中的不同点颜色表示两种不同的生物类型。
# 最后，我们可以看看表型。该命令将创建具有不同表型相似性的 PCA 图。
# 如果我们对表型的定义感兴趣，我们也可以检索原始表型矩阵
minePheno(eval)
pmat <- getPhenoMat(eval)
pmat[,which(colSums(pmat)>0)]

# 基于这些结果，我们可以看到，营养不良生物类型通常生长较慢，仅使用葡萄糖的发酵，而野生型可以在氧气的帮助下呼吸葡萄糖。

3.模拟群落生长

在例二中，实际上就是一个群落模拟生长的标准操作方法

此外，在BacArena中，作者也提供了一个模拟完成的分析结果，包含8株菌

data("sihumi_test")

# 使用print函数或summary函数查看数据集的基本信息
print(sihumi_test)

# 使用plotAbundance函数绘制不同物种随时间变化的丰度
plotAbundance(sihumi_test)

# 使用plotSpecActivity函数分析在模拟期间哪些代谢物质的吸收和产生模式变化最大
plotSpecActivity(sihumi_test)[[2]]

# 使用findFeeding3函数自动识别和可视化物种之间的交叉喂养关系
feeding <- findFeeding3(sihumi_test, time = 5, mets = c("代谢物质1", "代谢物质2")) # 由你自己设定
plot(feeding)

# 使用plotShadowCost函数分析饮食中缺失但可能对生长有贡献的成分
plotShadowCost(sihumi_test, spec_nr = 7)

# 使用plotSubUsage函数分析不同物种对饮食成分的消费情况
plotSubUsage(sihumi_test, subs = c("EX_sucr(e)", "EX_cellb(e)", "EX_ocdca(e)"))

4.其他功能与补充说明

• BacArena支持并行计算，可以使用simEnv_par函数进行并行计算，或者结合parLapply进行多个模拟以验证结果

    library(parallel)
  replicates <- 2  # 假设我们要进行2次重复
  cores <- ifelse(detectCores() >= 2, 2, 1)  # 根据可用核心数设置核心数
  cl <- makeCluster(cores, type="PSOCK")
  clusterExport(cl, "Ec_core")
  simlist <- parLapply(cl, 1:replicates, function(i) {
  # ... 重复模拟设置和执行步骤 ...
  })
  

• 如果不主动设置底物，那么模拟环境中将不会有任何底物

• 由于不同代谢重建流程参考的数据库差异，代谢产物的名称并不是标准化的，如果使用”错误“的代谢产物名称，则会导致菌株在模拟环境下无法摄取底物。因此如果想尝试手动设置培养基底物文件，一个简单的方法是查看模型中自带培养底物物质的名称： 例如执行simEnv()模拟后，模拟结果存储在bac对象中，那么我们可以通过bac@medium的方式来呼出底物名称，从而帮助编辑底物文件。 至于本案例中涉及到的产物，主要是以BIGG的命名规则，具体可参见下图

• BacArena默认配置的求解器是开源的glpkAPI，该求解器在小规模运行下表现良好，但在面对多菌株互作、大面积模拟环境下性能交叉。可以通过更换求解器来获得更好的性能预测：

    # 举例来说，在执行模拟前设置求解器：
  sybil::SYBIL_SETTINGS("SOLVER","cplexAPI")
  # 如果在配置R求解器有问题请参考后文的软件部署环境
  

测试抑菌效果

使用gapseq和BacArena进行交叉喂养

实际上，在gapseq的readdoc中，已经提供了一些基于bacarena的实际案例分享，包括胃肠道细菌的交叉喂养和酸奶生产中的乳酸菌。这两种案例实际上已经覆盖了较多的菌株代谢分析的场景，而作者们也试图在不断的扩展代谢模型的应用场景，例如gapseq最近正在尝试泛基因组范畴的代谢模型重建。

最近偶然想到，是否可以尝试使用BacArena进行抑菌实验的模拟，这个问题实际上早在19年就已经有人提出。实际上，基于代谢模型来判断抑菌效果已经有相关文献论证了其理论可行性，不过可惜的是，我并没看到这一主题有任何新的进展。

最近，我有幸对手头的一些进行过实验的菌株进行模拟测试。

用于模拟抑菌实验

首先声明，这是一个不严谨的比对实验。原因有2点，一是受限于实验菌株的代谢模型较多，无法一一调教到精准水平下。二是考虑到抑菌本身是一个复杂的过程，此处实验仅从代谢生长角度判断受试菌与致病菌之间的抑制效果，即定义受试菌最终的物种数量如果是致病菌的1倍以上，则认为能从代谢的角度上产生抑菌效果。

1.手头菌株的抑菌情况

目前，我手头有总共85株菌进行了抑菌实验，其中选择的抑菌对象为金黄色葡萄球菌，可以看到在受试菌株中，有27株显示出了一定的抑菌效果。

2.抑菌预测效果

通过宏基因组组装的方式，我们成功获取了其中部分菌株的MAG，并将其使用前文的方式进行建模和比对。结果可以看到，有49株菌的预测结果与实验结果一致，其余24株表现出不一致的情况。并且我也注意到，在不一致结果中，只存在假阴性结果，即实际判定为抑菌，但代谢结果表示为不抑制。

从整体预测角度来看，最终的预测准确率有67%，精确度达到1，受限于样本量，召回率只有4%。

👉 从这个简单的实验中，我们可以看出代谢模型在预测抑菌方向上潜在的应用潜力。

软件构建与使用

大家如果对软件的构建与使用有什么需求，请在同名wx公众号上获取相关信息。

本文由博客一文多发平台 OpenWrite 发布！