中國科技大學(xué)課件系列：《生物信息學(xué)》04講課講稿

1、中國科技大學(xué)課件系列：生物信息學(xué)04r第一節(jié)：數(shù)學(xué)基礎(chǔ)：概率及概率模型第一節(jié)：數(shù)學(xué)基礎(chǔ)：概率及概率模型r第二節(jié)：雙序列比對(duì)算法的介紹第二節(jié)：雙序列比對(duì)算法的介紹Dot matrix動(dòng)態(tài)規(guī)劃算法動(dòng)態(tài)規(guī)劃算法w(Needleman-Wunsch, Smith-Waterman算法算法) FASTA和和BLAST算法算法r第三節(jié)：打分矩陣及其含義第三節(jié)：打分矩陣及其含義r第四節(jié)：多序列比對(duì)第四節(jié)：多序列比對(duì)r71個(gè)蛋白質(zhì)家族的個(gè)蛋白質(zhì)家族的1572種變化；種變化；r序列相似性序列相似性 85%；r 功能同源的蛋白質(zhì)功能同源的蛋白質(zhì) 通過中性進(jìn)化，引入通過中性進(jìn)化，引入可接受的點(diǎn)突變；可接受的點(diǎn)突變；

2、r 進(jìn)化模型：進(jìn)化模型：A. 基本假設(shè)：中性進(jìn)化，基本假設(shè)：中性進(jìn)化，Kimura,1968;B. 進(jìn)化的對(duì)稱性進(jìn)化的對(duì)稱性: A-B = B-A;C. 擴(kuò)展性：通過對(duì)較短時(shí)間內(nèi)氨基酸替代關(guān)系擴(kuò)展性：通過對(duì)較短時(shí)間內(nèi)氨基酸替代關(guān)系的計(jì)算來計(jì)算較長時(shí)間的氨基酸替代關(guān)系；的計(jì)算來計(jì)算較長時(shí)間的氨基酸替代關(guān)系；r 兩個(gè)蛋白質(zhì)序列的兩個(gè)蛋白質(zhì)序列的1%氨基酸發(fā)生變化；氨基酸發(fā)生變化；r 定義進(jìn)化時(shí)間以氨基酸的變異比例為準(zhǔn)，定義進(jìn)化時(shí)間以氨基酸的變異比例為準(zhǔn)，而不是時(shí)間；因?yàn)楦鱾€(gè)蛋白質(zhì)家族進(jìn)化的速而不是時(shí)間；因?yàn)楦鱾€(gè)蛋白質(zhì)家族進(jìn)化的速度并不相等；度并不相等；r PAM2 = PAM1*PAM1 PAM3

3、 = (PAM1)3 PAM250= (PAM1)250r選取多個(gè)家族的相似性選取多個(gè)家族的相似性85%的保守序列；的保守序列；r根據(jù)匹配計(jì)分進(jìn)行多重比對(duì)根據(jù)匹配計(jì)分進(jìn)行多重比對(duì)(不含空位不含空位)；r以比對(duì)結(jié)果構(gòu)建進(jìn)化樹，反映氨基酸替換關(guān)以比對(duì)結(jié)果構(gòu)建進(jìn)化樹，反映氨基酸替換關(guān)系；系；r計(jì)算每種氨基酸轉(zhuǎn)換成其它氨基酸的次數(shù)；計(jì)算每種氨基酸轉(zhuǎn)換成其它氨基酸的次數(shù)；r計(jì)算每種氨基酸突變率；計(jì)算每種氨基酸突變率；r計(jì)算每對(duì)氨基酸突變率，得到突變概率矩陣計(jì)算每對(duì)氨基酸突變率，得到突變概率矩陣，將此矩陣自乘，將此矩陣自乘n次；次；r將突變概率矩陣轉(zhuǎn)化為將突變概率矩陣轉(zhuǎn)化為PAMn矩陣。矩陣。r 已知已知

