基因组规模代谢网络模型自动化修正
2017-11-01吴晓红石贵阳
吴晓红， 薛 卫， 张 梁*， 石贵阳
（1.粮食发酵工艺与技术国家工程实验室，江南大学，江苏 无锡 214122；2.南京农业大学 信息科学技术学院 江苏 南京 210095）
基因组规模代谢网络模型自动化修正
吴晓红1， 薛 卫2， 张 梁*1， 石贵阳1
（1.粮食发酵工艺与技术国家工程实验室，江南大学，江苏 无锡 214122；2.南京农业大学 信息科学技术学院 江苏 南京 210095）
基于KEGG在线数据库以及6个蛋白质区间预测数据库，对基因组规模代谢网络模型进行了自动化修正。作者提出了蛋白质区间预测结果的权重打分机制，同时利用图像处理算法确定可信度高的特异性反应。上述修正的研究均在Spathaspora passalidarum NRRL Y-27907基因组规模代谢网络精炼过程中得到运用实施，对于提高模型构建效率意义重大。
基因组规模；代谢网络；断点补齐；图像处理；区间预测
随着基因组高通量测序数据的涌现以及大量的生物学数据的产生，代谢网络模型构建成为研究生物信息学的热点之一。代谢网络构建是一个耗时费力的过程，因此许多自动化构建的工具随之应运而生。通常这些自动化工具侧重关注代谢网络粗模型的构建如 metaSHARK[1]和 AUTOGRAPH[2]，其次关注代谢网络模型的模拟过程，如CellNetAnalyzer[3]、OptFlux[4]和 COBRA Toolbox[5]，只有少量的自动化工具是针对代谢网络模型的精炼过程。目前能够提供代谢网络模型自动化精炼过程的工具有Model SEED、Pathway Tools、RAVEN 和 SuBliMinaL。
代谢网络的模型构建包括粗模型的构建、模型的精炼、数学模型的转换、模型的预测验证四个过程。一个高质量的代谢网络模型，应达到模型模拟结果和生物实际生长表型一致，否则要不断的重复精炼修正过程，直到模拟与表型一致。模型的精炼修正无疑是代谢网络模型构建过程中最耗时耗力的过程，现有模型精炼工具并不能真正实现真菌代谢网络模型精炼过程的自动化。模型的精炼过程必须包括漏洞代谢的填补、反应区间定位等。Model SEED[6]和Pathway Tools[7]只能提供原核生物的代谢网络模型的精炼自动化过程，不能提供反应区间的定位。RAVEN[8]和SuBliMinaL[9]是基于Wolf PSORT蛋白质区间预测数据库实现自动化定位区间的程序。但是Wolf PSORT[10]只是基于氨基酸组成特征的在线预测数据库。研究表明，基于氨基酸组成、二肽和物理化学三种综合特征的蛋白质区间定位预测结果更为准确[11]。
利用作者所在实验室自动化构建全基因组代谢网络模型的程序，自动构建了Spathasporapassalidarum NRRL Y-27907全基因组规模代谢的粗模型。以S.passalidarum NRRL Y-27907的基因组规模代谢网络模型的精炼过程为例，以简单、面向对象的Java语言为基础，对精炼过程中人工冗杂的断点补齐的方法进行了研究，提出了一种基于KEGG[12]在线数据库自动化填补漏洞反应的方法，并利用权重打分机制分析，6个真菌蛋白质定位数据库预测S.passalidarum NRRL Y-27907的结果，在保证模型中反应的物种特异性的同时，实现了真菌代谢网络模型精炼的自动化。自动化修正的流程见图1。图中进程g、进程n、进程o为一个小的流程循环。进程g中判断反应包含断点，则进入进程h，查找该反应在注释图谱中对应的坐标，并在进程i中读取此坐标，在进程j中判断此坐标是否为特异性反应，如果是，则在进程p中记录该反应。如果不是，则在进程l中判断此坐标是否为最后一个坐标，如果是最后一个坐标，则进入进程n，即进入进程g、进程n、进程o该流程循环。如果不是最后一个坐标，则进入进程m，读取下一个坐标，判断此坐标是否为特异性反应，重复此循环直至将所有的特异性反应都被找出，进入进程q，进行模型修正。在进程r中判断模型中是否已经包含此反应，若已经包含，则回到进程n，即进入进程g、进程n、进程o该流程循环，检查下一条反应。若不包含此反应，则进入进程s，将此反应加入到模型中。
图1 自动补齐断点流程Fig.1 Process of the auto-refinement of gap
1 自动填补网络漏洞
采用柴文平[13]等人的方法构建了S.passalidarum NRRL Y-27907代谢网络粗模型。构建的代谢网络粗模型需要进一步精细化与修正，最终完成一个高质量的基因组规模代谢网络模型。
1.1 代谢网络漏洞查找
模型导入到装有COBRA工具包和GLPK线性规划器的Matlab中，将模型转化为计算机可读的格式 （SBML）才能进行代谢网络漏洞查找。通过xls2model程序将模型Excel表读取为计量学S矩阵。S矩阵（828×984）表示该模型由828个代谢物和984个反应组成。同时通过GapFind程序完成代谢漏洞的查找，其中上游漏洞代谢物有为44个，下游漏洞代谢物有128个。
1.2 基于KEGG网络爬虫反应
KEGG是代谢网络构建常用数据库，含有多个在线子数据库，其中REACTION数据库包含迄今为止发现的所有生化反应。各个子数据库的网页数据格式比较统一明确，方便人们进行远程服务器访问。但是，KEGG数据库更新频繁，各个子数据库不能够免费下载，需要付费使用。而在基因组代谢网络断点补齐过程中，因为数据信息量浩大，频繁访问远程服务器比较耗时耗力。因此，实现一种批量在线获取并存取数据的方法意义重大。
1.2.1 方法概述 利用超文本转移协议和Java控件HttpClient相结合，实现对网页中特定信息的抓取KEGG提供物种特异性基因组信息以及所有反应式信息查询网页，通过一定的URL（Uniform Resource Locator，统一资源定位符）格式地址发送HTTP请求并获取网页中的基因信息。在漏洞填补的过程中需要访问大量不同的网络资源，获取相关的基因信息，由于数据量较大且人工操作比较繁琐，这里利用Java控件HttpClient实现爬虫技术，抓去符合特定条件的网络资源。HttpClient是Apache Jakarta Common下的子项目，可以用来提供高效的、最新的、功能丰富的支持HTTP协议的客户端编程工具包，并且它支持HTTP协议最新的版本和建议。利用HttpClient访问具体的URL地址，获取服务器端返回的获取html内容，html内容由标题、js代码、正文、相关链接、声明等区域组成，而有用信息只出现在正文中的各种html标签标记内，分析html标签并获取特定的网页信息。
1.2.2 漏洞填补算法实现
1）获取注释图谱：提交物种基因组蛋白质序列至KAAS自动注释服务器，获取注释信息，下载html和text格式。
2）查找包含断点的注释图谱：根据Matlab软件中GapFind程序返回的漏洞代谢物列表，在代谢网络模型Excel格式中确定代谢物的反应途径，依据KASS注释返回的途径图谱找到包含漏洞代谢物的所有反应。
注释返回的KEGG代谢途径为包含糖代谢等在内的110个途径。查找包含断点的代谢图谱的流程见图2。具体思路和伪代码步骤如下：
A：获取断点化合物所对应的Subsystem信息，记为sub。
B：向注释查询网页URL地址发送HTTP请求。
C：如果服务器端响应代码为HTTPStatus.SC_OK则正常响应，否则继续请求，获取html正文内容。
D：分析html内容，设i为行号，由第一行开始遍历标签对中的每一行，
For i from 1 to n
｛
if（该行中第二个标签中的内容与sub相等）
｛
提取对应的第一个标签中的内容，记为KO；
｝
else
忽略该行，遍历下一行；
｝
E：根据D中的KO号得到满足条件图谱的URL地址，向URL地址发送HTTP请求得到服务器端响应的网页图片记为T1，T1即为整个网络结构图，其中绿色酶号表示包含断点的特异性反应。
F：点击T1左上角途径方框，进去包含所有反应页面page1，网页中每一个EC号对应图谱中的一个具体反应，它的URL地址指向具体的反应方程式。
G：获取page1中所有EC号对应的反应，设ec_num为每一个EC号，从第一个开始
EC_K_Break.txt保存包含断点化合物的EC，K号的信息。3）查找EC_K_Break.txt中每个K对应的坐标根据K号获取其在T1中对应的坐标，判断特异性反应。
图2 断点代谢途径定位Fig.2 Orientation of gap metabolic pathway
1.3 判断特异性反应
KEGG所有的反应都包含在通路数据库（PATHWAY database）中，PATHWAY 图谱上有颜色标记的酶号是指这个物种特定的基因或酶，只有有颜色标记的酶号表示的反应才是具有该物种特异性的反应，也才能添加到代谢网络模型中。在代谢网络模型中添加非特异性的反应会改变整个代谢途径和代谢物流量，进而使模型模拟的结果偏离实验数据，影响模型的准确性和可信度。
构建代谢网络模型需要提取代谢途径中的特异性反应，图中特异性反应对应的酶号所在的方形框有颜色标记。因此通过网络爬虫技术获得方形框的位置列表，定位到某酶号所在的方形框后需要选取框内的像素点，读取其颜色值，如果颜色分量RGB均为0或255，则没有颜色标记，反之则有。代谢网络特异性反应获取流程见图3。
基本思路为：
根据得到的position坐标读取T1对应点的RGB色彩值。
Picture（Key：酶号；Value：代谢网络图中所有方形框的坐标向量集｛V1，V2，……，Vn｝）
For i from 1 to n
{
If（某酶号所在的方形框）
｛
沿方形框的长边内侧逐一选取像素点，读取其颜色值；
If颜色分量RGB均为0或255 then没有颜色标记
else有颜色标记；
If有颜色标记then该酶号对应的是特异性反应
do将反应加入菌的代谢网络模型中；
else舍弃该酶号对应的反应。
｝｝
反应式漏洞填补
遍历new_rec.TXT中每一个反应，查看模型中是否存在，存在则不处理，否则添加。
A：读取new_rec.TXT中每行反应记为new_rec，i为行号
For i from 1 to n
｛
if（模型中不包含 new_rec）
｛
将new_rec添加到模型中；
｝
else
忽略该反应，查找下一条反应；｝
图3 特异性反应获取流程Fig.3 Process of getting the pecificreaction
2 获取反应区间定位
细胞是生命活动的基本单位，它由执行不同机体功能的称为亚细胞的各部分组成，如细胞膜、细胞核、线粒体、高尔基体、内质网等。亚细胞功能是由位于其中的蛋白质执行的，蛋白质所在的亚细胞称为蛋白质的亚细胞位置[14]。蛋白质必须转运到其应在的亚细胞位置上才能正确行使其功能，否则就会出现机体功能紊乱，正确合理的蛋白区间定位是高质量模型构建的基础，见表1。
表1 真菌蛋白质亚细胞预测数据库Table 1 Database for subcellular localization of fungal proteins
确定一条蛋白质的亚细胞位置称为蛋白质亚细胞定位[15]。蛋白质亚细胞定位的传统方法是通过生物化学实验，如射线晶体衍射电子显微镜核磁共振等方法进行测定[16]。实验方法精确度高，但费时耗力代价昂贵，而且对难于结晶的蛋白质来说，实验方法不再有效。借助于先进高效的计算机自动化数据处理技术，出现了一些蛋白质定位预测网站。结合Spathasporapassalidarum NRRL Y-27907的生理生化性质和蛋白质特征提取方法、算法和准确性等，选取了6个真菌生物蛋白质区间预测网站，自动化提取分析网站的预测结果，在权重打分机制的基础上得到最佳的蛋白质定位区间。这6个网站是基于蛋白质的氨基酸组成、伪氨基酸组成、二肽、生物化学特征或是四种特征的综合。
2.1 区间定位算法实现
A：对每条反应获取对应的KO号。
B：将A中的KO号在KASS注释结果中查找基因号，并在本地下载Spathasporapassalidarum NRRL Y-27907蛋白质序列库提取其对应的蛋白质序列。
C：将蛋白质序列提交到对应网站的表单中，获取返回的定位信息。
D：获取定位区间的信息并填入反应式中。
在获取具体反应的区间信息过程中，需要将反应所对应的蛋白质序列提交到网页的表单中，提交后返回具体的区间定位信息，此处会遇到两个问题：1）表单提交过程中不支持大量蛋白质序列自动提交。由于模型中蛋白质序列数量较大，在有的网站中获取定位信息时不支持大量序列的一次性提交而只能分别提交单个序列获得定位信息，在提交过程中任务量大且人工耗费时间长。2）大量蛋白质序列提交耗费时间长，在网站中提交多个序列后等待服务器端反馈的定位信息耗费时间太长，甚至会发生无响应等问题，见图4。
图4 反应亚细胞定位流程Fig.4 Process of subcellular localization
HttpClient支持访问特定的URL地址，获取服务器端返回的html信息，并且能够分析html中form表单中的信息，实现内容的自动提交。由于涉及到的定位页面所有的表单提交方式都是POST提交，利用HttpClient中的PostMethod方法实现post提交。表单中的元素赋值过程：获取表单中需要赋值的元素标签，以蛋白质序列元素赋值标签为例，标签为