1、生物信息学中的算法问题 Version 1.0,卜东波 贺思敏 赵 屹中科院计算所生物信息学课题组华大曙光生物信息学联合实验室2002/12/03,主要内容,生物信息学中的算法问题我们的工作 (ICT & IBP & BGI),一、生物学 vs 信息科学,生物信息学的研究目标,特点：天然的形式化碱基：A,C,T,G四种常见氨基酸：20种目标：以DNA序列作为源头揭示“基因组信息结构的复杂性及遗传语言的根本规律”；之后进行蛋白质结构和功能预测。,生物信息学的两个挑战,高性能计算：海量的数据每14个月翻一番算法：海量的数据使得原有算法不适用新需求,生物信息学的研究流程,第一步：生物学问题的提出生物
2、学为主第二步：数学建模、算法设计信息科学为主第三步：结果解释、实验验证生物学,生物信息学脉络,生物信息学问题概览（1）,基因组时期：序列结构功能DNA测序和拼接比对进化树蛋白质质谱鉴定序列注释：基因预测、细胞定位结构预测：RNA结构预测、蛋白质折叠。,生物信息学问题概览（2）,后基因组时期：相互作用网络功能生物芯片（DNA芯片、蛋白质芯片）相互作用网络调控网络E-Cell药物设计。,1. 大规模测序和拼接,生物学问题：从DNA片段恢复原始序列,DNA整体,切成小段,小段和载体结合,结合后进行测序,全自动的测序仪器：MegaBace,需要拼接！,因为整个基因组太长（上M),而每次只能测得一个50
3、0的小片断(read)问题：如何根据read恢复原始顺序？类比：10本圣经，都从随机点起始剪成500个字母左右的小纸条，问：给你这么一堆小纸条，你能读出圣经来吗？,拼接问题的数学描述,数学问题：公共超串输入：设有字符串S，预先估计其长度大约为n，现在已知一个字符串集合R=R1,R2Rn，其中每个Ri都是S的一个子串。问：原始序列S是什么？算法：Hamiltion路径类Euler路径类Local Search类,2. 序列比对,生物学问题：序列的相似性同源性原始序列：S: acgctgT: catgt可行解：S: a c g c t g - T: - c a t g t S: a c g c t
4、 g - T: - c a t g t 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. 结构预测,结构大致决定功能 一级结构
5、(氨基酸序列)二级结构 (螺旋、片层、回环)超二级结构（aba)三级结构 (由二级结构组合成三维构像)四级结构：多个亚基实验测定方法：x-ray晶体衍射NMR核磁共振实验耗时、昂贵一个蛋白质结构测定需要$200K or more需数月或者更长有些蛋白质还无法测定,蛋白质结构（2）,理论上可计算的。能量最低原则变元：主干的psi/phi angles侧链的旋转优化问题，但是搜索空间极其巨大局部极值点,三种预测方法,ab initio 方法根据第一原理计算量极大，实际上不可行同源建模方法：基本假设：序列同源结构相似有效，但是必须具有同源的序列Threading方法：基本假设: 自然界中蛋白质主链模
6、式是有限的90% 新蛋白质和PDB某个已知蛋白质结构相似推论: 多个蛋白质会具有相同的主链模式预测问题识别问题能够处理序列上不相似，但是结构相似的情况,Threading方法,思路：将序列尽可能好地放入结构模板中；设计评价函数，对匹配情况进行打分； 关键:已知的结构模板库衡量匹配情况的打分函数寻求最优的算法；,序列： MTYKLILNGKTKGETTTEAVDAATAEKVFQYANDNGVDGEWTYTE模板库：,数学描述：,MTYKLILNGKTKGETTTEAVDAATAEKVFQYANDNGVDGEWTYTE,残基和结构的匹配： environment: E_s,两个残基相邻的衡量:
7、E_p,空位罚分: E_g,min （ E_p + E_s + E_g ）,Protein Threading by PROSPECT,prediction examples from CASP3 contest,actual,predicted,actual,actual,actual,predicted,predicted,predicted,t49,t68,t57,t70,5. 蛋白质DOMAIN识别,生物学观点：一个蛋白质结构可以包含多个DOMAIN: DOMAIN是蛋白质折叠、功能和演化的基本单位不同的蛋白质具有相同的DOMAIN识别DOMAIN有助于蛋白质折叠数学观点：最小割,ac
8、tin,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,最小割,节点：一个节点表示一个残基边
9、：残基残基之间的相互作用容量：根据生物学知识，比如相互作用的种类和强度确定边的容量,经典解法及其结果：,maximum flow:Edmonds-Karp algorithm (1972)enumeration of all minimum cuts:Picard-Queyranne algorithm (1980)complexity: O (n3)(n is the number of residues in structure),6. 序列注释,输入：DNA序列输出：各个功能位点：基因、启动子、ncRNA。可以利用的知识：生物学规律正例和反例当前最好的方法：HMM形式语言,7. 蛋白质质
10、谱鉴定,生物学问题：原有的Edman方法昂贵、耗时根据质谱来测定蛋白质氨基酸序列什么是质谱？,样品制备,一级质谱：片段分离,MassSpectrometry,短肽段离子化，经过电场，依据荷质比的不同分离,二级质谱：片段再打断,肽段可能形成的离子,MS/MS Spectrum,y-ions,b-ions,LGEYGFQNALLVR,DeNOVO: 从质谱猜原始序列,m/z,147,260,389,504,633,762,875,1022,1080,1166,1166,1020,907,778,663,534,405,292,145,88,% Relative Abundance,100,0,25
11、0,500,750,1000,我们的工作,组织结构：,生物组：4个硕士提生物学问题算法组：5人部分独立的算法研究，做技术储备生物信息算法设计与分析软件组：8人高性能算法、机器软件开发,合作：,生物物理所国家药物筛选中心华大基因中心李明老师,基于Local Search的序列拼接算法,现有算法：,Hamilton路径类算法Eulerian路径类算法,Hamilton路径类算法,包括：Phrap, CAP3, TIGR, GigAssembler生成图：结点：每个片段自成一个结点边：如果两个结点间有Overlap沿DNA序列从头走到尾，将经过每个结点一次且仅一次Hamilton Path,Hami
12、lton Path,拼接三步曲,STEP 1。OverlapSTEP 2。LayoutSTEP 3。Consensus/Mosaic,STEP1。Overlap,这一步对所有的Read进行两两比对，通常采用快速Smith-Waterman算法，以确定两个Read之间是否有Overlap。考虑到各个碱基的出错概率，常常对Overlap进行打分，衡量Overlap的可能性高低，一般采用LLR（Log Likehood Ratio）方法打分。,STEP 2。Layout,根据Read之间的重叠信息形成Contig，即将各个Read merge起来，形成一个逐次链接的链接体。这一步实际上是在求一条Ha
13、milton 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中的错误，获得ErrorFr
14、ee 的Read思路：使用Gk的近似将所有的Read切割成小片k-mersk-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，能够覆盖所有的Pa
15、th。变换方法，等价变换，使得成为寻找一条Eulerian Path的问题。,新方法：机器学习的方法,原始序列：概念C正例: 原始序列的所有子串反例：非子串目标：通过正例和反例，学习出原始序列的近似C。加强：不仅仅和已知的训练集一致，而且对于未知的样本，犯错误的概率很小,PAC模型,原始概念C, 正例集POS，反例集NEG，如果算法A对于给定的,，在多项式时间内结束，并输出概念C，满足：C和C的误差小于则称为PAC算法,拼接与最短公共超串,Blumer定理：对于概念C,如果使用m/2个正例,m/2个反例，并且学习得到概念C的size小于某个值时，那么A是一个PAC算法。利用最短公共超串，满足B
16、lumer定理的算法。,原始算法：,GroupMerge算法使用O(n logn logn )个片段以（1- )的概率保证：学习得到串C和原始串的误差小于,评价,优点：坚实的理论基础对结果的加强：不仅仅是包含所有子串而是犯错误的概率最小缺点：慢效果不好,Local Search算法：,可行解：片段的排列相邻关系：两个片段交换位置目标函数：超串长度,算法框架：,(0)计算并保存每两个子串之间的overlap数 Repeat: (1)随机选择排列P,得到相应的超串的长度t, (2)Repeat: 在P的邻域内进行搜索, 如果邻域内存在P对应的超串长度小于t, 则P-P; Until: P是邻域内的
17、最小点。 Until 终止条件满足,Heuristics:,1. 将排列环起来，搜索最短的环状超串2. 整个group的移动3. 寻找排列中使得overlap最小的地方，将环截断，形成多个线性的子超串contigs,结果：,找出的串可以近似看作最短超串，算法可独立用于求最短公共超串问题；重复的片段在contigs的两侧；通常按这种方法找出的contig是原始串的子串；拼接水稻基因组联盟1,10,11号染色体，效果得到了很高的评价用于华大基因中心水稻完全图项目,DeNOVO算法,质谱的简化模型：,考虑26个英文字母，每个字母有权重W(c)对字符串S，已知质谱问：字符串？,难点：,1。一次打断形成
18、两个离子，混杂在一起2。峰的类型不知道，是b离子,y离子还是混合体？3。离子的修饰：脱氨、脱水、同位素等,组合优化问题：,已知谱L，求序列S=A0A1A2An ,满足：同时 max f( spectrum(S), L )根据生物学知识确定合适的函数f,结果1：,答案:LVNELTEFAK结果：LVNELTEFAKLVNELTHLPKLVNELTYSPK,结果2：,答案:HPEYAVSVLLR结果：HPEYAVSVLLRHPEYAVWLLRHPEYAVPVFPK,结果3：,答案:HLVDEPQNLLK结果：HLVDEAGPNLLKHLVDEQPNLLKHLVDEPQNLLK,蛋白质相互作用网络分
19、析,复杂的蛋白质相互作用网络,酵母: 2617蛋白, 11855连接,任务1：根据拓扑关系聚类,目标：将蛋白质分类，使得类内蛋白质连接紧密类间蛋白质连接稀疏思路：第一步：先转化到欧式空间第二步：使用Ward法聚类,转换方法：,寻求方向Y，满足：每个结点都尽量得靠近其邻居即：min定理：Y是对称阵AT (L-A)A的最小特征值对应的特征向量其中L是Laplacian矩阵,聚类结果分析：,1。类内连接紧密，类间稀疏2。同一类蛋白质功能类似,聚类谱系图,任务2：蛋白质相互作用网络的谱分析,谱分析,1。正特征值相应的特征向量，绝对值较大分量相应蛋白质近似成团；2。正特征值相应的特征向量，绝对值较大分量相应蛋白质近似成二部图；,团和二部图,团：,分析：同一团的蛋白质生物学功能相似应用：预测未知蛋白质的功能结果：分析了48个团，预测了100个未知蛋白质的功能部分结果Nature 5月份文章实验验证,二部图：,分析：同一复合物的不同亚基重要蛋白质的备份一个完整生化流程，不同阶段的反应物应用：预测Pathway,目前工作：,1。系统生物学，ncRNA研究2。比较基因组学：猪、人、鼠非编码区比较3。生物信息专用计算机：快速的Blast及全基因组比对,,Thanks!,