4、3個(gè)蛋白質(zhì)家族若干保守序列片段：個(gè)蛋白質(zhì)家族若干保守序列片段：家族一：家族一：FKILK，F(xiàn)KIKK，F(xiàn)FILL，F(xiàn)FIKL家族二：家族二：IIFFF， IIFIF ， IKFFL ， IKFIL家族三：家族三： KIFKK，KIFLK，KLFKL，KLFLL按按Doyhoff方法構(gòu)建方法構(gòu)建PAM1與與PAM2矩陣矩陣r位置對(duì)齊，多重比對(duì)（不考慮空位）：位置對(duì)齊，多重比對(duì)（不考慮空位）：r統(tǒng)計(jì)每種氨基酸出現(xiàn)的頻率；統(tǒng)計(jì)每種氨基酸出現(xiàn)的頻率；fi = 氨基酸氨基酸i的數(shù)目的數(shù)目/總氨基酸數(shù)目總氨基酸數(shù)目fL = 12/60 = 0.2.家族一家族一家族二家族二家族三家族三F K I L KI

5、I F F FK I F K KF K I K KI I F I FK I F L KF F I L LI K F F LK L F K LF F I K LI K F I LK L F L Lr最大簡約法最大簡約法家族一家族一:wL和和K間相互轉(zhuǎn)換次數(shù)：間相互轉(zhuǎn)換次數(shù)：N(LK) = 3家族二，家族三家族二，家族三 FKILKFKIKKFKIKKFFIKLFFILLFFIKL(LK)(KF)(LK)(LK)r計(jì)算每種氨基酸轉(zhuǎn)換成其它氨基酸的次數(shù)。計(jì)算每種氨基酸轉(zhuǎn)換成其它氨基酸的次數(shù)。r假設(shè)兩種氨基酸間相互轉(zhuǎn)換一樣。假設(shè)兩種氨基酸間相互轉(zhuǎn)換一樣。e.g. N(LK)= 3 + 0 + 3 =

6、6KFILK116F121I121L611r每種氨基酸相對(duì)突變率每種氨基酸相對(duì)突變率miri：第：第i種氨基酸；種氨基酸；rfi ：每種氨基酸出現(xiàn)的頻率；：每種氨基酸出現(xiàn)的頻率；mK = 8/(122 fK 100) = 0.01251002iifim總替換數(shù)總共發(fā)生替換數(shù)氨基酸r氨基酸氨基酸i替換為替換為j的突變率的突變率mije.g.mKK = 1- mK = 0.9875mKF = mF 1/4 = 0.001389iiiiijmmjijjimmji1時(shí)，總共發(fā)生替換數(shù)氨基酸相互替換的次數(shù)與氨基酸時(shí)，r氨基酸突變概率氨基酸突變概率一步轉(zhuǎn)移概率矩陣一步轉(zhuǎn)移概率矩陣M1ij原氨基酸原氨基酸K

