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ự

Đại Học Quốc Gia Thành Phố Hồ Chí Minh Trường Đại Học Bách Khoa PHẠM MẠNH HÙNG CÁC KỸ THUẬT TOÁN HỌC CHO BÀI TOÁN SO SÁNH ĐA TRÌNH TỰ Chuyên ngành: Khoa học Máy tính LUẬN VĂN THẠC SĨ TP. HỒ CHÍ MINH, tháng 11 năm 2007 ĐẠI HỌC QUỐC GIA TP. HCM CỘNG HOÀ Xà HỘI CHỦ NGHIà VIỆT NAM TRƯỜNG ĐẠI HỌC BÁCH KHOA Độc Lập - Tự Do - Hạnh Phúc ---------------- ---oOo--- Tp. HCM, ngày . .05. . tháng . .11. . năm .2007. NHIỆM VỤ LUẬN VĂN THẠC SĨ Họ và tên học viên : Phạm Mạnh Hùng..............................Giới tính : Nam ;/ Nữ … Ngày, tháng, năm sinh : 26/2/1982....................................Nơi sinh : Phú Yên ................... Chuyên ngành : Khoa học Máy tính...................................................................................... Khoá : 2005 ......................................................................................................................... 1- TÊN ĐỀ TÀI : ................................................................................................................ CÁC KỸ THUẬT TOÁN HỌC CHO BÀI TOÁN SO SÁNH ĐA TRÌNH TỰ ........................................................................................................................................... ........................................................................................................................................... 2- NHIỆM VỤ LUẬN VĂN : .............................................................................................. ........................................................................................................................................... ........................................................................................................................................... ........................................................................................................................................... ........................................................................................................................................... 3- NGÀY GIAO NHIỆM VỤ : ........................................................................................... 4- NGÀY HOÀN THÀNH NHIỆM VỤ : .......................................................................... 5- HỌ VÀ TÊN CÁN BỘ HƯỚNG DẪN : TS. Nguyễn Văn Minh Mẫn .......................... Nội dung và đề cương Luận văn thạc sĩ đã được Hội Đồng Chuyên Ngành thông qua. CÁN BỘ HƯỚNG DẪN CHỦ NHIỆM BỘ MÔN (Họ tên và chữ ký) QUẢN LÝ CHUYÊN NGÀNH (Họ tên và chữ ký) TS. Nguyễn Văn Minh Mẫn TS. Đinh Đức Anh Vũ CÔNG TRÌNH ĐƯỢC HOÀN THÀNH TẠI TRƯỜNG ĐẠI HỌC BÁCH KHOA ĐẠI HỌC QUỐC GIA TP HỒ CHÍ MINH Cán bộ hướng dẫn khoa học : TS. Nguyễn Văn Minh Mẫn ........................................ Cán bộ chấm nhận xét 1 : ............................................................................................ Cán bộ chấm nhận xét 2 : ............................................................................................ Luận văn thạc sĩ được bảo vệ tại HỘI ĐỒNG CHẤM BẢO VỆ LUẬN VĂN THẠC SĨ TRƯỜNG ĐẠI HỌC BÁCH KHOA, ngày . . . . tháng . . . . năm . 2007 . 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 i LỜI CAM ĐOAN Tôi cam đoan rằng, ngoại trừ các kết quả tham khảo từ các công trình khác như đã ghi rõ trong luận văn, các công việc trình bày trong luận văn này là do chính tôi thực hiện và chưa có phần nội dung nào của luận văn này được nộp để lấy một bằng cấp ở trường này hoặc trường khác. Ngày 05 tháng 11 năm 2007 Phạm Mạnh Hùng 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 ii LỜI CẢM ƠN Tôi xin gởi lời cảm ơn chân thành nhất đến TS. Nguyễn Văn Minh Mẫn, người đã tận tình hướng dẫn, giúp đỡ tôi trong suốt quá trình thực hiện luận văn và tạo điều kiện để tôi có thể hoàn thành luận văn này. Xin cảm ơn gia đình và những người bạn đã dành cho tôi tình thương yêu và sự hỗ trợ tốt 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 iii TÓM TẮT LUẬN VĂN So sánh đa trình tự(Multiple Sequence Alignment-MSA) là một trong 10 bài toán lớn của Sinh tin học(Bioinformatics). MSA đóng vai trò quan trọng trong Sinh tin học nói chung và lĩnh vực tìm kiếm gene (Gene Finding) nói riêng. MSA là một bài toán NP, và hoàn toàn chưa có giải pháp trọn vẹn để tìm lời giải tối ưu của bài toán. Nhiều phương pháp sử dụng heuristic đã được đưa ra để giải quyết bài toán khi tập dữ liệu đầu vào lớn, các phương pháp này hướng tới việc tìm 1 lời giải cận tối ưu với thời gian tính toán và bộ nhớ sử dụng chấp nhận được. Progress Algorithm là một phương pháp tốt tiếp cận theo phương thức này. Đề tài này trình bày một giải thuật mới dựa trên Progressive Algorithm. Sử dụng lời giải của bài toán TSP để mô tả quá trình so sánh(align) các sequence. Để cung cấp một Progressive Algorithm có chất lượng, giải thuật đã tối ưu bài toán Pairwise Sequence Alignment(PSA) về độ chính xác và bộ nhớ sử dụng thông qua giải thuật ”chia để trị” kết hợp với việc sử dụng 3 ma trận đánh giá BLOSUM. Thông qua quá trình so sánh với CLUSTALW(một chương trình hiện thực Progressive Algorithm được đánh giá là cho kết quả tốt nhất), dựa trên kết quả kiểm thử với tập dữ liệu BAliBASE benchmark và một số nguồn dữ liệu từ NCBI(National Center for Biotechnology Information), chương trình hiện thực giải thuật đã cung cấp một lời giải có độ chính xác khá cao, tiết kiệm bộ nhớ và có thời gian tính toán chấp nhận được. Từ khoá: Algorithm, Sequence Alignment, Multiple Sequence Alignment, MSA, Pairwise Sequence Alignment, PSA, Progressive Algorithm, Dynamic Programming, Traveling Salesman Problem, TSP, CLUSTALW, BLOSUM, BAliBASE. 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 iv MỤC LỤC LỜI CAM ĐOAN ...........................................................................................................i LỜI CẢM ƠN ............................................................................................................... ii TÓM TẮT LUẬN VĂN .............................................................................................. iii DANH MỤC HÌNH ..................................................................................................... vi DANH MỤC BẢNG .................................................................................................. viii Chương 1. GIỚI THIỆU...............................................................................................1 1.1. Giới thiệu ..............................................................................................................1 1.2. Kết cấu của luận văn ...........................................................................................4 Chương 2. TỔNG QUAN VỀ KHÁI NIỆM SO SÁNH TRÌNH TỰ (SEQUENCE ALIGNMENT) ...........................................................................................6 2.1. So sánh trình tự....................................................................................................6 2.1.1. Định nghĩa So sánh trình tự(Sequence Alignment) ....................................................6 2.1.2. Phân loại .....................................................................................................................7 2.1.3. So sánh 2 trình tự (Pairwise Sequence Alignment-PSA)............................................8 2.1.4. So sánh nhiều trình tự (Multiple Sequence Alignment-MSA)....................................9 2.2. Các khái niệm khác ...........................................................................................10 2.2.1. Ma trận đánh giá(Scoring Matrix) ............................................................................12 2.2.2. Gap............................................................................................................................14 2.2.3. Phương pháp đánh giá(Scoring Method) ..................................................................15 2.3. Các phương pháp giải quyết bài toán so sánh trình tự ..................................18 2.3.1. Phương pháp Quy hoạch động(Dynamic Programming)..........................................19 2.3.2. Sử dụng các thiết bị phần cứng.................................................................................20 2.3.3. Phương pháp tìm kiếm cục bộ(Local Search)...........................................................21 2.3.4. Sử dụng giải thuật Di truyền(Genetic Algorithm) ....................................................21 2.3.5. Sử dụng Mô hình Markov ẩn(Hidden Markov Model-HMM). ................................21 Chương 3. CƠ SỞ LÝ THUYẾT VÀ PHƯƠNG PHÁP THỰC HIỆN .................24 3.1. Giới thiệu về Dynamic Programming ..............................................................24 3.2. Bài toán PSA và cách giải quyết bằng kỹ thuật quy hoạch động ..................24 3.2.1. Giải thuật quy hoạch động cho bài toán PSA ...........................................................25 3.2.2. Giải thuật Gotoh........................................................................................................28 3.2.3. Giải thuật cải tiến không gian nhớ ...........................................................................29 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 .........................................................................32 3.3.1. Giải thuật Center Star Alignment Algorithm............................................................33 3.3.2. Phương pháp Progressive Algorithm giải quyết bài toán MSA................................37 3.3.3. Feng-Doolittle Algorithm .........................................................................................38 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 .......................................................................42 4.1. Giải thuật sử dụng cho bài toán PSA...............................................................42 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 v 4.1.1. Giải thuật tính toán dựa theo kỹ thuật chia để trị......................................................43 4.2. Giải thuật hiện thực cho bài toán MSA ...........................................................49 4.2.1. Bài toán TSP(Travelling Salesman Problem-Bài toán người bán hàng). .................50 4.2.2. Giải thuật 1A.............................................................................................................51 4.2.3. Giải thuật 1B(Giải thuật cải tiến gom nhóm nhỏ nhất).............................................55 4.3. Giải thuật di truyền và bài toán TSP. ..............................................................57 4.3.1. Đặc điểm giải thuật di truyền....................................................................................57 4.3.2. Cấu trúc thuật giải di truyền tổng quát......................................................................59 4.4. Phần hiện thực giải thuật và chương trình: ....................................................60 Chương 5. KẾT QUẢ NHẬN XÉT............................................................................66 5.1. Một số kết quả chạy chương trình. ..................................................................66 5.2. BAliBASE (Benchmark Alignment Database)................................................68 5.3. So sánh kết quả ..................................................................................................69 5.3.1. Giới thiệu về các chương trình được sử dụng...........................................................70 5.3.2. So sánh độ chính xác của kết quả .............................................................................70 5.3.3. So sánh về mặt thời gian chạy, bộ nhớ .....................................................................77 Chương 6. KẾT LUẬN ...............................................................................................78 TÀI LIỆU THAM KHẢO...........................................................................................80 Phụ lục 1. Bảng đối chiếu Thuật ngữ Anh - Việt......................................................83 Phụ lục 2. Từ viết tắt ...................................................................................................87 Tham khảo Chỉ mục....................................................................................................88 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 vi DANH MỤC HÌNH Hình 2.1 Ví dụ về PSA............................................................................................................................7 Hình 2.2 Ví dụ về so sánh trình tự theo hướng toàn cục..........................................................................8 Hình 2.3 Ví dụ về so sánh trình tự theo hướng cục bộ.............................................................................8 Hình 2.4 Cấu trúc 1 PSA..........................................................................................................................8 Hình 2.5 Giới thiệu 1 MSA......................................................................................................................9 Hình 2.6 Giới thiệu các khái niệm của MSA .........................................................................................10 Hình 2.7 Quá trình biến đổi của 2 sequence...........................................................................................10 Hình 2.8 Ví dụ về các phép thay thế gap ..............................................................................................11 Hình 2.9 Ví dụ về Gap. ..........................................................................................................................15 Hình 2.10 Mối tương quan giữa các chương trình hiện thực cho các phương pháp. .............................19 Hình 2.11 Phương pháp tính toán chính xác bằng dynamic programming ............................................20 Hình 2.12 Mô hình Markov cho bài toán MSA. ....................................................................................22 Hình 3.1 Phương pháp quy hoạch động cho bài toán PSA ....................................................................25 Hình 3.2 Các ma trận S, D, I cho 2 chuỗi AGTAC and AAG. .............................................................31 Hình 3.3 Minh hoạ quá trình tìm 1 MSA tối ưu.....................................................................................33 Hình 3.4 Mô hình tiến hoá hình sao .......................................................................................................34 Hình 3.5 Minh họa Center Star Algorithm.............................................................................................35 Hình 3.6 Hình minh hoạ cho Progressive Algorithm. ............................................................................37 Hình 3.7 Minh họa Feng-Doolittle Algorithm .......................................................................................39 Hình 3.8 Ví dụ thực thi Feng-Doolittle Algorithm ................................................................................39 Hình 4.1 Mô hình quá trình thực hiện giải thuật PSA............................................................................43 Hình 4.2 Quá trình xây dựng ma trận của thuật giải cho bài toán PSA .................................................48 Hình 4.3 Quá trình align của Center Star Algorithm và phiên bản cải tiến ...........................................50 Hình 4.4 Bài toán TSP. ..........................................................................................................................50 Hình 4.5 Kết quả bài toán TSP...............................................................................................................51 Hình 4.6 Lưu đồ thuật giải 1A ...............................................................................................................52 Hình 4.7 Lưu đồ thuật giải 1B................................................................................................................55 Hình 4.8 Cấu trúc chương trình hiện thực..............................................................................................61 Hình 4.9 Module PSA ............................................................................................................................61 Hình 4.10 Sơ đồ các khối chức năng của Module MSA. .......................................................................62 Hình 4.11 Sơ đồ các khối chức năng của module TSP. .........................................................................63 Hình 5.1 Đồ thị tương quan về độ chính xác của MSAPR, CLUSTALW và MULTAL ......................72 Hình 5.2 Đồ thị tương quan về độ chính xác của MSAPR, CLUSTALW và HMMT...........................74 Hình 5.3 Đồ thị tương quan về độ chính xác của MSAPR, CLUSTALW và HMMT...........................75 Hình 5.4 Đồ thị tương quan về độ chính xác của MSAPR, CLUSTALW, SAGA................................75 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 vii Hình 5.5 Đồ thị tương quan về độ chính xác của MSAPR, CLUSTALW, SAGA................................76 Hình 5.6 So sánh thời gian thực thi của MSAPR và CLUSTALW .......................................................77 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 viii DANH MỤC BẢNG Bảng 2.1Ma trận BLOSUM62 lưu trữ hàm đánh giá độ tương đồng của tập 23 amino acid.................12 Bảng 2.2 Một phần ma trận Identity ......................................................................................................13 Bảng 3.1 Bảng kết quả giải thuật quy hoạch động cho bài toán PSA ....................................................26 Bảng 4.1 Định dạng của file dữ liệu đầu vào .........................................................................................63 Bảng 4.2 Định dạng của file dữ liệu đầu ra............................................................................................64 Bảng 4.3 Định dạng file dữ liệu đầu ra theo chuẩn MSF.......................................................................64 Bảng 4.4 Bảng tóm tắt các lớp của chương trình. ..................................................................................65 Bảng 5.1 TAT Protein HIV1..................................................................................................................66 Bảng 5.2 Kết quả Alignment của MSAPR và CLUSTALW với TAT HIV1 ........................................67 Bảng 5.3 Kết quả chạy chương trình với Nhóm 1 có chiều dài nhỏ ......................................................71 Bảng 5.4 Kết quả chạy chương trình với Nhóm 1 có chiều dài trung bình............................................71 Bảng 5.5 Kết quả chạy chương trình với Nhóm 1 có chiều dài lớn .......................................................72 Bảng 5.6 Kết quả chạy của các chương trình với các sequence của nhóm 2. ........................................73 Bảng 5.7 Kết quả chạy của các chương trình với các sequence của nhóm 3. ........................................74 Bảng 5.8 Kết quả chạy của các chương trình với các sequence của nhóm 4 .........................................75 Bảng 5.9 Kết quả chạy của các chương trình với các sequence của nhóm 5 .........................................76 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 1 Chương 1. GIỚI THIỆU 1.1. Giới thiệu Cùng với sự phát triển mang tính đột phá của Khoa học kỹ thuật, trong vài thập kỷ qua, sinh học phân tử đã có nhiều bước phát triển mạnh mẽ, một loạt các công cụ ứng dụng sinh học ra đời góp phần thúc đẩy quá trình giải mã một số lượng lớn trình tự bộ gen ở nhiều loài sinh vật. Cho đến nay, nhiều bộ gen vi khuẩn và các sinh vật bậc cao đã được giải mã gần như hoàn toàn. Dự án về bộ gen người được thành lập (1997), và quá trình giải trình tự tất cả 24 cặp nhiễm sắc thể của bộ gen người cũng đã hoàn thành từ cuối năm 2000, cũng như đã giải được khoảng 90% bộ gen người(2001). Lượng thông tin sinh học ngày càng trở nên phong phú và đa dạng. Ðể có thể xử lý và ứng dụng khối lượng thông tin đồ sộ như vậy , ngành Sinh tin học(hay Bioinformatics) ra đời, đó là sự kết hợp giữa công nghệ thông tin và sinh học, một cách đơn giản sinh tin học giải quyết các vấn đề của sinh học bằng cách sử dụng các kỹ thuật của khoa học máy tính. Các lĩnh vực lớn đang được Sinh tin học giải quyết hiện nay: ƒ Genomic: nghiên cứu cấu trúc và chức năng của gene. ƒ Proteinomics: Phân tích một tỉ lệ lớn các protein của một sinh vật ƒ Pharmacogenomics: phát triển các loại thuốc mới nhắm đến một căn bệnh xác định ƒ MicroArray: nghiên cứu về DNA chip, protein chip. Mục tiêu hàng đầu của sinh tin học gắn liền với quá trình phân tích các thông tin sinh học. Điều này được thể hiện qua các nghiên cứu về: ƒ Tìm kiếm các gene trên các trình tự DNA ở các sinh vật khác nhau. ƒ Phát triển các phương pháp nhằm dự đoán các trình tự RNA, cấu trúc và chức năng của các protein mới được phát hiện. ƒ Tập hợp các trình tự có sự tương đồng cao để đưa ra mô hình protein. ƒ So sánh các trình tự protein tương đồng và thành lập cây phả hệ mô tả mối quan hệ tiến hó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 2 Trong lĩnh vực nghiên cứu phân tích cấu trúc và chức năng của gene và protein, phân tích trình tự(chuỗi DNA, protein) đóng vai trò quan trọng. Để đơn giản cho việc nghiên cứu, trình tự DNA, protein sẽ được tuần tự hóa và nghiên cứu dưới dạng chuỗi các ký tự. Thông thường khi một gene được phát hiện, một trong những yêu cầu quan trọng là làm thế nào xác định được chức năng của gene này, yêu cầu tương tự cũng được đặt ra khi phát hiện ra protein mới. Một phương pháp tiếp cận phổ biến đó là chúng ta sẽ so sánh, đánh giá sự giống nhau(tương đồng) của chuỗi DNA, protein này với những chuỗi DNA, protein đã biết, từ đó có thể đưa ra dự đoán về chức năng cũng như cấu trúc của những gene mới phát hiện(Sequence Alignment). Quá trình tiến hóa của loài người là một quá trình biến đổi đa dạng, từ một gene(chuỗi DNA) tổ tiên dưới tác động của quá trình tiến hóa đã biến đổi tạo nên những khác biệt so với gene gốc ban đầu. Do đó việc nhận định sự giống nhau của các đoạn gene, trình tự là một vấn đề lớn của sinh tin học. Vấn đề được đặt ra (trong phân tích trình tự) đó là làm thế nào để có được phép so sánh tốt cho các trình tự DNA, khi mà số lượng tế bào trong cơ thể là khoảng 1014 và mỗi tế bào mang khoảng 3.109 ký tự trong đoạn DNA của chúng. Bài toán so sánh 2 trình tự(Pairwise Sequence Alignment-PSA) đã được giải quyết trọn vẹn bằng nhiều phương pháp khác nhau, đồng thời với việc giải quyết bài toán này, xuất hiện nhu cầu về việc so sánh nhiều trình tự, để so sánh nhiều đoạn gene hoặc tìm ra một phần tử đại diện cho một tập các gene nhằm đáp ứng nhu cầu ngày càng lớn của việc tìm kiếm dự đoán cấu trúc của các gene, protein, khi kho dữ liệu sinh học được tập hợp ngày càng lớn. Bài toán so sánh nhiều trình tự được đặt ra như vấn đề tất yếu. Không như bài toán so sánh 2 trình tự, bài toán so sánh nhiều trình tự(Multiple Sequence Alignment-MSA) là một bài toán NP mở, cho đến hiện tại (2007) vẫn chưa có một giải pháp nào có thể cung cấp một lời giải trọn vẹn cho bài toán, các lời giải thường tập trung vào việc tìm ra phép so sánh “gần” tốt nhất, và mỗi phương pháp tiếp cận sẽ chỉ cho những lời giải thực sự tốt tùy từng yêu cầu tiếp cận và bài toán cụ thể. Progressive Algorithm là một hướng giải quyết tốt cho bài toán so sánh nhiều trình tự. Đây là phương pháp kết hợp Qui hoạch động(Dynamic Programming) với heuristic. Phương pháp này sẽ tăng tốc độ tính toán, giảm độ phức tạp của giải thuật, có thể áp dụng cho các cơ sở dữ liệu gene lớn, phục vụ cho các dự án giải mã gene của các sinh vật bậc cao. 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 3 Từ khi được giới thiệu cho đến hiện nay, bài toán MSA đã và vẫn đang là một thách thức cho các nhà khoa học. Nghiên cứu và tìm ra một giải pháp cho bài toán vẫn là động lực thúc đẩy nhiều công trình khoa học về bài toán này. Xuất phát từ những đặc điểm của bài toán MSA đề tài này cố gắng tập trung vào giải quyết một số vấn đề sau: ƒ Khảo sát tổng quát các đặc điểm của bài toán MSA, các phương pháp giải quyết bài toán. ƒ Nghiên cứu về phương pháp dynamic programming, dynamic programming kết hợp với heuristic, Progressive Algorithm. ƒ Đề xuất một phương pháp giải quyết bài toán dựa trên Progressive Algorithm. ƒ Xây dựng chương trình hiện thực giải thuật được đề xuất và kiểm thử trên tập dữ liệu thực tế được lấy từ tổ chức NCBI(National Center for Biotechnology Information), và BAliBASE benchmark. Với những mục tiêu này đề tài đã thu được một số kết quả: ƒ Cung cấp cái nhìn tổng quan nhất về so sánh trình tự nói chung và bài toán MSA nói riêng. ƒ Phân loại các phương pháp giải quyết bài toán MSA, phân tích các ưu điểm và nhược điểm của từng phương pháp. ƒ Xây dựng giải thuật giải quyết bài toán MSA dựa trên việc cải thiện, tối ưu hoá bài toán PSA về độ chính xác cũng như bộ nhớ sử dụng, thông qua việc sử dụng 3 ma trận đánh giá BLOSUM, từ kết quả này của bài toán PSA sử dụng Progressive Algorithm kết hợp với lời giải của bài toán TSP để thực hiện quá trình so sánh nhiều trình tự, tìm ra lời giải cận tối ưu. ƒ Xây dựng thành công chương trình hiện thực giải thuật, cho phép tìm lời giải cho bài toán MSA với độ chính xác khá cao dựa trên kết quả kiểm thử trên các mẫu dữ liệu thực tế BAliBase benchmark và NCBI. Chương trình 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 4 cho phép tiết kiệm bộ nhớ sử dụng, cũng như thời gian tính toán chấp nhận được. 1.2. Kết cấu của luận văn Luận văn bao gồm 6 chương. Chuơng 1. GIỚI THIỆU Chương này trình bày về bối cảnh, mục tiêu cũng như kết quả thu được của luận văn. Chương 2. TỔNG QUAN VỀ KHÁI NIỆM SO SÁNH TRÌNH TỰ Chương này trình bày tổng quát về khái niệm so sánh trình tự, bài toán PSA, MSA, các phương pháp đánh giá chất lượng của MSA, các phương pháp giải quyết bài toán MSA. Chương 3. CƠ SỞ LÝ THUYẾT VÀ PHƯƠNG PHÁP THỰC HIỆN Chương này giới thiệu chung về phương pháp quy hoạch động(dynamic programming). Giới thiệu về phương pháp quy hoạch động giải quyết bài toán PSA, giải thuật tính giá trị PSA cải tiến về mặt bộ nhớ sử dụng. Phần tiếp theo của chương này trình bày về cách tiếp cận bài toán MSA hướng đến bài toán chính xác hoàn toàn bằng quy hoạch động thuần tuý, những khó khăn khi tiếp cận theo phương pháp này, giới thiệu một cách giải quyết bài toán MSA theo hướng gần đúng dựa trên kỹ thuật quy hoạch động kết hợp heuristic: Center Star Algorithm . Phần cuối chương này trình bày về 3 điểm chính, bao gồm giới thiệu Progressive Algorithm tổng quát, Progressive Algorithm phổ biến nhất, giải thuật Feng-Doolittle(Feng-Doolittle Algorithm) và một số chương trình hiện thực Progressive Algorithm trong thực tế. 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 Đây là chương dài nhất và cũng là chương giới thiệu những giải pháp mới của đề tài. Chương này trình bày về cách tiếp cận của luận văn để xây dựng giải thuật giải quyết bài toán MSA. Đầu chương giới thiệu về giải thuật tối ưu hoá tìm lời giải bài toán PSA dựa trên việc sử dụng kết hợp giải thuật tính giá trị PSA trình bày ở chương 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 5 3 và kỹ thuật chia để trị để tìm lời giải cho bài toán PSA. Phần này giới thiệu thêm việc sử dụng song song 3 ma trận BLOSUM làm hàm đánh giá để cải tiến độ chính xác, phù hợp với thực tế lời giải của bài toán PSA. Tiếp theo chương này đưa ra một giải pháp mới giải quyết bài toán MSA bằng cách kết hợp sử dụng giải thuật cho bài toán PSA vừa thu được, và giải thuật Feng-Doolittle để mô tả cách thức align các nhóm chuỗi trình tự(sequence) với nhau. Sử dụng kết quả bài toán TSP để tìm ra thứ tự align các nhóm sequence , lựa chọn điểm bắt đầu thực hiện quá trình align các sequence thông qua cách thức chọn điểm trung tâm của Center Star Algorithm. Sau nữa, chương này sẽ trình bày 1 cải tiến của giải thuật mới vừa nêu, nhằm nâng cao chất lượng của MSA bằng kỹ thuật gom nhóm theo khoảng cách ngắn nhất dựa trên thứ tự align thu được từ bài toán TSP. Gần cuối chương sẽ giới thiệu về phương pháp giải quyết bài toán TSP bằng giải thuật di truyền(Genetic Algorithm-GA), và cuối cùng sẽ giới thiệu về các module của chương trình hiện thực giải thuật vừa nêu. Chương 5. KẾT QUẢ, NHẬN XÉT Chương này giới thiệu về kết quả của chương trình hiện thực. Đánh giá kết quả này, so sánh với một số chương trình giải quyết bài toán MSA. Chương 6.KẾT LUẬN Chương này đề cập lại những việc đã thực hiện được của đề tài. Nêu lên hướng mở rộng và phát triển tiếp theo cho đề tài. 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 6 Chương 2. TỔNG QUAN VỀ KHÁI NIỆM SO SÁNH TRÌNH TỰ (SEQUENCE ALIGNMENT) 2.1. So sánh trình tự Như đã giới thiệu Gene Finding Problem là một hướng phát triển quan trọng của sinh tin học, Gene Finding dựa vào quá trình phân tích trình tự để có thể đưa ra được những kết quả. Một số các mục tiêu của quá trình phân tích trình tự: ƒ Xác định gene. ƒ Xác định chức năng của từng gene. Chúng ta giả thiết rằng các trình tự có cấu trúc tương tự nhau thì sẽ có chức năng tương tự nhau, do đó ta có thể tìm kiếm chức năng của một gene thông qua việc so sánh mức độ giống nhau của nó với một gene khác có chức năng đã được xác định. ƒ Xác định sự lặp lại của các trình tự ƒ Xác định protein dựa trên quy tắc sắp đặt của các biểu thức gene. ƒ Xác định các vùng chức năng khác nhau của DNA… Trong quá trình phân tích trình tự thì khái niệm so sánh trình tự đóng vai trò quan trọng. Đây là nền tảng cơ bản cho việc phân tích. So sánh trình tự sẽ giúp cho quá trình dự báo sự giống nhau về chức năng của các trình tự, dự báo cấu trúc bậc 3 của DNA, protein. Trong việc tìm hiểu một gene mới, chúng ta thường quan tâm đến việc xác định những đặc điểm để phân biệt gene đồng thời đưa ra những giả thuyết về chức năng của gene. Việc đưa ra những giả thuyết về chức năng của gene thường dựa vào những giải thuật đánh giá sự giống nhau, tương đồng giữa các trình tự. 2.1.1. Định nghĩa So sánh trình tự(Sequence Alignment) So sánh trình tự(một số tài liệu gọi là phép gióng hàng, gióng cột) là quá trình nghiên cứu sự giống nhau giữa các chuỗi trình tự(sequence), đo lường sự giống nhau giữa các trình tự. Là cách thức so sánh giữa 2 hay nhiều trình tự dựa trên việc so sánh một chuỗi các thành phần(ký tự) của trình tự để tìm ra những điểm tương đồng, giống nhau giữa các trình tự. Các trình tự được đề cập trong phần nghiên cứu này là các các chuỗi trình tự DNA, RNA hoặc các trình tự amino acid( protein). 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 7 Xét 2 chuỗi: A C G C T G và C A T G T. Chúng ta sẽ đo lường sự giống nhau giữa 2 chuỗi này. Bên dưới là một ví dụ về một trong những khả năng alignment của 2 chuỗi trên. Tiêu chí để đánh giá sự giống nhau sẽ dựa trên một hàm đánh giá (scoring function). Hình 2.1 Ví dụ về PSA Ký tự “–“ được gọi là 1 gap, gap thể hiện cho ý nghĩa trong sinh học đó là một phần của trình tự đã bị mất đi do, các hành vi của quá trình tiến hóa sinh vật như: đột biến, sự mất đi một thành phần trong chuỗi trình tự. Giá trị hàm đánh giá càng cao thì chúng ta có một kết quả alignment càng tốt. Xét một hàm đánh giá đơn giản như sau, nếu 2 thành phần trong chuỗi là giống nhau thì hàm đánh giá sẽ có kết quả +2, nếu 2 thành phần trong chuỗi khác nhau thì hàm đánh giá tại vị trí này sẽ có kết quả -1, như vậy kết quả so sánh ở trên sẽ có giá trị hàm đánh giá là: 3*(2)+(-1)*5=1 Ý nghĩa: Trên quan điểm sinh học, phép so sánh trình tự thể hiện quá trình biến đổi chọn lọc tự nhiên của các chuỗi trình tự, từ đó cho phép các nhà sinh học đưa ra kết luận về nguồn gốc của các đoạn gene, DNA, RNA, hay protein. 2.1.2. Phân loại Dựa trên phương pháp so sánh người ta chia ra làm 2 loại alignment: ƒ Phép alignment trình tự theo hướng toàn cục(Global Sequence Alignment): Phép toán alignment được áp dụng trên toàn bộ chuỗi trình tự. Thường được sử dụng khi các trình tự so sánh có kích thước gần tương đương và các trình tự này có độ tương đồng, giống nhau cao. ƒ Phép alignment trình tự theo hướng cục bộ(Local Sequence Alignment): Phép toán alignment được sử dụng trên một phần của chuỗi trình tự. Thường được sử dụng khi các trình tự có chiều dài lớn, độ tương đồng giống nhau không cao, chỉ có một số ít các gene giống nhau trên 2 trình tự, hoặc khi 2 trình tự có kích thước khác biệt lớn. A C – – G C T G – C A T G – 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 8 Ví dụ về so sánh trình tự theo hướng toàn cục: Hình 2.2 Ví dụ về so sánh trình tự theo hướng toàn cục Toàn bộ 2 chuỗi trình tự L G P S S K Q T G K G S − S R I W D N và L N − I T K S A G K G A I M R L G D A được so sánh Ví dụ về so sánh trình tự theo hướng cục bộ: Hình 2.3 Ví dụ về so sánh trình tự theo hướng cục bộ Chỉ một phần của 2 chuỗi được so sánh: TGKG và AGKG Tùy thuộc vào số lượng trình tự, bài toán so sánh trình tự được chia làm 2 mức độ: ƒ So sánh 2 trình tự ƒ So sánh nhiều trình tự. 2.1.3. So sánh 2 trình tự (Pairwise Sequence Alignment-PSA). Định nghĩa 2.1: Gọi S1 và S2 là 2 chuỗi, một sự alignment A giữa S1 và S2 sẽ tạo ra 2 chuỗi S’1 và S’2 bằng cách thêm vào các ký tự “-“ vào S1, S2, trong đó: ƒ |S’1|=|S’2| ƒ Nếu loại bỏ các ký tự “-“ khỏi S’1 và S’2 ta sẽ có S1 và S2 Với |S1|, |S2| lần lượt là chiều dài của S1 và S2. Hình 2.4 Cấu trúc 1 PSA L G P S S K Q T G K G S − S R I W D N L N − I T K S A G K G A I M R L G D A − − − − − − − T G K G − − − − − − − − − − − − − − − A G K G − − − − − − − − A C G C T G C A T G T A C – – G C T G – C A T G – T – (S1) (S2) (S’1) (S’2) 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 9 2.1.4. So sánh nhiều trình tự (Multiple Sequence Alignment-MSA) Trong mục này ta sẽ xem xét nguồn gốc sinh học của quá trình thực hiện Multiple Squence Alignment. Nguyên nhân chính cho sự ra đời của quá trình so sánh nhiều trình tự là việc so sánh sự tương đồng về trình tự của các protein với một tập protein đã có sẵn trong Cơ Sở Dữ Liệu (CSDL). Thông thường các protein lưu trong CSDL thường được tổ chức thành các nhóm chung(protein family), có sự tương đồng về cấu trúc, chức năng, và cấu trúc bậc 3. Khi cần khảo sát một protein mới, chúng ta hy vọng có thể đưa ra các giả thuyết về cấu trúc, chức năng và quá trình tiến hóa của protein này thông qua phép toán alignment. Tuy nhiên, chúng ta không thể thực hiện việc alignment chuỗi trình tự của protein này với từng trình tự của mỗi protein trong CSDL, điều này là không thể về mặt thời gian xử lý. Do đó cách tiếp cận tốt hơn là chúng ta sẽ so sánh trình tự của protein này với mỗi tập hợp protein trong CSDL, thông qua việc so sánh trình tự của protein này với một trình tự đại diện cho mỗi tập hợp protein. Vấn đề đặt ra là làm cách nào để tìm ra trình tự đại diện cho một tập hợp protein, câu trả lời sẽ được cung cấp thông qua việc thực hiện phép toán multiple alignment của tập hợp protein này, để tìm ra phần tử tương đồng nhất đại diện cho tập hợp protein. Định nghĩa 2.2: Cho k chuỗi S1, S2,…,Sk một phép toán Multiple Sequence Alignment(MSA) của k chuỗi này sẽ tạo ra k chuỗi mới S’1, S’2,…, S’k bằng cách thêm các ký tự “-” vào các chuỗi S1, S2,…,Sk trong đó: ƒ |S’1|=|S’2|=…=|S’k|=n ƒ Nếu bỏ đi các ký tự “-” khỏi S’i ta sẽ được lại chuỗi ban đầu Si (1≤i≤k) Hình 2.5 Giới thiệu 1 MSA k được gọi là số sequence của MSA n: là chiều dài của MSA. MSA có n cột, mỗi cột chứa các ký tự đại diện cho các sequence của MSA, các ký tự này có thể là các amino acid(nucleotide) hoặc phần tử ‘-‘. A G T − G T G A G T A G T G − G T C G T G − − T A G T G A G T G T G A G T A G T G G T C G T G T A G T G (S1) (S2) (S3) (S4) (S’1) (S’2) (S’3) (S’4) 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 10 Hình 2.6 Giới thiệu các khái niệm của MSA 2.2. Các khái niệm khác Trong quá trình tiến hóa của mình trên một đoạn gene, DNA, protein có thể xuất hiện các hành vi: đột biến(mutation), mất(delete), thêm hoặc giữ lại trạng thái di truyền của các nucleotide(DNA) hoặc amino acid với protein. Hình 2.7 Quá trình biến đổi của 2 sequence Như vậy kết quả của phép so sánh giữa ACTCGATT và AGCTAATC: Chúng ta xét 2 chuỗi A=a1a2…am , B=b1b2…bn. Việc biến đổi từ chuỗi trình tự A sang B là sự kết hợp của các quá trình: quá trình thay thế, sự xuất hiện các gap. Một sự thay thế trong quá trình biến đổi từ A sang B là sự thay thế của 1 phần tử của A bằng 1 phần tử của B. Sự thay thế có thể là 1 trong 2 quá 1 trình: quá trình đột biến hoặc quá trình giữ lại trạng thái di truyền. A G T − G T G A G T A G T G − G T C G T G − − T A G T G Cột 7 của MSA k=4 sequence Chiều dài của MSA n=7 ACTCGATT ACTCGATT ACGATT Mất TC tại vị trí thứ 3, 4 AGCATC ACTCGATT Đột biến thay C bằng G tại vị trí 2 Đột biến thay G bằng C tại vị trí 3 Đột biến thay T bằng C tại vị trí 6 AGCTAATC Thêm TA vào vị trí 4 A C T C G - - A T T A G - - C T A A T C 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 11 Một gap bao gồm các phần tử của 1 chuỗi mà mỗi phần tử này tương ứng với các phần tử có ký hiệu là “-“ của chuỗi còn lại. Gap gồm có 2 loại: deletion gap và insertion gap, tương ứng với quá trình thêm vào hoặc mất đi các phần tử di truyền Hình 2.8 Ví dụ về các phép thay thế gap Như vậy thông qua kết quả alignment của 2 chuỗi có thể biết được quá trình biến đổi, tiến hóa từ chuỗi này sang chuỗi khác. Chúng ta có thể quan sát một ví dụ thực tế sau: GCGCTCCGGGACGCCTTCCGCCGTCGGGAGCCCTACAACTACCTGCAGAGGGCCTATTAC -------------------------||||||| ||||||||||||||||||||||| ||| GGGAGCCTTACAACTACCTGCAGAGGGCCTACTAC CAGGTGGGGAGCGGGCCGGGCAG TAG |||||| ||---||||||| |||------------------------------------- CAGGTGCGG GGGCCGGCCAGGGTGCTACCCCAAGCCTACTGACTGTCTTACTGG CCTTCCCCAGAGCCCCCTAGCCGCAGGCACCAGAGGGTCCAAGACAAGACTGGAAGGGCA -----------------------|| || ||| | ||||| || || |||| | | | CAAGCTTCAGCGAGTCCAGGAGAAAGCTGGGAAGCCC CCTCGGGTTCGG GAGGAGCTGTGAGTGGCT | ||||| |||------||||| |||||| |||||------------------------ CGCCGGGTCCGGGTCCGAGAGGAACTGTGAATGGCTGAGCCTGCTTCTCGAGGATCAGGC Mỗi một alignment có một giá trị thông qua việc đánh giá các phần tử tạo thành nó. Giá trị này phản ánh quá trình biến đổi chuỗi A thành chuỗi B và ngược lại. Định nghĩa giá trị này một cách tối ưu sẽ cho phép ta tìm được một quá trình biến đổi tối ưu từ A thành B(và ngược lại), điều này cũng đồng nghĩa với việc tìm phép alignment tối ưu của A và B. Việc tính toán giá trị cho một alignment phụ thuộc vào: ƒ Ma trận đánh giá ƒ Giá trị của Gap ƒ Phương pháp đánh giá A C T C G - - A T T A G - - C T A AT C Sự thay thế Deletion gap Insertion gap 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 12 2.2.1. Ma trận đánh giá(Scoring Matrix) Kết quả của việc tính giá trị cho mỗi alignment phụ thuộc nhiều vào kết quả của hàm đánh giá sự tương đồng của mỗi cặp amino acid(nucleotide) ( , )a bσ . Trong phần ví dụ về alignment ta có đưa ra một ví dụ đơn giản cho ( , )a bσ tuy nhiên trong thực tế sinh học khả năng xuất hiện của mỗi cặp amino acid(cũng như nucleotide) là khác nhau, xác suất xuất hiện cùng lúc của cặp amino acid này có thể cao trong khi xác suất xuất hiện của cặp amino acid kia có thể thấp. Vì thế độ tương đồng của các cặp amino acid thường được lưu trữ dưới dạng một ma trận 2 chiều gọi là ma trận đánh giá. Một ví dụ: Bảng 2.1Ma trận BLOSUM62 lưu trữ hàm đánh giá độ tương đồng của tập 23 amino acid Xét trên phương diện toán, ma trận đánh giá là 1 ánh xạ được định nghĩa như sau: 2: ( ')σ ∑ →\ trong đó ' {' '}∑ = ∑∪ − và ∑ là tập các amino acid hoặc nucleotide Có nhiều hình thức ma trận đánh giá khác nhau dựa trên quá trình nghiên cứu, thống kê thực tế sinh học. Hiện tại có 4 hình thức ma trận đánh giá: identity matrix, genetic code matrix, chemical similarity matrix và substitution matrix. A R N D C Q E G H I L K M F P S T W Y V B Z X * A 6 -2 -2 -3 -1 -1 -1 0 -2 -2 -2 -1 -1 -3 -1 2 0 -4 -3 0 -2 -1 -1 -6 R -2 8 -1 -2 -5 1 0 -3 0 -4 -3 3 -2 -4 -3 -1 -2 -4 -3 -4 -2 0 -2 -6 N -2 -1 8 2 -4 0 0 -1 1 -5 -5 0 -3 -4 -3 1 0 -6 -3 -4 5 0 -2 -6 D -3 -2 2 9 -5 0 2 -2 -2 -5 -5 -1 -5 -5 -2 0 -2 -6 -5 -5 6 1 -2 -6 C -1 -5 -4 -5 13 -4 -5 -4 -4 -2 -2 -5 -2 -4 -4 -1 -1 -3 -4 -1 -5 -5 -3 -6 Q -1 1 0 0 -4 8 3 -3 1 -4 -3 2 -1 -5 -2 0 -1 -3 -2 -3 0 5 -1 -6 E -1 0 0 2 -5 3 7 -3 0 -5 -4 1 -3 -5 -2 0 -1 -4 -3 -4 1 6 -1 -6 G 0 -3 -1 -2 -4 -3 -3 8 -3 -6 -5 -2 -4 -5 -3 0 -2 -4 -5 -5 -1 -3 -2 -6 H -2 0 1 -2 -4 1 0 -3 11 -5 -4 -1 -2 -2 -3 -1 -3 -4 3 -5 -1 0 -2 -6 I -2 -4 -5 -5 -2 -4 -5 -6 -5 6 2 -4 2 0 -4 -4 -1 -4 -2 4 -5 -5 -2 -6 L -2 -3 -5 -5 -2 -3 -4 -5 -4 2 6 -4 3 1 -4 -4 -2 -2 -2 1 -5 -4 -2 -6 K -1 3 0 -1 -5 2 1 -2 -1 -4 -4 7 -2 -5 -2 0 -1 -4 -3 -3 -1 1 -1 -6 M -1 -2 -3 -5 -2 -1 -3 -4 -2 2 3 -2 8 0 -4 -2 -1 -2 -1 1 -4 -2 -1 -6 F -3 -4 -4 -5 -4 -5 -5 -5 -2 0 1 -5 0 9 -5 -4 -3 1 4 -1 -5 -5 -2 -6 P -1 -3 -3 -2 -4 -2 -2 -3 -3 -4 -4 -2 -4 -5 11 -1 -2 -5 -4 -4 -3 -2 -2 -6 S 2 -1 1 0 -1 0 0 0 -1 -4 -4 0 -2 -4 -1 6 2 -4 -3 -2 0 0 -1 -6 T 0 -2 0 -2 -1 -1 -1 -2 -3 -1 -2 -1 -1 -3 -2 2 7 -4 -2 0 -1 -1 -1 -6 W -4 -4 -6 -6 -3 -3 -4 -4 -4 -4 -2 -4 -2 1 -5 -4 -4 16 3 -4 -6 -4 -3 -6 Y -3 -3 -3 -5 -4 -2 -3 -5 3 -2 -2 -3 -1 4 -4 -3 -2 3 10 -2 -4 -3 -2 -6 V 0 -4 -4 -5 -1 -3 -4 -5 -5 4 1 -3 1 -1 -4 -2 0 -4 -2 6 -5 -4 -1 -6 B -2 -2 5 6 -5 0 1 -1 -1 -5 -5 -1 -4 -5 -3 0 -1 -6 -4 -5 5 0 -2 -6 Z -1 0 0 1 -5 5 6 -3 0 -5 -4 1 -2 -5 -2 0 -1 -4 -3 -4 0 5 -1 -6 X -1 -2 -2 -2 -3 -1 -1 -2 -2 -2 -2 -1 -1 -2 -2 -1 -1 -3 -2 -1 -2 -1 -2 -6 * -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 -6 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 13 Identity matrix: đây là cơ chế đánh giá độ tương đồng đơn giản nhất, trong ma trận này các cặp amino acid giống nhau sẽ có giá trị của phần tử tương ứng trong ma trận là 1, các cặp amino acid còn lại sẽ nhận giá trị 0. Ví dụ Bảng 2.2 Một phần ma trận Identity Ma trận mã di truyền(Genetic code matrix): trong ma trận này hàm đánh giá của mỗi cặp amino acid dựa trên độ tương đồng về mã di truyền. Ngày nay ma trận này hiếm khi được sử dụng trong việc alignment các chuỗi amino acid. Ma trận tương đồng hóa học(chemical similarity matrix): trong ma trận này, các amino acid có cấu trúc tương đồng về cấu trúc vật lý cũng như thuộc tính hóa học như kích thước, hình dạng, khả năng phân cực,… thì phần tử tương ứng trong ma trận sẽ nhận giá trị lớn hơn so với các cặp còn lại. Ma trận thay thế(substitution matrix): Ma trận này được tính toán và xây dựng dựa trên các quan sát thống kê về tần số thay đổi của các amino acid trong việc alignment các chuỗi trình tự. Ma trận thay thế được đánh giá là tốt hơn so với 3 hình thức ma trận trên và hiện nay cũng được sử dụng phổ biến nhất. Trong phần này, xin giới thiệu 1 hình thức ma trận thay thế hay được sử dụng trong các công trình nghiên cứu cũng như trong các công cụ phục vụ cho việc tính toán sinh học: ma trận BLOSUM. Ma trận BLOSUM(Block Subtitutation Matrix): Khái niệm ma trận BLOSUM được xây dựng dựa trên 1 MSA của tập các protein có sự khác biệt về quá trình tiến hóa. Hàm đánh giá của các cặp amino acid trong ma trận được tính dựa vào tần số quan sát sự thay đổi trong các khối(block) của các local alignment của các họ protein trong cơ sở dữ liệu BLOCK(Henikoff và Henikoff 1992)[14]. Cơ sở dữ liệu này gồm 3000 block của các protein đặc trưng, đại diện cho vài trăm nhóm protein. Các ma trận BLOSUM được thiết kế nhằm giải quyết cho các A R N D C Q A 1 0 0 0 0 0 R 0 1 0 0 0 0 N 0 0 1 0 0 0 D 0 0 0 1 0 0 C 0 0 0 0 1 0 Q 0 0 0 0 0 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 14 trường hợp align các sequence có sự biến đổi, tiến hóa về mặt di truyền trong một khoảng thời gian dài. Ma trận BLOSUM bao gồm nhiều cấp độ, ký hiệu BLOSUMn. Ma trận BLOSUMn(1 100n≤ ≤ ) cho biết độ tương đồng của các chuỗi được dùng để tính ra chúng. Ví dụ, chúng ta xét ma trận BLOSUM62, giá trị của các phần tử trong ma trận được tính từ tập các protein có độ tương đồng không lớn hơn 62%. Trong tập các ma trận BLOSUMn, các ma trận có chỉ số n nhỏ thường được sử dụng trong việc align các sequence có độ khác biệt cao(độ tương đồng thấp), và các ma trận có chỉ số n lớn thường được sử dụng cho các sequence có độ tương đồng cao. Ví dụ: ma trận BLOSUM62 thường được sử dụng nhất cho việc align các sequence khi chưa xác định độ tương đồng của chúng, ma trận BLOSUM45 thường được sử dụng cho các sequence có sự khác biệt cao, ma trận BLOSUM100 thường được sử dụng cho các ma trận đó độ tương đồng cao. Việc tính toán tập các ma trận BLOSUM dựa trên công thức xác suất biến đổi: ( , ) ( , ) log( )O E Pscore a b a b k P σ= = trong đó OP là xác suất chuyển đổi từ amino acid a sang amino acid b trong tập quan sát EP : là xác suất xuất hiện của amino acid b trong tập quan sát. k là hệ số làm tròn. Thông thường k=10, và giá trị hàm đánh giá được làm tròn thành số nguyên. ( , ) 0score a b > : cho biết sự thay thế giữa amino acid a và b có khả năng xảy ra cao hơn so với sự thay đổi một cách ngẫu nhiên ( , ) 0score a b < : cho biết sự thay thế giữa amino acid a và b có khả năng xảy ra thấp hơn so với sự thay đổi một cách ngẫu nhiên. ( , ) 0score a b = : cho biết sự thay thế giữa amino acid a và b tương đương với việc thay thế 2 amino acid một cách ngẫu nhiên. 2.2.2. Gap Việc tính giá trị của mỗi gap phụ thuộc vào bản chất tự nhiên của các chuỗi trình tự, có thể tính giá trị của gap theo các hàm tuyến tính hoặc đa thức. Trong thực 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 15 đơn giản hầu hết các phương pháp đề xuất tính giá trị của gap dựa trên một hàm tuyến tính theo chiều dài của gap. Chúng ta sẽ xem xét “Phương pháp tính giá trị của gap dựa trên hàm tuyến tính”. Thông thường khả năng xuất hiện phần tử đầu tiên của gap(khả năng mở gap) thường khác với khả năng xuất hiện của các phần tử theo sau, gọi q(q>0) là giá trị xác định khả năng mở một gap, và khả năng xuất hiện của mỗi phần tử trong gap(khả năng mở rộng của gap) là r(r>0). Gọi )(kγ là hàm giá trị của gap ta định nghĩa: ( ) ( * )γ = − +k q r k Ví dụ : CAGGTGGGGAGCGGGCCGGGCAG |||||| ||---||||||| ||| CAGGTGCGG GGGCCGGCCAG Hình 2.9 Ví dụ về Gap. Trong ví dụ trên: Khả năng xuất hiện của cặp đầu tiên A − là q+r . Khả năng xuất hiện của các phần tử tiếp theo trong gap lần lượt là r. Giá trị của gap: (3) ( *3)γ = − +q r . 2.2.3. Phương pháp đánh giá(Scoring Method) Phương pháp đánh giá cho phép đánh giá sự giống nhau, tương đồng giữa các trình tự dựa trên một số tiêu chí nào đó. Việc so sánh giữa 2 hay nhiều chuỗi trình tự sẽ cho ra nhiều kết quả so sánh khác nhau từ một tập chuỗi trình tự ban đầu. Cơ sở để đánh giá sự giống nhau, tương đồng giữa các trình tự sau phép alignment thường căn cứ vào một hàm đánh giá cụ thể. Việc xây dựng hàm đánh giá tốt sẽ cho phép xác định được kết quả nào của phép so sánh là tối ưu. Hàm đánh giá chính là cốt lõi của một phương pháp đánh giá. Đối với PSA, phương pháp đánh giá phổ biến nhất là dựa vào tổng giá trị của các cặp ký tự đại diện không phải là gap, và giá trị của các gap trong PSA. Gọi K1 là tập các phần tử không phải là gap của PSA, K2 là tập các gap trong PSA. Hàm đánh giá PSA của 2 chuỗi S1, S2 với kết quả S’1, S’2 có dạng: 1 2( , )Score S S = Score(K1)+Score( K2) 1 1 1 2( ) ( ' [ ], ' [ ])σ ∈ = ∑ i K Score K S i S i , với S’1[i], S’2[i] là các phần tử thuộc K1. 2 2( ) ( )i i K Score K kγ ∈ = ∑ với γ(ki) là các hàm tính cho các gap thuộc K2. 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 16 Đối với một MSA do bản chất phức tạp của dữ liệu sinh học nên tất cả các phương thức đánh giá đều có những hạn chế, và không có một tiêu chuẩn tổng quát nào trong việc đo lường chất lượng của nó. Trong phần này xin được giới thiệu một số phương pháp xây dựng hàm đánh giá phổ biến cho một MSA [1]. ƒ Sum-of-Pair(SP): Đây là phương pháp được sử dụng phổ biến nhất. Nội dung của phương pháp này là đánh giá MSA của k sequence dựa trên tổng kết quả alignment của tất cả ( k2 ) cặp sequence có trong MSA. Theo phương pháp này giá trị của mỗi cột của MSA sẽ được tính bằng tổng tất cả các hàm đánh giá độ tương đồng của các cặp phần tử trong cột này. Gọi ci,j là ký tự tại dòng i, cột j trong MSA và ( , )a bσ là hàm đánh giá sự tương đồng của cặp amino acid a, b. Ta sẽ có giá trị hàm đánh giá độ tương đồng SP của cột i như sau: , ,( ) ( , ) 1 , ,1x i y i x y SP i c c x y k i nσ < = ≤ ≤ ≤ ≤∑ Trong đó n là chiều dài của MSA. Gọi 1( ,..., )kSPScore S S là hàm đánh giá độ tương đồng của MSA theo SP. 1( ,..., )kSPScore S S sẽ bằng tổng tất cả các hàm đánh giá độ tương đồng của các cột trong MSA: 1 1 1 ( ,..., ) ( ,..., ) ( ) k k k i Score S S SPScore S S SP i = = =∑ MSA tối ưu là MSA có SP tối ưu. Công thức này cũng có thể được viết lại dưới dạng sau: 1 1 1 1 ( ,..., ) ( ,..., ) ( , ) k k k k i j i j i Score S S SPScore S S SPScore S S = = + = =∑ ∑ trong đó ( , )i jSPScore S S là hàm đánh giá theo SP của 2 chuỗi Si và Sj Ví dụ: Xét một MSA của 3 sequence: ACCTG , CCTGT, AGCTAT (S1) A C - C T G - (S2) - C - C T G T (S3) A - G C T A T Nếu chúng ta định nghĩa ( , )a bσ như sau: 0 ( , ) 1 a b a b a b σ ≠⎧= ⎨ =⎩ Khi đó SP(2)= ( , ) ( , ) ( , )C C C Cσ σ σ+ − + − =1+0+0=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 17 SP=SP(1)+SP(2)+SP(3)+SP(4)+SP(5)+SP(6)+SP(7)=1+1+1+3+3+1+1=11 Kết quả đánh giá độ tương đồng của 3 chuỗi này cũng có thể được biểu diễn: SP(S1,S2,S3)=SP(S1,S2)+SP(S2,S3)+SP(S3,S1)=5+3+3=11 ƒ Phương pháp #LOG#: Phương pháp này đánh giá độ tương đồng của MSA dựa trên các công thức sau: Hàm đánh giá độ tương đồng của MSA tại cột thứ i: # ( ) # ( ) log ( )i i a LOG i c a c a ∈∑ = ∑ Công thức đánh giá độ tương đồng của toàn bộ MSA: 1 1 1 ( ,..., ) # ( ,..., ) # ( ) log ( ) n k k i i i a Score S S LOG S S c a c a = ∈∑ = =∑∑ Trong đó a∈∑ với ∑ là tập các amino acid hoặc nucleotide ( )ic a là số lần xuất hiện của a tại cột thứ i. ƒ Phương pháp Trung bình thông tin(Average Information Content): Trong phương pháp này hàm đánh giá độ tương đồng của MSA tại cột thứ i được xây dựng dựa trên công thức sau: ( )1( ) ( ) log ( ) i i a f aIC i f a n p a∈∑ = ∑ fi(a) là tần số xuất hiện của a trong cột thứ i của MSA. ( ( )( ) ii c af a k = ) p(a) là tần số xuất hiện của a trong toàn bộ MSA. Khi đó công thức tính độ tương đồng của MSA: 1 1 1 ( )1( ,..., ) ( ,..., ) ( ) log ( )= ∈Σ = = ∑∑n ik k i i a f aScore S S IC S S f a n p a ƒ Phương pháp Entropy: Phương pháp này đánh giá độ tương đồng của MSA dựa trên lý thuyết về Entropy Hàm đánh giá độ tương đồng của cột thứ i trong MSA: ( ) ( ) log ( )i i a entropy i c a f a ∈∑ = −∑ Hàm đánh giá độ tương đồng của toàn bộ MSA: 1 1 1 ( ,..., ) ( ) ( ) log ( ) n n k i i i i a score S S entropy i c a f 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 18 Trong phạm vi của luận văn này phương pháp được sử dụng để đánh giá là phương pháp Sum-of-Pair. 2.3. Các phương pháp giải quyết bài toán so sánh trình tự Hiện nay đã có nhiều cách tiếp cận khác nhau của các nhà nghiên cứu về vấn đề này, Pairwise Sequence Alignment-PSA đã được giải quyết khá triệt để, tuy nhiên bài toán Multiple Sequence Alignment-MSA dựa trên phương pháp đánh giá Sum-of-Pair vẫn còn là một vấn đề để ngỏ. Wang và Jiang [32] đã chứng minh rằng việc tìm kiếm phép Alignment tối ưu cho bài toán Multiple Alignment là bài toán NP. Đã có nhiều mô hình, thuật giải được áp dụng để giải quyết bài toán này: Kỹ thuật quy hoạch động(Dynamic Programming), sử dụng mô hình Markov ẩn (Hidden Markov Model), các kỹ thuật tìm kiếm cục bộ(Local Search)…. Hầu hết các phương pháp giải quyết rất tốt bài toán PSA. Tuy nhiên vì MSA là bài toán NP nên các kỹ thuật đưa ra đều không thể đánh giá chính xác về mức độ tốt xấu của từng giải thuật, mỗi giải thuật sẽ có ưu điểm trong từng trường hợp cụ thể [6], [28], [35]. Việc giải quyết bài toán được chia thành 2 hướng: ƒ Phương pháp tính toán chính xác, tìm ra MSA tối ưu ƒ Phương pháp tính toán gần đúng tìm ra MSA cận tối ưu. Phương pháp tính toán chính xác: bao gồm phương pháp quy hoạch động tổng quát(Needleman-Wuns)[22], và một số giải pháp phần cứng. Phương pháp tính toán cận tối ưu: Phương pháp này sử dụng các tiếp cận heuristic để giải quyết bài toán trên tập dữ liệu lớn theo hướng tìm ra các lời giải cận tối ưu, có thể chấp nhận được về độ chính xác cũng như thời gian và không gian bộ nhớ sử dụng. Phương pháp này chia thành 2 nhóm: ƒ Progressive Algorithm ƒ Iterative Algorithm Progressive Algorithm: về cơ bản phương pháp này vẫn dựa trên nền tảng của dynamic programming. Ý tưởng của Progressive Algorithm là tìm MSA gần tối ưu thông qua việc align các sequence với nhau, 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. Ưu điểm của phương pháp này giúp tìm lời giải nhanh hơn. Chi tiết về phương pháp này sẽ đề cập trong chương 3. Một số chương trình hiện thực theo 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 19 phương pháp này: CLUSTALW, MULTALIGN, PILEUP, BLAST, FASTA. MULTAL, DIALIGN. Iterative Algorithm:Ý tưởng của các giải thuật theo phương pháp này là : đầu tiên sẽ xây dựng một MSA mà không quan tâm đến độ tốt xấu, sau đó sẽ cải thiện chất lượng của MSA này theo thời gian thông qua các bước lặp. Một số các giải thuật theo phương pháp này: Giải thuật di truyền, Mô hình Markov ẩn. So sánh Progressive Algorithm và Iterative Algorithm: Cả 2 phương pháp đều có ưu điểm và khuyết điểm, tùy theo điều kiện cũng như yêu cầu cụ thể mà chúng ta có thể lựa chọn phương pháp hợp lý nhất. Iterative Algorithm thông thường cho lời giải có độ chính xác cao hơn tuy nhiên trong một số trường hợp, phương pháp này có thể cho ra các lời giải không tốt. Iterative Algorithm đòi hỏi thời gian tính toán cao, một số trường hợp giải thuật có thể cho thời gian tính toán cao hơn so với cách tiếp cận theo phương pháp tính toán chính xác của Needlman-Wuns(Nicholas, 2002)[23]. Đây là một trong những nhược điểm lớn của phương pháp này. Hình 2.10 Mối tương quan giữa các chương trình hiện thực cho các phương pháp. 2.3.1. Phương pháp Quy hoạch động(Dynamic Programming) Đây là phương pháp phổ biến và cơ bản nhất để tiếp cận bài toán so sánh trình tự. Kỹ thuật dynamic programming giải quyết bài toán dựa trên lý thuyết về Edit distance được đưa ra bởi Vladimir Levenshtein, khoảng cách này định nghĩa số các thao tác nhỏ nhất để có thể chuyển đổi một chuỗi trình tự này thành một chuỗi trình tự khác. Edit distance có thể biểu diễn cho quá trình tiến hóa của một đoạn gene, DNA, protein UPGMA Multalign, Pileup NJ ClustalW GA Saga HMM Hmmt Iterative Progressive Mutal Dialign 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 20 theo thời gian. Điều này hoàn toàn hữu dụng cho lý thuyết về so sánh trình tự. Kỹ thuật dynamic programming có thể tính được edit distance. Bằng việc xây dựng các ma trận đánh giá k chiều với mỗi chiều của ma trận tương ứng với một trong k trình tự của bài toán, phương pháp dynamic programming có thể tìm ra được giá trị hàm đánh giá tối ưu cho ma trận này, từ đó sử dụng kỹ thuật lưu vết để tìm ra được phép so sánh nhiều trình tự tối ưu nhất. Needleman-Wuns đã xây dựng phương pháp tìm MSA tối ưu của k sequence dựa trên việc tổng quát hóa giải thuật quy hoạch động cho 2 sequence. Trong giải thuật này các tác giả đã thay thế việc tính toán điền các giá trị cho ma trận 2 chiều bằng ma trận n chiều. Giải thuật phụ thuộc vào số sequence, chiều dài cũng như sự khác biệt của các sequence. Giải thuật này đòi hỏi độ phức tạp là hàm số mũ theo chiều dài và số trình tự ( )knθ . Trong thực tế việc hiện thực giải thuật này đòi hỏi chi phí tiêu tốn rất cao về mặt thời gian cũng như bộ nhớ sử dụng. Khi n và k đủ lớn việc thực thi giải thuật là không thực tế. Một mini-supercomputer với 4GB bộ nhớ có thể align 20 sequence thuộc nhóm phospholipase A2, mỗi sequence có khoảng 130 amino acid[23]. Đã có một số phương pháp cải thiện giải thuật của Needleman-Wuns, bằng cách áp dụng heuristic vào giải thuật này, tuy nhiên vẫn có hạn chế về giá trị của n, k. 0 0 0 0 0 0000000 G A A T T C G G A T C 2 54310-1 13532-1 -3-11340 -2-2-2-112 -1-1-1-10 Hình 2.11 Phương pháp tính toán chính xác bằng dynamic programming 2.3.2. Sử dụng các thiết bị phần cứng Phương pháp sử dụng thiết bị phần cứng để hiện thực hóa kỹ thuật dynamic programming, hoặc sử dụng các thiết bị phần cứng để song song hóa, chia nhỏ bài toán. Các phương pháp này sẽ giúp thực thi bài toán nhanh hơn tuy nhiên cả hai phương pháp này đều đòi hỏi chi phí quá cao. 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 21 2.3.3. Phương pháp tìm kiếm cục bộ(Local Search) Các kỹ thuật tìm kiếm cục bộ sử dụng các kỹ thuật lặp với các tham số đầu vào xác định, giải quyết bài toán so sánh nhiều trình tự theo hướng tìm ra lời giải gần tối ưu nhất dựa trên các bước di chuyển tốt của của mỗi trạng thái của lời giải so với các trạng thái xung quanh nó trong mỗi vòng lặp xác định. Chủ yếu được sử dụng cho việc giải quyết bài toán so sánh theo hướng toàn cục. Phương pháp này có ưu điểm thừa kế được các ưu điểm của các giải thuật tìm kiếm cục bộ: giảm không gian tìm kiếm của lời giải, tốc độ tìm kiếm nhanh…, tuy nhiên do phương pháp tìm kiếm cục bộ chỉ đảm bảo tính đúng(sound) chứ không đảm bảo tính đầy đủ(complete), và do việc lựa chọn các bước đi để thoát khỏi tình trạng tối ưu cục bộ là ngẫu nhiên nên có thể thuật toán sẽ không tìm ra được lời giải tốt trong một số lần chạy(phụ thuộc vào việc ấn định các tham số đầu vào). Một số giải thuật : Tabu Search, Simulated Annealing(giải thuật luyện kim). 2.3.4. Sử dụng giải thuật Di truyền(Genetic Algorithm) Phương pháp sử dụng thuật giải di truyền cho bài toán so sánh nhiều trình tự là một phương pháp phù hợp với bản chất sinh học [34]. Ý tưởng tạo ra một quần thể các alignment sau đó chọn ra các lời giải tốt nhất có thể. Phương pháp này tạo dựng một quần thể các phép so sánh trình tự từ các phép toán cơ bản của sinh học, lai tạo, đột biến và cố gắng cải thiện hàm thích nghi của quần thể các phép so sánh trình tự. Phương pháp này cho kết quả tốt, tuy nhiên hạn chế của phương pháp này là thời gian tính toán lớn(đặc điểm chung của giải thuật di truyền). Một số phần mềm: SAGA(Notredame và Higgins, 1996)[24]. 2.3.5. Sử dụng Mô hình Markov ẩn(Hidden Markov Model-HMM). Áp dụng mô hình Markov ẩn để mô hình hóa bài toán so sánh nhiều trình tự. Ý tưởng của phương pháp này là sử dụng mô hình Markov ẩn để biểu diễn MSA(được giới thiệu bởi Eddy(1995) [8]), sau đó tối ưu khả năng mà một mô hình HMM có thể biểu diễn cho các sequence đã được align. Trong mô hình này các nucleotide(A,C,T,G) hoặc 23 amino acid sẽ là tập các ký tự. Các trạng thái của mô hình sẽ thuộc 3 loại trạng thái: match, insert, delete. Mỗi ký tự sẽ có 1 xác suất xuất hiện nhất định tại mỗi trạng thái. Giữa các trạng thái sẽ có xác suất chuyển đổi từ trạng thái này sang trạng thái khác. Dựa trên mô hình này mỗi chuỗi trình tự bất kỳ trong 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 22 sinh học sẽ được sinh ra bằng một con đường tập các trạng thái. Tập hợp các trạng thái của các chuỗi trình tự trong mô hình HMM sẽ là 1 kết quả của bài toán so sánh nhiều trình tự. Và như vậy bài toán so sánh nhiều trình tự sẽ trở thành bài toán tìm xác suất điều kiện cực đại của các chuỗi trình tự khi biết mô hình. Hình 2.12 Mô hình Markov cho bài toán MSA. Hình trên minh họa quá trình sinh ra một kết quả của phép so sánh giữa 2 trình tự từ mô hình Markov. Trạng thái match được biểu diễn bằng các hình chữ nhật, trạng thái insert được biểu diễn bởi các hình thoi, trạng thái delete được biểu diễn bởi các hình tròn. Bài toán tìm xác suất cực đại của các chuỗi trình tự trong mô hình có thể được giải quyết bằng nhiều kỹ thuật: giải thuật Viterbi, Forward-Backward Algorithm, Baum-Welch Algorithm… Khác với các phương pháp ở trên phụ thuộc vào cơ chế đánh giá(hàm đánh giá). Mô hình HMM sẽ xây dựng bài toán so sánh nhiều trình tự dựa trên chính bản thân thông tin của các trình tự, sự biến đổi của các thông số đánh giá sẽ do chính bản thân các trình tự quyết định. Phần mềm hiện thực theo phương pháp này: HMMT. Ngoài các phương pháp được trình bày ở trên, vẫn còn một số phương pháp tiếp cận khác để giải quyết bài toán như mạng Neuron … Tựu trung lại các phương pháp đều có điểm mạnh và điểm khuyết của chúng trong việc giải bài toán so sánh nhiều trình tự. Tuy nhiên như tôi đã trình bày không một phương pháp nào có thể gọi là hoàn hảo và chiếm ưu thế tuyệt đối, tất cả các phương pháp đều hướng đến việc tìm ra một lời giải gần tốt nhất với khoảng thời gian A1 A2 A3 − A4 − A5 − − B1 B2 B3 B4 B5 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 23 hợp lý, thay vì hướng đến việc tìm ra lời giải tốt nhất với một khoảng thời gian của một bài toán NP. Trong các phương pháp trên phương pháp dynamic programming là phương pháp ra đời sớm nhất, rất hiệu quả với bài toán so sánh 2 trình tự(PSA), tuy nhiên lại không thật sự tốt với bài toán MSA có dữ liệu lớn, điều này có thể được khắc phục nếu tiếp cận thuật giải theo hướng Progressive Algorithm, xây dựng một sự cải tiến bằng cách kết hợp các kỹ thuật heuristic trong quá trình xử lý, giảm thiểu không gian tìm kiếm, loại bỏ cách tiếp cận dynamic programming truyền thống(xây dựng ma trận đánh giá k chiều), đưa về bài toán xử lý dựa trên việc đánh giá tất cả các bài toán so sánh 2 trình tự được sinh ra từ các trình tự của bài toán MSA. Phương pháp dựa trên Progressive Algorithm này sẽ cho phép giảm độ phức tạp của giải thuật, tăng tốc độ tính toán, giúp tìm kiếm lời giải trong khoảng thời gian ngắn, mà vẫn đáp ứng được lời giải gần tốt nhất trong hầu hết các trường hợp. Trong phạm vi nghiên cứu của mình tôi xin đề xuất phương pháp thực hiện này cho việc giải quyết bài toán MSA. Giải pháp được đưa ra chỉ giới hạn giải quyết bài toán MSA dưới hình thức global alignment. 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 24 Chương 3. CƠ SỞ LÝ THUYẾT VÀ PHƯƠNG PHÁP THỰC HIỆN Chương này xin được giới thiệu về phương pháp quy hoạch động, phương pháp Progressive Algorithm sử dụng heuristic dựa trên quy hoạch động để giải bài toán. 3.1. Giới thiệu về Dynamic Programming Dynamic Programming [7] là phương pháp giải bài toán bằng cách kết hợp các lời giải của các bài toán con. Đặc điểm của giải thuật dynamic programming: ƒ Các bài toán con không độc lập với nhau: chúng có chung các bài toán con nhỏ hơn. ƒ Giải mỗi bài toán con chỉ một lần, và ghi nhớ lời giải đó trong một bảng để truy cập khi cần đến. Một giải thuật dynamic programming được xây dựng qua bốn bước: 1. Xác định cấu trúc của một lời giải tối ưu. 2. Định nghĩa đệ quy cho giá trị của một lời giải tối ưu. 3. Tính giá trị của một lời giải tối ưu từ dưới lên (“bottom-up”). 4. Xây dựng lời giải tối ưu từ các thông tin đã tính. Trong phần này, xin được trình bày tổng quan cách tiếp cận bằng kỹ thuật dynamic programming cho vấn đề so sánh trình tự, từ bài toán so sánh 2 trình tự PSA đến bài toán so sánh nhiều trình tự MSA [25], [27], [30]. Phần trình bày này chỉ xét đối với phép toán Global Alignment. Vì bài toán PSA chính là nền tảng cho kỹ thuật dynamic programming giúp giải quyết bài toán MSA, nên trong phần đầu tôi sẽ trình bày về bài toán cơ bản PSA. 3.2. Bài toán PSA và cách giải quyết bằng kỹ thuật quy hoạch động Phần 2.1.3 đã giới thiệu về bài toán PSA. Phần này sẽ giới thiệu hướng giải quyết bài toán này bằng quy hoạch động. Kỹ thuật này sẽ cho phép chúng ta xây dựng một phép alignment tối ưu. Xét 2 chuỗi trình tự S1 và S2, |S1|=n, |S2|=m, mục đích của chúng ta là tìm kiếm một phép alignment tối ưu cho S1 và S2. 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 25 3.2.1. Giải thuật quy hoạch động cho bài toán PSA Định nghĩa 3.1: Xét định nghĩa về 1 cơ chế đánh giá như sau ƒ Nếu a và b là 2 ký tự đại diện cho các amino acid(nucleotide) trong 2 trình tự. Gọi ( , )a bσ là hàm đánh giá của cặp a, b trong trong ma trận đánh giá ƒ Hàm giá trị của gap phụ thuộc vào hằng số r>0 và chiều dài k ( ( ) *γ = −k r k ). ƒ T là một chuỗi, ta định nghĩa |T| là chiều dài của chuỗi T. Giá trị của phép alignment A của 2 chuỗi S1, S2 với kết quả S’1, S’2: 1 1 2 2( '[ ], '[ ]) *σ ∈ −∑ i K S i S i l r với K1 là tập các phần tử không là gap, |K1|=l1, l2 là tổng chiều dài của gap, ' '1 2 1 2| | | |l l l S S+ = = = ƒ Phép alignment tối ưu là phép alignment có giá trị lớn nhất có thể có của 2 chuỗi. Định nghĩa 3.2: Gọi S(i,j) là giá trị của một phép alignment tối ưu của chuỗi S1i(S1[1],.., S1[i]) và S2j (S2[1],..,S2[j])( 1≤ ≤i n ,1≤ ≤j m ) như vậy S(n,m) sẽ là giá trị của một phép alignment tối ưu của S1 và S2. S(i, j) được định nghĩa như sau: S(0, 0)=0 S(i,0)=S(i-1,0) -r, i>0 S(0,j)=S(0,j-1) -r, j>0 Từ định nghĩa của S(i, j) ta sẽ có công thức tính S(i, j) như sau: S(i,j)= max { S(i-1, j-1) + 1 2( [ ], [ ])S i S jσ , S(i-1, j) - r , S(i, j-1) - r } i>0, j>0 Trong đó S(i-1, j-1) chính là giá trị của một phép alignment tối ưu của 2 chuỗi S1i(S1[1],.., S1[i-1]) và S2j(S2[1],…,S2[j-1]) Hình 3.1 Phương pháp quy hoạch động cho bài toán PSA 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 26 Công thức S(i,j) ở trên cho phép ta dễ dàng vận dụng kỹ thuật quy hoạch động. Ví dụ: Xét 2 chuỗi “ACBCDB” và “CADBD” Hàm đánh giá : 2 ( , ) 1 a b a b a b σ =⎧= ⎨− ≠⎩ Cho r=1 Ta thu được ma trận giá trị S(i,j) của 2 chuỗi trình tự như hình dưới j 0 1 2 3 4 5 i C A D B D 0 0 -1 -2 -3 -4 -5 1 A -1 -1 1 0 -1 -2 2 C -2 1 0 0 -1 -2 3 B -3 0 0 -1 2 1 4 C -4 -1 -1 -1 1 1 5 D -5 -2 -2 1 0 3 6 B -6 -3 -3 0 3 2 Từ ma trận kết quả này ta có S(6,5)=2. Như vậy giá trị của phép alignment tối ưu của 2 chuỗi trên là 2. Từ kết quả này ta sẽ đi tìm các phép alignment tối ưu mà cho S(6,5)=2. Sử dụng kỹ thuật lưu vết ta sẽ có được 3 trong số các kết quả, như hình dưới. j 0 1 2 3 4 5 i C A D B D 0 0 -1 -2 -3 -4 -5 1 A -1 -1 1 0 -1 -2 2 C -2 1 0 0 -1 -2 3 B -3 0 0 -1 2 1 4 C -4 -1 -1 -1 1 1 5 D -5 -2 -2 1 0 3 6 B -6 -3 -3 0 3 2 Bảng 3.1 Bảng kết quả giải thuật quy hoạch động cho bài toán PSA 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 27 Dựa vào các con đường được sinh ra do kỹ thuật lưu vết từ S(m,n) đến S(0,0), các alignment sẽ được sinh ra dựa trên nguyên tắc: ƒ Nếu con đường đi theo hướng đường chéo từ S(i,j) đến S(i-1,j-1) thì 2 ký tự đại diện cho S1[i] và S2[j] sẽ được ghi vào kết quả(đã được align) ƒ Nếu con đường đi theo hướng từ S(i,j) đến S(i-1,j) thì ký tự đại diện cho S1[i] và ‘-‘ sẽ được ghi vào kết quả(1 insertion gap được sinh ra) ƒ Nếu con đường đi theo hướng từ S(i, j) đến S(i, j-1) thì ký tự đại diện cho S2[j] và ‘-‘ sẽ được ghi vào kết quả(1 deletion gap được sinh ra) Theo nguyên tắc này, trong ví dụ trên ta thu được 3 kết quả tương ứng với 3 con đường truy hồi: Như vậy cách tiếp cận quy hoạch động gồm 3 bước: ƒ Khởi tạo ma trận từ 2 chuỗi sequence. ƒ Tính toán, điền giá trị cho ma trận. ƒ Sử dụng kỹ thuật lưu vết để tìm ra kết quả. Giải thuật tính S(m,n) và lưu vết S(0,0)=0; for i=1 to n do S(i,0)= S(i-1,0) - r; for j=1 to m do S( j,0)= S(0, j-1) - r; for i=1 to |S1| do for j=1 to |S2| do { S(i, j)= max ( S(i-1, j-1) + 1 2( [ ], [ ])S i S jσ , S(i-1, j) - r , S(i, j-1) - r );(1) BackPTR(S(i, j))=(imax, jmax);//lưu trữ vết vào mảng BackPTR. } – A C B C D B C A D B – D – A C B C D B – – C A – D B D A C B C D B – – C – A D B D 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 28 Với imax, jmax là giá trị tương ứng của i, j trong biểu thức (1). Xây dựng phương thức tìm alignment theo nguyên tắc nêu ở trên, dựa vào mảng lưu trữ vết BackPTR, ta sẽ thu được kết quả. Giải thuật tìm S(m,n) và lưu vết có độ phức tạp (( 1) ( 1))n mΟ + + + , quá trình tìm alignment từ tập vết BackPTR có độ phức tạp (( 1) ( 1))n mΟ + + + , vì thế độ phức tạp của giải thuật quy hoạch động cho bài toán là: (( 1)( 1))n mΟ + + + (( 1) ( 1))n mΟ + + + = ( )nmΟ . 3.2.2. Giải thuật Gotoh Trong giải thuật trên, chúng ta xét trường hợp gap đơn giản, bỏ qua các tham số khả năng mở gap và khả năng xuất hiện phần tử tiếp theo trong gap . Khi giá trị của gap phụ thuộc vào các tham số này, công thức tính S(i,j) ( i>0, j>0) sẽ có dạng: 1 2( , ) max{ ( 1, 1) ( [ ], [ ]), max{ ( , 1) ( 1)}, max{ ( 1, ), ( 1)}} 2 , 2 S i j S i j S i S j S i x j x S i j y y x i y j σ γ γ = − − + − − + − − − − ≤ ≤ ≤ ≤ Gotoh[13] đưa ra 1 giải thuật cho phép tính khi gap là 1 hàm tuyến tính dạng : ( ) ( * )k q r kγ = − + Để tính toán S chúng ta đưa vào thêm các định nghĩa sau: Gọi D(i, j) là giá trị lớn nhất của alignment giữa S1i và S2j mà được kết thúc bằng một deletion gap. I(i, j) là giá trị lớn nhất của một alignment giữa S1i và S2j mà được kết thúc bằng một insertion gap. Chúng ta sẽ có công thức tính cho S(i, j), D(i, j), I(i, j) S(0, 0) = 0, S(i, 0) = max{D(i, 0) } for i > 0, S(0, j ) = max{I (0, j ) } for j > 0, S(i, j ) = 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) ( , [0, ]i j n∈ )là một điểm nằm trên con đường sinh ra alignment này. Khi đó, SP bằng tổng giá trị alignment tối ưu của Ai với Bj và giá trị alignment tối ưu của TiA với TjB : ( , ) ( , )T Ti j i jSP score A B score A B= + (3) Chứng minh: Thật vậy ta có ( , ) ( , )= = T Ti i j jSP score A B score A A B B (1’) Do M(i,j) thuộc con đường sinh ra phép alignment tối ưu nên ta có thể tách thành 2 bài toán con hay ( , ) ( , ) ( , )= +T T T Ti i j j i j i jscore A A B B score A B score A B (2’) Từ (1’), (2’) ta có điều phải chứng minh Để tìm M*, phương pháp được sử dụng là cố định i* và tìm j*. Ta chọn i* = [m/2]. Chúng ta sẽ tìm * [0, ]j n∈ thỏa mãn yêu cầu đặt ra. Xét 1 giá trị j bất kỳ [0, ]j n∈ . Gọi SS(j) là giá trị của phép alignment tối ưu của Ai* và Bj , DD(j) là giá trị của phép alignment tối ưu của Ai* với Bj mà kết thúc bằng một gap. Gọi RR(j) là giá trị của phép alignment tối ưu của rev(A)m-i* với rev(B)j, và CC(j) là giá trị của phép alignment tối ưu của rev(A)m-i* với rev(B)j kết thúc bằng một gap. 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 45 SS, DD, RR, CC là các vector gồm n phần tử. Định lý 4.2 Giá trị j* cần tìm sẽ thỏa đẳng thức sau: j*=max j∈[0,n] {max (SS( j) + RR(n − j) , DD( j) + CC(n − j) + q) } Chứng minh: Thật vậy: Do rev(rev(A)m-i* )= * T iA và rev(rev(B)n-j )= TjB nên ta có: RR(n-j) là giá trị của alignment tối ưu của * T iA với TjB CC(n-j) là giá trị của alignment tối ưu của * T iA với TjB bắt đầu với một gap Như vậy theo định l ý 4.1, ứng với mỗi giá trị j, giá trị của toàn bộ alignment ứng với kiểu (1): SP=SS(j) +RR(n-j). Giá trị của toàn bộ alignment ứng với kiểu (2) là SP=DD(j)+CC(n-j)+q. Do đó việc tìm giá trị tối ưu của việc chuyển đổi từ A sang B sẽ là việc tìm j* sao cho giá trị max (SS( j*) + RR(n − j*) , DD( j*) + CC(n − j*) + q) là lớn nhất: j*=max j∈[0,n] {max (SS( j) + RR(n − j) , DD( j) + CC(n − j) + q) } M*(i*, j*) chính là cặp điểm giữa tối ưu cần tìm, định lý được chứng minh. Để ý rằng SS, DD có thể tính được bằng cách áp dụng giải thuật Gotoh cải tiến(xem 3.2.3) cho Ai* và B; RR, CC có thể tính được bằng cách áp dụng giải thuật Gotoh cải tiến cho rev(A)m-i* với rev(B). Như vậy quá trình tìm M*(i*, j*) được chia làm 2 phần: ƒ Tìm SS(j), DD(j) bằng cách áp dụng giải thuật Gotoh cải tiến cho Ai* với B.(hướng tới). ƒ Tìm RR, CC bằng cách áp dụng giải thuật Gotoh cải tiến cho rev(A)m-i* với rev(B)(hướng lui). Giá trị tối ưu của phép alignment có thể tìm được bằng cách gọi đệ qui để tìm giá trị tối ưu của Ai* với Bj* , đệ qui để tính giá trị tối ưu của TiA *với T jB * và kết nối 2 phần này lại. Bài toán sẽ chia thành các bài toán con nhỏ hơn(chia để trị) như hình trên(đường nối các điểm giữa trên hình biểu diễn toàn bộ con đường tính giá trị tối ưu). Điều kiện dừng của giải thuật đệ quy là khi i* nhỏ hơn 1 giá trị iM. Khi đó chúng ta sẽ thực hiện giải thuật quy hoạch động với ma trận M gồm iM dòng và jM cột(jM là giá trị tương ứng của j khi i*=iM). Áp dụng kỹ thuật lưu vết trong quy hoạch động kế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 46 hợp với việc sử dụng một cấu trúc dữ liệu lưu trữ kết quả alignment từ các ma trận M, ta sẽ dễ dàng tìm được kết quả alignment của 2 chuỗi A, B. Giải thuật shared vectors SS[0..m] , DD[0..m] , RR[0..n] , CC[0..n] procedure DIFF(A, B, m, n) { diff(A, B, m, n) } recursive procedure diff(A, B, m, n){ if (n<=iM) { Áp dụng kỹ thuật quy hoạch động cho A,B. Tìm vết của giải thuật quy hoạch động. Lưu trữ kết quả. } else{ i*← [m/ 2