生物信息学中的算法问题 Version 1.0 卜东波 贺思敏 赵 屹 中科院计算所生物信息学课题组 华大－曙光生物信息学联合实验室 2002/12/03

主要内容 生物信息学中的算法问题 我们的工作 (ICT & IBP & BGI)

一、生物学 vs 信息科学

生物信息学的研究目标 特点：天然的形式化 碱基： A,C,T,G 四种 常见氨基酸： 20 种 目标： 以 DNA 序列作为源头 揭示 “ 基因组信息结构的复杂性及遗传语言 的根本规律 ” ； 之后进行蛋白质结构和功能预测。

生物信息学的两个挑战 高性能计算： 海量的数据 每 14 个月翻一番 算法： 海量的数据使得原有算法不适用 新需求

生物信息学的研究流程 第一步：生物学问题的提出 生物学为主 第二步：数学建模、算法设计 信息科学为主 第三步：结果解释、实验验证 生物学

生物信息学脉络

生物信息学问题概览（ 1 ） 基因组时期：序列－结构－功能 DNA 测序和拼接 比对 进化树 蛋白质质谱鉴定 序列注释：基因预测、细胞定位 结构预测： RNA 结构预测、蛋白质折叠 。。。。。。

生物信息学问题概览（ 2 ） 后基因组时期：相互作用－网络－功能 生物芯片（ DNA 芯片、蛋白质芯片） 相互作用网络 调控网络 E-Cell 药物设计 。。。。。。

1. 大规模测序和拼接 生物学问题： 从 DNA 片段恢复原始序列

DNA 整体 切成 小段 小段和载体结合 结合后进行测序

全自动的测序仪器： MegaBace

需要拼接！ 因为整个基因组太长（上 M), 而每次只能 测得一个 500 的小片断 (read) 问题：如何根据 read 恢复原始顺序？ 类比： 10 本圣经，都从随机点起始剪成 500 个字母左右的小纸条，问：给你这么 一堆小纸条，你能读出圣经来吗？

拼接问题的数学描述 数学问题： 公共超串 输入：设有字符串 S ，预先估计其长度大约为 n ，现在 已知一个字符串集合 R={R 1,R 2 …R n } ，其中每个 R i 都是 S 的一个子串。问：原始序列 S 是什么？ 算法： Hamiltion 路径类 Euler 路径类 Local Search 类

2. 序列比对 生物学问题： 序列的相似性－ > 同源性 原始序列： S: acgctg T: catgt 可行解： 1. S: a c g c t g - T: - c – a t g t 2. S: a c g c t g - T: - c a – t g t 3. S: - a c g c t g T: c a t g - t -

序列联配： 两序列联配： 全局联配 (Global Alignment) 局部联配 (Local Alignment) 空位处罚 (Gap Penalty) 多序列联配 全基因组比对

Open Problems: 快速的多序列比对算法 快速的全基因组比对算法

3. 进化树 生物学问题： 根据形态、 DNA 、行为学特征 推导种群进化关系树

进化树问题的数学描述 输入： N 个物种的特征（ DNA 、形态。。。） 输出： 以这 N 个特征为叶节点的一颗树 距离法： 聚类谱系树 简约法： 最小突变树

4. 结构预测 结构大致决定功能 一级结构 ( 氨基酸序列 ) 二级结构 ( 螺旋、片层、回环 ) 超二级结构（ aba … ) 三级结构 ( 由二级结构组合成三维构像 ) 四级结构：多个亚基 实验测定方法： x-ray 晶体衍射 NMR 核磁共振 实验耗时、昂贵 一个蛋白质结构测定需要 $200K or more 需数月或者更长 有些蛋白质还无法测定

蛋白质结构（ 2 ） 理论上可计算的。 能量最低原则 变元： 主干的 psi/phi angles 侧链的旋转 优化问题，但是 搜索空间极其巨大 局部极值点

三种预测方法 ab initio 方法 根据第一原理 计算量极大，实际上不可行 同源建模方法： 基本假设：序列同源－ > 结构相似 有效，但是必须具有同源的序列 Threading 方法： 基本假设 : 自然界中蛋白质主链模式是有限的 ~90% 新蛋白质和 PDB 某个已知蛋白质结构相似 推论 : 多个蛋白质会具有相同的主链模式 预测问题－ > 识别问题 能够处理序列上不相似，但是结构相似的情况

Threading 方法 思路： 将序列尽可能好地放入结构模板中； 设计评价函数，对匹配情况进行打分； 关键 : 已知的结构模板库 衡量匹配情况的打分函数 寻求最优的算法； 序列： MTYKLILNGKTKGETTTEAVDAATAEKVFQYANDNGVDGEWTYTE 模板库：

数学描述： MTYKLILNGKTKGETTTEAVDAATAEKVFQYANDNGVDGEWTYTE 残基和结构的匹配： environment: E_s 两个残基相邻的衡量 : E_p 空位罚分 : E_g min （ E_p + E_s + E_g ）