7、FIL替換氨替換氨基酸基酸K0.98750.0015630.0015630.009375F0.0013890.9944440.0027780.001389I0.0017860.0035710.9928570.001786L0.01250.0020830.0020830.983333r由突變率由突變率mij計(jì)算計(jì)分矩陣中的分值計(jì)算計(jì)分矩陣中的分值rij：r將將rij = rji取平均值，再取整數(shù)；取平均值，再取整數(shù)；（按先前假設(shè)，（按先前假設(shè)， rij = rji） rKK = 10lg(mkk/ fk) = 5.6857 6 (rKF + rFK )/2 = -22.833 -23 )/lg(

8、10iijijfmr r三個(gè)家族序列片段得到的三個(gè)家族序列片段得到的PAM1計(jì)分矩陣：計(jì)分矩陣：KFILK6F-235I-22-196L-13-22-207r將氨基酸突變概率矩陣自乘一次，得到兩將氨基酸突變概率矩陣自乘一次，得到兩步轉(zhuǎn)移概率矩陣步轉(zhuǎn)移概率矩陣M2ij M2ij = M1ij M1ijr三個(gè)家族序列片段得到的三個(gè)家族序列片段得到的PAM2計(jì)分矩陣：計(jì)分矩陣：KFILK6F-205I-19-166L-10-19-187r PAM250: 250%期望的突變；期望的突變；r 蛋白質(zhì)序列仍然有蛋白質(zhì)序列仍然有15-30%左右的相似性；左右的相似性；rPAM250: 15-30%的序列相

9、似性的序列相似性；rPAM120: 40%的序列相似性；的序列相似性；rPAM80: 50%rPAM60: 60%r如何選擇最合適的矩陣？如何選擇最合適的矩陣？r 多種嘗試多種嘗試r1. PAM系列矩陣存在的問題：系列矩陣存在的問題：A. 氨基酸的打分矩陣，不關(guān)心核酸；氨基酸的打分矩陣，不關(guān)心核酸；B. 進(jìn)化模型的構(gòu)建需要系統(tǒng)發(fā)育樹的分析，因進(jìn)化模型的構(gòu)建需要系統(tǒng)發(fā)育樹的分析，因此，成為一個(gè)循環(huán)論證的問題：序列比對(duì)此，成為一個(gè)循環(huán)論證的問題：序列比對(duì)矩矩陣構(gòu)建陣構(gòu)建打分打分進(jìn)行新的序列比對(duì)；進(jìn)行新的序列比對(duì)；C. 數(shù)據(jù)集很小；數(shù)據(jù)集很??；r2. 打分矩陣的改進(jìn)打分矩陣的改進(jìn)A. 選用大量的序列

10、數(shù)據(jù)，構(gòu)建選用大量的序列數(shù)據(jù)，構(gòu)建PAM矩陣；矩陣；B. BLOSUM系列矩陣系列矩陣;C. 核酸的打分矩陣核酸的打分矩陣;r最被廣泛使用的氨基酸打分矩陣最被廣泛使用的氨基酸打分矩陣;r根據(jù)蛋白質(zhì)模塊數(shù)據(jù)庫根據(jù)蛋白質(zhì)模塊數(shù)據(jù)庫BLOCKS中蛋白質(zhì)中蛋白質(zhì)序列的高度保守部分的比對(duì)而得到的，最常序列的高度保守部分的比對(duì)而得到的，最常用的是用的是BLOSUM62;rBLOCK: 蛋白質(zhì)家族保守的一段氨基酸，無蛋白質(zhì)家族保守的一段氨基酸，無gap，一般幾個(gè)至上百個(gè)氨基酸；，一般幾個(gè)至上百個(gè)氨基酸；rProsite家族：至少有一個(gè)家族：至少有一個(gè)BLOCK存在于該存在于該家族的所有蛋白質(zhì)序列中；家族的所

11、有蛋白質(zhì)序列中；rBLOSUM62: 序列的平均相似性為序列的平均相似性為62%的的BLOCK構(gòu)建的打分矩陣；構(gòu)建的打分矩陣；r提取提取Prosite數(shù)據(jù)庫中數(shù)據(jù)庫中504個(gè)家族的個(gè)家族的2萬多蛋萬多蛋白質(zhì)序列，合并其中相似性白質(zhì)序列，合并其中相似性62%的序列；的序列；r統(tǒng)計(jì)各統(tǒng)計(jì)各BLOCK的氨基酸對(duì)數(shù)量的氨基酸對(duì)數(shù)量f；r計(jì)算氨基酸對(duì)的出現(xiàn)頻率計(jì)算氨基酸對(duì)的出現(xiàn)頻率q；r計(jì)算每種氨基酸的期望頻率計(jì)算每種氨基酸的期望頻率p；r計(jì)算氨基酸對(duì)出現(xiàn)的期望頻率計(jì)算氨基酸對(duì)出現(xiàn)的期望頻率e；r計(jì)算計(jì)算BLOSUM62矩陣分量矩陣分量rij)/(lg22eqrijr序列相似性與序列相似性與PAM及及B

12、LOSUM矩陣的大致矩陣的大致對(duì)應(yīng)關(guān)系：對(duì)應(yīng)關(guān)系：序列相似性序列相似性 %999080706050403020PAM數(shù)值數(shù)值11123385680112159 246BLOSUM數(shù)值數(shù)值908062-45r 不同物種中，許多基因的功能保守，序列相不同物種中，許多基因的功能保守，序列相似性較高，通過多條序列的比較，發(fā)現(xiàn)保守似性較高，通過多條序列的比較，發(fā)現(xiàn)保守與變異的部分；與變異的部分；r 可構(gòu)建可構(gòu)建HMM模型，搜索更多的同源序列；模型，搜索更多的同源序列；r 構(gòu)建進(jìn)化的樹的必須步驟；構(gòu)建進(jìn)化的樹的必須步驟；r 比較基因組學(xué)研究；比較基因組學(xué)研究；r 兩類：全局或局部的多序列比對(duì)；兩類：全局或

13、局部的多序列比對(duì)；Made by GENEDOCGapVDSCYGap0-11-22-33-44-55V-114-7-18-29-40E-22-76-5-16-27S-33-18-510-1-12L-44-29-16-19-3C-55-40-27-1287Y-66-51-38-23-31542時(shí)間復(fù)雜度：時(shí)間復(fù)雜度：O(n2)三條序列：時(shí)間復(fù)雜度：三條序列：時(shí)間復(fù)雜度：O(lmn) = O(n3)四條序列：時(shí)間復(fù)雜度：四條序列：時(shí)間復(fù)雜度：O(n4)，非多項(xiàng)式時(shí)間！，非多項(xiàng)式時(shí)間！多項(xiàng)式時(shí)間復(fù)雜度要求：O(n3)m條序列：時(shí)間復(fù)雜度：條序列：時(shí)間復(fù)雜度：O(nm)，NPC問題問題！Sequen

14、ce ASequence BSequence C 搜索有限空間，類似于搜索有限空間，類似于BLAST算法算法r 最優(yōu)的多序列比對(duì)，其兩兩序列之間的比對(duì)最優(yōu)的多序列比對(duì)，其兩兩序列之間的比對(duì)不一定最優(yōu)。不一定最優(yōu)。 最優(yōu)的多序列比對(duì)最優(yōu)的多序列比對(duì)非最優(yōu)的雙序列比對(duì)非最優(yōu)的雙序列比對(duì)rMSA - Multiple Sequence AlignmentrDavid Lipman等，等，1989年初始開發(fā)；年初始開發(fā)；r應(yīng)用多維動(dòng)態(tài)規(guī)劃算法，得到最優(yōu)的全局應(yīng)用多維動(dòng)態(tài)規(guī)劃算法，得到最優(yōu)的全局比對(duì)。比對(duì)。r工具資源：工具資源：/CBBresearch

15、/Schaffer/msa.html/general/software/packages/msa/manual/manual.phpr1. 漸進(jìn)方法：漸進(jìn)方法：progressive methods代表：代表：ClustalW/X, T-Coffeer2. 迭代方法：迭代方法：iterative methods 代表代表: PRRP, DIALIGNr3. 部分有向圖算法：部分有向圖算法：Partial Order Algorithm (POA)r4. 全局多序列比對(duì)的隱馬爾科夫模型全局多序列比對(duì)的隱馬爾科夫模型profile HMMr5. 整合算法：整合算法

16、： MUSCLEr(1) ClustalW/Xr(2) T-Coffeer1. Clustal: 1988年開發(fā)；年開發(fā)；r2. ClustalW: 1994年，年，Julie D. Thompson等人改進(jìn)、發(fā)展；等人改進(jìn)、發(fā)展；r3. ClustalX: 1997年，圖形化軟件；年，圖形化軟件；r1. 將所有序列兩兩比對(duì)，計(jì)算距離矩陣；將所有序列兩兩比對(duì)，計(jì)算距離矩陣；r2. 構(gòu)建鄰接進(jìn)化樹構(gòu)建鄰接進(jìn)化樹(neighbor-joining tree)/指導(dǎo)樹指導(dǎo)樹(guide tree)；r3. 將距離將距離最近最近的兩條序列用動(dòng)態(tài)規(guī)劃的算法的兩條序列用動(dòng)態(tài)規(guī)劃的算法進(jìn)行比對(duì)；進(jìn)行比對(duì)；r

17、4. “漸進(jìn)漸進(jìn)”的加上其他的序列。的加上其他的序列。兩兩比對(duì)，構(gòu)建距離矩陣指導(dǎo)樹的構(gòu)建指導(dǎo)樹的構(gòu)建漸進(jìn)比對(duì)漸進(jìn)比對(duì)每條序列的權(quán)值每條序列的權(quán)值Score:BLOSUM62的分?jǐn)?shù)的分?jǐn)?shù)r1. FASTA序列格式，多序列：序列格式，多序列：r BioEdit, GeneDoc等軟件等軟件GeneDocGeneDoc軟件，導(dǎo)入軟件，導(dǎo)入.aln.aln文件文件r1. 采用采用Clustal程序計(jì)算兩兩序列之間的全程序計(jì)算兩兩序列之間的全局最優(yōu)比對(duì)結(jié)果；局最優(yōu)比對(duì)結(jié)果；r2. 采用采用LALIGN程序計(jì)算兩兩序列之間的局程序計(jì)算兩兩序列之間的局部最優(yōu)比對(duì)的結(jié)果；部最優(yōu)比對(duì)的結(jié)果；r3. 設(shè)計(jì)加權(quán)系統(tǒng)

