Khóa luận Các phương pháp sắp hàng đa chuỗi nhanh

ng công thức sau để tính khoảng cách giữa từng “hàng xóm” với nút mới. Ở đây: f, g là 2 hàng xóm và u là nút mới được tạo thành: Khi một nút mới được tạo thành ta dùng công thức sau để tính khoảng cách của nó với các nút cũ: Ở đây, u là nút mới, k là nút cũ, f và g là 2 phần tử tạo nên nút mới u. 2.1.3 Progressive alignment Dựa vào cây hướng dẫn (guide tree) được tạo ra từ bước 2, chúng ta sử dụng sắp hàng các chuỗi từ nút lá cho đến gốc của cây. Mỗi bước sẽ là quá trình sắp hàng 2 nhóm chuỗi đã được sắp hàng trước đó sử dụng thuật toán quy hoạch động [9] [10]. Gap ở những nhóm chuỗi này được giữ nguyên trong kết quả tạo thành. Việc này lặp đi lặp lại cho đến khi gặp gốc của cây hướng dẫn. Đó là kết quả cuối cùng. 2.2. MUSCLE 2.2.1 Các loại khoảng cách và các cách xây dựng cây hướng dẫn MUSCLE[6] sử dụng hai cách để xác định khoảng cách giữa các chuỗi đó là khoảng cách K-mer[11] cho những cặp chuỗi chưa được sắp hàng và khoảng cách Kimura[12] cho những cặp đã được sắp hàng (có độ dài bằng nhau). K-mer[11] được định nghĩa là một chuỗi có độ dài bằng K của các amino acid đứng liền kề nhau. Thuật toán sử dụng K-mer không đòi hỏi cặp chuỗi cần phải sắp hàng mà ưu điểm lớn của nó là kết quả thu được với tốc độ khá cao (độ phức tạp thuật toán là O(L) với L là độ dài của chuỗi) [11]. Hình 1: Ví dụ về k-mer [6] Hình 1 thể hiện một ví dụ của K-mer, với K = 3, ở các chuỗi trên dễ dàng nhận thấy K-mer là 6 và tại các chuỗi ở phần dưới K-mer là 13. Tương tự như vậy với K = 4 các chuỗi bên trên K-mer là 4 và các chuỗi dưới K-mer là 9. Khoảng cách Kimura[12] giữa 2 chuỗi đã được sắp hàng (có độ dài bằng nhau) được tính theo công thức sau: DK2p= -1/2 ln(1-2P-Q) – ¼ ln(1/2Q) Ở đây, P là số lượng transition và Q là số lượng transversion trong hai chuỗi. Transition là một phép thay thế khi thay A bởi G, C bởi T hoặc ngược lại. Trong khi đó Transversion là phép thay thế A bởi C hoặc T hay các trường hợp tương tự. Sau khi xây dựng xong ma trận khoảng cách, MUSCLE[6] sử dụng thuật toán UPGMA[13] để nhóm các chuỗi lại với nhau. UPGMA[13] là một thuật toán xây dựng cây hướng dẫn cho việc sắp hàng từng profile. Đầu vào của UPGMA[13] là một ma trận khoảng cách. Ở mỗi bước 2 nút gần nhất được nhóm lại với nhau tạo thành một nút mới và ma trận khoảng cách được tính lại theo công thức sau: Ở đây, chúng ta tính khoảng cách giữa 2 nút A và B; x và y lần lượt là các nút ban đầu trong A và B; d(x, y) là khoảng cách giữa 2 nút x và y. 2.2.2 Profile alignment Để áp dụng bắt cặp sóng đôi vào các profiles, điều cần thiết là xác định một hàm tính điểm cho một cách sắp hàng. Ta coi i, j là 2 amino acid, pi là xác xuất i xuất hiện, pij là xác xuất i và j được align với nhau, fxi là xác suất của i xuất hiện tại cột x của profile thứ nhất, fxG là xác suất xuất hiện một gap tại cột x trong các profile. Khi đó MUSCLE[6] đưa ra hàm tính Log-Expectation theo công thức sau: LExy = (1 - fxG)(1- fyG) log ∑i∑j fxi fyi pij/pipj MUSCLE[6] sử dụng các tham số pi và pij là các tham số của ma trận 240 PAM VTML [14]. 2.2.3 Thuật toán Thuật toán MUSCLE[6] là một loạt các bước sử dụng các khái niệm đã trình bày ở trên. Tuy nhiên tổng quát lại nó bao gồm 3 phần chính là: Phần 1: draft progessive. Phần 2: improved progessive. Phần 3: Refinement. Mỗi phần làm một nhiệm vụ riêng biệt, nhưng kết nối chặt chẽ với nhau bởi đầu ra của phần này là đầu vào của phần tiếp theo. Hình 2: Các bước thực hiện của MUSCLE [6] Phần 1: draft progessive. Mục tiêu chính của phần này là tạo ra một đa chuỗi thẳng hàng mà vấn đề chính là tốc độ chứ không phải là sự chính xác. Sử dụng khoảng cách K-mer[8] để xác định khoảng cách giữa mỗi cặp của chuỗi dầu vào. Tạo ra ma trận khoảng cách D1. Sử dụng thuật toán UPGMA[13] để xây dựng cây hướng dẫn TREE1. Ở mỗi lá của cây Tree1, một profile được tạo ra từ một chuỗi đầu vào. Các nút trong cây được thăm theo thứ tự tiền tố (nghĩa là các nút con được thăm trước cha mẹ). Ở mỗi nút trong (internal node) phương pháp một bắt cặp sóng đôi (pairwise alignment - một đa chuỗi thẳng hàng nhưng chỉ có được xây dựng từ 2 profile) được xây dựng từ 2 nút con. Việc này lặp đi lặp lại cho đến khi gặp gốc của cây. Đó chính là đa chuỗi thẳng hàng – mục tiêu của bước này. Phần 2: improved progessive. Hầu hết các lỗi trong phần một đều do việc sử dụng khoảng cách K-mer[11], phần 2 sẽ tạo lại cây bằng cách sử dụng khoảng cách Kimura[12]. Phương pháp này cho kết quả tốt hơn nhưng đòi hỏi các chuỗi phải được sắp hàng (có độ dài bằng nhau). Với mỗi cặp chuỗi trong MSA1, chúng ta tính khoảng cách cho chúng sử dụng khoảng cách Kimura[12]. Kết quả ta có ma trận khoảng cách D2. Từ ma trận khoảng cách D2, sử dụng phương pháp UPGMA[13] ta tạo nên cây Tree2. Làm tương tự bước 1.3 với cây Tree2 ta được đa chuỗi thẳng hàng MSA2. Phần 3: Refinement. Tìm phương án tốt nhất. - Bước 3.1 Một cạnh được chọn từ Tree2 (cạnh được chọn bằng cách giảm dần khoảng cách tới gốc) - Bước 3.2 Chia cây Tree2 thành 2 phần bằng cách bỏ cạnh vừa chọn, sau đó tính lại các profile trên mỗi phần đó. - Bước 3.3 Tạo ra một đa chuỗi thẳng hàng (sequences alignment) mới bằng cách sắp hàng một lần nữa 2 profile vừa được tạo ra. - Bước 3.4 Nếu điểm SP (sum of pair)[23] được cải thiện thì cách sắp hàng mới được giữ lại, ngược lại ta bỏ đi. - Lặp lại các bước 3.1-3.4 cho tới khi hội tụ hoặc đi tới giới hạn do người sử dụng định nghĩa[15]. 2.3 MAFFT 2.3.1 Bắt cặp nhóm sử dụng FFT Biểu diễn một amino acid dưới dạng một vector Tần suất của sự thay thế các amino acid phụ thuộc mạnh mẽ vào các thuộc tính lí hóa của chúng, đặc biệt là các thuộc tính khối lượng (volume) và độ phân cực (polarity)[16]. Do đó để biểu diễn một amino acid a dưới dạng vector thì ta cần một vector trong đó các thành phần của vector này là v(a) – thể hiện thuộc tính khối lượng và p(a) – thể hiện thuộc tính độ phân cực[17]. Ở đây, MAFFT[4] đã sử dụng các giá trị v(a) và p(a) đã được chuẩn hóa để thuận lợi cho việc tính toán sau này. Công thức xác định các giá trị chuẩn hóa đó là: Trong dó các giá trị và là các giá trị trung bình, các giá trị là độ lệch tiêu chuẩn. Khi đó một chuỗi amino acid sẽ được chuyển thành một chuỗi các vector. Tính toán mối quan hệ giữa hai chuỗi amino acid Mối tương quan về mặt khối lượng giữa 2 chuỗi amino acid được tính theo công thức sau Trong đó N và M là độ dài của 2 amino acid. là giá trị của thuộc tính khối lượng của 2 amino acid tại vị trí thứ n. k là độ trễ trong phép tính. Việc định nghĩa k sẽ được nêu ở phần sau. Bằng một cách hoàn toàn tương tự ta sẽ tính được giá trị . Khi đó hàm quan hệ giữa 2 chuỗi amino acid này được định nghĩa theo công thức c(k) = cv(k) + cp(k). Nếu ta coi N = M trong trường hợp trên, khi đó độ phức tạp của thuật toán để tính mối quan hệ này là O(N2). Tuy nhiên, dựa vào biến đổi Fourier ta có thể tối ưu bước này và giảm độ phức tạp thuật toán xuống còn O(N logN) [18]. Xét với trường hợp tính toán cv(k), ta coi V1(m) và V2(m) là dạng biến đổi của v1(n) và v2(n) ta sẽ có Ta coi dấu dùng để biểu diễn các cặp khi sử dụng Fourier và dấu * thể hiện sự chuyển vị ma trận. Khi đó giá trị của hàm phụ thuộc sẽ được tính theo công thức Tìm kiếm đoạn tương đồng Độ trễ k giữa hai chuỗi là độ chậm của chuỗi thứ nhất với chuỗi thứ 2. Ở phần B hình dưới đây biểu diễn độ trễ k giữa 2 chuỗi 1 và 2 Hình 3: Ví dụ về độ trễ [4] Bằng thực nghiệm ta thấy, với các giá trị c(k) lớn, đồng nghĩa với việc tại độ trễ đó có thể xuất hiện các chuỗi tương đồng (như trên hình vẽ). Khi đó thuật toán sử dụng một bảng trượt với độ lớn của bảng là 30 amino acid (trong hình ví dụ thì độ lớn của bảng là 4). Sau đó lấy 20 giá trị độ trễ k có c(k) lớn nhất. Sử dụng bảng trượt trên hai chuỗi này, nếu tỉ lệ tương đồng vượt quá giới hạn (giá trị giới hạn mặc định là 0.7) thì coi như ta đã tìm được một vùng tương đồng, nếu có nhiều vùng tương đồng liên tiếp với nhau thì nối chúng lại với nhau. Tuy nhiên độ lớn của vùng tương đồng này không vượt quá 150 amino acid, nếu lớn hơn 150 thì phải tách chúng ra làm 2 đoạn tương đồng. Tạo ma trận tương đồng Sau khi đã tìm ra được các cặp đoạn tương đồng giữa 2 chuỗi. Khi đó ta định nghĩa một ma trận hai chiều là ma trận tương đồng như sau: Ma trận S(i,j) (i, j = 1 … n) – n là số đoạn tương đồng giữa h chuỗi. Tại S(i, j) nếu đoạn i của chuỗi một tương đồng với đoạn j của chuỗi hai thì điểm S(i,j) sẽ là điểm được tính như trên. Trường hợp còn lại S(i,j) = 0. Hình 4: Ví dụ về việc tạo ma trận tương đồng [4] Hình trên là một ví dụ của việc tạo ma trận tương đồng với 5 đoạn tương đồng giữa chuỗi 1 và chuỗi 2. Việc chọn đường đi tối ưu phụ thuộc vào S(2,3) và S(3,2). Nếu S(2,3) có giá trị lớn hơn thì đường đi sẽ theo đường in đậm như hình vẽ. Mở rộng từ bắt cặp giữa 2 sequence thành bắt cặp nhóm Ta có thể dễ dàng mở rộng tính toán bắt cặp giữa 2 chuỗi thành bắt cặp giữa 2 nhóm bằng cách thay công thức tính cv(k) và cp(k) với vi(n) và pi(n) bằng: Trong đó wi là trọng số của chuỗi i trong group1 Đặc biệt, trong trường hợp sử dụng cho các nucleotide chúng ta có thể áp dụng công thức sau để tính toán hàm quan hệ: c(k) = cA(k) + ct(k) + cc(k) + cg(k). 2.3.2 Hệ thống tính điểm Ma trận điểm: Để nâng cao hiệu năng của việc bắt cặp, một hệ thống tính điểm (ma trận tương đồng và điểm phạt gap) đã được sửa đổi [19]. Ma trận điểm cho thuật toán được định nghĩa theo công thức : Trong đó: a, b là hai amino acid các giá trị average được tính Mab là ma trận ban đầu, mặc định ở đây là ma trận PAM 200 [15]. fa là tần số xuất hiện của amino acid a [15]. Sa là tham số thêm vào. Đối với các acid nucleic ta có các giá trị fa = 0.25 và Sa = 0.06 Điểm phạt Công thức tính điểm phạt được định nghĩa theo công thức: Trong đó: Sop là điểm phạt khi tạo ra một gap gstart(j) là số gap xuất hiện bắt đầu từ vị trí thứ j gend(i) là số gap xuất hiện từ vị trí bắt đầu cho tới i zm(i) = 1 và am(i) = 0 nếu tại vị trí thứ i của chuỗi có gap ngược lại zm(i)= 0 và am(i) = 1. 2.4 PROBCONS PROBCONS[5] sử dụng mô hình Markov ẩn [20] cho thuật toán progressive. Điểm khác biệt chính giữa PROBCONS[5] và các phương án tiếp cận khác đó chính là việc sử dụng ước lượng chuẩn xác cực đại (maximum expected accuracy) chứ không phải là cách sử dụng mô hình Viterbi[21], ngoài ra PROBCONS[5] còn sử dụng ước lượng các phép biến đổi để bảo toàn thông tin của các chuỗi trong khi bắt cặp. Mặt khác PROBCONS[5] còn sử dụng ma trận chuyển đổi chuẩn BLOSUM62[22] và điểm phạt phát sinh khi thêm dấu cách trong việc bắt cặp cũng được huấn luyện bởi ước lượng cực đại. Các bước thực hiện của thuật toán PROBCONS[5] như sau: Cho m chuỗi, S = {S(1) ,… , S(m)} Bước 1: Tính các ma trận xác suất: Với mỗi cặp chuỗi x và y của S, và với mỗi và . Ta tính ma trận Pxy (i, j) = P(xi ~ yj  a*| x, y). Đây là xác suất ký tự xi và yj được bắt cặp với nhau trong a* - pairwise alignment được sinh ra từ mô hình HMM. Bước 2: Tính toán kỳ vọng của độ chính xác. Định nghĩa kỳ vọng của độ chính xác của một bắt cặp sóng đôi (pairwise alignment) a được tạo bởi hai chuỗi x và y được tính bằng số lượng liên kết chính xác của các ký tự chia cho chiều dài của chuỗi ngắn hơn. Ea* (accuracy(a, a*)|x, y) = Với mỗi cặp chuỗi x và y ta tính toán cực đại của kỳ vọng của độ chính xác bằng phương pháp quy hoạch động và đặt E(x, y) = Ea* (accuracy(a, a*)|x, y) Bước 3: Ước lượng tính vững chắc của các bước chuyển đổi. Ở bước này ta ước lượng lại P(xi ~ yj  a*| x, y) ở bước trên. Ở dạng ma trận chuyển đổi ta có công thức: P’xy Bước 4: Xây dựng cây hướng dẫn (guide tree). Xác định ma trận tương đồng giữa 2 chuỗi x và y bằng giá trị E(x, y) được tính ở bước 2. Từ đó ta xây dựng cây hướng dẫn. Bước 5: Progressive alignment Dựa vào cây hướng dẫn được dựng từ bước 4, ta tiến hành bắt cặp. Qua các thực nghiệm với các bộ dữ liệu chuẩn, kết quả của PROBCONS[5] là rất cao (cao nhất với rất nhiều bộ dữ liệu chuẩn). Tuy nhiên, một hạn chế cực lớn của PROBCONS[5] là tốc độ, về mặt này PROBCONS[5] không thể so sánh với các thuật toán trên. Chương 3. Cây quyết định 3.1 Cách giải quyết của Chuong B. Do và Kazutaka Katoh Như đã trình bày ở trên, hiện nay có khá nhiều phương pháp sắp hàng đa chuỗi, nhưng mỗi phương pháp lại có một đặc điểm riêng kèm theo đó là những ưu khuyết điểm riêng. Đôi khi một phương pháp cho kết quả tốt với bộ dữ liệu này, lại không phù hợp với bộ dữ liệu khác. Một phương pháp cho kết quả rất cao nhưng tốc độ lại quá chậm, hoặc không thể xử lý những dữ liệu quá lớn. Qua đó có thể thấy việc xây dựng cây quyết định để giải quyết vấn đề chọn phương pháp tối ưu nhất cho mỗi loại dữ liệu đầu vào là vô cùng quan trọng. Hai nhà khoa học Chuong B. Do và Kazutaka Katoh đã đề ra một giải pháp[26] là nghiên cứu về từng phương pháp và từng loại dữ liệu, qua đó có thể đưa ra được những phương pháp phù hợp với từng bộ dữ liệu cả về mặt điểm chuẩn lẫn thời gian xử lý. Hai tác giả đã chia các loại dữ liệu ra thành 3 phần riêng biệt cần phải xem xét là: Dữ liệu yêu cầu tìm thành phần lặp. Dữ liệu đầu vào có số lượng chuỗi lớn ( > 200 chuỗi ). Dữ liệu đầu vào có chuỗi có độ dài lớn ( > 2000 amino acid). Đối với loại dữ liệu thứ nhất, chúng ta không xem xét nó trong phạm vi của khóa luận này. Đối với dữ liệu đầu vào có số lượng chuỗi lớn ( > 200 chuỗi ). Hai tác giả đã chỉ ra độ phức tạp thuật toán trong việc tính toán ma trận khoảng cách là nhân tố chủ yếu trong việc dẫn đến thời gian thực hiện quá lâu của các phương pháp. Cho nên những phương pháp có cách xây dựng ma trận khoảng cách với độ phức tạp thấp sẽ được chọn. Ở đây, 2 phương pháp MUSCLE và FFT-NS-2 đã được chọn. Với dữ liệu có chuỗi có độ dài lớn ( > 2000 amino acid), thì độ phức tạp không gian của thuật toán là nguyên nhân chính dẫn đến việc thuật toán có xử lý được loại dữ liệu này không. Hầu hết các phương pháp có độ phức tạp không gian là O(L2) với L là độ dài trung bình của các chuỗi. Đối với loại dữ liệu này, những phương pháp có độ phức tạp không gian tuyến tính ( CLUSTALW, FFT-NS-1 và FFT-NS-2 ) được sử dụng. Tuy nhiên cách chia của Katoh và Chuong B Do còn chưa được rõ ràng, chưa chỉ rõ đối với từng khoảng nhỏ dữ liệu. Do đó tôi sẽ phát triển tiếp phương pháp của 2 tác giả Chuong B Do và Kazukata Katoh trong khóa luận này. Trong khóa luận này, ta tập trung nghiên cứu về 4 chương trình sắp hàng đa chuỗi tốt nhất hiện nay là: CLUSTALW, MUSCLE, PROBCONS, MAFFT (bao gồm L-INS-i, E-INS-i, G-INS-i, FFT-NS-1, FFT-NS-2). Ở đây, chúng ta tập trung vào 2 vấn đề tốc độ và điểm chuẩn (benchmark) để đưa ra 2 cây quyết định cho 2 yêu cầu về tốc độ và benchmark. 3.2 Vấn đề tốc độ Với các phương pháp kể trên, chúng đều có một giới hạn về dữ liệu đầu vào khác nhau. Vấn đề ở đây, là kiểm tra giới hạn của từng phương pháp. Để kiểm tra, chúng ta sử dụng 3 bộ dữ liệu là BAliBASE2[23], BAliBASE3[24] và Pfam-A[25]. Chúng ta chia thành 3 tình huống riêng biệt: - Dữ liệu với số lượng chuỗi lớn ( > 200 chuỗi). - Dữ liệu với số lượng chuỗi nhỏ, tổng số amino acid nhỏ. - Dữ liệu với độ dài của chuỗi quá lớn ( > 2000 amini acid). 3.2.1 Dữ liệu với số lượng chuỗi lớn ( > 200 chuỗi) Trong trường hợp này, tốc độ đóng một vai trò vô cùng quan trọng. Trong trường hợp này ta có các phương pháp có thể chạy được là: MUSCLE, FFT-NS-1, FFT-NS-2. Những phương pháp này cho kết quả tương đối thấp, nhưng có tốc độ chạy khá cao. Do đó ta sẽ kiểm tra giới hạn của 3 phương pháp này, các test được trích từ bộ dữ liệu Pfam-A ( là bộ dữ liệu chỉ là một tập hợp các chuỗi protein). Với các test có số lượng chuỗi từ 200 đến 500 chuỗi. Bảng 3: Kiểm tra các MUSCLE, FFT-NS2, FFT-NS1 với các test có số lượng chuỗi từ 200 đến 500 chuỗi. Số lượng chuỗi MUSCLE FFT-NS-2 FFT-NS-1 200 – 250 28 / 128 / 68.8 5 / 12 / 9.4 3 / 11 / 8 250 - 300 63 / 284 / 113.8 9 / 17 / 12.6 9 / 18 / 11.8 300 – 350 47 / 183 / 136.8 8 / 19 / 13.2 7 / 20 / 12.4 350 – 400 74 / 339 / 167.6 12 / 31 / 18 8 / 21 / 12.8 400 – 450 102 / 604 / 257.4 15 / 56 / 26.2 13 / 48 / 22.8 450 – 500 145 / 738 / 372.2 18 / 60 / 34 16 / 49 / 27.6 Kết quả chỉ ra thời gian chạy nhanh nhất, lâu nhất của một test, và thời gian trung bình của các test đó. Từ những số liệu trên ta có thể thấy. MUSCLE chỉ nên chạy với các test dưới 400 chuỗi. Tiếp tục với các phương pháp FFT-NS2 và FFT-NS1, ta có: Bảng 4: Kiểm tra FFT-NS2 với các dữ liệu có số lượng chuỗi lớn hơn 400 Số lượng chuỗi FFT-NS-2 500 – 1000 20 / 306 / 90.4167 1000 – 2000 97 / 454 / 219.455 2000 – 3000 250 / 535 / 397.727 3000 – 4000 350 / 631 / 486 4000 – 5000 497 / 824 / 651.4 Dựa vào các số liệu trên, ta có thể thấy FFT-NS-2 chỉ nên giới hạn chạy với các dữ liệu dưới 4000 chuỗi. Trong 2 phương pháp FFT-NS-2 và FFT-NS-1 thì FFT-NS-2 có bước xử lý thô là FFT-NS-1. Do đó FFT-NS-2 có tốc độ chậm hơn nhưng cho kết quả cao hơn FFT-NS-1. Trong các phương pháp được xét, FFT-NS-1 là phương pháp có tốc độ cao nhất, nhưng cho kết quả tồi nhất. Đây là phương pháp chỉ nên sử dụng khi các phương pháp khác đã không thể chạy được. 3.2.2 Dữ liệu với số lượng sequence nhỏ, tổng số amino axit nhỏ Với trường hợp này, hầu hết các phương pháp đều có thể chạy được, khi đó các điểm đánh giá phải được đặt nên hàng đầu. Chúng ta nên xử dụng các phương pháp có độ chính xác cao như: PROBCONS, L-INS-i, G-INS-i, E-INS-i. Với PROBCONS, kiểm tra với 2 bộ dữ liệu BAliBASE2[23] và BAliBASE3[24], ta có các thông số như sau: Bảng 5: Thời gian chạy của PROBCONS theo tống số amino acid Tổng số amino acid của dữ liệu Thời gian chạy 0 – 1000 2.3 1000 – 2000 8 2000 – 3000 18.1 3000 – 4000 43.3 4000 – 5000 71.5 5000 – 6000 103.5 6000 – 7000 140.6 7000 – 8000 190.8 8000 – 9000 250.25 9000 – 10000 301 Từ những số liệu trên có thể thấy, PROBCONS chỉ nên dùng với những bộ dữ liệu nhỏ (có tổng số amino acid nhỏ). Một cách chính xác, khi cần chạy nhanh và chạy chính xác nên giới hạn tổng số amino acid của dữ liệu lần lượt là 7000 và 9000. Đối với 3 phương pháp L-INS-i, G-INS-i, E-INS-i, nên giới hạn dữ liệu đầu vào ở mức 200 sequences (Theo Katoh [26]). 3.2.3 Dữ liệu với độ dài của chuỗi quá lớn ( > 2000 amino acids) Đối với các dữ liệu loại này, thì độ phức tạp không gian là một vấn đề quan trọng cần phải xtôi xét đầu tiên trong việc lựa chọn các phương pháp sắp hàng đa chuỗi. Hiện nay, hầu hết các phương pháp sắp hàng đa chuỗi đều hướng đến việc sử dụng thuật toán quy hoạch động với việc sử dụng độ phức tạp không gian là O(L2) (Ở đây, L là độ dài trung bình của các chuỗi). Đối với các chuỗi đặc biệt dài (> 2000 amino acids) các phương pháp có độ phức tạp không gian tuyến tính O(L) là sự lựa chọn tối ưu để giải quyết vấn đề này. Với các phương pháp đang được xem xét thì CLUSTALW, FFT-NS-2 và FFT-NS-1 là những phương pháp như thế. 3.3 Vấn đề điểm chuẩn (benchmark) 3.3.1 Với các chuỗi có độ tương đồng cao Độ tương tự (identity), là thuật ngữ chỉ việc mức độ giống nhau của các chuỗi đầu vào. Theo Katoh và Chuong B. Do, với dữ liệu đầu vào có mức độ tương đồng cao ( > 35 % ), thì việc chạy bất cứ chương trình nào cũng không ảnh hưởng quá lớn đến kết quả cuối cùng [26]. Còn việc kiểm tra mức độ tương tự, tôi đã sử dụng một chương trình sắp hàng đa chuỗi có tốc độ cao (cụ thể ở khóa luận này là FFT-NS-1) để tạo ra các chuỗi sắp hàng (có độ dài bằng nhau), sau đó kiểm tra mức độ tương tự với độ phức tạp tuyến tính (O(L) với L là độ dài của chuỗi sau khi sắp hàng). 3.3.2 Với các chuỗi có độ tương đồng thấp Với các chuỗi có mức độ tương tự thấp (<= 35 % ), các hệ thống tính điểm chuẩn khác nhau đã thống nhất xác định PROBCONS và L-INS-i là các phương pháp cho kết quả cao nhất hiện nay. Tuy nhiên phương pháp PROBCONS khi chạy với dữ liệu là các chuỗi DNA luôn rất chậm. Do đó, với các dữ liệu là chuỗi DNA ta không nên sử dụng phương pháp PROBCONS. Nói chung, sắp hàng các chuỗi có độ tương tự thấp được hiểu là 1 trong 3 trường hợp sau: Trường hợp 1: global homology – tương đồng (homology) trên toàn chiều dài của chuỗi protein Hình 5: Ví dụ về global homology [4] Ở đây, X chỉ ra phần có thể được align, o là các phần không được align và – là gap. Theo hình trên ta có thể thấy, toàn chiều dài các chuỗi là các phần có thể được align. Đây là trường hợp đơn giản nhất và phương pháp PROBCONS và G-INS-i là 2 phương pháp cho kết quả tốt nhất trong các phương pháp đang xét. Trường hợp 2: local homology – tương đồng (homology) được bao quanh bởi các miền không tương đồng Hình 6: Ví dụ về local homology [4] Hình trên chỉ ra một tập các chuỗi có chứa trong nó một miền có thể align và xung quanh nó là các phần không tương đồng. Khi đó, L-INS-i là phương pháp tối ưu. Trường hợp 3: Các đoạn gap nội khối dài - các khoảng tương đồng (homology) ngắn chia tách bởi các đoạn gap nội khối Hình 7: Ví dụ về các đoạn gap nội khối [4] Trong trường hợp này, có nhiều vùng có thể align, nhưng hầu hết chúng khá rời rạc và được tách bởi những đoạn gap rất dài. Khi đó E-INS-i là phương pháp cho kết quả tốt nhất trong các phương pháp kể trên. Tuy nhiên, trong hầu hết các hệ thống tính điểm chuẩn, PROBCONS và L-INS-i là hai phương pháp cho kết quả tốt nhất. 3.4 Cây quyết định Có hai yêu cầu cần phải giải quyết là: tốc độ và benchmark, cho nên ta sẽ tạo hai cây quyết định dựa trên những lý thuyết đã trình bày ở trên. 3.4.1 Cây quyết định cho yêu cầu tốc độ xử lý cao Hình 8: Cây quyết định với yêu cầu xử lý tốc độ cao 3.4.2 Cây quyết định cho yêu cầu tốc điểm chuẩn cao Hình 9: Cây quyết định với yêu cầu xử lý với điểm chuẩn cao Trong phương pháp được đề xuất ở đây, tôi chưa xử lý việc tìm cách phát hiện kiểu của các dữ liệu đầu vào có độ tương đồng thấp. Ở đây ta mặc định với các chuỗi có độ tương tự nhỏ (<= 35 %), chúng ta chỉ sử dụng 2 phương pháp L-INS-i và PROBCONS (là 2 phương án cho kết quả tốt nhất hiện nay). Tuy nhiên kết quả cuối cùng khi chưa xử lý vấn đề này cũng rất khả quan. Chương 4: Kết quả thực nghiệm và bình luận Một đánh giá toàn diện và so sánh được các chương trình sắp hàng đa chuỗi đòi hỏi một số lượng lớn các sự dữ liệu được sắp xếp chính xác mà có thể được sử dụng như các bộ kiểm thử. Các bộ dữ liệu này có thể chỉ ra được hiệu suất của các chương trình sắp hàng đa chuỗi phụ thuộc vào số lượng các chuỗi, mức độ giống nhau giữa các chuỗi và số lượng các phép chèn thêm vào liên kết này. Các yếu tố khác cũng có thể ảnh hưởng đến chất lượng liên kết chẳng hạn như độ dài của chuỗi, … BAliBASE là bộ dữ liệu đáp ứng được đầy đủ các yêu cầu như thế. Do đó trong khóa luận này tôi sẽ sử dụng BAliBASE để kiểm tra hiệu năng của hai phương pháp sử dụng cây quyết định (cả về tốc độ lần điểm chuẩn) và so sánh nó với các chương trình sắp hàng đa chuỗi khác. 4.1 Giới thiệu về BAliBASE BAliBASE - Benchmark Alignment dataBASE là một bộ dữ liệu được xây dựng bởi các nhà khoa học Julie D. Thompson, Olovier Poch và một số nhà khoa học khác. Việc xây dựng bộ dữ liệu BAliBASE hoàn toàn dựa trên những kết quả đã được kiểm chứng trước đó đồng thời bắt cặp dựa trên kinh nghiệm của chính những nhà khoa học này. BAliBASE là một bộ dữ liệu mở, được thiết kế để phục vụ cho mục đích đánh giá các chương trình sắp hàng đa chuỗi. Nó đặt ra tất cả các trường hợp gặp phải trong quá trình sắp hàng. Cơ sở dữ liệu của BAliBASE được làm một cách thủ công với các chú thích chi tiết. 4.1.1 BAliBASE 2 BAliBASE 2 có tất cả 8 reference, nhưng chỉ thường sử dụng 5 reference đầu tiên. Các file dữ liệu được cung cấp dưới định dạng RSF hoặc MSF. 4.1.2 BAliBASE 3 BAliBASE 3 bao gồm 5 reference. Mỗi reference bao gồm một số lượng file. Các file có tên là BBnnnnn bao gồm các chuỗi full-length, trong khi các file có tên là BBSnnnnn là các chuỗi chỉ chứa các vùng tương đồng (homologous). 4.1.3 Cách đánh giá của BAliBASE BAliBASE sử dụng hai hệ số điểm là sum of pair (SP) và total colum (TC) để kiểm tra tính chính xác của một đa chuỗi thẳng hàng so với kết quả mà các nhà khoa học bắt cặp một cách thủ công. Điểm SP được tính theo thuật toán sau: - Đặt s(x, y). Ở đây x và y là hai amino axit, s(x, y) là điểm khi bắt cặp x với y trong đa chuỗi thẳng hàng. khi đó ta sẽ có giá trị của s(x, y) tương ứng là: s(x, y) = 1 nếu x và y đều là một amino axit. s(x, y) = -1 nếu x và y là hai amino acid khác nhau. s(x, y) = -2 nếu x là gap, y khác gap và ngược lại. s(x, y) = 0 nếu x và y đều là gap. - Giả sử SP(mi) là giá trị của điểm “sum of pair” ở cột thứ i của đa chuỗi thẳng hàng mà ta cần tính điểm ( phân biệt với đa chuỗi thẳng hàng mà các nhà khoa học đã sắp hàng bằng tay để làm kết quả so sánh ), giá trị SP(mi) được tính bằng cách lấy tổng của các s(x, y) trong đó x và y là các amino axit được lấy từ cột thứ i của đa chuỗi thẳng hàng. Sau đó điểm SP của cả đa chuỗi sẽ được tính bằng cách lấy tổng tất cả các điểm SP(mi) của tất cả các cột trong đa chuỗi. Một ví dụ đơn giản cho việc tính toán SP(mi) như sau: Bảng 6: Tính toán SP(mi) m1 m2 m3 m4 m5 m6 m7 m8 m9 seq1 G T T C C T G - T seq2 - T G C - T G - T seq3 G T G C - T T - T Score -3 3 -1 3 -4 3 -1 0 3 Ví dụ trên thể hiện một đa chuỗi thẳng hàng với 3