Protein Threading by PROSPECT prediction examples from CASP3 contest actualpredicted actual predicted t49 t68 t57 t70

5. 蛋白质 DOMAIN 识别 生物学观点： 一个蛋白质结构可以包含多个 DOMAIN: DOMAIN 是蛋白质折叠、功能和演化的 基本单位 不同的蛋白质具有相同的 DOMAIN 识别 DOMAIN 有助于蛋白质折叠 数学观点： 最小割 actin

DOMAIN 识别 生物学的不严格表述： DOMAIN 连接紧致，接近球状 DOMAIN 之间作用相对较弱 可操作的定义： DOMAIN 内部残基相互作用较强 DOMAIN 之间残基相互作用较弱 现有识别方法不实用 SCOP 数据库靠手工来维护

DOMAIN 识别与最小割 bottleneck interface

Network Flow Problem source sink capacity edge node Ford-Fulkerson Theorem: the minimum cut of a network is equal to its maximum flow

最小割 节点：一个节点表示一个残基 边：残基－残基之间的相互作用 容量：根据生物学知识，比如相互作用 的种类和强度确定边的容量

经典解法及其结果： maximum flow: Edmonds-Karp algorithm (1972) enumeration of all minimum cuts: Picard-Queyranne algorithm (1980) complexity: O (n 3 ) (n is the number of residues in structure)

6. 序列注释 输入： DNA 序列 输出： 各个功能位点：基因、启动子、 ncRNA 。。。 可以利用的知识： 生物学规律 正例和反例 当前最好的方法： HMM 形式语言

7. 蛋白质质谱鉴定 生物学问题： 原有的 Edman 方法昂贵、耗时 根据质谱来测定蛋白质氨基酸序列 什么是质谱？

样品制备 酶切 或者 物理方法切割

一级质谱：片段分离 Mass Spectrometry 短肽段离子化，经过电场，依据荷质比的不同分离

二级质谱：片段再打断

肽段可能形成的离子

MS/MS Spectrum y-ions b-ions LGEYGFQNALLVR

DeNOVO: 从质谱猜原始序列 m/z KLEDEELFGS % Relative Abundance

我们的工作

组织结构： 生物组： 4 个硕士 提生物学问题 算法组： 5 人 部分独立的算法研究，做技术储备 生物信息算法设计与分析 软件组： 8 人 高性能算法、机器 软件开发

合作： 生物物理所 国家药物筛选中心 华大基因中心 李明老师

基于 Local Search 的序列拼接算法

现有算法： Hamilton 路径类算法 Eulerian 路径类算法

Hamilton 路径类算法 包括： Phrap, CAP3, TIGR, GigAssembler 生成图： 结点：每个片段自成一个结点 边：如果两个结点间有 Overlap 沿 DNA 序列从头走到尾，将经过每个结点一 次且仅一次 Hamilton Path

拼接三步曲 STEP 1 。 Overlap STEP 2 。 Layout STEP 3 。 Consensus/Mosaic

STEP1 。 Overlap 这一步对所有的 Read 进行两两比对，通 常采用快速 Smith-Waterman 算法，以确 定两个 Read 之间是否有 Overlap 。 考虑到各个碱基的出错概率，常常对 Overlap 进行打分，衡量 Overlap 的可能 性高低，一般采用 LLR （ Log Likehood Ratio ）方法打分。

STEP 2 。 Layout 根据 Read 之间的重叠信息形成 Contig ， 即将各个 Read merge 起来，形成一个逐 次链接的链接体。 这一步实际上是在求一条 Hamilton Path ， 通常采用的是贪心法。

STEP 3 。 Consensus 对于每个 Contig ，按照投票或者其他的 原则计算出一个 Sequence 。

评价 判定是否存在 Hamilton Path 是 NP 完全问 题 现在采用的是贪心法 不存在有效算法 出错

出错例示

EULER Path 类算法 转换成图论问题 de Bruijn 图 结点： l-mers 边：两个 l-mers 重叠（ l-1) 个单元 片段：图上一条路径 沿着 DNA 从头走到尾，将经过每条边一 次且仅一次。 Euler Path 附加条件：经过每条路径的特殊 Euler Path.

STEP 1 。 Consensus 目的：排除 Read 中的错误，获得 Error － Free 的 Read 思路：使用 Gk 的近似－将所有的 Read 切 割成小片 k-mers k-mers: 是 Solid 如果出现在至少 M 个 Read 中；否则称为 Weak 。

STEP 2 。 de Bruijn 图的构造 结点：每个 k-mers 是结点 边： de Bruijn 边的构造方法， v-u 当前仅当 v 的 尾巴和 w 的头相同。 每个 Read 表示成一条 Path ， 每个 Repeat 表示成一个多入口、多出口的单一 链，但是不知道出入口之间的对应关系。如果 没有 Read 来覆盖这条单一链，则称为 Tangle 。