18、，綜合考慮以上兩類結(jié)果的設(shè)計(jì)加權(quán)系統(tǒng)，綜合考慮以上兩類結(jié)果的因素，構(gòu)建指導(dǎo)庫；因素，構(gòu)建指導(dǎo)庫；r4. 最后，采用漸進(jìn)式比對(duì)算法，得到最終的最后，采用漸進(jìn)式比對(duì)算法，得到最終的結(jié)果。結(jié)果。同時(shí)進(jìn)行全局和局部的同時(shí)進(jìn)行全局和局部的雙序列比對(duì)雙序列比對(duì)對(duì)以上打分的結(jié)果設(shè)計(jì)對(duì)以上打分的結(jié)果設(shè)計(jì)權(quán)重系統(tǒng)，找到序列中權(quán)重系統(tǒng)，找到序列中最保守的部分最保守的部分漸進(jìn)方法的比對(duì)，基于上述漸進(jìn)方法的比對(duì)，基于上述計(jì)算的計(jì)算的primary libraryr1. 距離最近的，有兩組序列距離最近的，有兩組序列AB和和CD，哪組，哪組最先比對(duì)？兩種方案：最先比對(duì)？兩種方案：A. 分別、同時(shí)比對(duì)。但是，是以分別、同時(shí)

19、比對(duì)。但是，是以AB為準(zhǔn)，加入為準(zhǔn)，加入CD，然后再加上其他序列，還是，然后再加上其他序列，還是CD為準(zhǔn)？結(jié)為準(zhǔn)？結(jié)果可能出入很大果可能出入很大B. 隨機(jī)挑選一組作為基準(zhǔn)隨機(jī)挑選一組作為基準(zhǔn)r2. 當(dāng)序列差異較大時(shí)，上述問題更加明顯。當(dāng)序列差異較大時(shí)，上述問題更加明顯。r1. 三條序列：三條序列：r2.若若Seq1,2先比對(duì)先比對(duì)，再加入，再加入Seq3：r3. Seq1,3先比對(duì)，先比對(duì)，再加入再加入Seq2:r4. Seq2,3先比對(duì)，先比對(duì)，再加入再加入Seq1:Seq1: ARKCVSeq2: ARCVSeq3: AKCVARKCVAR-CVA-KCVARKCVA-RCVA-KCVAR

