代谢网络重构与通量平衡分析,这两个词放在一起,乍一听像是实验室里才有的专业黑话,但它们解决的其实是一个很朴素的问题:给你一条基因组序列,你能不能算出这个生物在特定条件下怎么生长、能产什么、想让它多产某个东西应该改哪里。
基因组规模代谢模型(Genome-scale Metabolic Model,GEM)就是这份答卷。它把所有代谢基因、对应的酶催化反应、涉及的代谢物全部串成一张完整的网络地图,再用数学约束描述这张网络里的物质流动。有了这张地图,过去靠"试错"完成的很多实验,现在可以先在电脑上跑一遍模拟,说句大实话,能省下大量摇瓶子和看发酵罐的时间。
这篇文章写给三类人:刚接触代谢建模的研究生、想用模型辅助菌株改造的发酵工程师、以及想彻底搞懂FBA到底在算什么的分析爱好者。我会从代谢网络重构的基础逻辑讲起,把构建流程、通量平衡分析的数学原理、程序实现、应用案例和常见坑全部串起来,每一步尽量给出能直接上手的操作细节,不绕弯子、不堆术语。
1. 项目核心思路:为什么需要基因组规模代谢模型
1.1 代谢网络重构到底在重构什么
细胞里同时发生着上千个化学反应,这些反应不是孤立的,而是串成路径、编成网络。比如葡萄糖进入细胞之后,先被磷酸化,再裂解成两分子三碳化合物,然后氧化脱羧进入三羧酸循环,产出还原力和ATP,同时各种前体物又被分流去合成氨基酸、核酸和脂类。这个网络中的每一个基因,通常对应一个或几个催化反应;每一个反应,又连接着若干代谢物。
传统代谢工程的做法,是在纸上手画几条主路径,然后试错性地敲基因、加酶、调启动子。但细胞内的事情远比一张A4纸复杂:你敲掉一个基因,底物可能流向支路;你增强一个酶的表达,产物积累的反馈抑制又会压制上游反应;换一种碳源,整个网络走的路径完全不一样。这就是为什么很多人在做菌株改造时,总感觉"改一个基因就冒出一个新问题"。
基因组规模代谢模型要做的事,就是把这一千多个反应和代谢物放进一个计算机可读、可算的数学框架里。它不是画示意图,而是在严格定义:网络中有哪些节点(代谢物)、哪些管道(反应)、每条管道的大小限制(反应速率上下限),以及物质在管道之间流动必须遵循的守恒定律(质量平衡)。
这种重构的价值在于"全局视角"。以前你看到的是局部路径,现在你能看到整个网络在约束条件下的所有可行流动方式。某个基因敲除之后产量是升是降,不用靠猜,先在模型里算一算,方向基本就有数了。
1.2 通量平衡分析解决的三个核心问题
模型建好之后,一个核心问题马上摆到面前:怎么算?
网络里有数百甚至上千个未知通量,而细胞只在特定条件下进入稳定生长。通量平衡分析(Flux Balance Analysis,FBA)的核心假设是:细胞处于稳态,即细胞内所有代谢物的浓度不随时间变化,此时每个代谢物的生成速率等于消耗速率,这可以用一个巨大的线性方程组来约束所有反应通量。
光有约束还不够解出唯一值,因为未知数比方程数多得多。所以FBA又引入一个目标函数,比如最大化生物量合成速率,也就是让细胞长得最快,然后在这个目标下求最优解。这样一来,FBA就能回答三个核心问题:
第一,特定底物下细胞的生长速率上限是多少?这对应实验中的最大比生长速率,可以直接拿来验证模型预测跟实测是否一致。
第二,每一个反应的通量是多少?通量较高的反应往往是产能或供前体的主路径,通量几乎为零的反应则是冗余旁路。这些信息可以直接指导基因操作方向。
第三,某个基因被敲掉、某个反应被关闭之后,最优生长和产物合成会怎么变?这种模拟叫基因敲除预测,是FBA在代谢工程里最常用的场景。
这三个问题,本质上把"生物实验"变成了"数学优化问题",也把代谢工程从盲目试错推向理性设计模式。理解了这一点,你就能明白为什么代谢网络重构在过去二十年里能成为系统生物学最活跃的领域之一。
2. 模型构建全流程实操
2.1 第一步:基因组注释与代谢基因鉴定
构建模型的原料,首先是基因组序列。无论你研究的是大肠杆菌、酿酒酵母、谷氨酸棒杆菌还是乳酸菌,第一步都是拿到一份高质量基因组注释结果,把每条基因的生物功能搞清楚。
高质量注释是模型质量的绝对地基。我见过太多人拿一份粗糙的注释文件直接开始建模型,结果后面花了几周时间在排查一个低质量注释引入的"幽灵反应"。我推荐的操作流程是:先用Prokka或RAST这类自动化工具跑一遍快速功能注释,再用EggNOG-mapper对每个基因做精细的直系同源簇注释,把EC号(酶学委员会编号)和KO号(京都基因与基因组百科全书编号)标出来。
光有EC号还不够。建模型时你需要知道的不是"这个基因是某种酶",而是"这个酶催化哪个反应,反应物和产物是什么,需要哪些辅因子,反应方向如何"。EC号到反应方程的映射,通常用BIGG、BiCG或MetaCyc这类数据库来完成。比如EC 2.7.1.1(己糖激酶),在BIGG里对应的是HEX1反应:ATP + D-Glucose → ADP + D-Glucose 6-phosphate,而且反应方程已经帮你做好了原子平衡和电荷平衡。
这一步实操时最容易踩的坑,是EC号版本差异和底物特异性注释过粗。同一个EC号下面可能包含多个同工酶,底物谱完全不同,必须结合序列同源性和文献做手工修正。基因组注释不是纯粹的全自动流程,它是自动化与人工打磨反复迭代的过程,不要奢望一步到位。
2.2 第二步:从生化反应到代谢网络
拿到基因-反应映射之后,下一步是组装代谢网络。组装过程中必须解决四个关键问题。
第一个是反应的唯一性和方向性。一个反应在你的物种里到底存不存在?它在生理条件下是可逆还是不可逆?这些信息必须从文献、数据库和热力学估算中综合判断。像鞭毛组装、蛋白质分泌这类不在代谢范围内的过程,要果断过滤掉。
第二个是辅因子和能量代谢的守恒。NADH/NADPH、ATP/ADP、乙酰辅酶A/辅酶A这些辅因子的生成与消耗,必须和主代谢路径严格平衡。很多自动构建工具会在这一步漏掉部分辅因子循环,导致模型出现"凭空产能"这种不合理现象,这种模型跑出来的结果没有一点点参考价值。
第三个是隔室问题。原核生物相对简单,大多只有一个胞质隔室;但真核生物有细胞质、线粒体、过氧化物酶体、内质网等多个隔室,同一个反应在不同隔室中对应不同的代谢物节点,而同一代谢物在不同隔室之间需要专门的转运反应来"搬运"。真核模型的隔室标注一旦出错,后面的通量分析很可能全盘皆错。
第四个是反应方程的完整性和原子守恒。每个反应方程里的C、H、O、N、P、S原子数必须左右平衡,如果某个反应少了一个水分子,或者多了一个质子,最后模型的质量平衡就会出问题。现在很多自动化工具比如CarveMe、GapFill和RAVEN Toolbox都能产出初步模型,但原子平衡的精细修补依然需要人眼来核对。
我的经验是:自动化工具产出的初版模型只能当"粗坯",后续必须经过至少两轮人工打磨。第一轮用脚本把所有反应方程过一遍质量平衡,把明显错误抓出来;第二轮逐个检查"死胡同代谢物"——即只能生成不能被消耗的代谢物——这些代谢物通常意味着网络中还缺反应,是模型补全的重点。
2.3 第三步:模型格式与基本属性校验
最终模型应该导出成标准格式。目前最通用的格式是SBML,它是系统生物学界的"通用语言",几乎所有建模工具和数据库都支持。另一个值得留意的格式是COBRA的MATLAB结构或COBRApy的JSON格式,这两个更贴近后续计算环境。
模型里有几个基本对象绕不开:代谢物(Metabolite)、反应(Reaction)、基因(Gene)、以及基因-蛋白-反应(GPR)关联规则。GPR规则尤其重要,它刻画的是一个反应被执行需要哪些基因参与:如果是"且"逻辑,说明多个基因共同编码一个酶复合体,缺失任何一个都不行;如果是"或"逻辑,说明多个同工酶都可以催化同一反应,敲掉其中一个还有替补。
模型建好后,第一轮校验一定要看四个指标:所有反应是否满足原子守恒;能否在标准培养基下模拟出非零生长;交换反应的上下限是否跟实际培养基组分一致;敲除一个必需基因后模拟生长是否变为零。这四个指标分别对应质量平衡、网络完整性、环境约束正确性和GPR规则合理性。
我做完这些之后还会多跑一个"回归测试":对每个基因做一遍单基因敲除模拟,计算出预测必需基因集合,然后跟文献里的必需基因列表对比重合率。重合率低于八成时,基本可以断定模型里有些反应方向或者辅因子定义有问题,需要回头修改。这个过程比较枯燥,但它是模型可信度的试金石。
3. 通量平衡分析原理与实现
3.1 FBA的数学本质
通量平衡分析听起来高深,但数学内核一句话就能讲清楚:它是一个带约束的线性规划问题。
先定义变量v_j为第j个反应的通量,单位通常是mmol/gDW/h,即每小时每克细胞干重消耗或生成的毫摩尔数。细胞内每种代谢物i都要满足稳态条件,于是有矩阵方程:
S · v = 0
这里S是化学计量矩阵,S_ij表示代谢物i在反应j中的化学计量系数,生成取正、消耗取负。这个方程的意思很朴素:稳态下每个代谢物的净生成速率等于净消耗速率,所以它不会在细胞内凭空积累或消失。
与此同时,每个反应还受到上下界限约束:
lb_j ≤ v_j ≤ ub_j
培养基里没有的底物,对应摄取反应的上限就应该是0;氧气浓度有限,摄氧量就存在上限;不可逆反应的下限是0。这些上下限给求解加入了"物理现实"的边界。
然后需要设定一个目标函数Z,最常用的是最大化生物量反应通量v_biomass,写成:
Maximize: Z = c^T · v
其中c是一个向量,只有生物量反应那一项为1,其余为0。求解这个线性规划,就能找到在满足所有约束的前提下让目标函数最大化的通量向量v,作为预测的代谢通量分布。
这个框架之所以成立,是因为它抓住了细胞生长最基本的"权衡"本质:碳架、能量、还原力都是有限资源,流向生物量的多了,流向产物的就少,反之亦然。线性规划还有一个很好的理论保证——只要可行域不空,就一定有全局最优解,而且最优解在可行域的顶点取得。这个顶点往往正好对应着通量分布中大量反应为零、少数反应处于饱和状态的极端网络状态,跟实际发酵中"少数主路径承担绝大多数通量"的观察非常吻合。
3.2 目标函数与约束条件的设定技巧
目标函数的选择是整个FBA流程里最主观、也最影响结果的一环。
最基本的选项是最大化生物量合成速率。此时"生物量"不是真实反应,而是人为构造的一个"汇"反应,它按一定化学计量系数一次性消耗所有必需前体物,包括氨基酸、核苷酸、脂质、辅因子、细胞壁成分,合成1克干细胞。这个反应的化学计量系数来自实验测定的细胞组成数据,不同菌种、不同培养条件下差异很大。所以你在文献里看到的生物量反应,一定是从特定菌株生长数据来的,不能随手套用。
如果你关心的是目标产物产量,可以把目标函数改成最大化产物净输出通量,或者最大化产物通量与生物量通量的加权组合。更精细的做法是双目标优化:先最大化生长,再在生长速率不低于某个阈值(比如最大值的95%)的约束下最大化产物通量。这更贴近真实发酵过程"先长菌、后产东西"的两阶段模式。
约束条件的设定同样需要分步走。最简单的约束只有几个交换反应上下限和不可逆反应下限;但随着分析深入,往往需要叠加更多约束:把实测的葡萄糖摄取速率、氧气摄取速率、CO2排出速率当作已知数据收紧相应交换反应边界,让模型预测更贴近实际发酵曲线;加入转录组数据作为转录水平约束;加入非生长相关维持能(NGAM)这一项固定ATP消耗;以及加入热力学可行性约束,排除方向违背热力学规律的反应。
NGAM这一点经常被新手忽略。细胞即使不生长,也在做离子平衡维护、蛋白周转,这部分ATP消耗不是零。如果模型里没有这一项,FBA往往预测出"所有碳源全部用来生长"的过于乐观结果。常规做法是设置一个ATP水解反应ATPM,通量下界取实测经验值。比如大肠杆菌经典模型里NGAM一般取3到8 mmol/gDW/h,但这个值不同文献差异挺大,最好用自己批次培养的底物消耗数据来校准。
3.3 程序实现:用Python完成一次FBA
理解了原理,写代码就很简单了。现在的标准工具链是Python生态下的COBRApy。这里给一个最小可行的示例,以公开获取的大肠杆菌GEM模型为例:
import cobra from cobra.io import read_sbml_model # 加载SBML格式的模型 model = read_sbml_model("ecoli_gem.xml") # 设定培养基:葡萄糖为碳源,氧气充足,无机盐不限制 model.medium = { "EX_glc__D_e": 10.0, # 葡萄糖摄取上限 10 mmol/gDW/h "EX_o2_e": 18.5, # 氧气摄取上限 18.5 mmol/gDW/h "EX_h2o_e": 1000.0, # 水不限制 "EX_pi_e": 1000.0, # 无机磷不限制 "EX_nh4_e": 1000.0, # 氨态氮不限制 "EX_fe3_e": 1000.0, # 铁离子不限制 } # 设置目标函数为生物量反应 model.objective = "R_biomass" # 求解 solution = model.optimize() print("生长速率:", solution.objective_value) print("葡萄糖摄取通量:", model.reactions.EX_glc__D_e.flux) print("乙酸分泌通量:", model.reactions.EX_ac_e.flux)这段代码跑完之后,输出会告诉你在给定培养基下这个菌株的最优生长速率、葡萄糖摄取速率和副产物分泌通量。拿实验实测值做对比,如果差异不大,说明模型骨架合理;如果差异很大,优先怀疑生物量反应系数和NGAM的取值。
基因敲除模拟同样简单:
# 模拟敲除某个基因:该基因关联的所有反应通量强制设为0 with model: model.genes['b0720'].knock_out() ko_solution = model.optimize() print("敲除后的生长速率:", ko_solution.objective_value)这个with上下文管理器非常方便,它保证所有修改只在代码块内生效,退出后自动恢复原始模型,不用手动备份。
如果要做全基因组范围的敲除扫描,比如找出所有能提升产物得率的单基因敲除靶点,写个循环逐一把非必需基因敲一遍、记录产物通量变化、排序取前几十名就行。实操时建议先做一次"基因必需性预筛",把敲除后生长直接降为零的必需基因跳过,能省下不少计算时间。
3.4 结果解读与通量可视化
算出一堆通量数字之后,不能直接把原始数值扔给实验人员看。对FBA结果,我觉得至少有三件事值得做。
第一件事是通量分布的整体检查。最优解里哪些反应通量特别高、哪些是零,那些高通量反应尤其要注意跟实验上已知的代谢流分析结果做对比。比如模型预测三羧酸循环完全关闭,但实验数据表明该菌有氧培养时TCA循环活跃,那大概率是模型约束出问题了。
第二件事是优化解的鲁棒性检查。一次FBA解只是一个可行顶点,这个顶点对参数变化是否敏感,目标函数稍微扰动会不会导致通量分布剧烈变化,这些需要做参数扰动或者通量变异性分析(FVA)来确认。
第三件事是可视化。把通量画到代谢网络图上是最直观的解读方式。COBRApy自带的绘图接口功能比较基础,我更多是把通量结果导出成表格,在Cytoscape或者画图脚本里叠加到代谢路径图上,用线条粗细表示通量大小,用颜色表示反应方向变化。
FVA这里的价值值得多说两句。FBA只给你一个最优通量值,但很多反应在最优解上的自由度其实很大,只要目标函数不降太多,它们的通量可以在一个区间内浮动。FVA做的事情就是:在约束目标函数不低于最优值一定比例(比如99%)的前提下,分别求解每个反应通量的最小值和最大值。这样你得到的是一段"通量可行区间"而不是单纯一个点。区间越窄,说明该反应在最优方案里的角色越被锁定;区间很宽,说明它是可替代的柔性环节。做基因改造靶点筛选时,大多数人当然更倾向于优先动那些区间窄、角色明确、影响大的反应。
4. 模型应用:从预测到改造
4.1 基因敲除靶点预测:从单敲到组合设计
GEM模型在代谢工程里最经典的应用,就是给菌株设计"在哪里动刀"。
最简单的方案是单基因敲除扫描:每个基因敲一遍,模拟目标产物通量的变化,按效果排序。但实操中很快会发现问题:很多基因单敲之后产物通量变化并不显著。这时就需要换一种视角,用"生长偶联生产"来判断。敲除某个基因之后,如果目标产物的通量与生物量反应形成正相关——细胞每多长一点就必须多产一点目标产物——那这个突变体就是极有实验价值的候选菌株,因为它把产物合成和生存需求绑在了一起,不容易在传代中退化。
组合敲除的设计又是另一番考量。单敲效果不好时,往往需要两三个基因协同操作。这时候先别急着做组合穷举,应该先用最小化代谢调整算法(MOMA)或调控优化算法(ROOM)做一轮预筛,把候选基因集合缩小到几十个,再在这个范围内做两两组合模拟敲除。组合搜索的计算量是平方级别增长的,几百个基因做全组合就是上万次FBA求解,虽然单机也能跑,但没有必要每次都跑全量。
我印象很深的一个案例是:某个课题组想提高一株大肠杆菌的氨基酸产量,模型扫描发现单敲某个转运蛋白基因后,产物通量反而轻微下降;但模拟结果表明,把这个转运蛋白基因和另一个副产物合成路径的关键基因同时敲掉,产物通量能提升接近两倍。这种非直觉的组合,恰恰是传统试错法很难发现的。
4.2 培养基与底物优化
模型的另一个高频应用方向,是培养基组分和碳源选择优化。
在菌株固定的前提下,培养基里用哪种碳源产物得率更高?要不要补加某种氨基酸或维生素?某个转运环节是不是瓶颈?这些问题都能在模型里用"通量扫描"或"培养基组分优化"给出初步答案。
具体做法也不复杂:在一个包含多种候选碳源、氮源、无机盐摄取反应的模型基础上,把目标函数设为产物分泌通量最大化,然后对每一种碳源与氮源组合跑一次FBA,比较产物通量。对筛选出来的最优培养基,再做一轮"营养限制分析":把每种组分的摄取上限逐一设为零,看目标产物通量降多少,降得最多的组分就是潜在的瓶颈营养物质。
这一步需要清醒认知:模型预测的是"这个基因型在最优调控策略下的理论表现上限",真实菌株的转录调控和翻译调控未必能完全达到模型设想的最优策略。所以培养基优化的结果更适合作为实验设计方向,而不是最终答案。我习惯把模型排出来的前几个组合方案配合小规模发酵矩阵做验证,经常有意外收获,但也遇到过模型说行、实验完全不长的情况,那通常说明建模时忽略了某些严格的调控机制或底物毒性效应。
4.3 与组学数据整合的扩展玩法
把转录组、蛋白组、代谢组数据整合进GEM模型,是当前系统生物学非常活跃的扩展方向。
最有代表性的一类方法是"基于数据构建约束",典型的如GIMME算法:把转录组数据映射到基因-蛋白-反应关联规则上,低表达的基因对应反应会被施加额外的通量上限约束,然后模型在"通量分布尽量与转录组低表达情况一致"的约束条件下重新求最优解。这样得到的通量分布比纯FBA更贴近真实细胞状态,尤其适合分析处在复杂调控条件下的细胞。
另一类方法是"代谢物覆盖度分析"。把代谢组检测到的代谢物集合映射到模型上,标记哪些代谢物在网络里是死胡同节点——能合成但不能被代谢。死胡同代谢物往往意味着某条下游路径缺失,而这可能恰好跟目标产物积累或者副产物毒性相关。这种整合分析能帮你把代谢组里的"异常峰"定位到具体的代谢路径上。
还有一个容易上手又很实用的做法:把通量变异性分析结果和转录组数据叠加比对。如果一个反应在模型里可行区间很窄、通量很大,但其编码基因在转录组里表达量却很低,这就产生了"模型预测与实验数据冲突"的信号。冲突点要么说明模型约束错了,要么说明体内实际情况存在调控抑制。这类冲突点集合起来,往往是论文里最有价值、也最能推动下一步实验方向的内容。
5. 实操中的常见问题与排查技巧
5.1 模型不可行:长不出菌怎么办
建好模型跑FBA,第一个经常遇到的结果就是求解器报"不可行",也就是不管怎么分配通量都满足不了所有约束,模型里压根不存在可行解。
遇到这种情况,第一反应不是怀疑求解器,而是系统性排查约束条件。最常见的原因有四个。
第一个是交换反应上下限设错。比如你本意是让模型只用葡萄糖,却忘了把其他碳源(乙酸、甘油、乳糖等)的摄取上限设成0,模型可能走捷径靠别的碳源生长。但反过来,如果你把某个生物量合成必需的无机盐摄取上限误设为0,模型就会直接不可行。
第二个是死胡同代谢物累积。某个代谢物只能被反应生成、却没有反应消耗它,稳态方程"生成等于消耗"就没法满足。此时要么用模型补全算法(GapFill)自动添加缺失反应,要么人工检查死胡同代谢物,定位缺失的下游路径。
第三个是无效循环导致的数值问题。某些反应组合能形成一个"凭空循环",不消耗净底物就能持续生成代谢物,这让求解过程中的数值性质变得病态。遇到时可以运行CobraPy的反应循环检查工具,把这部分反应的下限做微调。
第四个是隔室方向定义错误。真核模型的跨膜转运方向一旦搞反,某个隔室里的代谢物无法被消耗、不断累积,稳态方程同样会产生矛盾。
排查顺序建议是:先查交换反应,再查生物量反应的组分完整性——你是否把模型需要的关键代谢物都放进了生物量方程——最后才用GapFill自动补全。不要一上来就开GapFill,那样补出来的反应往往纯粹是为了让模型可行而存在,生物学意义缺失,会给后续所有分析埋雷。
5.2 通量分布符合最优但不符合常识
另一种高频问题是:FBA能顺利求解,生长速率也不低,但通量分布看起来跟常识对不上。比如明明是有氧培养,模型预测TCA循环完全关闭;或者某些反应方向跟预想完全相反。
这种情况大多数是"多目标性"造成的。真实细胞的代谢调控并不只追求最大化生长一个目标,还有应激响应、氧化还原平衡、蛋白表达成本等众多考量。FBA把这一切简化成线性约束加一个线性目标函数,给出的自然是一个"偏科"的极端解。
化解办法分三个层次。第一层是调整目标函数,比如改用双目标优化,漫帮产品通量和生长通量的加权组合;或者采用"最少总通量"策略,即在满足最优生长的前提下让全网络总通量最小化,这样得到的分布更接近实验测量值。第二层是增加额外生理约束,比如设定呼吸速率下限、维持NADPH供需平衡。第三层是放弃静态稳态假设,改用动态通量平衡分析(dFBA),把时间维度加进模型。
还有一个容易被忽略的问题:代谢环路。模型里如果存在一个能够净产生ATP但不消耗任何底物的无效循环,求解器有可能在某些约束组合下利用这个循环"空转产生能量"。真实细胞里这类无效循环会被严格调控抑制。排查方式是对每个循环计算净ATP产率和碳流走向,必要时在循环中加入热力学约束,把不合理方向堵死。
5.3 交换反应通量异常与培养基模拟偏差
交换反应是整个模型里最容易出错的部分。它代表外部环境与细胞内部物质交换的边界,比如葡萄糖摄取对应的EX_glc__D_e、乙酸分泌对应的EX_ac_e。因为培养基组分几乎全由这些反应的上下限定义,一旦出错,所有下游分析都会受影响。
常见的问题包括:培养基里加了某种物质,却忘了把对应摄取上限设为正值;单位搞混,把mmol/L当成mmol/gDW/h用,直接差出数量级;或者把同一个物质既设为限制性底物又设为非限制性底物,约束自相矛盾。
给我的建议是:开始任何FBA分析之前,先把所有交换反应上下限列成一张表,与实际发酵培养基配方逐项比对。我做项目的第一步永远是跑一遍"模型现状打印",把交换反应全部输出成表格,跟实验记录的培养基成分一字一句核对。这个步骤很枯燥,但能阻止后续大量返工。
另一个高频事故是"敏感性分析"时上下限没清理干净。比如想观察葡萄糖摄取速率从1到10 mmol/gDW/h对生长速率的影响,做法是直接改EX_glc__D_e上限并重复求解,但一定同时要把这个反应的下限设为0。否则模型可能给你算出一个"自动反向分泌葡萄糖"的荒谬结局。很多新手拿到一个怪异的敏感性曲线,第一反应是怀疑求解器,其实九成情况只是下限没设干净。
5.4 模型质控清单
最后整理一份我个人实战中每次拿到新模型都会逐项过一遍的质控清单:
| 检查项 | 方法 | 参考标准 |
|---|---|---|
| 原子质量平衡 | 检查每个反应方程原子数平衡 | 无未平衡反应 |
| 电荷平衡 | 检查反应方程两侧净电荷差 | 无显著失衡 |
| 基础可行性 | 富培养基下FBA求解 | 有非零生长解 |
| 生物量反应完整性 | 核对生物量方程包含的代谢物是否都在网络内 | 无缺失节点 |
| 必需基因验证 | 单基因敲除模拟并对比文献 | 重合率大于80% |
| 交换反应设置 | 输出交换反应与培养基配方逐项比对 | 无多余或缺失 |
| 死胡同代谢物 | 脚本列出只能生成或只能消耗的代谢物 | 数量越少越好 |
| 通量分布合理性 | 对比中心碳代谢通量已知范围 | 数量级一致 |
这份清单正常速度跑完需要半天,遇到问题多可能要两三天。花这个时间绝对值得,因为模型产出的所有预测结论,其可信度都建立在模型自身自洽性之上。地基不稳的话,后面算出再漂亮的结果也没有意义。
越做模型越发现,基因组规模代谢模型不是"做完就完事"的静态产物,而是一个需要跟实验数据持续互相校准的动态工具。我自己的习惯是:每拿到一批新的发酵数据,就把实测的底物摄取速率、产物分泌速率、生长速率塞回模型里,重新校准生物量方程和ATP维持项,看看哪些旧参数需要修正。这样经过三五个迭代,模型对菌株的预测能力会有肉眼可见的提升。
还有一个小心得:FBA算出来的通量分布永远不是"真相",它更接近"这个基因型在给定约束下的最优可能性空间"。真正能说服同事和审稿人的,永远是模型预测与实验验证之间的交叠部分,以及那些存在于冲突缝隙里的新假说。所以做模型的人一定不要闭门造车,多跑发酵罐,多跟实验端的伙伴对齐细节,哪怕只是确认一句他们用的培养基成分表,"翻车"概率都能显著下降。
如果这篇分享能帮你把第一个GEM模型建起来、把第一轮FBA跑通,那就算达到目的了。后续想深入的话,可以往动态建模、多菌共培养模型、或者酶约束模型(比如GECKO这类框架)方向探索,这些都是在GEM基础上长出来的新分支。路还很长,但起点就是今天讲的这些基础逻辑。