STEP 3. Eulerian Super Path 在 de Bruijn 图中寻找一条 Eulerian Path ， 能够覆盖所有的 Path 。 变换方法，等价变换，使得成为寻找一 条 Eulerian Path 的问题。

新方法：机器学习的方法 原始序列：概念 C 正例 : 原始序列的所有子串 反例：非子串 目标：通过正例和反例，学习出原始序 列的近似 C` 。 加强：不仅仅和已知的训练集一致，而 且对于未知的样本，犯错误的概率很小

PAC 模型 原始概念 C, 正例集 POS ，反例集 NEG ， 如果算法 A 对于给定的 ε,δ ，在多项式时 间内结束，并输出概念 C ’ ， 满足： C ’ 和 C 的误差小于 ε 则称为 PAC 算法

拼接与最短公共超串 Blumer 定理： 对于概念 C, 如果使用 m/2 个正例,m/2 个反 例，并且学习得到概念 C ’ 的 size 小于某个 值时，那么 A 是一个 PAC 算法。 利用最短公共超串，满足 Blumer 定理的 算法。

原始算法： Group － Merge 算法 使用 O(n logn logn ) 个片段 以（ 1- δ) 的概率保证： 学习得到串 C ’ 和原始串的误差小于 ε

评价 优点： 坚实的理论基础 对结果的加强： 不仅仅是包含所有子串 而是犯错误的概率最小 缺点： 慢 效果不好

Local Search 算法： 可行解： 片段的排列 相邻关系： 两个片段交换位置 目标函数： 超串长度

算法框架： (0) 计算并保存每两个子串之间的 overlap 数 Repeat: (1) 随机选择排列 P, 得到相应的超串的长度 t, (2)Repeat: 在 P 的邻域内进行搜索, 如果邻域内存在 P ’ 对应的超串长度小于 t, 则 P<-P ’ ; Until: P 是邻域内的最小点。 Until 终止条件满足

Heuristics:  1. 将排列环起来，搜索最短的环状超串  2. 整个 group 的移动  3. 寻找排列中使得 overlap 最小的地方，将 环截断，形成多个线性的子超串 contigs

结果：  找出的串可以近似看作最短超串，算法 可独立用于求最短公共超串问题；  重复的片段在 contigs 的两侧；  通常按这种方法找出的 contig 是原始串的 子串；  拼接水稻基因组联盟 1,10,11 号染色体， 效果得到了很高的评价  用于华大基因中心水稻完全图项目

DeNOVO 算法

质谱的简化模型： 考虑 26 个英文字母， 每个字母有权重 W(c) 对字符串 S ， 已知质谱＝ 问：字符串＝？

难点： 1 。一次打断形成两个离子，混杂在一起 2 。峰的类型不知道，是 b 离子,y 离子还 是混合体？ 3 。离子的修饰：脱氨、脱水、同位素等

组合优化问题： 已知谱 L ， 求序列 S=A 0 A 1 A 2 … A n, 满足： 同时 max f( spectrum(S), L ) 根据生物学知识确定合适的函数 f

结果 1 ： 答案 : LVNELTEFA K 结果： LVNELTEFA K LVNELTHLP K LVNELTYSP K

结果 2 ： 答案 : HPEYAVSVL LR 结果： HPEYAVSVL LR HPEYAVWLL R HPEYAVPVF PK

结果 3 ： 答案 : HLVDEPQNL LK 结果： HLVDEAGP NLLK HLVDEQPNL LK HLVDEPQNL LK

蛋白质相互作用网络分析

复杂的蛋白质相互作用网络

酵母 : 2617 蛋白, 连接

任务 1 ：根据拓扑关系聚类 目标：将蛋白质分类，使得 类内蛋白质连接紧密 类间蛋白质连接稀疏 思路： 第一步：先转化到欧式空间 第二步：使用 Ward 法聚类

转换方法： 寻求方向 Y ，满足： 每个结点都尽量得靠近其邻居 即： min 定理： Y 是对称阵 A T (L-A)A 的最小特征值 对应的特征向量 其中 L 是 Laplacian 矩阵

聚类结果分析： 1 。类内连接紧密，类间稀疏 2 。同一类蛋白质功能类似

聚类谱系图

任务 2 ：蛋白质相互作用网络的谱分析

谱分析 1 。正特征值相应的特征向量，绝对值较 大分量相应蛋白质近似成团； 2 。正特征值相应的特征向量，绝对值较 大分量相应蛋白质近似成二部图；

团和二部图

团： 分析：同一团的蛋白质生物学功能相似 应用：预测未知蛋白质的功能 结果： 分析了 48 个团， 预测了 100 个未知蛋白质的功能 部分结果 Nature 5 月份文章实验验证

二部图： 分析： 同一复合物的不同亚基 重要蛋白质的备份 一个完整生化流程，不同阶段的反应物 应用： 预测 Pathway

目前工作： 1 。系统生物学， ncRNA 研究 2 。比较基因组学： 猪、人、鼠非编码区比较 3 。生物信息专用计算机： 快速的 Blast 及全基因组比对

Thanks!