20、KCVAR-CVAK-CVr1. 部分解決漸進(jìn)算法存在的問題部分解決漸進(jìn)算法存在的問題,主要是主要是ClustalW/X存在的問題；存在的問題；r2. PRRPr3. DIALIGN1. 1. 先用先用“漸進(jìn)漸進(jìn)”算法進(jìn)行算法進(jìn)行多序列比對(duì)多序列比對(duì); ;2. 2. 基于多序列比對(duì)的結(jié)果基于多序列比對(duì)的結(jié)果構(gòu)建進(jìn)化樹；構(gòu)建進(jìn)化樹；3. 3. 重新計(jì)算序列之間的距重新計(jì)算序列之間的距離，再用離，再用“漸進(jìn)漸進(jìn)”算法進(jìn)行算法進(jìn)行多序列比對(duì)；多序列比對(duì)；4. 4. 重復(fù)上述步驟，直到結(jié)重復(fù)上述步驟，直到結(jié)果不再發(fā)生改變?yōu)橹埂９辉侔l(fā)生改變?yōu)橹埂1. 對(duì)所有序列進(jìn)行兩兩之間的局部最優(yōu)化的對(duì)所有序列進(jìn)

21、行兩兩之間的局部最優(yōu)化的比對(duì)；比對(duì)；r2. 找到所有能夠匹配的部分找到所有能夠匹配的部分M1；將重疊的；將重疊的、前后連續(xù)、前后連續(xù)(consistency)的匹配部分連接的匹配部分連接起來起來(diagonals)，為，為M2；r3. 將剩下的未比對(duì)的序列重新比對(duì)，再發(fā)現(xiàn)將剩下的未比對(duì)的序列重新比對(duì)，再發(fā)現(xiàn)能夠匹配的部分，構(gòu)成新能夠匹配的部分，構(gòu)成新M1，將，將consistency部分構(gòu)成部分構(gòu)成M2；r4. 重復(fù)上述步驟，直到結(jié)果收斂。重復(fù)上述步驟，直到結(jié)果收斂。r主要改進(jìn)：主要改進(jìn)：1. 所有序列的兩兩比對(duì)，通過所有序列的兩兩比對(duì)，通過profile HMM的的方法進(jìn)行雙序列比對(duì)；方法進(jìn)行雙序列比對(duì)；2. 將漸進(jìn)算法與迭代算法整合；將漸進(jìn)算法與迭代算法整合；3. 目前，性能最優(yōu)。目前，性能最優(yōu)。r算法分為三個(gè)部分，每個(gè)部分相對(duì)獨(dú)立；算法分為三個(gè)部分，每個(gè)部分相對(duì)獨(dú)立；r1. Draft progressive: (1) 對(duì)兩條序列，計(jì)算距離采用對(duì)兩條序列，計(jì)算距離采用k-mer的思想；的思想；(2) 用用UPGMA算法構(gòu)建引導(dǎo)樹；算法構(gòu)建引導(dǎo)樹；(3) 使用漸進(jìn)算法進(jìn)行多序列比對(duì)；使用漸進(jìn)算法進(jìn)行多序列比對(duì)；r優(yōu)點(diǎn)：兩條序列之間的距離不采用動(dòng)態(tài)規(guī)劃優(yōu)點(diǎn)：兩條序列之間的距離不采用動(dòng)態(tài)規(guī)劃算法進(jìn)行比對(duì)，節(jié)省