Khóa luận Sắp hàng hoàn chỉnh hai hệ genome - Hà Tuấn Cường

ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC CÔNG NGHỆ Hà Tuấn Cường SẮP HÀNG HOÀN CHỈNH HAI HỆ GENOME KHOÁ LUẬN TỐT NGHIỆP ĐẠI HỌC HỆ CHÍNH QUY Ngành: Công Nghệ Thông Tin HÀ NỘI – 2010 HÀ NỘI - 20 (chữ hoa, 12pt, đậm, căn giữa) ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC CÔNG NGHỆ Hà Tuấn Cường SẮP HÀNG HOÀN CHỈNH HAI HỆ GENOME KHOÁ LUẬN TỐT NGHIỆP ĐẠI HỌC HỆ CHÍNH QUY Ngành: Công Nghệ Thông Tin GV hướng dẫn: TS. Lê Sỹ Vinh HÀ NỘI – 2010 HÀ NỘI - 20 (chữ hoa, 12pt, đậm, căn giữa) Lời cảm ơn Lời đầu tiên, em xin gửi lời cảm ơn sâu sắc nhất đến thầy giáo TS. Lê Sỹ Vinh người đã không quản vất vả tận tình hướng dẫn em trong suốt thời gian làm khóa luận tốt nghiệp vừa qua. Em cũng xin bày tỏ lòng biết ơn tới các thầy, cô giáo trong trường Đại học Công nghệ - Đại học Quốc gia Hà Nội. Các thầy cô đã dạy bảo, chỉ dẫn chúng em và luôn tạo điều kiện tốt nhất cho chúng em học tập trong suốt quá trình học đại học. Em cũng xin gửi lời cảm ơn tới thầy giáo PGS.TS. Từ Minh Phương, người đã cho em những lời khuyên bổ ích trong quá trình làm khóa luận. Tôi cũng xin cảm ơn những người bạn của mình, các bạn đã luôn ở bên tôi, giúp đỡ và cho tôi những ý kiến đóng góp quý báu trong học tập cũng như trong cuộc sống. Cuối cùng con xin gửi tới bố mẹ và toàn thể gia đình lòng biết ơn và tình cảm yêu thương nhất. Con xin dành tặng bố mẹ kết quả mà con đã đạt được trong suốt bốn năm học đại học. Con cám ơn bố mẹ và chị nhiều. Khóa luận được tài trợ một phần bởi đề tài nghiên cứu QC.09.09 thuộc Đại học Quốc Gia Hà Nội. Hà Nội, tháng 5 năm 2010 Hà Tuấn Cường Tóm tắt Sự phát triển của công nghệ giải mã trình tự đã giúp giải mã ngày càng nhiều các hệ gen, đặc biệt là những hệ gen có kích thước vừa và nhỏ như vi rút hay vi khuẩn (hơn 7000 bộ gen của vi rút và vi khuẩn đã được giải mã). Bên cạnh đó hệ gen của những sinh vật bậc cao cũng đã được giải mã hoàn chỉnh như người, chó, chuột. Điều đó dẫn đến một nhu cầu cấp thiết là phải nghiên cứu các phương pháp và xây dựng một chương trình so sánh và bắt cặp trình tự cho hai hệ gen. Trong khóa luận này, em xin được trình bày phương pháp và xây dựng một chương trình so sánh bắt cặp trình tự hoàn chỉnh cho hai hệ gen. Chương trình cho phép bắt cặp toàn bộ các ADN trên cả hai hệ gen, xác định được cả những biến đổi của tửng nucleotide và các biến đổi ở mức độ gen. Chương trình được xây dựng dựa trên cở sở kết hợp và cải tiến các phương pháp đã có như “Pairwise Alignment with Rearrangement” [23], BLASTZ [18] và “Optimal Alignment with Linear space” [9]. Qua đó khắc phục những hạn chế và lựa chọn những ưu điểm của chúng để tạo thành một chương trình sắp hàng hệ gen hoàn chỉnh. Chương trình đã được thực nghiệm kết quả trên các dữ liệu mô phỏng và các dữ liệu thật được lấy từ Gen Bank tại NCBI và thu được những kết quả khả quan. Đối với các dữ mô phỏng, kết quả sắp hàng của chương trinh cho thấy đã xác định được các đoạn gen có độ tương đồng rất cao, tỷ lể sắp hàng giữa các nucleotide giống nhau đạt mức trên 97%. Khi thực nghiệm với dữ liệu thật và so sánh độ tương đồng với giá trị bắt cặp thu được khi chạy phương thức Hungarian[8] với các hệ gen được chia sẵn bằng cách sử dụng các đoạn gen cung cấp tại Gen Bank cũng cho kết quả tương đương thậm chí tốt hơn trong hầu hết các trường hợp. Mục lục Danh sách hình vẽ Hình 1. Ví dụ về một trình tự. 8 Hình 2: Các biến đổi ở mức độ gen giữa Người và Khỉ 13 Hình 3:Ví dụ về phép biến đổi trong “Simulaneous Character Swapping". 20 Hình 4: Single Swap (trái) và Couple Swap (phải) 22 Hình 5:Bắt cặp trình tự với Ukkonen Barrier 29 Hình 6: Giao diện chương trình 31 Hình 7: Kết quả chương trình 32 Danh sách các bảng Bảng 1: Ma trận trọng số của BLASTZ 26 Bảng 2: Kết quả Test với số Inversion – Move là 0 34 Bảng 3: Kết quả Test với số Inversion – Move là 1 34 Bảng 4: Kết quả Test với số Inversion – Move là 2 34 Bảng 5: Kết quả Test với số Inversion – Move là 3 35 Bảng 6: Kết quả chạy dữ liệu thật 37 Lời mở đầu Năm 1854, Charles Darwin cho xuất bản quyển sách “Nguồn gốc của các loài sinh vật”, một công trình nghiên cứu sinh học nổi tiếng và đặt nền tảng cho thuyết tiến hóa của ông. Trong đó có viết “tất cả các động vật tương tự nhau phải tiến hóa từ một tổ tiên chung và tất cả các sinh vật phải tiến hóa từ một vài hoặc một tổ tiên chung đã sống cách đây nhiều triệu năm.” [7] Bộ gen của sinh vật là một trình tự ADN, theo thuyết tiến hóa thì chúng cùng được biến đổi và phát triển từ một tổ tiên chung. Trải qua hàng triệu năm tiến hóa và phát triển, một số đoạn gen được mất đi cũng như bị di chuyển vị trí so với ban đầu, hình thành lên những hệ gen khác nhau đại diện cho hàng tỷ sinh vật trên trái đất. Một trong những nhiệm vụ cần thiết là phải tìm ra mối quan hệ về mặt cấu trúc giữa các hệ gen sinh vật, qua đó xây dựng lên một bức tranh toàn cảnh về sự tương tự và tiến hóa giữa các sinh vật trên hành tinh. Với sự phát triển của công nghệ giải mã trình tự, ngày càng nhiều các hệ gen đã được giải mã hoàn chỉnh và được lưu trữ trong các ngân hàng cơ sở dữ liệu gen. Việc so sánh và tìm ra sự tương đồng giữa các hệ gen một cách thủ công là không thể thực hiện được. Do đó dẫn đến một nhu cầu cấp thiết phải nghiên cứu phương pháp và xây dựng chương trình để so sánh và bắt cặp trình tự cho hai hệ gen. Mặc dù một số phương pháp đã được nghiên cứu và phát triển, chúng mới chỉ tập trung vào xác định và bắt cặp cho những vùng ADN có độ tương đồng cao giữa hai hệ gen. Tức là, một phần lớn trong hệ gen có thể không được bắt cặp và so sánh khi ta tiến hành với các loài sinh vật có hệ gen khác nhau nhiều. Vì vậy cần phải xây dưng một hệ thống có khả năng bắt cặp được toàn bộ các ADN trên hai hệ gen. Chương 1. Giới thiệu Chương này giới thiệu về những kiến thức cơ bản về tin sinh học, bài toán bắt cặp trình tự và bắt cặp trình tự theo hệ gen. Nội dung giới thiệu được dựa một phần trên bài giảng của Viện Đại học Ohio State, Hoa Kỳ [13] Trình tự Một hệ gen của một sinh vật được thể hiện là một trình tự của các ADN. Trình tự là một dãy tuyến tính các phần tử được sặp xếp theo thứ tự. Như vậy một trình tự chứa hai loại thông tin: thông tin về phần tử và thông tin định vị - thông tin về vị trí tương đối của từng phần tử so với các phần tử khác. Các thông tin định vị có thể được xác định theo nhiều cách như theo trục, theo thời gian, vị trí của nhiễm sắc thể hoặc trong 1 vòng protein. Hình 1. Ví dụ về một trình tự. Hình trên cùng: 1 đoạn 18S rDNA của sâu bọ khác cánh. Hình giữa trên: Tổng quát cơ thể động vật chân dốt. Hình giữa dưới: Orthopteran stridulation. Hình dưới cùng: Đoạn gen mtDNA [13] Các loài sinh vật được tiến hóa từ một tổ tiên chung, biến đổi qua các dạng hình thái trong quá trình tiến hóa và phát triển. Khi đề cập đến trình tự, có ba vấn đề chúng ta cần phải nói đến là hệ thống các ký tự trong trình tự, các phép biến đổi trình tự và hàm khoảng cách đánh giá sự thay đổi của trình tự. 1.1.1. Hệ thống ký tự Tập hợp các phần tử có thể xuất hiện trong trình tự được gọi là một hệ thống các ký tự ( ∑ ) , trong ADN, người ta sử dụng một thệ thống kí tự ∑ = {A, C, G, T, λ ) trong đó A, C, G, T đại diện cho 4 nucleotides : adenine (A),  cytosine (C),  guanine (G) và thymine (T).  λ là ký tự đặc biệt đại diện cho 1 khoảng trống là 1 vị trí mà không có nucleotide thực tế. Trong hầu hết các trường hợp, ký tự gap (‘-‘) có thể được sử dụng để thay thế cho λ. Bất kỳ một trình tự nào cũng là một sự thể hiện bởi các phần tử có thể xuất hiện trong trình tự và được định nghĩa trong ∑. 1.1.2. Các phép biến đổi Trong quá trình tiến hóa, có 4 phương thức chính để biến đổi một trình tự là phép thay thế, phép chèn – xóa, đảo ngược và dịch chuyển. Biến đổi phức tạp xảy ra là sự kết hợp của 2 phép đảo ngược và dịch chuyển, sự kết hợp này gây ra sự khó khăn trong việc dựng lại mối quan hệ giữa các trình tự trong những hệ thống phân tích phức tạp. Phép thay thế: Phép chèn – xóa: Phép đảo ngược: Phép dịch chuyển: 1.1.3. Khoảng cách Một điều quan trọng và cần thiết là xây dựng một hàm mục tiêu đánh giá khoảng cách giữa hai trình tự, qua đó đánh giá sự tương đồng, mối quan hệ giữa hai hệ gen. Khoảng cách này có thể được tính toán theo một số hàm như thay thế, chèn, xóa làm biến đổi một trình tự này thành một trình tự khác. Khoảng cách giữa hai trình tự có thể chỉ được tính đơn giản chỉ là chi phí thay thế (Hamming Hamming [10]) trong những trình tự có độ dài bằng nhau hay phức tạp hơn bao gồm cả chi phí chèn – xóa và dịch chuyển. Sắp hàng trình tự Sắp hàng trình tự là một thủ tục cực kỳ quan trọng trong Tin sinh học, nó được xem là nền tảng cho tất cả các thủ tục khác. Vấn đề đặt ra là tạo ra những sự sắp hàng giữa các nucleotide thông qua việc chèn các ký tự gap, làm cho khoảng cách giữa hai trình tự tức chi phí sửa chữa (là tổng chi phí cho các sự kiện chèn – xóa, thay thế các nucleotide) giữa hai trình tự là nhỏ nhất (hoặc lớn nhất). Đầu vào là 2 trình tự X = (x1, x2, …xp) và Y = (y1, y2, …yq), sắp hàng trình tự X và Y là cách chèn các kí tự trống vào hai trình tự X và Y sao cho chúng có độ dài bằng nhau và khoảng cách (chi phí sửa chữa) giữa hai trình tự là nhỏ nhất (hoặc lớn nhất). Các thuật toán quy hoạch động đầu tiên cho việc sắp hàng giữa các chuỗi ký tự được trình bày bởi Levenshtein [14], với độ phức tạp về thời gian và bộ nhớ là O(n2). Needleman và Wunsch [16] lần đầu tiên áp dụng thuật toán này vào lĩnh vực Tin sinh học năm 1970. Yêu cầu bộ nhớ giảm xuống còn O(n) bởi Hirschberg[12] trong khi thời gian chạy vẫn là O(n2). Những cải tiến của Ukkonen [21,22] với những cặp trình tự có khoảng cách độ dài là d, thuật toán yêu cầu thời gian O(nd) cho trường hợp xấu nhất và độ phức tạp thời gian trung bình là O(d2+n). Thuật toán Quy hoạch động tính toán chi phí bắt cặp theo công thức sau: (1) Cost[i][j] là chi phí bắt tới vị trí i của trình tự 1 và vị trí j của trình tự 2, σi,j là chi phí thay thế nucleotide ở vị trí thứ i của trình tự 1 và ở vị trí j của trình tự 2, σindex là chi phí chèn- xóa một nucleotide. Pairwise Alignment by Needleman and Wunsch Cost[0][0] ← 0 {Khởi tạo cột đầu tiên} for i = 0 to |X| do Cost[i][0] ← Cost[i-1][0] + σindex {Khởi tạo hàng đầu tiên} for j = 0 to |Y| do Cost[0][i] ← Cost[0][j-21] + σindex {Quy hoạch động} for i = 1 to |X| do for j = 1 to |Y| do ins ← Cost[i-1][j] + σindex del ← Cost[i][j-1] + σindex sub ← Cost[i-1][j-1] + σi,j Cost[i-1][j-1] ← min(ins, del, sub) end for end for return Cost[|X|][|Y|] Thuật toán Needleman – Wunsch hoạt động khi mà chi phí cho việc chèn – xóa các nucleotide là một trọng số cố định. Tức chi phí cho việc chèn một đoạn gap có độ dài k là wk = kw1. Trên thực tế, việc tính toán chi phí chèn – xóa thường phức tạp hơn, bao gồm chi phí cho việc bắt đầu và mở rộng các đoạn gap. Waterman [25], tiến hành thực nghiệm trên một khối lượng lớn các trình tự với trọng số cho việc chèn gap wk ≤ kw1 với độ phức tạp thời gian là O(n3). Lý do của việc tăng độ phức tạp về thời gian là do việc bổ sung thêm việc tính toán chi phí chèn – xóa gap trong các trường hợp. Công thức được đưa ra: (2) Trong đó P[i][j] và Q[i][j] là chi phí chèn và xóa ở vị trí ( i , j) Trong các trường hợp đặc biệt, chi phí chèn gap là một hàm tuyến tính wk = uk +v trong đó v được gọi là chi phí bắt đầu một đoạn gap và v là chi phí mở rộng đoạn gap. Gotoh (1982) [9] đã đưa ra một công thức tính toán tối ưu hóa việc tính toán ma trận P và Q giảm độ phức tạp thời gian xuống còn O(n2). Công thức mà Gotoh đưa ra là : (3) Sắp hàng trình tự hệ gen Trong quá trình tiến hóa của các sinh vật, bên cạnh những biến đổi ở mức độ điểm (sự thay thế chèn – xóa của các nucleotide) còn có những sự biến đổi ở mức độ gen. Có 3 phép biến đổi chính ở mức độ gen là phép chèn gen, xóa gen và dịch chuyển gen. Hình 2 mô tả một ví dụ về sự biến đổi ở mức độ gen giữa Người và Khỉ. Ta thấy gen số 1 đã bị dịch chuyển, nó nằm ở đầu của hệ gen Người nhưng lại nằm ở cuối ở hệ gen của Khỉ. Ngoài ra, gen số 2 tồn tại ở Khỉ nhưng không tồn tại ở Người. Tức là hoặc nó bị xóa khỏi hệ gen của Người hoặc nó được chèn thêm vào hệ gen của Khỉ. Do ta không phân biệt được phép chèn gen, và xóa gen, ta gọi chung là phép chèn/xóa gen. Trải qua hàng triệu năm tiến hóa, với sự biến đổi ở mức độ gen, hệ gen của các sinh vật ngày nay đã có sự khác nhau rất lớn về kích thước, số lượng gen, thứ tự các gen cũng như về nội dung của các gen. Hình 2: Các biến đổi ở mức độ gen giữa Người và Khỉ Sắp hàng trình tự hệ gen là một trường hợp riêng của sắp hàng trình tự, trong đó đầu vào là toàn bộ trình tự ADN của một hệ gen sinh vật. Sắp hàng trình tự hệ gen giúp xây dựng bức tranh toàn cảnh về sự tương tự và tiến hóa giữa các sinh vật, là cơ sở cho hướng nghiên cứu Comparative genomics [4], cho phép nâng cao độ chính xác dự đoán gen. Về mặt tính toán, bắt cặp hệ gen đặt ra nhiều vấn đề cần giải quyết như kích thước trình tự lớn, thứ tự các phần tương đồng trên các hệ gen thường thay đổi. Do tính quan trọng cũng như đặc thù phương pháp, vấn đề so sánh và sắp hàng trình tự hệ gen được trình bày thành một phần riêng, tách khỏi sắp hàng trình tự nói chung. Các thuật toán sắp hàng trình tự thông thường mới chỉ xác định được các biến đổi ở mức độ điểm (sự biến đổi của các nucleotide) cũng như chỉ làm việc được với các dữ liệu nhỏ. Khi nghiên cứu về việc sắp hàng trình tự theo hệ gen, chúng ta phải tính toán cả những biến đổi ở mức độ điểm lẫn mức độ gen. Đặc biệt thời gian thực thi cũng là một vấn đề hết sức quan trọng do kích thước rất lớn của các hệ gen. Ví dụ kích thước của hệ gen người lên tới 3 tỉ ADN. Một trong những hệ thống sắp hàng hệ gen đầu tiên là BLASTZ [18] được phát triển bới nhóm của Webb Miller vào đầu những năm 2000 tại đại học Pennsylvania để sắp hàng hệ gen của người và chuột. Cũng như các phương pháp sắp hàng hệ gen khác, Phương pháp BLASTZ được phát triển từ tư tưởng thuật toán tìm kiếm BLAST [2] (thuật toán xác định những đoạn giống nhau cao giữa hai chuỗi). Tư tưởng chung của thuật toán gồm ba bước: Bước 1: Tìm kiếm những cặp đoạn ADN ngắn rất giống nhau ở cả hai hệ gen được gọi là hạt giống (seed). Những đoạn này có độ dài vào khoảng 7 đến 13 ADN và được gọi là seed. Để thực hiện việc tìm kiếm này, có thể sử dụng nhiều kỹ thuật khác nhau như bảng băm, cây hậu tố (suffix tree). Bước 2: Mở rộng các hạt giống về cả hai phía sao cho trong quá trình mở rộng chi phí không vượt qua một ngưỡng cho trước. Quá trình mở rộng này không cho phép chèn gap Bước 3: Tiến hành nối các cặp ADN được mở rộng ở bước 2 lại với nhau để tạo thành những cặp ADN lớn hơn, bước này được phép chèn thêm gap. Sau khi nối, các cặp ADN này sẽ được đánh giá độ tương đồng. Các nghiên cứu hiện tại tập trung vào cải tiến bước thứ 1 và bước thứ 3. Nổi bật là các nghiên cứu của Aaron Darling và đồng nghiệp tại đại học Wisconsin–Madison trong việc cải tiến cách xác định các hạt giống ở bước 1. Họ định nghĩa hạt giống là những cặp ADN giống nhau và xuất hiện duy nhất trên cả hệ gen. Nhóm tác giả đã xây dựng hệ thống MAUVE để sắp hàng đa hệ gen và thu được những kết quả có độ chính xác cao trên những hệ gen có độ tương đồng cao [1]. Bên cạnh đó, nhóm tác giả Michael Brudno tại đại học Standford tập trung vào cải tiến bước 3 để kết nối các đoạn ADN và phát triển hệ thống SLAGAN [5]. Nhóm tác giả áp dụng phương pháp quy hoạch động để tìm ra cách kết nối các đoạn ADN tốt nhất, trong đó cho phép các đoạn ADN được phép dịch chuyển và đảo chiều. Kết quả so sánh hai hệ thống MAUVE và SLAGAN cho thấy MAUVE tốt hơn SLAGAN trên những tập dữ liệu có độ tương đồng cao, còn SLAGAN cho kết quả tốt hơn MAUVE trên những tập dữ liệu tồn tại nhiều phép thay thế ADN ở mức độ điểm và ít phép đảo chiều đoạn ADN ở mức độ gen. Mặc dù một số phương pháp đã được nghiên cứu và phát triển, chúng mới chỉ tập trung vào xác định và bắt cặp cho những vùng ADN có độ tương đồng cao giữa hai hệ gen. Tức là, một phần lớn trong hệ gen có thể không được bắt cặp và so sánh khi ta tiến hành với các loài sinh vật có hệ gen khác nhau nhiều. Để giải quyết vấn đề trên, những nghiên cứu đầu tiên của TS. Lê Sỹ Vinh và đồng nghiệp tại Bảo Tàng Lịch Sử Tự Nhiên Hoa Kỳ, và tại trường Đại Học Công Nghệ nhằm so sánh và sắp hàng toàn bộ hệ gen đã được tiến hành và cho kết quả thử nghiệm khả quan [23,24]. Nhóm nghiên cứu định nghĩa việc sắp hàng toàn bộ hệ gen phải thỏa mãn ba điều kiện chính sau: Xác định được các phép biến đổi ở mức độ gen (chèn, xóa, dịch chuyển vị trí). Xác định được các phép biến đổi ở mức độ điểm (thay thế, chèn, xóa). Bắt cặp toàn bộ các ADN trên hệ gen. Hệ thống bắt cặp thỏa mãn ba điều kiện trên sẽ cho phép bắt cặp các gen với các mức độ tương đồng khác nhau. Để đáp ứng được ba yêu cầu trên, Vinh và các đồng nghiệp đã nghiên cứu cách kết hợp điểm phạt cho các phép biến ở mức độ điểm, và các phép biến đổi ở mức độ gen vào thành một hệ thống tính điểm phạt chung. Điều này cho phép chúng ta xây dựng hàm tục tiêu rõ ràng để tìm ra cách bắt cặp toàn bộ hệ gen tốt nhất. Kết quả thí nghiệm với 760 bộ gen ty thể của các loài động vật cho thấy hệ thống tính điểm cho kết quả tốt [23]. Sử dụng phương pháp bắt cặp toàn bộ hệ gen, nhóm tác giả đã xây dựng quá trình tiến hóa của 11 Corona vi rút và tái khẳng định lại kết luận vi rút Corona gây ra dịch bệnh hô hấp cấp (SARs) có chung nguồn gốc với vi rút Corona ở loài dơi chứ không phải là loài chồn hôi (canivor) [24]. Chương 2. Bài toán sắp hàng hoàn chỉnh hai hệ gen 2.1. Tổng quan Theo những nghiên cứu của TS. Lê Sỹ Vinh và đồng nghiệp [23,24], một hệ thống sắp hàng hoàn chỉnh hai hệ gen phải thỏa mãn ba điều kiện chính : Xác định được các phép biến đổi ở mức độ gen (chèn, xóa, dịch chuyển vị trí). Xác định được các phép biến đổi ở mức độ điểm (thay thế, chèn, xóa). Bắt cặp toàn bộ các ADN trên hệ gen. Các phương pháp bắt cặp hệ gen tiêu biểu như BLASTZ[18], SLAGAN[5], MAUVE[1] mới chỉ dừng lại ở mức phát hiện và sắp hàng các đoạn ADN tương đồng. Như vậy sẽ có một phần lớn trong hệ gen có thể không được bắt cặp và so sánh khi ta tiến hành với các loài sinh vật có hệ gen khác nhau nhiều. Giải quyết vấn đề này, Lê Sỹ Vinh và các đồng nghiệp đã giới thiệu một phương pháp có khả năng sắp hàng hoàn toàn được hệ gen của hai sinh vật bất kỳ “Pairwise Alignment with Rearrangement” [23]. 2.2 Pairwise Alignment with Rearrangement “Pairwise Align with rearrangement” là thuật toán sắp hàng trình tự trong đó cho phép có sự sắp xếp lại của các ký tự trong trình tự do Lê Sỹ Vinh và các đồng nghiệp tại Bảo tàng lịch sử tự nhiên Hoa Kỳ đưa ra vào năm 2006[23]. Ưu điểm của thuật toán này so với các thuật toán bắt cặp trình tự trước đây ở chỗ nó cho phép có sự di chuyển vị trí của các ký tự. Đầu vào hai trình tự, “Pairwise Alignment with Rearrangement” sẽ đưa ra cách sắp hàng sao cho tổng 3 chi phí : Chi phí thay đổi vị trí (Break Cost), chi phí chèn – xóa và chi phí thay thế các ký tự là nhỏ nhất. Thực nghiệm cho thấy, việc sắp hàng trình tự có sự đổi chỗ của các ký tự cho kết quả tối hơn so với cắt bắt cặp trình tự thông thường không có sự đổi chỗ. 2.2.1. Cơ sở lý thuyết Trong “Pairwise Alignment with rearrangement”, ta xem hai hệ gen như hai chuỗi ký tự, tức là mỗi đoạn gen sẽ được xem như là một ký tự trong chuỗi đầu vào. Có X = (x1, x2, …, xp) là một chuỗi gồm p ký tự, Y = (y1, y2, …, yq ) là một chuỗi gồm q ký tự. C(xi, yj) là chi phí để thay thế ký tự xi thành ký tự yj với i = 1 … p, j = 1 ... q. C(xi, y0) và C(x0, yj) là chi phí chèn/xóa kí tự xi và yj tương ứng. Khi thực hiện với hai hệ gen, ta có xi và yj là một chuỗi các nucleotide, khi đó chi phí C(xi, yj) là chi phí nhỏ nhất để biến đổi xi thành yj. Gọi R(Y, YR) là hàm chi phí chuyển đổi giữa Y và một hoán vị YR của nó. Thông thường R(Y, YR) được tính là khoảng cách breakpoint[11] hoặc khoảng cách inversion[17]. Một chuỗi X’ = (x1’, x2’, …, xk’) được gọi là một chuỗi phát triển (edited sequence) từ X khi và chỉ khi X thu được từ X’ sau khi xóa hết các ký tự gap. Ví dụ X’ = (‘-‘, 1, 2, ‘-‘, 3, 4) là một chuỗi phát triển từ X = (1,2,3,4) . Một cặp A(X,Y) = A(X’,Y’) của hai chuỗi X’ = (x1’, x2’, …, xk’) phát triển từ X và Y’ = (y1’, y2’, …, yk’) phát triển từ Y được gọi là một bắt cặp hoàn chỉnh của X và Y. Chi phí C(A) của một bắt cặp A là tổng chi phí thay thế và chèn/xóa của các ký tự trong X và Y. (4) Chi phí tối ưu để bắt cặp A*(X,Y) = argminA(X,Y){C(A)} có thể được tính với độ phức tạp thời gian là O(pq) sử dụng kỹ thuật quy hoạch động (Xem phần 1.2). Một cách bắt cặp được xây dựng bằng cách chèn thêm ký tự gap ở cả hai chuỗi sao cho thứ tự các ký tự trong chuỗi phải được giữ nguyên. Một chuỗi XR’ = (x1’, x2’, …, xk’) được gọi là một chuỗi phát triển có sắp xếp lại (edited rearrangement sequence) từ X nếu sau khi loại bỏ gap ở XR’ ta thu được XR là một hoán vị của X. Ví dụ với XR’ = (‘-‘, 1, 4, 2, ‘-‘, 3) là một chuỗi phát triển có sắp xếp lại từ X = (1,2,3,4). Một cặp A = (XR’, YR’) của hai chuỗi phát triển có sắp xếp lại XR’ = (x1’, x2’, …, xk’) và YR’ = (y1’, y2’, …, yk’) được gọi là một bắt cặp trình tự có sắp xếp lại (PAR) của hai chuỗi X và Y. Chi phí CR(AR) của PAR AR là tổng các chi phí thay thế giữa các ký tự, chi phí chèn - xóa gap và chi phí sắp xếp lại giữa các ký tự. Ta có công thức: (5) Mục đích của bài toàn là tìm một PAR AR* có chi phí bé nhất. Tức là: (6) 2.2.2. Thuật toán Thuật toán “Pairwise Alignment with Rearrangement” sử dụng chiến lược leo đồi để tìm ra cặp PAR AR*. Chiến lược này gồm 2 bước. Bước đầu tiên một PAR AR xuất phát sẽ được tạo ra. Sau đó ở bước thứ 2 chúng ta tìm ra PAR AR* bằng cách lần lượt tìm ra các cặp AR tối ưu hơn. Đầu tiên, một PAR AR xuất phát sẽ được khởi tạo bằng cách sử dụng thuật toán “Stepwise addition”. Đầu tiên chúng ta khởi đầu với một PAR chưa đầy đủ là AR = (XR’ = X, YR’ = ). Sau đó ta lần lượt chèn các ký tự yj Y, j = 1 … |Y| vào YR’ để tạo thành một PAR hoàn chỉnh. Ký tự yj được thêm vào vị trí sao cho chi phí của AR mới là thấp nhất. Thuật toán được mô tả như sau : Stepwise Addition Method YR ← Ø for j = 1 to |Y| BestCost ← +∞, BestPosition ← nil for each position p in YR do Insert yj into YR at position p AR ← A*(X, YR) if CR(AR) < BestCost then BestPosition ← p, BestCost ← CR(AR) Remove yj from YR end for Insert yj into YR at BestPosition end for Return AR = A*(X, YR) Bước thứ 2, chúng ta tìm cách từng bước một tối ưu hóa PAR AR hiện có, tức là tìm cách giảm chi phí CR(AR). Để thay đổi một PAR, ta tiến hành đổi vị trí giữa các ký tự trong chúng. Với 1 PAR AR(XR’, YR’) ta định nghĩa một phép biến đổi M(i, j, t | i ≤ j ≤ t-1) của chuỗi YR = (y1, y2, …, yp) thu được sau khi loại bỏ ký tự gap từ YR là phép dịch chuyển đoạn (yi,…,yj) trong YR đến vị trí t để thu được một chuỗi không có ký tự gap mới YR. Một phép biến đổi M(i, j, t) được gọi là có thể thực hiện được (possible move) nếu chi phí CR của PAR mới tạo được từ XR và YR tốt hơn chi phí của PAR AR cũ. Thuật toán được mô tả như sau: Character Moving Method Build an initial PAR AR = (XR’, YR’) by stepwise addition method iteration ← 0 repeat positionMove ← false foreach triple positions (i, j, t | i ≤ j < t-1) if M(i, j, t) is a possible move then Move character(yi, …, yj) in YR to position t possibleMove ← true end if end for iteration ← iteration + 1 until possibleMove = false return A*(XR, YR) Để tối ưu hóa một PAR AR, một thuật toán khác phức tạp hơn được để xuất, có tên là “Simulaneous Character Swapping". Trong thuật toán này, 1 biến đổi S(k, t) được định nghĩa là phép đổi vị trí yi và yj của chuỗi YR thu được sau khi loại bỏ ký tự gap từ YR’. Một phép biến đổi được gọi có thể thực hiện được nếu chi phí của AR mới giảm đi. Xét hai phép biến đổi S1(k1, t1) và S2(k2, t2) với k1 t1 hoặc t2 < t1. Thực nghiệm cho thấy nếu ta đồng thời thực hiện hai phép biến đổi S1 và S2 độc lập với nhau có thể sẽ tạo ra 1 PAR tốt hơn. Ngược lại nếu S1 và S2 không độc lập, chỉ thực hiện phép biến đổi cho chi phí thấp hơn. Hình 3:Hình trái: ví dụ 2 phép biến đổi S1(k1, t1) và S2(k2, t2) là độc lập với nhau. Hình phải:Đổi chỗ đồng thời 2 phép biến đổi độc lập. Thuật toán “Simulaneous Character Swapping” được mô tả như sau: Simulaneous Character Swapping Method Build an initial PAR AR = (XR’, YR’) by stepwise addition method iteration ← 0 bestCost ← 0 repeat positionSwap ← false Find independent possible swap L and sort L increasing ly if |L| > 0 then r ← L repeat swap simultaneous L1, …, LR in YR to get ŸR ÄR ← A*(XR, ŸR) if CR(ÄR) < bestCost then possibleSwap ← true bestCost ← CR(ÄR), YR ← ŸR else r = r/2 until possibleMove = false iteration ← iteration + 1 until possibleSwap = true return A*(XR, YR) 2.2.3. Độ phức tạp của thuật toán Như ta đã biết thuật toán “Pairwise Alignment with Rearrangement” bao gồm hai bước. Bước một là tìm một PAR xuất phát bằng cách sử dụng phương thức “Stepwise addtion”. Bước hai là từng bước tối huy hóa PAR hiện có bằng cách một trong hai cách : “Character moving” hoặc “Simulaneous Character Swapping”. Ở bước 1, nhìn vào phương thức “Stepwise addtion”, ta có thể thấy từ dòng 5 đến dòng 9, độ phức tạp thời gian yêu cầu là O(pq) với p, q là độ dài của hai chuối X và Y, độ phức tạp thời gian yêu cầu từ dòng 4 cho đến dòng 10 sẽ là O(pq2). Do vậy yêu cầu về độ phức tạp về thời gian trong bước này là O(pq3). Ở bước 2, nếu ta sử dụng phương thức “Character Moving”, ta thấy yêu cầu vê độ phức tạp thời gian từ dòng 5 đến dòng 10 sẽ là O(pq4). Vì vậy yêu cầu về độ phức tạp thời gian nên sử dụng phương thứ này sẽ là O(pq4 x iteration) với iteration là số vòng lặp. Ngược lại nếu ta sử dụng theo phương thức “Simulaneous character swapping “, yêu cầu độ phức tạp thời gian trong trường hợp xấu nhất từ dòng 6 đến dòng 16 chỉ còn là O(pq3), độ phức tạp thời gian yêu cầu cho cả bài toán là O(pq3) x iteration với iteration là số vòng lặp. 2.3. Bắt cặp với những trình tự lớn Như đã trình bày ở trên, thuật toán “Pairwise Alignment with Rearrangement” yêu cầu về độ phức tạp thời gian tối thiếu là O(pq3) trong mỗi vòng lặp. Điều này vẫn đến thuật toán gặp khó khăn khi áp dụng với những trình tự lớn. Nghiên cứu của Đinh Ngọc Thắng đã đề xuất một phương thức mới : “Fast Swapping”, giúp giảm độ phức tạp thời gian yêu cầu trong mỗi vòng lặp xuống còn O(pq) [20]. Từ công thức (5) ta có chi phí CR(AR) bao gồm hai phần. Chi phí sắp hàng giữa các ký tự (chèn, xóa, sửa) và chi phí sắp xêp lại thứ tự. Chi phí sắp xếp lại thức tự R(X, XR) + R(Y, YR) là một giá trị không âm và chỉ bằng 0 khi và chỉ khi sau khi loại bỏ các ký tự gap ở XR và YR ta thu được X và Y. Chi phí sắp hàng giữa các ký tự C(X’, Y’) = có thể được tính dễ dàng theo phương thức quy hoạch động với độ phức tạp về thời gian là O(pq) (Xem phần 1.2). Xét một phép đổi chỗ yi’ và yj’, chi phí sắp hàng có thể dễ dàng được tính lại với độ phức tạp thời gian yêu cầu là O(1) là : C(X’, Y’) – [C() + C()] + [C() + C()] (7) Trong “Fast Swapping” ta sử dụng 2 phép đổi chỗ để tối ưu hóa chi phí CR(AR), trong đó cách đầu tiên (Single Swap) được dùng để tối ưu chi phí sửa chữa, cách thứ 2 (Couple Swap) được sử dụng để tối ưu hóa chi phí sắp xếp lại. Hình 4: Single Swap (trái) và Couple Swap (phải) Single Swap là cách đổi chỗ hai ký tự bất kì của X’ hoặc Y’ trong khi Couple Swap đổi chỗ các ký tự ở cả hai chuỗi cùng 1 lúc. Tức là Couple Swap chính là thực hiện 2 phép đổi chỗ Single Swap cùng lúc trên cả hai chuỗi. Tuy nhiên tác dụng của chúng là hoàn toán khác nhau. Trong khi Single Swap chủ yếu làm thay đổi chi phí sửa chữa giữa các ký tự thì Couple Swap thay đổi chi phí sắp xếp lại, trong khi chi phí sửa chữa giữa các ký tự vẫn được giữ nguyên. Một phép đổi chỗ được gọi là có thể thực hiện nếu nó làm cho chi phí CR(AR) giảm xuống. Trong thủ tục dưới đây, đầu vào sẽ là một PAR xuất phát, chừng nào còn tìm được một cách đổi chỗ có thể thực hiện thì ta sẽ tiến hành đổi chỗ để tìm ra PAR mới tối uy hơn, cập nhật lại PAR. Quá trình này dừng lại khi không có một phép đổi chỗ có thể thực hiện nào được tìm thấy. Fast Swapping Method iteration ← 0 while (true) do iteration ← iteration + 1 for i = 1 to |XR’| do for j = 0 to |YR’| do if  have possible swap S in {SX’(i,j), SY’(i,j), CS(i,j)} then Perform S if no swap performed then break out of while loop XR ← gapless(XR’), YR ← gapless(XR’) (XR’, YR’) ← PairwiseAlignment(XR, YR) end while return A*(XR, YR) Ta có một phép đổi chỗ chỉ yêu cầu thời gian là O(1). Như vậy nhìn vào chương trình ta thấy từ dòng 4 đến dòng 7 độ phức tạp thời gian yêu cầu chỉ là O(pq). Dòng thứ 11 khởi động thuật toán quy hoạch động để sắp hàng hai chuỗi XR, YR cũng chỉ cần yêu cầu thời gian là O(pq) (Xem phần 1.2) . Như vậy độ phức tạp thời gian áp dụng phương thức này sẽ là O(pq x iteration) với iteration là số lần chạy vòng while. Chương 3. Full Genome Alignment 3.1. Xây dựng hệ thống Thuật toán “Pairwise Alignment with Rearrangement” [23] tuy đã sắp hàng được hoàn toàn hai hệ gen, tuy nhiên nhược điểm của nó chỉ có thể bắt cặp với những hệ gen đã được chia sẵn. Vấn đề đặt ra là khi đưa vào một hệ gen mới bất kỳ, chúng ta phải tiến hành xác định những đoạn gen trong đó để tiến hành sắp hàng các đoạn gen. Điều này đòi hỏi phải khi đưa vào hai hệ gen của hai sinh vật bất kỳ, phải tìm cách chia các hệ gen thành nhiều đoạn gen con liên tiếp sao cho khi sắp hàng với nhau chúng tạo được những cặp gen có độ tương đồng cao, tức là những cặp gen có nhiều khả năng là được biến đổi và tiến hoá từ cùng một đoạn gen trong hệ gen tổ tiên chung xa xưa của chúng. Như đã trình bày ở phần đầu, một số chương trình sắp hàng hệ gen hiện có tập trung vào việc tìm kiếm những vùng gen tương đồng trên hai hệ gen như MUAVE[1], SLAGAN[5] và nổi bật là BLASTZ[18]. Với những cải tiến độc đáo của mình, BLASTZ có khả năng tìm kiếm những vùng gen có độ tương đồng cao một cách tương đối tốt và với thời gian chấp nhận. Trong khóa luận này, em xây dựng lên một hệ thống bắt cặp hệ gen hoàn chỉnh, dựa trên ý tưởng sử dụng BLASTZ như một bước tiền xử lý trước khi áp dụng thuật toán “Pairwise Alignment with Rearrangement”, qua đó khắc phục được nhược điểm hiện có của phương pháp này. Cụ thể chúng ta sẽ sử dụng những vùng tương đồng mà BLASTZ nhận dạng được để tiến hành chia cắt hệ gen thành những đoạn gen ngắn liên tiếp, tạo thành đầu vào cho chương trình “Pairwise Alignment with Rearrangement” với “Fast Swapping” [20], trong quá trình này vấn đề cần lưu ý là phải tiến hành loại bỏ các đoạn trùng lặp, lựa chọn và giữ lại các cặp gen có độ tương đồng cao hơn. Để xác định các biến đổi về điểm cũng như tính toán khoảng cách giữa các đoạn gen. Thuật toán “Pairwise Alignment with Rearrangement”, sử dụng phương pháp bắt cặp trình tự theo thuật toán quy hoạch động của Needleman – Wunsch [16]( Xem phần 1.2). Trong hệ thống này, em sẽ sử dụng thay thế bằng một phương pháp khác sử dụng thuật toán quy hoạch động được đưa ra bởi Gotoh “Optimal Alignment with Linear space” [9] và có sự cải tiến dựa trên nhận xét của Ukkonen[22]. Phương pháp mới này cho phép tính toán khoảng cách giữa các đoạn gen chính xác và hợp lý hơn so với phương pháp cũ của Needleman – Wunsch. Chương trình gồm những bước cụ thể như sau: Đầu vào là hai hệ gen hoàn chỉnh bất kỳ. Sử dụng chương trình BLASTZ để xác định và bắt cặp những vùng ADN tương đồng. Tiến hành tách từng hệ gen thành một dãy các đoạn ADN thành phần nhỏ liên tục dựa vào các vùng có độ tương đồng cao xác định được bởi BLASTZ. Dựa trên thuật toán Gotoh và nhận xét của Ukkonen, xây dựng chương chình bắt cặp trình tự. Áp dụng để sắp hàng từng cặp ADN thành phần trong hai hệ gen, xác định các phép biến đổi ở mức độ điểm (thay thế, chèn, xóa). Sắp hàng toàn bộ hệ gen, xác định các biến đổi ở mức độ gen bằng thuật toán Pairwise Alignment with Rearrangement với Fast Swapping.. Đầu ra đưa ra danh sách nhưng cặp gen đã được sắp hàng, trong đó chỉ rõ những sự biến đổi ở mức độ điểm ở từng cặp gen. Cho biết thông tin về các đoạn gen đã bị dịch chuyển, bị đảo ngược, tồn tại ở hệ gen này nhưng không tồn tại ở hệ gen kia. Sắp hàng hoàn chỉnh hai hệ gen. 3.2. Giới thiệu về BLASTZ BLASTZ là phương pháp tìm kiếm các đoạn ADN tương đồng được phát triển bởi nhóm Miller thuộc trường Đại học Pennsylvania Hoa Kỳ. Nó được áp dụng thành công trong việc bắt cặp giữa hai hệ Gen Người và Chuột [18] 3.2.1. Tính năng của BLASTZ BLASTZ sử dụng 3 chiến lược đã được sử dụng bởi Gapped BLAST [3] đó là: Tìm kiếm những cặp đoạn ngắn ADN rất giống nhau ở cả hai hệ gen được gọi là hạt giống (seed). Mở rộng các hạt giống về cả hai phía sao cho trong quá trình mở rộng chi phí không vượt qua một ngưỡng cho trước. Quá trình mở rộng này không cho phép chèn gap Tiến hành tiếp tục mở rộng các cặp ADN ở bước 2 lại với nhau để tạo ra những cặp ADN lớn hơn bằng cách cho phép chèn thêm gap. Việc mở rộng đảm bảo chi phí không vượt quá một ngưỡng nhất định. Tuy nhiên so với Gapped BLAST và các chương trình sắp hàng hệ gen khác, BLASTZ có những ba sự cải tiến quan trọng. Trước tiên, BLASTZ sử dụng một cách tính điểm bắt cặp được đánh giá bởi Chiromonte [6]. Theo đó thay vì chi phí sắp hàng gồm chi phí thay thế và chi phí sắp hàng chính xác chỉ là một giá trị chung với tất cả các nucleotide thì trong BLASTZ chi phí sắp hàng các nucleotide được cho bởi ma trận sau : A C G T A 91 -114 -31 -123 C -114 100 -125 -31 G -31 -125 100 -114 T -123 -31 -114 91 Bảng 1: Ma trận trọng số của BLASTZ Chi phí chèn – xóa các ký tự gap được cho bởi một hàm tuyến tính. Việc chèn – xóa k ký tự gap liên tiếp sẽ phải chụi một điểm phạt là 400 + 30k. Hai thay đổi tiếp theo giúp cải tiến đáng kể tốc độ thực hiện và độ nhay của BLASTZ trong việc bắt cặp toàn bộ bộ gen. Thứ nhất là việc loại bỏ những đoạn trùng lặp. Ví dụ khi chương trình nhận ra rằng nhiều khu vực trong một bộ gen của chuột được sắp hàng với cùng một phân khúc trong bộ gen của người, chương trình sẽ tự động đánh dấu để nó được bỏ qua trong các bước sau của quá trình bắt cặp. Cải tiến này giúp BLASTZ không bắt cặp những đoạn ADN trùng nhau – những đoạn ADN có thể được nhân lên trong quá trình biến đổi và tiến hóa. Thứ hai BLASTZ áp dụng một ý tưởng thông minh của Ma [15] trong việc xác định các đoạn ngắn gần giống nhau ban đầu (seed). Ma đề xuất việc tìm kiếm trong 19 nucleotide liên tiếp, trong đó 12 nucleotide được chỉ định bằng 1 trong chuỗi 1110100110010101111 là giống hệt nhau. Để tăng độ nhạy, BLASTZ còn cho phép 1 vị trí bất kỳ trong 12 vị trí ở trên được phép có một sự thay thế giữa các cặp nucleotide tương đồng (A – G, G – A, C – T, T – C) 3.2.2. Chương trình BLASTZ Chương trình BLASTZ được cài đặt theo năm bước chính. Bước 1 thực hiện tìm kiếm và đánh dấu các đoạn giống nhau bị lặp lại ở cả hai trình tự (Repeat Finding) Bước 2 tiến hành tìm kiếm các cặp hạt giống (seed) có độ tương đồng cao ở cả hai hệ gen, trong bước này BLAST sử dụng độ tương đồng 12of19 được đề xuất bởi Ma[15]. Ngoài ra cho phép có 1 sự thay thế giữa các cặp nucleotide ở 1 vị trí bất kỳ trong 12 vị trí ở trên. Bước 3 BLASTZ sẽ tiến hành mở rộng các cặp seed tìm được. Quá trình mở rộng này được thực hiện như sau: Lần lượt mở rộng các cặp seed về hai phía, không cho phép chèn gap. Quá trình mở rộng dừng lại khi chi phí phạt vượt quá một ngưỡng cho phép (X). Nếu điểm số sắp hàng cặp ADN thu được đạt mức K cho trước. Tiếp tục tiến hành mở rộng cặp ADN lần lượt về 2 phía với việc cho phép chèn gap. Quá trình mở rộng tiến hành cho đến khi chi phí phạt vẫn chưa vượt quá ngưỡng cho phép (Y) Giữ lại các đoạn ADN có điểm số sắp hàng đạt một mức cho trước (L). Bước 4 là thực hiện ghi nhận lại vị trí các đoạn tương đông tìm được sao cho phù hợp với 2 hệ gen ban đầu. Bước 5 là tiến hành điều chỉnh lại kết quả cho phù hợp với các tùy chọn. Cuồi cùng, BLASTZ có 1 tùy chọn cho phép việc đảo ngược một hệ gen rồi tiến hành sắp hành lại với hệ gen còn lại. Việc này cho phép tìm kiếm các đoạn tương đồng trong trường hợp 1 đoạn gen bị đảo ngược. BLASTZ xây dưng một hệ thống tùy chọn, cho phép người dùng thay đổi các tham số của chương trình phù hợp với mục đích của người sử dụng. 3.3. Optimal Alignment with Linear space “Pairwise Alignment with Rearrangement” sử dụng thuật toán quy hoạch động của Needleman – Wunsch [16] (Xem phần 1.2) để tính toán trọng số khoảng cách giữa các đoạn gen đồng thời xác định những biến đổi về điểm. Thuật toán “Pairwise Alignment” của Needleman – Wunsch có những hạn chế nhất định khi nó chỉ có thể làm việc khi mà chi phí chèn – xóa gap là một trọng số cố định. Trên thực tế trong quá trình tiến hóa khi xóa bỏ những đoạn ADN liên tiếp, việc xóa nucleotide đầu tiên bao giờ cũng khó khăn hơn rất nhiều so với các nucleotide tiếp theo. Do đó hàm chi phí cho việc xóa – chèn những đoạn gap liên tiếp bao giờ cũng là một hàm tuyết tính w(k) = a + bk cho việc xóa k nucleotide liên tiếp trong đó w(k) < kw(1) Vì vậy, trong chương trình mới, chúng ta tiến hành thay thế thuật toán “Pairwise Alignment” đơn giản của Needleman – Wunsch bằng thuật toán “Optimal Alignment with Linear space” của Gotoh[9]. Trong thuật toán này Gotod đưa ra các định nghĩa sau: dAB( Ai , Bj ) là chi phí bắt cặp giữa hai đoạn Ai và Bj trong đó Ai được sắp hàng với Bj. dA- ( Ai , Bj ) là chi phí bắt cặp giữa hai đoạn Ai và Bj trong đó Ai được sắp hàng với kí tự gap d-B ( Ai , Bj ) là chi phí bắt cặp giữa hai đoạn Ai và Bj trong đó Bj được sắp hàng với kí tự gap Với w(Ai, Bj) là chi phí khi sắp hàng ký tự Ai và ký tự Bj. w(k) = a +bk là chi phi chèn – xóa k ký tự ta có công thức quy hoạch động: dAB( Ai, Bj ) = min (dAB( Ai-1, Bj-1 ), dA-( Ai-1, Bj-1 ), d-B( Ai-1, Bj-1 ) ) + + w( Ai, Bj ) dA- ( Ai, Bj ) = min(dAB( Ai-1, Bj ) + a, dA-(Ai-1, Bj ), d-B( Ai-1, Bj) + a)+ b d-B ( Ai,Bj ) = min(dAB( Ai, Bj-1 ) + a, dA-(Ai, Bj-1) + a, d-B( Ai, Bj-1 ) )+ b Chi phí tối ưu để bắt cặp hai trình tự A và B là giá trị nhỏ nhất trong ba giá trí dAB(A, B), dA-(A, B) và d-B(A, B). Thuật toán Gotoh có độ phức tạp về thời gian là O(n2) và yêu cầu không gian bộ nhớ là O(n2). Do vậy có tồn tại một số khó khăn khi làm việc với những chuỗi có độ dài lớn. Trong chương trình của mình, em đưa thêm một số cải tiến để giải quyết vấn đề này. Thứ nhất, áp dụng nhận xét của Ukkonen cho việc bắt cặp những trình tự gần giống nhau, Ukkonen chỉ ra rằng khi bắt cặp những trình tự gần giống nhau, chỉ có khu vực quanh đường chéo chính là được sử dụng. [22] Do vậy chúng ta sử dụng thêm 1 barrier trong quá trình quy hoạch động để giảm thời gian thực hiện chương trình. Áp dụng nhận xét này cho thuật toán Gotoh, độ phức tạp giời gian giảm xuống O(dm) với d là khoảng cách độ dài giữa hai trình tự, m là độ dài của trình tự ngắn hơn. Hình 5:Sắp hàng trình tự với Ukkonen Barrier [13] Thứ hai, trong quá trình quy hoạch động theo thuật toán Gotod, giá trị của 1 hàng chỉ được tính dựa vào 1 hàng trước nó, do vậy ta có thể sử dụng 2 mảng một chiều để thay thế cho 1 hàng hai chiều. Như vậy có thể giảm không gian bộ nhớ xuống chỉ còn O(n) Chương trình cụ thể như sau : Optimal Alignment with Linear space 1 barrier = |X| - |Y| + pairwiseBarrier 2 for i=1 to |Y| do 3 d2n[0] = maxInt 4 d1n[0] = a + bi 5 d0n[0] = maxInt 7 u = max(1, i - barrier) 8 v = min( |X|, i + barrier) 9 if u > 1 then 10 d0n[u - 1] = maxInt 11 d1n[u - 1] = maxInt 12 d2n[u - 1] = maxInt 13 end if 14 if v<|X| then 15 d0n[v+1] = maxInt 16 d1n[v+1] = maxInt 17 d2n[v+1] = maxInt 18 end if 19 for j=u to v do 20 d0n[j] = min(d0p[j-1], d1p[j-1], d2p[j-1])+ w[Y[i]][X[j]] 21 d1n[j] = min(d0p[j]+a, d1p[j], d2p[j]+a)+b 22 d2n[j] = min(d0n[j-1]+a, d1n[j-1]+a, d2n[j-1])+b 23 end for 24 d0p = d0n 25 d1p = d1n 26 d2p = d2n 27 end for 28 return min(d0p[|X|, d1p[|X|, d2p[|X|) Chương 4. Kết quả 4.1. Chương trình Chương trình được cài đặt và sử dụng trên cả môi trường Windows và Linux. Đầu vào là toàn bộ trình tự ADN của một hệ gen của hai sinh vật. Dữ liệu được cho dưới định dạng chuẩn FASTA được lấy từ file hoặc dán thẳng vào chương trình. Kết quả được xuất dưới dạng file Text trong đó có thông tin về các cặp gen đã được bắt cặp toàn bộ, các biến đổi ở mức điểm (chèn – xóa, thay thế) cũng như các biến đổi ở mức độ gen (di chuyển, đảo ngược, xóa). Chương trình sử dụng một phần source code của chương trình BLASTZ được cung cấp tại Hình 6: Giao diên chương trình Các tùy chọn tham số được đưa vào gồm có : Ma trận trong số thay thế giữa các nucleotide (A, C, G, T) Chi phí chèn – xóa các đoạn Gap (Chi phí bắt đầu, mở rộng 1 đoạn Gap) Các tham số của BLASTZ : Giá trị ngưỡng K, L, X, Y (Xem phần 3.1.3. ) Chi phí Breakpoint Ngưỡng Barrier Kết quả thống kê số đoạn gen tách được ở hai hệ gen 1 và 2. Cho biết thông tin về các biến đổi ở mức độ gen: số đoạn gen bị đảo ngược, dịch chuyển, chèn – xóa, … Tiếp theo cho biết thông tin chi biết về các biến đổi ở mức độ điểm ở từng cặp gen: Số nucleotide bị thay thế, chèn, xóa, độ dài các đoạn gen, … Hình 7: Kết quả chương trình 4.2. Kiểm thử Tiến hành kiểm thử trên hai loại dữ liệu là dữ liệu mô phỏng và dữ liệu thật. Do chưa có một chương trình tương tự được công bố có khả năng so sánh và sắp hàng toàn bộ hai hệ gen nên ở phần này, việc dánh giá kết quả được thực hiện ở mức kiểm tra và so sánh độ tương đồng của hai hệ gen sau khi sắp hàng. 4.2.1. Dữ liệu mô phỏng Một Mitochondrial Genome bao gồm 2 sợi chính. 1 sợi giàu guanine thường bao gồm 13 đoạn polypeptide-encoding genes cấu thành 13 protein chính của sinh vật. Sợi thứ 2 giàu cytosine gồm 2 gene rRNA và 22 gen tRNA. [19] Thực hiện tạo dữ liệu mô phỏng bằng cách lấy 13 đoạn polypeptide-encoding gen ở sợi thứ nhất tạo thành một đoạn gen lớn. Tiến hành sinh các test khác nhau bằng cách chèn – xóa và thay thế ngẫu nhiên một số nucleotide, sau đó đảo ngược – di chuyển một số đoạn gen nhỏ trong đoạn gen ở trên ta được một đoạn gen mới. Trong quá trình này chúng ta lưu lại vị trí của từng nucleotide so với vị trí bạn đầu. Sử dụng chương trình tiến hành bắt cặp giữa hai đoạn gen rồi so sánh và đánh giá kết quả thu được từ chương trình. Kết quả được so sánh theo 2 tiêu chí, thứ nhất tỷ lệ số nucleotide sau khi sắp hàng được bắt cặp chính xác của với vị trí của chúng trong gen ban đầu. Thứ hai là tỷ lệ số nucleotide được sắp hàng chính xác so với nucleotide ở hệ gen ban đầu, tức là những cặp sắp hàng dạng : A – A, C – C, G – G, T – T Với mỗi cách chèn, xóa, đảo ngược ta thực hiện ngẫu nhiên 20 test, tổng cộng số test được thực hiện là 20 x 4 x 4 = 320 test. Tham số chương trình để mặc định. Máy tính CPU P8400 2x2.2GHz. Ta có kết quả thu được như sau: Inversion 0 Deletion 2% Deletion 5% Substitution 2% Match (by position) : 98.306% Match(by nucleotide): 99.023% Time : 5.2s Match (by position) : 98.203% Match(by nucleotide) : 98.983% Time : 6.65s Substitution 5% Match (by position) : 95.841% Match(by nucleotide): 97.470% Time : 5.6s Match (by position) : 95.548% Match(by nucleotide): 97.554% Time : 6.65s Bảng 2: Kết quả Test với số Inversion – Move là 0 Inversion 1 Deletion 2% Deletion 5% Substitution 2% Match (by position) : 91.733% Match(by nucleotide): 99.108% Time : 4.2s Match (by position) : 89.413% Match(by nucleotide): 98.920% Time : 4.25s Substitution 5% Match (by position) : 84.983% Match(by nucleotide): 97.769% Time : 3.85s Match (by position) : 90.884% Match(by nucleotide): 97.293% Time : 3.7s Bảng 3: Kết quả Test với số Inversion – Move là 1 Inversion 2 Deletion 2% Deletion 5% Substitution 2% Match (by position) : 92.231% Match(by nucleotide) : 98.942% Time : 3.95s Match (by position) : 79.552% Match(by nucleotide): 98.887 % Time : 4.2s Substitution 5% Match (by position) : 80.456% Match(by nucleotide) : 97.194% Time : 3.8s Match (by position) : 86.830% Match(by nucleotide): 97.379% Time : 3.9s Bảng 4: Kết quả Test với số Inversion – Move là 2 Inversion 3 Deletion 2% Deletion 5% Substitution 2% Match (by position) : 75.970% Match(by nucleotide) : 98.842% Time : 2.9s Match (by position) : 70.106% Match(by nucleotide) : 98.771% Time : 2.9s Substitution 5% Match (by position) : 77.388% Match(by nucleotide) : 97.273% Time : 3.05s Match (by position) : 70.414% Match(by nucleotide) : 97.206% Time : 2.9s Bảng 5: Kết quả Test với số Inversion – Move là 3 Nhìn vào kết quả thực nghiệm ta có thể thấy chương trình luôn xác định ra những đoạn gen có tỷ lệ tương đồng cao, tỷ lệ sắp hàng chính xác của theo từng nucleotide vẫn luôn đạt ở mức cao, lớn hơn 97%. Đối với những bộ test càng có nhiều sự biến đổi phức tạo, tức là có sự kiện đảo ngược kết hợp với di chuyển các đoạn gen, tỷ lệ sắp hàng đúng của từng nucleotide với vị trí cũ trong hệ gen gốc giảm dần. Điều này được lý giải do khi có những biến đổi phức tạp, cấu trúc của hệ gen bị thay đổi, cộng thêm với việc một số nucleotide trong hệ gen bị loại bỏ hoặc thay thế khỏi hệ gen. Trong quá trình tìm kiếm các đoạn gen tương đồng, nếu phát hiện nhiều cặp gen có độ tương đồng giống nhau, chương trình sẽ ưa tiên lựa chọn căp gen tìm thấy đầu tiên. Do đó dẫn tới việc vị trí của từng nucleotide bị sai lệch so với vị trí ban đầu mặc dù chúng vấn đước sắp hàng với cùng một loại nucleotide tương ứng. 4.2.2. Dữ liệu thật Sử dụng bộ dữ liệu Metazoa Mitochondira được lấy từ Gen Bank tại NCBI ( Lựa chọn 5 hệ Gen đặc trưng của các chủng người Việt (chủng da vàng – 2 bộ), người Uganda (chủng da đen), người Đức (chủng da trắng), người thổ dân Mỹ (chủng da đỏ) và 10 hệ Gen từ các loài sinh vật khác trên trái đất Khỉ, Cá, Vịt xiêm, Tôm, Gấu ngựa, Hải cẩu, Ếch, Kỳ nhông, Voi châu Á và rùa Biển. Mội hệ gen bao gồm 13 đoạn polypeptide-encoding genes, 2 đoạn rRNA gene là 12S và 16S rRNA genes và 22 đoạn tRNA genes. Dựa trên 37 đoạn gen chuẩn này ta tiến hành chia hệ gen thành nhiều đoạn con liên tiếp. Trung bình một hệ gen được chia thành khoảng 55 đoạn và có độ dài từ 16000 – 20000 nucleotide. Để kiểm nghiệm ta so sánh kết quả về điểm số đương đồng thu được từ chương trình sau khi đã chạy với hai hệ gen gốc (là tổng chi phí để sắp hàng các đoạn gen bao gồm chi phí sửa chữa ở mức độ điểm và chi phí sửa chữa ở mức độ gen) với kết quả khi bắt cặp hệ gen được chia sẵn các đoạn gen được chia dựa trên các đoạn gen nhỏ được chia sẵn trong Gen Bank sử dụng phương thức Kuhn-Munkres (Hungarian)[8]. Thuật toán (Hungarian) là thuật toán có độ phức tạp O(n3) được dùng để bắt cặp hai bộ dữ liệu có số dữ liệu bằng nhau và trọng số thay đổi vị trí bằng 0. Với 15 hệ gen đặc trưng của các chủng loài thu thập từ Gen Bank, ta thực hiện 105 test để bắt cặp đôi một cho từng hệ gen gốc, sau đó tính tỷ lệ giữa kết quả thu độ tương đồng thu được được từ chương trình với kết quả độ tương đồng thu được từ chương trình Hungarian chạy trên bộ gen đã chia sẵn. Tỷ lệ này càng cao càng thể hiện chương trình có khả năng tìm ra những độ tương đồng rất cao. Chương trình chạy với các tham số mặc định. Máy tính CPU P8400 2.2 GHz. Kết quả thu được như sau : Thời gian rung bình : 4.905s Tỷ lệ độ tương đồng thu : 100.698% Kết quả cụ thể trên một số test Gen 1 Gen 2 Thời gian Tỷ lệ Tôm chân trắng Người Uganda 8 111.36% Gấu ngựa Voi châu Á 4 105.04% Cá Người Đức 5 104.89% Hải cẩu Người thổ dân Mỹ 3 104.53% Khỉ sóc Người Việt 3 103.94% Vịt xiêm Rùa biển 4 100.42% Khỉ sóc Ếch 7 99.79% Ếch Người Việt 7 99.59% Vịt xiêm Ếch 6 98.86% Cá Tôm chân trắng 8 91.67% Tôm chân trắng Kỳ nhông 6 91.24% Bảng 6: Kết quả chạy dữ liệu thật Với kết quả thực nghiệm ta có thể thấy trong đa số trường hợp chương trình cho chi phí bắt cặp tương đương và tốt hơn so với Hungarian. Điều này là do chương trình xác định những cặp gen con tương đồng để chia nhỏ các hệ gen ra trước khi tiến hành bắt cặp một cách linh hoạt và hiệu quả hơn nhiều so với việc chia hệ gen cố định có sẵn. Trên thực tế với một số test, nếu thay đổi các ngưỡng các BLASTZ, số đoạn tương đồng nhận ra có thể tăng lên và các hệ gen được chia nhỏ thành nhiều đoạn con hơn, có thể nâng cao tỷ lệ độ tương đồng so với kết quả thu được từ thuật toán Hungarian khi chạy với các hệ gen được chia sắn bằng các đoạn gen con tại Gen Bank. Chương 5. Kết luận Với sự phát triển và tiến bộ nhanh chóng của công nghệ sinh học, số lượng gen được giải mã ngày càng nhiều và xây dưng một chương trình sắp hàng toàn bộ hệ gen nhanh chóng và chính xác là một điều hết sức cần thiết. Với kết quả của việc sắp hàng hai hệ gen, các nhà khoa học có thể dựa vào những những đoạn gen tương đồng để tìm ra chủng loại, phỏng đoán tính chất của những hệ gen mới. Điều này rất quan trọng trong việc phát hiện và nghiên cứu những chủng loại virut mới hay xác định nguồn gốc, mối liên hệ giữa các sinh vật trên hành tinh, cũng như trong các lĩnh vực như khảo cổ học hay công nghệ dược. Trong khóa luận này, em đã trình bày phương pháp và xây dựng một chương trình bắt cặp toàn bộ hệ gen, trên cơ sở kết hợp và cải tiến những phương pháp riêng lẻ đã từng được triển khai. Ba phương pháp chính được kết hợp trong chương trình này là BLASTZ[18], Optimal Pairwise Alignment with linear space[9] và Pairwise Alignment with Rearrangement[23]. Về cơ bản chương trình đã đạt được đầy đủ những yêu cầu của một hệ thống bắt cặp hệ gen cần có: Bắt cặp được toàn bộ hệ gen, chỉ ra được những biến đổi ở mức độ điểm (chèn – xóa, thay thế) cũng như những biến đổi ở mức độ gen (đảo ngược, dịch chuyển, chèn – xóa). Kết quả thực nghiệm trên hai loại dữ liệu mổ phỏng và dữ liệu thật cho kết quả khả quan. Chương trình có khả năng tìm kiếm và sắp hàng tạo nên những cặp gen có độ tương đồng rất cao, những cặp gen có khả năng cùng tiến hóa từ một đoạn gen trên hệ gen tổ tiên của chúng. Với những dữ liệu mô phỏng tạo ra ngẫu nhiên từ 13 đoạn gen con đặc trưng, kết quả thực nghiệm cho thấy mặc dù có những sự sai lệch về của từng nucleotide so với vị trí ban đầu tuy nhiên tỷ lệ sắp hàng chính xác của từng nucleotide với nucleotide tương ứng giữa hai hệ gen đạt mức trên 97%, cho thấy tỷ lệ tương đồng xác định được giữa hai hệ gen đạt được ở mức rất cao. Với các dữ liệu thật là các hệ gen được thu thập từ Gen Bank tại NCBI ( khi so sánh với kết quả thu được khi tính toán bằng phương thức Hungarian[8] chạy với các hệ gen được chia sẵn, thực nghiệm cho thấy độ tương đồng xác định được bằng chương trình thường tương đương thậm chí tốt hơn. Tài liệu tham khảo [1] Aaron C E, Darling, Bob Mau, Frederick R. Blatter, Nicole T Perna. Mauve: multiple alignment of conserved genomic sequence with rearrangements. 2004. [2] Altschul S F, Gish W, Miller W, Myers E W, Lipman D J. Basic local alignment search tool. 1990. [3] Altschul S.F Madden T L, Schaffer A, Zhang J, Zhang Z, Miller W, Lipman D J. Gapped BLAST and PSI-BLAST—a new generation of protein database search programs. 1997. [4] Bergman N.H. Comparative Genomics. 2007 [5] Brudno M, Do C, Cooper G, Kim M F, Davydov E, Green E D, Sidow A, Batzoglou S. LAGAN and Multi-LAGAN: Efficient tools for large- scale multiple alignment of genomic DNA. 2003. [6] Chiaromonte F, Yap V B, Miller W. Scoring pairwise genomic sequence alignments. Pacific Symp. Biocomput. 2002. [7] Darwin C. On the Origin of Species. John Murray, London, 6th edn. 1872. [8] Frank A. On Kuhn’s Hungarian Method - A tribute from Hungary. 2004. [9] Gotoh O. An improved algorithm for matching biological sequences. 1982. [10] Hamming R W. Error-detecting and error-correcting codes. 1950. [11] Hannenhalli S, Pevzner P. Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals. 1995. [12] Hirschberg D. A linear space algorithm for computing maximal common subsequences, 1975. [13] Lecture Notes in Systematics. The Ohio State University. Chapter 8. [14] Levenshtein V I. Binary codes capable of correcting deletions, insertions and reversals. 1996. [15] Ma B Tromp J, Li M. PatternHunter: Faster and more sensitive homology search.  2002. [16] Needleman S B  Wunsch C D. A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology. 1970.  [17] Sankoff D, Blanchette M. Multiple genome rearrangement and breakpoint phylogeny [Article]. - 1998. [18] Schwartz S, Kent WJ, Smit A, Zhang Z, Baertsch R, Hardison RC, Haussler D, Miller W. Human-Mouse Alignment with BLASTZ. 2000. [19] Shoffner JM, Wallace DC Oxidative phosphorylation diseases. The Metabolic and Molecular Bases of Inherited Disease. 1995. [20] Dinh Ngoc Thang, Le Sy Vinh, Hoang Xuan Huan. A Fast Algorithm for Genome Pairwise Alignment with Rearrangements, In Proceedings of KLLBI Workshop in the 10th Pacific Rim International Conference on Artificial Intelligence. Hanoi. 2008. [21] Ukkonen E On approximate string matching. Foundations of Computation Theory. 1993.  [22] Ukkonen E Algorithms for approximate string matching. Information and Control.  1985. [23] Le Sy Vinh, Varon A, Ward A W. Pairwise Alignment with Rearrangements. Genome Informatics. 2006. [24] Le Sy Vinh, Varon A, Janies D, Wheeler C W. Towards phylogenomic reconstruction. Proceeding of The 2007 International Conference on Bioinformatics and Computational Biology. 2007. [25] Waterman M S, Smith T F, Beyer W A Some biological sequence metrics. Advances in Mathematics. 1976.