Luận văn Các kỹthuật toán học cho bài toán so sánh đa trình tự

max{S(i − 1, j − 1) + σ(S1[i] , S2[j] ),D(i, j ), I (i, j) } for i > 0, j > 0. D(0, j ) = S(0, j ) − q for j ≥ 0, D(i, 0) = D(i − 1, 0) − r for i > 0, D(i, j ) = max{D(i − 1, j ) − r, S(i − 1, j ) − q − r }for i > 0, j > 0. I (i, 0) = S(i, 0) − q for i ≥ 0, I (0, j ) = I (0, j − 1) − r for j > 0, I (i, j ) = max{I (i, j − 1) − r, S(i, j − 1) − q − r }for i > 0, j > 0. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 29 Phần tiếp theo xin đề cập về cách tính giá trị của phép chuyển đổi từ chuỗi S1 sang S2 theo các khái niệm được giới thiệu ở trên. Chúng ta sẽ sử dụng 3 mảng 2 chiều để lưu trữ kết quả cho S, D,I. Gotoh’s algorithm arrays S[0..m, 0..n] , D[0..m, 0..n] , I[0..m, 0..n] S(0, 0)← 0 for j ← 1 to n do { S(0, j) ← −( q+r*j) D(0, j) ← S(0,j) - q } for i ← 1 to m do { S(i, 0) ← −(q + r*i) I(i, 0) ← S(i,0)−q for j ← 1 to n do { I(i, j) ← max {I(i, j - 1) , S(i, j - 1) −q} - r D(i, j) ← max {D(i - 1, j) , S(i - 1, j) - q} - r S(i, j) ← max {D(i, j) , I(i, j) , S(i − 1, j − 1) + 1 2( [ ], [ ])S i S jσ } } } write "cost is" S(m, n) 3.2.3. Giải thuật cải tiến không gian nhớ Từ giải thuật Gotoh ta có nhận xét sau: ƒ Các giá trị của các phần tử trong dòng thứ i của ma trận S và D chỉ phụ thuộc vào giá trị của các phần tử trong dòng thứ i-1 và dòng thứ i. ƒ Giá trị của các phần tử của dòng thứ i của ma trận I chỉ phụ thuộc vào các phần tử trước đó cũng của dòng này. Từ 2 nhận xét trên ta thấy rằng để tính giá trị S(m,n) chúng ta không cần phải lưu trữ toàn bộ 3 ma trận như giải thuật Gotoh. Ta có thể giảm không gian lưu trữ khi thực thi giải thuật xuống còn (n)θ thông qua, kỹ thuật đổi chỗ(swap). Bằng kỹ thuật này ta chỉ cần một không gian nhớ n phần tử, lưu trữ giá trị các phần tử của dòng thứ i-1 khi bắt đầu tính giá trị các phần tử của dòng thứ i. Giá trị tính được của dòng thứ i sẽ lại được ghi đè lên không gian nhớ này. Như vậy lúc nào ta cũng chỉ cần (n)θ không gian Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 30 nhớ. Tuy nhiên giải thuật chỉ tính giá trị của S(m,n), không lưu vết và tìm kết quả. Việc tìm kết quả alignment dựa trên giải thuật này sẽ được trình bày trong chương 4. Gọi SS và DD là 2 mảng có chiều dài n sẽ dùng để lưu trữ giá trị của các phần tử thuộc dòng thứ i của các ma trận S và D, e là giá trị của I . Sau vòng lặp thứ i, vector SS sẽ lưu trữ giá trị của hàng thứ i của ma trận S. Tương tự như vậy sau vòng lặp thứ i của giải thuật, DD sẽ lưu trữ giá trị của hàng thứ i của ma trận D. Giải thuật: vectors SS[0..n] , DD[0..n] Linear(S1,S2,SS,DD,-q) Procedure Linear(S1,S2,SS,DD,p) Var e, c, s, t SS(0) ← 0 for j ← 1 to n do/*Khởi tạo*/ { SS( j) ← −(q + r*j) DD( j) ← SS(j) − q } t←p for i ← 1 to m do/*Tính SS, DD*/ { s ← SS(0) SS(0) ← c← t← t−r /∗ −(q + r*i)*/ e← SS(0) − q for j← 1 to n do { e← max {e, c −q} − r DD( j) ← max {DD( j) , SS( j) − q} − r c← max {DD( j) , e, s + 1 2( [ ], [ ])S i S jσ } s ← SS( j) SS( j) ← c } } write "cost is" SS(n) Như vậy sau khi kết thúc giải thuật SS sẽ lưu trữ giá trị của hàng thứ i của ma trận S. Khi đó SS(n)=S(m,n). Tại mỗi thời điểm ứng với giá trị của i>0, j>0 giá trị của các thông số SS(k), DD(k), s, c, e (k>0)sẽ được hiểu như sau: SS(k) = ⎪⎩ ⎪⎨ ⎧ < ≥− jkkiS jkkiS ),( ),1( Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 31 DD(k) = ⎪⎩ ⎪⎨ ⎧ < ≥− jkkiD jkkiD ),( ),1( e = I(i, j −1) c = S(i, j − 1) s = S(i − 1, j − 1) Ví dụ: Cho ),( baσ = -1 nếu a≠ b, ),( aaσ = 0, and kk 5.02)( −−=γ . Hai chuỗi S1, S2 lần lượt là AGTAC và AAG. Quá trình chuyển đổi từ chuỗi S1 sang S2 bao gồm các bước delete GT, thay thế C bằng G, và giá trị của phép chuyển này là -4. Áp dụng giải thuật Gotoh ta sẽ có các bảng như sau.‘*’ là các giá trị không xác định. Hình 3.2 Các ma trận S, D, I cho 2 chuỗi AGTAC and AAG. Với giải thuật cải tiến, ta chỉ cần sử dụng 2 vector có 4 phần tử S và D và các giá trị e,c, s tương ứng với I(i, j − 1), S(i, j − 1), và S(i − 1, j − 1). Hình trên biểu diễn giá trị của SS, DD, e, c, s khi i= 5 và j = 2, giá trị của SS và DD được biểu diễn trong ô chữ nhật, giá trị của e, c, và s được biểu diễn trong các ô hình elip. A G T A C A − − Α G Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 32 3.3. Giải thuật tính toán phép Alignment tối ưu cho bài toán Multiple Alignment sử dụng kỹ thuật dynamic programming Phương pháp dynamic programming truyền thống tiếp cận bài toán MSA bằng cách xây dựng hàm đánh giá theo phương pháp Sum-of-Pair(sử dụng tổng các kết quả của phép alignment từng cặp trình tự(pairwise alignment)). Với phương pháp này thay vì sử dụng hàm σ ta sẽ định nghĩa một hàm δ (x,y) để đo lường khoảng cách giữa 2 ký tự x, y. Sự khác nhau càng nhiều giữa 2 chuỗi trình tự thì δ càng lớn. Như vậy với từng cặp trình tự S, T mục tiêu của phương pháp này là sẽ tìm giá trị nhỏ nhất của biểu thức: d*(S,T)=min d(S,T)=min(∑ = l i iTiS 1 ])['],['(δ ) với l=|S’|=|T’| Định nghĩa 3.3: Giá trị của một MSA của k chuỗi S1, S2,…,Sk theo cách tiếp cận sử dụng tổng kết quả của các phép alignment từng cặp trình tự(Sum of pair-SP) là tổng của tất cả ( k2 ) phép alignment từng cặp trình tự từ k chuỗi này. Ví dụ: Xét một MSA của 3 chuỗi trình tự: (S1) A C - C D B - (S2) - C - A D B D (S3) A - B C D A D Nếu chúng ta định nghĩa ),( yxδ như sau: 0, ( , ) 1, x y x y x y δ =⎧= ⎨ ≠⎩ Giá trị của MSA trên sẽ là d(S1,S2)+d(S1,S3)+d(S2,S3)=3+4+5=12 Định nghĩa 3.4: Alignment tối ưu của k chuỗi S1, S2,…,Sk theo phương pháp sử dụng tổng kết quả của các phép alignment từng cặp trình tự, là một alignment có giá trị của tổng kết quả của các phép alignment từng cặp trình tự, mà giá trị này là nhỏ nhất. (minimize biểu thức * 1 1 ( ) ( , ) = = + = ∑ ∑k k i j i j i f M d S S ) Xét bài toán MSA với k chuỗi trình tự S1,S2,…,Sk có chiều dài |S1|,|S2|,…,|Sk|. Gọi min 1min{| |,...,| |)kn S S= , max 1max{| |,...,| |)kn S S= Rõ ràng theo cách tiếp cận sử dụng cho bài toán với 2 chuỗi trình tự thì ta sẽ cần một ma trận k chiều: 1 2(| |)(| |)...(| |)k k S S S  . Ma trận này có 1 2(| |)(| |)...(| |)k k S S S  phần tử. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 33 Và mỗi phần tử sẽ phụ thuộc vào 2k -1 phần tử xung quanh nó. Như vậy quá trình tìm phép alignment tối ưu như đã tiếp cận với bài toán PSA sẽ cho độ phức tạp của giải thuật là 1 2 min min((| |)(| |)...(| |)) (2 ) ( ). (2 ) ((2 ) )k k k kkS S S n nΟ Ο ≥Ο Ο = Ο . Rõ ràng với nmin lớn khoảng vài trăm phần tử và k≥3 thì giải thuật là không hiệu quả về mặt thời gian. Hình 3.3 Minh hoạ quá trình tìm 1 MSA tối ưu Để cải thiện thời gian tính toán (cũng là giảm độ phức tạp của giải thuật ), ta có thể thay đổi cách tiếp cận, thay vì tìm lời giải tối ưu ta có thể chỉ cần tìm lời giải gần tối ưu nhưng có độ phức tạp nhỏ hơn nhiều. Điều này tương đương với việc tìm ra một phép alignment có giá trị ∑ ∑k k i j i=1 j=i+1 d (S ,S ) xấp xỉ với giá trị của phép alignment tối ưu ∑ ∑k k * i j i=1 j=i+1 d (S ,S ) , và quá trình tìm ra phép alignment có độ phức tạp chấp nhận được. Một trong những giải pháp được đưa ra đó là sử dụng heuristic dựa trên phương pháp quy hoạch động truyền thống. Phần tiếp theo xin được giới thiệu về giải thuật Center Star Alignment, một giải thuật heuristic kinh điển hiện thực theo hướng này. 3.3.1. Giải thuật Center Star Alignment Algorithm [30] Giải thuật được xây dựng dựa trên ý tưởng về sự tiến hóa trong tự nhiên của các chuỗi trình tự đại diện cho các loài. Quá trình tiến hóa trong tự nhiên của các loài có thể được mô tả theo một đồ thị. Giải thuật Center Star Alignment dựa trên một trong những hình thức của đồ thị để mô tả quá trình alignment. Hình thức đồ thị ở đây là đồ thị có kiến trúc hình sao. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 34 Hình 3.4 Mô hình tiến hoá hình sao Nội dung cơ bản của giải thuật là chọn ra một phần tử trung tâm và thực hiện quá trình align các phần tử còn lại với phần tử này để tìm ra một kết quả MSA. Giải thuật có độ phức tạp nhỏ hơn 2 2max( )n kΟ . Ta đưa ra một số giả thuyết (assumption) về hàm đánh giá khoảng cách như sau(làm bổ đề cho quá trình đánh giá sai số): 1. 0),( =xxδ 2.Bất đẳng thức tam giác: ),(),(),( zyyxzx δδδ +≤ Định nghĩa 3.5: Xét 2 chuỗi S và T. Gọi D(S,T) là giá trị hàm đánh giá khoảng cách của phép alignment giữa S và T(theo hướng toàn cục) mà hàm đánh giá khoảng cách của nó là nhỏ nhất. Định nghĩa 3.6: K là tập lưu các thể hiện của các trình tự sau mỗi bước của phép multiple alignment. Sau mỗi bước của quá trình alignment, K sẽ chứa các thể hiện mới của các chuỗi trình tự ban đầu. Giải thuật: Xét tập I gồm k chuỗi trình tự. Tìm chuỗi SC là một chuỗi trong Ι sao cho: ∑ −∈ }{ ),( CSIS C SSD đạt giá trị nhỏ nhất(Điều này có thể thực hiện dựa trên việc áp dụng giải thuật giải quyết bài toán PSA cho )(2k cặp trình tự của I) Đưa SC vào K: K }{ CS← Xét k-1 trình tự còn lại S1, …, SC-1, SC+1,…, Sk, lần lượt thêm vào trong tập K, thông qua phương thức Add(Si,K) . … …………………. …………………. …………………. …………………. …………………. …………………. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 35 Giải thuật kết thúc khi tất cả các sequence được đưa vào K. Ta sẽ định nghĩa phương thức Add(Si,K). Add(Si,K) có thể được giải thích như sau: Giả sử trong K có chứa các chuỗi trình tự S’C, S’1,..,S’i-1 là các thể hiện của SC, S1,..,Si-1,. Phương thức Add(Si, K) sẽ thực hiện tìm kết quả alignment của S’C và Si. Kết quả thu về sẽ là 2 chuỗi trình tự S’’C và S’i. Thay thế S’C bằng S’’C, thêm S’i vào K. Điều chỉnh lại các trình tự S’1,…S’i-1 theo sự thay đổi của S’C →S’’C : Nếu cột nào của S’’C xuất hiện gap so với S’C, thì cột tương ứng của S’1,…,S’i-1 cũng sẽ được thêm gap (ký tự ‘-‘). Giải thuật sẽ cho kết quả gần đúng so với lời giải tối ưu có thể có của phép toán multiple alignment. Độ phức tạp của giải thuật: Ta có thể chứng minh giải thuật có độ phức tạp nhỏ hơn 2 2max( )n kΟ . Thật vậy giải thuật cho mỗi bài toán PSA có độ phức tạp tối đa là 2 max( )nΟ , do có )(2k cặp trình tự từ I nên việc tìm SC đòi hỏi độ phức tạp nhỏ hơn 2 max 2( ( )) knΟ = 2 2max( )n kΟ . Mỗi thao tác thêm Si vào K sẽ tạo ra một thể hiện mới S’’C của SC trong K, chuỗi S’’C có chiều dài tối đa là ' ' max| | *CS n i≤ . Mỗi lần thêm một phần tử mới Si+1 sẽ đòi hỏi phải thực hiện một lần giải thuật cho bài toán PSA với chuỗi Si+1 ( 1 max| |iS n+ ≤ ), S’’C( ' ' max| | *CS n i≤ ). Do đó độ phức tạp cho giải thuật Add(Si,K) là '' 1 (| |* | |) k C i i i C S S =≠ Ο∑ thỏa: 1'' 2 2max max max 1 1 (| |* | |) ( * * ) ( ) k k C i i i i C S S i n n n k − = =≠ Ο ≤ Ο = Ο∑ ∑ Vậy độ phức tạp của giải thuật sẽ là 2 2 2 2 2 2max max max( ) ( ) ( )T n k n k n k≤ Ο + Ο = Ο Hình 3.5 Minh họa Center Star Algorithm s2 s1 s3 s4 s1: MPE s2: MKE s3: MSKE s4: SKE MPE | | MKE MSKE - || MKE SKE || MKE MPE MKE -MPE -MKE MSKE -MPE -MKE MSKE -SKE Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 36 Đánh giá: Gọi MC là một kết quả của giải thuật. Gọi f(MC) là giá trị hàm đánh giá của phép alignment này. Và M* là một phép alignment tối ưu của k chuỗi. f(M*) là giá trị hàm đánh giá của phép alignmnet tối ưu. Ta có: ( ) 2( 1) 2 ( *) Cf M k f M k −≤ < . Thật vậy ta chứng minh điều này như sau. Gọi V(MC)= ∑ ∑ = ≠= k i k ji j ji SSd 1 1 ),( và V(M*)= ∑ ∑ = ≠= k i k ji j ji SSd 1 1 ),(* Rõ ràng V(MC)=2f(MC) , V(M*)=2f(M*) Ta có V(MC)= ∑ ∑ = ≠= k i k ji j ji SSd 1 1 ),( ≤ ∑∑ = ≠= +k i k ji j jCCi SSdSSd 1 1 )],(),([ 1 1 [ ( , ) ( , )] k k i C C j i j i D S S D S S j = = ≠ ≤ +∑∑ (do MC là 1 lời giải của thuật toán nên d(Si, SC)=D(Si,SC), d(SC,Sj)=D(SC,Sj)) Dễ thấy ∑ ∑ = ≠= +k i k ji j jCCi SSDSSD 1 1 )],(),([ = ∑ ≠= − k Cj j Cj SSDk 1 ),()1(2 Suy ra V(MC)≤ ∑ ≠= − k Cj j Cj SSDk 1 ),()1(2 (1) V(M*)= ∑ ∑ = ≠= k i k ji j ji SSd 1 1 ),(* ≥ ∑ ∑ = ≠= k i k ji j ji SSD 1 1 ),( = ∑ ≠= k Cj j Cj SSDk 1 ),( (2) (1)(2) suy ra 2)1(2 *)( )( *)( )( <−≤= k k MV MV Mf Mf cc Như vậy với tỉ lệ giữa kết quả hàm đánh giá của giải thuật với hàm đánh giá của phép alignment tối ưu lúc nào cũng thuộc [1,2). Giải thuật sẽ cho độ chính xác cao khi chuỗi trình tự đưa vào có quá trình tiến hoá trong tự nhiên tuân theo mô hình đồ thị hình sao. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 37 Center Star Alignment là một trường hợp đặc biệt của phương pháp Progressive Algorithm được trình bày trong chương 2. Trong phần tiếp theo ta sẽ đề cập chi tiết về phương pháp đã được giới thiệu này. 3.3.2. Phương pháp Progressive Algorithm giải quyết bài toán MSA. Ý tưởng của Progressive Algorithm là align các sequence dựa trên việc áp dụng bài toán PSA lặp đi lặp lại nhiều lần. Progressive Algorithm có nhiều phiên bản hiện thực khác nhau, tuy nhiên về cơ bản một Progressive Alogrithm gồm có 3 bước: 1. Chọn ra 2 sequence từ tập các sequence, align 2 sequence này với nhau. 2. Chọn một sequence khác và align với nhóm sequence đang align. 3. Lặp lại bước 2 cho đến khi tất cả các sequence đã được align. Hình 3.6 Hình minh hoạ cho Progressive Algorithm. Từ ý tưởng ban đầu là align mỗi sequence với 1 nhóm sequence, một sự mở rộng của giải thuật là phân hoạch tập các sequence thành các nhóm sequence, sau đó tiến hành align các nhóm sequence này cho đến khi tất cả các nhóm sequence được align thành 1 nhóm chung duy nhất. Về bản chất, việc align 1 sequence với 1 nhóm các sequence cũng chính là việc align 2 nhóm sequence với nhau mà trong đó 1 nhóm chỉ có 1 phần tử. Điểm khác biệt chính trong các phiên bản hiện thực của Progressive Algorithm là phép phân hoạch tập sequence, cũng như trật tự để align các nhóm sequence với nhau. Kỹ thuật cơ bản thường được sử dụng trong việc align 2 nhóm sequence là giải thuật quy hoạch động. Phiên bản Progressive Algorithm cơ bản nhất, … Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 38 thường được thừa kế, cải tiến sử dụng trong các phiên bản Progressive Algorithm sau này là giải thuật Feng-Doolittle. 3.3.3. Feng-Doolittle Algorithm Feng-Doolittle Algorithm [10], [11] triển khai ý tưởng phân hoạch tập sequence ban đầu thành các nhóm sequence, sau đó thực hiện việc align các nhóm sequence với nhau cho đến khi chỉ còn 1 nhóm duy nhất. Để áp dụng bài toán PSA cho các nhóm sequence, giải thuật xem mỗi nhóm sequence như 1 sequence. Sau đó thực hiện việc tính toán giá trị hàm đánh giá cho các cặp vị trí của 2 nhóm. Từ ma trận hàm đánh giá các cặp vị trí của 2 nhóm sequence, kết quả alignment của 2 nhóm sequence sẽ được xây dựng dựa trên quá trình thêm gap vào chúng. Xét 2 nhóm sequence A có k1 sequence chiều dài n1, B có k2 sequence chiều dài n2 . Công thức tính toán hàm đánh giá cột thứ i của nhóm A và cột thứ j của nhóm B: 1 2 , , 1 1 1( , ) ( , ) 1* 2 k k x i y j x y S i j a b k k σ = = = ∑∑ Trong đó 11 i n≤ ≤ là cột thứ i của nhóm sequence A 21 j n≤ ≤ là cột thứ j của nhóm sequence B. ,x ia : ký tự đại diện cho amino acid tại dòng x, cột i của nhóm A. ,y jb : ký tự đại diện cho amino acid tại dòng y, cột j của nhóm B. Ví dụ chúng ta xem xét việc align 2 nhóm A, B như sau: Giá trị hàm đánh giá cho vị trí cột thứ 4 của A và vị trí cột thứ 5 của B là: 1(4,5) ( ( , ) ( , ) ( , ) ( , ) ( , ) ( , )) 2*3 S R V R X K V K X G V G Xσ σ σ σ σ σ= + + + + + Quá trình thêm gap vào 2 nhóm sequence được thực hiện theo nguyên tắc: nếu xuất hiện gap tại 1 vị trí của nhóm sequence, thì tất cả các sequence của nhóm sẽ phải thêm gap vào vị trí này DGM R EA DGM K GA DNM G GS EAM SVGS ENM I XGC Nhóm A Nhóm B Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 39 Hình 3.7 Minh họa Feng-Doolittle Algorithm Hình 3.8 Ví dụ thực thi Feng-Doolittle Algorithm Quá trình align các nhóm sequence đã được phân hoạch theo một trình tự nào đó, sẽ tạo nên một cây nhị phân. MSA tối ưu sẽ tương ứng với 1 cây tối ưu. Việc tìm kiếm MSA tối ưu của Progressive Algorithm sẽ tương đương với việc tìm ra cây tối ưu. Khoảng cách giữa các sequence: Trong phần nói về Center Star Algorithm chúng ta có đề cập đến khái niệm khoảng cách giữa 2 sequence, ở mục này sẽ trình bày chi tiết hơn về khái niệm này, và cũng giải thích lý do khái niệm này thường được sử dụng trong phương pháp Progressive Algorithm. Giải thuật quy hoạch động đề cập đến MSA tối ưu là MSA có giá trị hàm đánh giá lớn nhất(dựa trên các ma trận hàm đánh giá), giá trị hàm đánh giá càng lớn nghĩa là các sequence tham gia có độ tương đồng càng cao, đồng nghĩa với việc là về mặt tiến hóa khoảng cách di truyền giữa các sequence là ngắn. Trong thực tế sinh học, người ta thường sử dụng khái niệm khoảng cách giữa các sequence để mô tả quá trình tiến hóa của các sequence, cũng như biểu hiện khả năng tương đồng của chúng. Vì phù hợp với bản chất phát triển tự nhiên của Thêm gap A B Align 2 nhóm A ,B DGMREA DGMKGA DNMGGS EAMSGVS ENMI GXC DGMR E - A DGMK G - A DNMG G - S EAMS G V S ENMI G X C Thêm gap vào vị trí thứ 6 của nhóm A Nhóm A Nhóm B Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 40 sinh học, cách tiếp cận hàm đánh giá theo phương pháp khoảng cách di truyền được sử dụng nhiều trong các bài toán sinh học. Định nghĩa 3.7 Khoảng cách giữa 2 sequence(Pairwise Distance), được tính bằng phương pháp Sum-of-Pair, là một hàm theo giá trị hàm đánh giá SPScore của 2 chuỗi này. ( , ) ( ( , ))i j i jDist S S f SPScore S S= với :f →\ \ Trong phần này chúng ta sẽ xem xét các giải thuật heuristic dựa trên Progressive Algorithm: CLUSTALW, MULTIALIGN,... Các thuật giải heuristic hầu hết đều phát triển dựa trên nền tảng của Feng-Doolittle Algorithm, trong đó phần lớn tập trung vào việc xây dựng cách thức phân hoạch các nhóm sequence, cũng như mô tả trình tự align của các nhóm này. Các thuật giải này về cơ bản đều gồm những bước chính sau: 1. Tính toán khoảng cách giữa các cặp sequence(có 2kC khoảng cách) thông qua kỹ thuật quy hoạch động 2. Xây dựng cây nhị phân mô tả trật tự align của các sequence dựa trên các cặp khoảng cách giữa các sequence, sao cho cây này là tối ưu 3. Align các sequence dựa trên cây mô tả. Một số chương trình trong thực tế xây dựng cây mô tả để thực hiện quá trình align MULTALIGN(Barton và Sternberg,1987) và PILEUP(Genetic Computer Group) sử dụng giải thuật UPGMA-Unweighted Pair-Group Medthod using arithmetic Averages(Sneath và Sokakk, 1973) để xây dựng cây mô tả. Chương trình sẽ chọn ra 2 sequence có khoảng cách ngắn nhất trong tất cả các cặp sequence, align chúng với nhau tạo thành một nhóm sequence. Tính lại tập khoảng cách của n-2 phần tử còn lại và nhóm sequence vừa được tạo. Tiếp tục tìm ra 1 cặp phần tử có độ tương đồng cao nhất và align chúng với nhau. Chương trình kết thúc khi tất cả các phần tử đã được align. Việc tính toán khoảng cách giữa nhóm sequence và 1 sequence giống như quá trình tính toán khoảng cách giữa 2 nhóm sequence của Feng-Doolittle Algorithm, trong đó một nhóm được rút gọn thành 1 sequence. CLUSTALW(Thompson 1994)[9]: Sử dụng giải thuật Neighbor-Joinning(Saitou và Nei,1987) để xây dựng cây mô tả. Đầu tiên các cây không gốc sẽ được tạo thông qua việc tính toán dựa vào khoảng cách của các cặp sequence, mỗi nút lá đại diện cho Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 41 một sequence. Từ các cây không gốc tính toán để xây dựng một nút gốc nối liền các cây này. Các nhánh cây có trọng số khác nhau, các trọng số này được tính toán và dẫn xuất từ khoảng cách của các cặp sequence ban đầu. CLUSTALW là một chương trình được biết đến và sử dụng nhiều nhất trong các chương trình giải quyết bài toán MSA. Đây cũng là chương trình mà kết quả cho độ chính xác tốt nhất so với các chương trình được hiện thực theo phương pháp Progressive Algorithm. Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 42 Chương 4. THIẾT KẾ GIẢI THUẬT VÀ HIỆN THỰC PHƯƠNG PHÁP GIẢI QUYẾT BÀI TOÁN MSA Trên cơ sở lý thuyết đã trình bày về các phương pháp giải quyết bài toán MSA tại thời điểm hiện tại. Trong phần này xin được trình bày giải pháp tiếp cận cũng như việc hiện thực chương trình để giải quyết bài toán MSA. Nội dung bao gồm: ƒ Thiết kế giải thuật giải quyết bài toán PSA. ƒ Thiết kế phương pháp phân hoạch, giải pháp gom nhóm cũng như trình tự align các sequence. ƒ Xây dựng giải thuật mới gọi là giải thuật 1A giải quyết bài toán MSA, và giải thuật cải tiến của 1A, giải thuật 1B. ƒ Phương pháp hiện thực các giải thuật. 4.1. Giải thuật sử dụng cho bài toán PSA Như đã trình bày trong phần cơ sở lý thuyết, bài toán PSA có thể được hiện thực một cách dễ dàng bằng phương pháp quy hoạch động, tuy nhiên phương pháp này có một số hạn chế : ƒ Việc xây dựng ma trận để phục vụ cho việc tính toán đòi hỏi vấn đề cấp phát bộ nhớ cũng như khối lượng tính toán lớn. Khi chiều dài của các sequence tham gia đủ lớn(n>1000), đây sẽ là một vấn đề trong việc hiện thực. ƒ Vấn đề tìm lại con đường thông qua các vết được lưu. ƒ Độ chính xác của phép alignment phụ thuộc rất nhiều vào ma trận đánh giá, tùy vào tính tương đồng mà mỗi cặp sequence sẽ phù hợp nhất với một ma trận đánh giá. Giải thuật quy hoạch động thuần túy chỉ sử dụng một ma trận đánh giá, do đó không cho phép thay đổi ma trận đánh giá phù hợp với cặp sequence đang xét, điều này hạn chế khả năng linh động của giải thuật. Trong phạm vi của luận văn này tôi xin trình bày một giải pháp mới để giải quyết bài toán PSA có thể hạn chế được những khuyết điểm của phương pháp quy hoạch Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 43 động thuần túy, cho phép tối ưu hóa vấn đề không gian nhớ sử dụng, thay đổi linh hoạt các ma trận đánh giá, nâng cao độ chính xác của giải thuật. Giải thuật quy hoạch động có độ phức tạp ( )nmΟ và đòi hỏi không gian bộ nhớ ( )Ο nm . Giải thuật được xây dựng chỉ đòi hỏi không gian bộ nhớ ( lg )Ο + +n m m và thời gian chạy vẫn là ( )nmΟ . 4.1.1. Giải thuật tính toán dựa theo kỹ thuật chia để trị Trong chương 2 chúng ta đã đề cập đến phương pháp cải tiến giải thuật Gotoh cho phép tính toán S(m,n) với độ phức tạp ( )nmΟ và không gian nhớ sử dụng là ( )Ο +m n . Trong phần này chúng ta sẽ phát triển giải thuật cho phép tìm ra kết quả alignment dựa trên giải thuật Gotoh đã cải tiến. Thuật giải để tìm kiếm phép alignment được xây dựng dựa trên ý tưởng của Hirschberg khi giải quyết bài toán tìm chuỗi con lớn nhất của 2 chuỗi. Hirschberg (1975)[15] đã giới thiệu một kỹ thuật chia để trị nhằm giải quyết bài toán tìm chuỗi con dài nhất của 2 chuỗi với độ phức tạp ( )nmΟ và cần ( lg )Ο +n m không gian nhớ. Nội dung chính của giải thuật được xây dựng trong phần này sẽ là tìm kiếm điểm giữa (midpoint) của một alignment tối ưu, sau đó xem điểm vừa tìm được là điểm biên sử dụng kỹ thuật đệ quy để tiếp tục tìm phép alignment tối ưu của phần trước và sau của nó. Ý tưởng này được minh họa như hình dưới đây. Hình 4.1 Mô hình quá trình thực hiện giải thuật PSA Xét 2 chuỗi A, B có chiều dài m, n(1 ,1m n≤ ≤ ) A(a1a2…am ),B(b1b2…bn). Ta đưa ra một số qui ước sau: i* j* Điểm giữa Các kỹ thuật toán học cho bài toán so sánh đa trình tự Phạm Mạnh Hùng Trang 44 Ai là chuỗi con của A bắt đầu từ phần tử đầu tiên đến phần tử thứ i, Ai(a1 a2…ai). T iA : là phần sau của chuỗi A từ vị trí thứ i+1. T iA (ai+1…am) rev(A)m: chuỗi ngược của chuỗi A( amam-1…a1). rev(A)m-i: chuỗi ngược của chuỗi ai+1ai+2…am. Chúng ta gọi PSA là tối ưu nếu giá trị hàm đánh giá của nó là lớn nhất. Giải thuật sẽ tìm kiếm điểm M*(i*,j*) ( *, * [0, ]i j n∈ ) trên ma trận sao cho M* nằm trên con đường sinh ra PSA tối ưu(đồng nghĩa với M* chia ma trận thành 2 vùng mà cách chia của nó là tối ưu), với M* ta có thể thực hiện quá trình chia nhỏ phục vụ cho quá trình gọi đệ quy. M* có thể thuộc 1 trong 2 hình thức: ƒ M* là điểm kết thúc của một alignment của Ai* với Bj* và là điểm bắt đầu của 1 alignment của TiA *với * T jB (1). ƒ M* là điểm kết thúc của một alignment của Ai* với Bj* kết thúc bằng gap, và là điểm bắt đầu của 1 alignment của TiA *với * T jB , bắt đầu với 1 gap (2). Định lý 4.1 Gọi SP là giá trị của phép alignment bất kỳ của A và B. M(i,j)