690 likes | 917 Views
Lecture 14:. 基因发现算法的比较和分析. 何淼 PhD 改编自 Prof. Chun-Ting Zhang Center of Bioinformatics Tianjin University http://tubic.tju.edu.cn. 内容. 基因发现算法简介 不同算法比较 算法分析 基于 Z-curve 算法的基因发现软件 结语. 1984 : 30 亿美元测定一个人类基因组. 1994 : 3 亿美元测定一个人类基因组 (1:10). 2004 : 3 千万美元测定一个人类基因组 (1:100).
E N D
Lecture 14: 基因发现算法的比较和分析 何淼PhD 改编自 Prof. Chun-Ting Zhang Center of Bioinformatics Tianjin University http://tubic.tju.edu.cn
内容 • 基因发现算法简介 • 不同算法比较 • 算法分析 • 基于Z-curve算法的基因发现软件 • 结语
1984: 30亿美元测定一个人类基因组 1994: 3亿美元测定一个人类基因组(1:10) 2004:3千万美元测定一个人类基因组(1:100) 2006:1百50万美元测定一个人类基因组(1:2000) 2008:50万美元测定一个人类基因组(1:6000) 未来目标:1000/100 美元测定一个人类基因组
人类基因组研究正在推向纵深不可或缺的三个基本主线人类基因组研究正在推向纵深不可或缺的三个基本主线 时间 我国从参与到独立启动 比较基因组计划 人类基因组计划 1 00 脊椎动物 哺乳类 灵长类 脊索动物 多细胞 人类基因组单倍体型图计划 10 单细胞 05 癌症基因组计划 特定群体基因组 单倍型体图 病原体基因组计划 10 多发疾病基因组计划 100 常见疾病基因组计划 个体化基因组单倍体图 环境疾病基因组计划 15 老龄疾病基因组计划
1.介绍 什么是计算机辅助基因发现? Mathematical algorithms Core; Computer Tool; Finding genes in un-annotated genomic sequence 什么类型的基因需要发现? 两类: protein coding genes ,RNA genes 在完整的蛋白质编码基因发现过程中包含哪些任务? 两个任务: 找到CDS和基因的其他部分,包括promoter等
1.介绍 有多少种 CDS发现工作 ? 两种: 原核生物ORFs (Open Reading Frames) 真核生物introns,exons, 剪切位点 基因发现研究的三种基本途径是什么? 组成 信号 相似性的原理
基于数据库和字统计的算法 2.算法I 1.密码子选择 编码区中的密码子使用频率表 非编码区中的“密码子”使用频率 S = AGGACGGGATCA … PS(AGG) = 0.01209; PS(ACG) = 0.00680, … PS = PS(AGG) PS(ACG) … P0(AGG) = 1/64 = 0.056; P0(ACG) =1/64 = 0.056 P0 = P0(AGG) P0(ACG) … 定义似然率的算法:LP = ln (PS/P0) 如果 LP > 0 编码区 如果 LP < 0 非编码区 滑移窗口 (具有重叠) 技术 例如: 窗口长度 = 120 bp 重叠长度 = 10 bp 参考文献Staden, R. & McLachlan, A. (1982) Nucleic Acids Res. 10, 141-156. 似然率(likelihood ratio,LR)是反映真实性的一种指标,属于同时反映灵敏度和特异度的复合指标。即有病者中得出某一筛检试验结果的概率与无病者得出这一概率的比值。该指标全面反映筛检试验的诊断价值,且非常稳定。似然比的计算只涉及到灵敏度与特异度,不受患病率的影响。
人类密码子使用频率表 2.算法I Gly GGG 17.08 0.23 Arg AGG 12.09 0.22 Trp TGG 14.74 1.00 Arg CGG 10.40 0.19 Gly GGA 19.31 0.26 Arg AGA 11.73 0.21 End TGA 2.64 0.61 Arg CGA 5.63 0.10 Gly GGT 13.66 0.18 Ser AGT 10.18 0.14 Cys TGT 9.99 0.42 Arg CGT 5.16 0.09 Gly GGC 24.94 0.33 Ser AGC 18.54 0.25 CyS TGC 13.86 0.58 Arg CGC 10.82 0.19 Glu GAG 38.82 0.59 Lys AAG 33.79 0.60 End TAG 0.73 0.17 Gln CAG 32.95 0.73 Glu GAA 27.51 0.41 Lys AAA 22.32 0.40 End TAA 0.95 0.22 Gln CAA 11.94 0.27 Asp GAT 21.45 0.44 Asn AAT 16.43 0.44 Tyr TAT 11.80 0.42 His CAT 9.56 0.41 Asp GAC 27.06 0.56 Asn AAC 21.30 0.56 Tyr TAC 16.48 0.58 His CAC 14.00 0.59 Val GTG 28.60 0.48 Met ATG 21.86 1.00 Leu TTG 11.43 0.12 Leu CTG 39.93 0.43 Val GTA 6.09 0.10 Ile ATA 6.05 0.14 Leu TTA 5.55 0.06 Leu CTA 6.42 0.07 Val GTT 10.30 0.17 Ile ATT 15.03 0.35 Phe TTT 15.36 0.43 Leu CTT 11.24 0.12 Val GTC 15.01 0.25 Ile ATC 22.47 0.52 Phe TTC 20.72 0.57 Leu CTC 19.14 0.20 Ala GCG 7.27 0.10 Thr ACG 6.80 0.12 Ser TCG 4.38 0.06 Pro CCG 7.02 0.11 Ala GCA 15.50 0.22 Thr ACA 15.04 0.27 Ser TCA 10.96 0.15 Pro CCA 17.11 0.27 Ala GCT 20.23 0.28 Thr ACT 13.24 0.23 Ser TCT 13.51 0.18 Pro CCT 18.03 0.29 Ala GCC 28.43 0.40 Thr ACC 21.52 0.38 Ser TCC 17.37 0.23 Pro CCC 20.51 0.33
基于数据库和字统计的算法 2.算法I 2. 氨基酸使用频率 在CDS中和在非编码区中分别发生的平均氨基酸频率 S = AGGACGGGATCA … AGG Arg, ACG Thr, … PS(Arg) = 0.056; PS(Thr) = 0.075, … PS = PS(Arg) PS(Thr) … P0(Arg) = Nd/61 = 6/61; P0(Thr) =Nd/64 = 4/61, ... P0 = P0(AGG) P0(ACG) … 定义似然率算法为: LP = ln (PS/P0) 如果 LP > 0 coding; 如果 LP < 0 non-coding 参考文献MaCaldon, P. & Argos, P. (1988) Proteins 4, 99-122.
基于数据库和字统计的算法 2.算法I 3. 密码子偏好性 基于同义密码子的相对使用频率. 例如, PS(GGG) = 0.23; PS(GGA) = 0.26; PS(GGT) = 0.18; PS(GGC) = 0.33. P0 = 1/Nd = 1/4 (for Gly) 参考文献Glibskov, M. et al. (1984) Nucleic Acids Res. 12, 539. 4.六聚物使用频率(6-mers= 6-nucleotides) 444444 = 4096 frequencies of 6-mers 1 6, 7 12, ….. P0 = 1/4096 参考文献Claverie, J.-M., et al. (1990) Meth. Enzymol. 183, 237.
基于数据库和字统计的算法 2.算法I 5.密码子原型 Shepherd (1981) —— RYN code Trifinov (1987) —— G-nonG-N code Zhang (1991-1994) —— R-nonG-S/W code 在第一密码子位置的R-dominant原则, 适用于所有物种 在第二密码子位置的G-least dominant原则, 适用于所有物种 在所有物种中存在S/W-相关性 参考文献 Shepherd, J.C. (1981) PNAS, USA 78,1596-1600. Trifonov, E. (1987) J. Mol. Biol. 194, 643-692. Zhang, C.T. & Zhang, R. (1991) Nucleic Acids Res. 19, 6313-6317. Zhang, C.T. & Chou, K.C. (1994) J. Mol. Biol. 238, 1-8.
基于数据库和字统计的算法 2.算法I S = AGGACGGGATCA … PS(AGG) = p(A1)p(G2)P(G3) = 0.270.200.29; PS(ACG) = P(A1)P(C2)p(G3) = 0.270.240.29; … PS = PS(AGG) PS(ACG) … P0(codon) = (1/4)(1/4)(1/4) = 1/64; P0 = P0(AGG) P0(ACG) … 定义似然率算法为: LP = ln (PS/P0) 如果 LP > 0 coding; 如果 LP < 0 non-coding
基于数据库和字统计的算法 2.算法I 6.Markov 模型 (1-8阶) 关键词:马尔可夫链(Markov chain); 阶(order); 转移概率(transition probabilities); 转移矩阵(transition matrix) S = AGGACGGGATCA … PS = p(A1)p1(GA)p2(GG)p3(AG)p1(CA)p2(GC)… p0 = 0.250.250.250.250.250.25…. 定义似然率算法为: LP = ln (PS/P0) 如果 LP > 0 coding; 如果 LP < 0 non-coding 参数的数目: 1阶: 44 + 3(44) = 16 + 48 = 64; 2阶: 416 + 3(416) = 256; … … 5阶: 4096 (34096 = 12288, GeneMark) 参考文献Zhang, C.T. (1990) J. Theor. Biol. 142, 281-284. Borodovsky, M. & McIninch, J. (1993) Comp. Chem. 17, 123-130. (GeneMark) Salzberg, S. L. et al. (1998) Nucleic Acids Res. 26, 544-548. (GLIMMER)
基于数据库和字统计的算法 2.算法I 当临近密码子位置核苷酸已知条件下,四个核苷酸在不同密码子位置的概率 A C G T Codon position 1 Codon position 3 A 0.22 0.33 0.24 0.13 C 0.21 0.29 0.27 0.21 G 0.44 0.15 0.37 0.53 T 0.13 0.22 0.12 0.13 Codon position 2 Codon position 1 A 0.36 0.27 0.35 0.18 C 0.21 0.23 0.24 0.27 G 0.19 0.14 0.23 0.23 T 0.24 0.35 0.19 0.31 Codon position 3 Codon position 2 A 0.16 0.19 0.15 0.07 C 0.28 0.44 0.41 0.33 G 0.40 0.12 0.27 0.45 T 0.16 0.25 0.17 0.16
基于公理和字统计的算法 2.算法II 1.位置的不对称性 pS (base)= (pS(base-1)+pS(base-2)+pS(base-3))/3 asym(base) = i:1-3[pS (base-i) - pS(base)]2 PAS = asym(A) + asym(C) + asym(G) + asym(T) 参考文献Fickett, J.W. (1982) Nucleic Acids Res. 10, 5303-5308. Staden, R. (1984) Nucleic Acids Res. 12, 551-567. 2. 周期不对称性指数(Periodic Asymmetry Index) 考虑一个重要的统计量: Nij(k) S = AGGACGGGATCA NGA(1) = 2; NAT(0) = 1; NGG(0) = 3; NAA(7) = 2, …… 相同核苷酸(AA, CC, GG, TT)的概率 at k = 2, 5, 8, … p1 k = 0, 3, 6, … p2 k = 1, 4, 7, … p3 PAIS = max (p1, p2, p3)/min (p1, p2, p3) 参考文献Konopka, A.K. (1994) In Structure and Methods: VI. Human Genome Initiative and DNA Recombination. (eds. R. H. Sarma, et al.) Adenine Press, pp. 113-125.
基于公理和字统计的算法 2.算法II 3.平均相互信息(Average Mutual Information) Mutual information I (k) (Shannon, 1948) I(k) = Sigma Pijlog[Pij/(pipj)], i, j A, C, G, T I(k) at k = 2, 5, 8, … Iin I(k) at k = 4, 7, 10, … Iout AMIS = (Iin + Iout + Iout)/3 参考文献Herzel, H. & Grosse, I. (1995) Physica A 216, 518-542. 4.傅立叶能谱(Fourier Power Spectrum) 基本事实: 3-base 周期 参考文献Tiwari, S., et al. (1997) Comp. Appl. Biosci. (CABIOS) 13, 263-270. (GeneScan) Zhang, C.T. et al. (1998) Bioinformatics 14, 685-690.
基于Z curve原理的算法 2.算法III Z Curve 令An , Cn , Gn ,Tn为研究的DNA序列中子序列由第一个碱基开始至第n个碱基结束碱基A, C, G , 和T 的累积数. DNA序列的 Z 变换为 xn = (An + Gn) – (Cn + Tn) = Rn – Yn , yn = (An + Cn) – (Gn + Tn) = Mn – Kn , zn = (An + Tn) – (Gn + Cn) = Wn – Sn , xn , yn , zn [-N, N]. 则(xn , yn , zn), 其中n = 0, 1, 2, 3, … 被称为用于研究DNA序列的 Z curve. Z curves用于基因组研究的例子 Z curves的简化 利用Z curve 终止位点的简化 利用直线的简化 利用高阶多项式的简化 利用剪切功能的简化
基于Z curve原理的算法 2.算法III 基于Z Curve理论的基因发现方法 特定相位(Phase-specific)Z curve Bases at 1, 4, 7, phase-1 curve, xn(1) , yn(1) , zn(1); Bases at 2, 5, 8, phase-2 curve, xn(2) , yn(2) , zn(2); Bases at 3, 6, 9, phase-3 curve, xn(3) , yn(3) , zn(3). 用直线拟合9 个组分的曲线: xn(1) = u1n, yn(1) = u2n, zn(1) = u3n ; xn(2) = u4n, yn(2) = u5n, zn(2) = u6n ; xn(3) = u7n, yn(3) = u8n, zn(3) = u9n . 用 (u1, u2, …, u9)生成了一个 9-D 空间。 Fisher 线性判别算法 结果: 对于原核生物和一些真核生物的基因组来说,两者的敏感度和特异性 指标介于: 95.0 % - 99.8 % 参考文献 Zhang, C.T. & Wang, J. (2000) Nucleic Acids Res. 28, 2804-2814. Guo, F.B., Ou, H.Y. & Zhang, C.T. (2003) Nucleic Acids Res. 31, 1780-1780
(给位) (受位) 不同算法的比较 Gene FindingApproaches • 信号搜索 • (SIGNAL) • 组成识别 • (CONTENT) • 序列比对 • (HOMOLOGY)
基因识别算法的评估与比较 研究背景 • 基因识别算法的核心是建立一个衡量序列编码能力的标准尺度(coding measures),通常是借助于对任意窗口的序列计算出的一个数或矢量来衡量序列的编码能力。常见的例子包括密码子使用矢量,碱基组成矢量和对序列进行的某些类型的傅立叶变换。 • Fickett和Tung在1992年发表的综述文章中列出了当时关于基因识别的二十种算法,并认为六核苷酸频率法(hexamer usage) 是描述编码区的最好方法。 • 自Fickett和Tung在1992年发表综述文章以来的十余年间,又出现了许多新的基因识别算法,有必要进行重新的评价和比较。 参考文献Fickett, J.W. and Tung, C.S. (1992) Nucleic Acids Res., 20, 6441-6450.
评估基因识别算法 Z曲线方法,密码子使用,六核苷酸频率法,码的偏好,氨基酸使用,码的原型,延长打乱快速傅立叶变换 ,马尔科夫模型 评价指标:灵敏度、特异度和准确率 数据库 数据库1:人类基因编码/非编码序列数据库 (Length = 192, 162, 129, 108, 87, 63 and 42 bp) 数据库2:通过已知mRNA验证的人类基因数据库 (Length = 192, 162, 129, 108, 87, 63 and 42 bp) 参考文献F. Gao and C.T. Zhang (2004). Bioinformatics, 20(5), 673-681.
TP(true positive): 实际编码区的核酸中被成功预测的核酸数目; TN(true negative): 实际非编码区的核酸中被成功预测的核酸数目; FN(false negative): 实际编码区的核酸中被误测为非编码的核酸数目; FP(false positive): 实际非编码区的核酸中被误测为编码的核酸数目。 REALITY coding noncoding coding PREDICTION TP FP TP+FP noncoding FN TN FN+TN TP+FN FP+TN
基于TP、TN、FP、FN,主要引进四个参数:Sn、Sp、CC、AC。基于TP、TN、FP、FN,主要引进四个参数:Sn、Sp、CC、AC。 敏感性(sensitivity,Sn): 特异性(specificity,Sp): Sn:实际编码区核酸序列中被成功预测的比例; Sp:预测编码核酸序列中被成功预测的比例。
DNA序列的Z 曲线理论 The Z curve theory of DNA sequence 描述编码区特征所采用的Z曲线变量 • 用Z曲线表示相位特异单核苷酸频率 (9) (1) • 用Z曲线表示相位特异双核苷酸频率 (36) (2)
用Z曲线表示相位特异三核苷酸频率 (144) (3) • 用Z曲线表示非相位特异双核苷酸频率 (12) (4) • 用Z曲线表示非相位特异三核苷酸频率 (48) (5)
Fisher判别 • 计算任一DNA序列片断的Z曲线的n个参数,将其映射为n维空间V中一点,利用Fisher判别函数对正负样本进行判定。 • Fisher判别的判别方程代表n维空间V中的一个超平面,用一个有n个分量分别为c1, c2, … ,cn的矢量c表示。利用训练集 (包括正负样本) 中的数据,确定判别编码/非编码的阈值c0。 • 判断任一DNA序列片断编码/非编码的准则如下: • Z(u)>0 / Z(u)<0 • 这里Z(u)=cu-c0,称之为Z分数或Z指数。
Z曲线方法对数据库I进行十重交叉检验的平均灵敏度、特异度和准确率Z曲线方法对数据库I进行十重交叉检验的平均灵敏度、特异度和准确率
密码子使用,六核苷酸频率法,码的偏好,氨基酸使用,码的原型,延长打乱密码子使用,六核苷酸频率法,码的偏好,氨基酸使用,码的原型,延长打乱 快速傅立叶变换,对数据库1进行十重交叉检验平均灵敏度、特异度和准确率
马尔科夫模型对数据库1进行十重交叉检验的平均灵敏度、特异度和准确率马尔科夫模型对数据库1进行十重交叉检验的平均灵敏度、特异度和准确率
Z曲线方法对数据库2进行十重交叉检验的平均灵敏度、特异度和准确率Z曲线方法对数据库2进行十重交叉检验的平均灵敏度、特异度和准确率
其他各种算法对数据库2进行十重交叉检验的平均灵敏度、特异度和准确率其他各种算法对数据库2进行十重交叉检验的平均灵敏度、特异度和准确率
利用Z曲线方法识别人类β球蛋白基因中的外显子利用Z曲线方法识别人类β球蛋白基因中的外显子
缩略词: • CU, Codon Usage; • HU, Hexamer Usage; • Cpre, Codon Preference; • AAU, Amino Acid Usage; • Cpro,Codon Prototype; • MM-1, Markov Model k = 1; • MM-2, Markov Model k = 2; • MM-5, Markov Model k = 5; • LS-FFT, lengthen-shuffling FFT; • ZC-9, ZC-69 and ZC-189, the Z curve methods with 9, 69 and 189 parameters, respectively.
算法分析 • 基因识别算法初步的分析 • Local <=> global • Statistical <=> non-statistical
Local ==Global GeneMark ZCURVE Glimmer (Markov chain model)
Statistical == Non-statistical GeneMark ORPHEUS Glimmer CRITICA ZCURVE