Luận án Xây dựng đồ thị tái tổ hợp di truyền cho dữ liệu hệ gen

112 trang thiennha21 14/04/2022 5650

Download

Bạn đang xem 20 trang mẫu của tài liệu "Luận án Xây dựng đồ thị tái tổ hợp di truyền cho dữ liệu hệ gen", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

luan_an_xay_dung_do_thi_tai_to_hop_di_truyen_cho_du_lieu_he.pdf

Nội dung text: Luận án Xây dựng đồ thị tái tổ hợp di truyền cho dữ liệu hệ gen

ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC CÔNG NGHỆ Nguyễn Thị Phương Thảo XÂY DỰNG ĐỒ THỊ TÁI TỔ HỢP DI TRUYỀN CHO DỮ LIỆU HỆ GEN LUẬN ÁN TIẾN SĨ CÔNG NGHỆ THÔNG TIN Hà Nội – 2020
ĐẠI HỌC QUỐC GIA HÀ NỘI TRƯỜNG ĐẠI HỌC CÔNG NGHỆ Nguyễn Thị Phương Thảo XÂY DỰNG ĐỒ THỊ TÁI TỔ HỢP DI TRUYỀN CHO DỮ LIỆU HỆ GEN Chuyên ngành: Khoa học Máy tính Mã số: 9480101.01 LUẬN ÁN TIẾN SĨ CÔNG NGHỆ THÔNG TIN NGƯỜI HƯỚNG DẪN KHOA HỌC: 1.PGS.TS. Lê Sỹ Vinh 2.PGS.TS. Lương Chi Mai Hà Nội – 2020
Lời cam đoan Tôi xin cam đoan đây là công trình nghiên cứu của riêng tôi. Các kết quả được viết chung với các tác giả khác đều được sự đồng ý của các đồng tác giả trước khi đưa vào luận án. Các kết quả nêu trong luận án là trung thực và chưa từng được ai công bố trong các công trình nào khác. Tác giả Nguyễn Thị Phương Thảo
Lời cảm ơn Luận án được thực hiện tại Trường Đại học Công nghệ, Đại học Quốc gia Hà Nội, dưới sự hướng dẫn của PGS. TS. Lê Sỹ Vinh và PGS. TS. Lương Chi Mai. Tôi xin bày tỏ lòng biết ơn sâu sắc tới PGS. TS. Lê Sỹ Vinh, PGS. TS. Lương Chi Mai và TS. Lê Sĩ Quang, những người đã có những định hướng giúp tôi thành công trong việc nghiên cứu của mình. Các Thầy Cô cũng đã động viên và khích lệ tinh thần, giúp tôi vượt qua những khó khăn để tôi hoàn thành được luận án này. Tôi cũng chân thành cảm ơn thầy Hồ Tú Bảo, Thầy đã cho tôi nhiều kiến thức quý báu về nghiên cứu khoa học. Những sự chỉ bảo quý giá của các Thầy Cô đã giúp tôi hoàn thành tốt luận án này. Tôi cũng xin cảm ơn tới các Thầy, Cô thuộc Khoa Công nghệ Thông tin, Trường Đại học Công nghệ, Đại học Quốc gia Hà Nội đã tạo mọi điều kiện thuận lợi giúp tôi trong quá trình làm nghiên cứu sinh. Tôi xin chân thành cảm ơn các đồng nghiệp trong phòng Nhận dạng và Công nghệ Tri thức, Viện Công nghệ Thông tin, Viện Hàn lâm Khoa học và Công nghệ Việt Nam đã luôn động viên, tạo điều kiện thuận lợi, bố trí thời gian tốt nhất cho tôi trong suốt quá trình làm nghiên cứu sinh. Cuối cùng, tôi xin gửi lời cảm ơn sâu sắc tới gia đình và bạn bè, những người đã cho tôi điểm tựa vững chắc để tôi có được thành công như ngày hôm nay. 2
MỤC LỤC Lời cam đoan 1 Lời cảm ơn 2 MỤC LỤC 3 Danh mục các ký hiệu và chữ viết tắt 6 Danh mục các bảng 7 Danh mục các hình vẽ, đồ thị 8 Danh mục các thuật toán 12 MỞ ĐẦU 13 Chương 1. GIỚI THIỆU 16 1.1. Giới thiệu chung 16 1.1.1. Hệ gen người 16 1.1.2. Mạng phát sinh loài 21 1.2. Xây dựng đồ thị tái tổ hợp di truyền 23 1.2.1. Sự kiện tái tổ hợp 23 1.2.2. Đồ thị tái tổ hợp di truyền 25 1.2.3. Bài toán xây dựng đồ thị ARG 32 1.3. Các phương pháp xây dựng đồ thị ARG 35 1.3.1. Các phương pháp xây dựng đồ thị ARG tối thiểu 35 1.3.2. Các phương pháp xây dựng đồ thị ARG hợp lý 39 1.3.3. Tổng hợp các phần mềm xây dựng đồ thị ARG 41 1.4. Ứng dụng ARG trong nghiên cứu tương quan toàn hệ gen 42 3
1.5. Kết luận chương 45 Chương 2. THUẬT TOÁN ARG4WG XÂY DỰNG ĐỒ THỊ TÁI TỔ HỢP DI TRUYỀN HỢP LÝ CHO DỮ LIỆU HỆ GEN 47 2.1. Giới thiệu 47 2.1.1. Các định nghĩa 47 2.1.2. Thuật toán Margarita xây dựng đồ thị ARG 48 2.2. Thuật toán ARG4WG 51 2.2.1. Chiến lược tìm đoạn đầu chung dài nhất 51 2.2.2. Thuật toán ARG4WG 54 2.3. Kết quả thực nghiệm 61 2.3.1. Các kết quả trên dữ liệu thật 61 2.3.2. Các kết quả trên dữ liệu mô phỏng 65 2.4. Kết quả ứng dụng ARG4WG vào bài toán tìm vùng gen liên quan đến bệnh sốt rét ở Châu Phi 67 2.5. Kết luận chương 72 Chương 3. PHƯƠNG PHÁP TỐI ƯU HÓA SỐ SỰ KIỆN TÁI TỔ HỢP TRONG QUÁ TRÌNH XÂY DỰNG ĐỒ THỊ ARG 75 3.1. Giới thiệu 75 3.2. Một số định nghĩa và khái niệm sử dụng trong các thuật toán 76 3.3. Hạn chế của thuật toán ARG4WG 78 3.4. Thuật toán REARG 79 3.4.1. Động cơ nghiên cứu 79 3.4.2. Thuật toán REARG 80 4
3.5. Thuật toán GAMARG 83 3.5.1. Động cơ nghiên cứu 83 3.5.2. Thuật toán GAMARG 83 3.6. Kết quả thực nghiệm 88 3.6.1. Kết quả trên các tập dữ liệu nhỏ 89 3.6.2. Các kết quả trên các tập dữ liệu từ dự án 1kGP 90 3.7. Kết luận chương 98 KẾT LUẬN 100 DANH MỤC CÁC CÔNG TRÌNH KHOA HỌC CỦA TÁC GIẢ LIÊN QUAN ĐẾN LUẬN ÁN 102 TÀI LIỆU THAM KHẢO 103 5
Danh mục các ký hiệu và chữ viết tắt D Tập các trình tự N Số lượng trình tự trong một tập các trình tự m độ dài của trình tự Sx Trình tự thứ x trong một tập các trình tự Sx[i] Giá trị của Sx tại vị trí thứ i ARG Đồ thị tái tổ hợp di truyền 1KGP Dự án 1000 hệ gen GWAS Nghiên cứu tương quan toàn hệ gen SNP Đa hình đơn nucleotit MRCA Tổ tiên chung gần nhất CwR Mô hình kết hợp và tái tổ hợp STT Số thứ tự RF Khoảng cách Robinson-Fould 6
Danh mục các bảng Bảng 1.1: Các phần mềm xây dựng đồ thị ARG tiêu biểu. 41 Bảng 2.1: Tập dữ liệu trích xuất từ dự án 1000 hệ gen người. 62 Bảng 3.1: Tập dữ liệu từ dự án 1kGP. 89 Bảng 3.2: Các kết quả của các thuật toán khác nhau trên các tập dữ liệu nhỏ. 89 Bảng 3.3: Số sự kiện tái tổ hợp ít nhất được tìm thấy bởi 5 thuật toán cho 100 trình tự của (a) DS1, (b) DS2 và (c) DS3. 91 Bảng 3.4: Số sự kiện tái tổ hợp ít nhất được tìm thấy bởi 5 thuật toán cho 200 trình tự của (a) DS1, (b) DS2 và (c) DS3. 92 Bảng 3.5: Trung bình thời gian chạy (giây) của 5 thuật toán cho 100 trình tự của các tập dữ liệu (a) DS1, (b) DS2, và (c) DS3. 95 Bảng 3.6: Trung bình thời gian chạy (giây) của 5 thuật toán cho 200 trình tự của các tập dữ liệu (a) DS1, (b) DS2, và (c) DS3. 97 7
Danh mục các hình vẽ, đồ thị Hình 1.1: Cấu trúc hệ gen người. Hệ gen người gồm 23 cặp nhiễm sắc thể, có khoảng 3 tỉ phân tử DNA, khoảng 20.000 đến 25.000 gen. Nguồn hình: 16 Hình 1.2: Các kiểu biến thể trình tự: (a) Thay thế một cặp bazơ đơn. Trong ví dụ, biến thể xuất hiện ở 2 vị trí so với trình tự tham chiếu, đó là thay thế nucleotit T↔A và G↔A. (b) Chuỗi GCA được chèn vào so với trình tự tham chiếu. (c) Chuỗi CG bị xóa so với trình tự tham chiếu. 17 Hình 1.3: Các loại biến thể cấu trúc: xóa, thêm, lặp, đảo hay lặp nhiều lần 1 đoạn DNA. Đoạn đột biến cấu trúc có kích thước lớn hơn 1kb. 18 Hình 1.4: Ví dụ dữ liệu SNP chứa biến thể 2 alen và nhiều alen. Có 8 vị trí SNP đều là 2 alen, gồm alen tham chiếu và 1 alen biến thể, ví dụ như A và G ở vị trí 1; T và C ở vị trí 2. Chỉ có vị trí 7 là 3 alen: alen tham chiếu (G) và 2 alen biến thể C, T. 19 Hình 1.5: Ví dụ 4 haplotype của 4 cá thể trên một vùng gen. Một haplotype được tạo thành từ sự kết hợp của các SNP được di truyền cùng nhau trong các đoạn DNA. 19 Hình 1.6: Cây phân loài biểu diễn mối quan hệ tiến hóa của một số loài linh trưởng. Đười ươi và Khỉ đột rẽ nhánh sớm hơn các loài linh trưởng khác. Con người rẽ ra một nhánh riêng và nhánh còn lại cho ra Tinh tinh và vượn Bonobo. 21 Hình 1.7: Khái quát hóa các mạng phát sinh loài điển hình [36]. 23 Hình 1.8: Hai hiện tượng tái tổ hợp phổ biến của người: (a) trao đổi chéo và (b) chuyển đổi gen. 24 Hình 1.9: Biến đổi dữ liệu SNP thành dạng nhị phân. Vị trí có giá trị giống với tham chiếu là 0, giá trị khác tham chiếu là 1. 28 8
Hình 1.10: Đồ thị ARG cho tập dữ liệu M gồm 7 trình tự độ dài 5 [26]. Trình tự tổ tiên là “00000”; 5 sự kiện đột biến tại các vị trí tương ứng (1,2,3,4,5) được ghi trên các cạnh xảy ra đột biến của đồ thị; 2 sự kiện tái tổ hợp xảy ra tại vị trí 3 và 4. 29 Hình 1.11: Điểm cắt tái tổ hợp. 30 Hình 1.12: Một ví dụ đồ thị ARG cho 4 trình tự với các ký hiệu: ■: trạng thái di truyền, ◘: trạng thái di truyền đột biến, □: trạng thái không xác định. 31 Hình 1.13: Các cây thành phần (đường đậm nét) của đồ thị ARG trong Hình 1.12. Nguồn hình [43]. 33 Hình 1.14: (a) Ví dụ cặp vị trí tương thích: cặp vị trí này chỉ chứa 3 loại giao tử và có thể có được từ 1 tổ tiên chung thông qua 2 sự kiện đột biến. (b) Cặp vị trí không tương thích: cặp vị trí chứa 4 loại giao tử và trong trường hợp này phải có ít nhất 1 sự kiện tái tổ hợp xảy ra dưới giả định các vị trí vô hạn (kí hiệu * biểu thị vị trí không có thông tin). 36 Hình 1.15: Một cây có nốt sùi cho tập trình tự giống với tập trong Hình 1.10 với 2 nốt sùi tương ứng với 2 chu trình tái tổ hợp không chung nút với nhau [27]. 38 Hình 1.16: (a) Đồ thị ARG cho tập 4 trình tự, trong đó trình tự s1, s2 là từ 2 cá thể khỏe mạnh, trình từ s3, s4 là từ 2 cá thể bị bệnh. (b) Đột biến 3 (vùng khoanh tròn) trên cây biên tại vị trí 3 của đồ thị ARG trong (a) cho ra sự phân biệt rõ nhất giữa các trình tự bệnh và trình tự không bệnh. 44 Hình 2.1: Lưu đồ thuật toán Margarita. 49 Hình 2.2: Vấn đề trong việc thực hiện sự kiện tái tổ hợp của Margarita. Hai trình tự S1 và S2 với đoạn chung dài nhất giữa hai trình tự được biểu diễn bằng đoạn màu đen. Thuật toán thực hiện lần lượt 2 sự kiện tái tổ hợp R1 và R2 trên trình tự S1 để sinh ra 3 trình tự con S11, S12 và S13. Sau đó, trình tự con chứa đoạn chung dài nhất S13 sẽ được kết hợp với S2. Vì vậy, khi đoạn chung dài nhất được tìm thấy bên trong 9
trình tự, thuật toán phải thực hiện 2 sự kiện tái tổ hợp trên một trình tự và từ 2 trình tự ban đầu (S1 và S2) sẽ thành 3 trình tự ở thế hệ tiếp theo (S11, S12 và S' (S' = S2)). 50 Hình 2.3: Tất cả các trình tự con từ phía bên trái của s mà có thể kết hợp với một trình tự trong D là một tập con của đoạn bên trái dài nhất của s ( sl ) 52 Hình 2.4: Phân tách s bằng cách chọn các đoạn chung dài nhất trong s để kết hợp với các trình tự trong D có thể không dẫn tới số cực tiểu sự kiện tái tổ hợp. 53 Hình 2.5: Sự kiện tái tổ hợp được biểu thị trong thuật toán ARG4WG. (a) Xét 2 trình tự S1 và S2, các đoạn đầu chung của 2 trình tự từ phía bên trái (hình lượn sóng) và từ phía bên phải (màu đen) được xác định. (b) Với 1 tập 3 trình tự S1, S2 và S3, các đoạn đầu chung của mỗi cặp được tính toán (hình lượn sóng) và đoạn đầu chung dài nhất được xác định được mô tả bằng đoạn màu đen giữa trình tự S1 và S2. (c) Một sự kiện tái tổ hợp được thực hiện trên trình tự S1 để sinh ra 2 trình tự con S11 và S12. S12 chứa đoạn đầu chung dài nhất sau đó sẽ được kết hợp với S2. Như vậy, ARG4WG luôn thực hiện 1 tái tổ hợp trên 1 trình tự và từ 2 trình tự ban đầu (S1, S2) sẽ thành 2 trình tự ở thế hệ tiếp theo (S11, S’), trong đó S’ = S2 và S11 có ít vật liệu di truyền hơn S1. 55 Hình 2.6: Trung bình thời gian chạy của Margarita, Margarita1.0 và ARG4WG cho: (a) 500 haplotype; (b) 1000 haplotype; và (c) 2000 haplotype. 63 Hình 2.7: Trung bình số sự kiện tái tổ hợp của Margarita, Margarita1.0 và ARG4WG cho: (a) 500 haplotype; (b) 1000 haplotype; và (c) 2000 haplotype. 65 Hình 2.8: Khoảng cách RF của các cây được tạo ra bởi thuật toán Margarita và ARG4WG so với các cây đúng tương ứng trên các khoảng tỉ lệ đột biến và tái tổ hợp khác nhau. 67 Hình 2.9: Sự tương quan đến bệnh từ 106 kiểm định hoán vị trên: (A) 10 ARG xây dựng trên toàn bộ NST 11; (B) 30 ARG xây dựng trên vùng 5000 SNP quanh gen 10
HBB; và (C) Tổng hợp kết quả cho các thực nghiệm trên vùng 1000 SNP quanh gen HBB. 70 Hình 2.10: Sự tương quan với bệnh khi sử dụng thuật toán Margarita trên vùng 4M- 6M quanh gen HBB. 72 Hình 3.1: Một ví dụ đồ thị ARG tối thiểu cho tập dữ liệu D(5) gồm 5 trình tự độ dài 5. Xét ngược chiều thời gian, thứ tự thực hiện các sự kiện đột biến, kết hợp hay tái tổ hợp để xây dựng đồ thị ARG được đánh số trong hình tròn. Trong ví dụ này, sự kiện tái tổ hợp được thực hiện đầu tiên trên trình tự “01010” sinh ra 2 trình tự “01 ” và “ 010”. Tiếp theo là sự kiện kết hợp trình tự “ 010” và “00010” thành trình tự “00010”. Sự kiện đột biến được thực hiện sau đó biến đổi trình tự “00010” thành trình tự “00000”. Quá trình xây dựng đồ thị ARG được tiếp tục thực hiện cho tới khi tổ tiên chung “10001” được tìm thấy. 77 Hình 3.2: Quá trình xây dựng đồ thị ARG cho tập dữ liệu D={S1,S2,S3,S4,S5} của thuật toán ARG4WG. 푅푖, 푗 biểu thị một sự kiện tái tổ hợp giữa vị trí i và vị trí j; biểu thị sự kiện kết hợp thứ x; 푖 biểu thị một sự kiện đột biến tại vị trí i. 79 Hình 3.3: Thuật toán ARG4WG xác định được 3 cặp ứng cử viên có cùng đoạn đầu chung dài nhất cho tập 5 trình tự như trong các khung hình chữ nhật. Một trong 3 cặp sẽ được chọn ngẫu nhiên để thực hiện tái tổ hợp. 80 Hình 3.4: Cho tập dữ liệu D={S1,S2,S3,S4,S5}, lựa chọn thực hiện tái tổ hợp trên trình tự S4 giữa vị trí 1 và vị trí 2 dẫn đến việc phải thực hiện thêm 1 sự kiện tái tổ hợp nữa để phá vỡ cặp vị trí không tương thích (1,2). 84 Hình 3.5: Lưu đồ thuật toán GAMARG. 85 Hình 3.6: Số sự kiện tái tổ hợp ít nhất được tìm thấy bởi 3 thuật toán ARG4WG, REARG và GAMARG cho 100 và 200 trình tự với 2000, 5000, và 10000 SNP của tập DS1, DS2, và DS3. 94 11
Danh mục các thuật toán Thuật toán 2.1: Thuật toán ARG4WG xây dựng một ARG từ một tập trình tự D cho trước. 60 Thuật toán 3.1: Thuật toán REARG. 81 Thuật toán 3.2: Thuật toán GAMARG. 87 12
MỞ ĐẦU Những thành tựu gần đây trong công nghệ giải trình tự hệ gen thế hệ mới (Next Generation Sequencing - NGS) đã giảm đáng kể chi phí giải trình tự toàn bộ hệ gen và dẫn đến sự gia tăng nhanh chóng về số lượng chuỗi DNA và cả hệ gen. Do đó, việc phát triển các mô hình và phương pháp tính toán mới là vấn đề thời sự đang được các nhà tin sinh học quan tâm để có thể quản lý, phân tích và khai thác bộ dữ liệu lớn các trình tự này một cách hiệu quả. Những dữ liệu này đại diện cho một nguồn thông tin rất hữu ích và đặt ra các vấn đề tính toán mới trong các nghiên cứu trên toàn hệ gen, điển hình là các nghiên cứu về phân bố của các biến thể di truyền trong một quần thể hay xác định các vùng gen có tác động và có ý nghĩa về mặt sinh học đối với các đặc điểm quan trọng mà ta quan tâm. Để giải quyết những bài toán này đòi hỏi nhiều phương pháp mới, trong đó có những hướng đi mới sử dụng lý thuyết đồ thị và thuật toán để mô hình hóa và tính toán các mô hình tiến hóa trong quần thể. Đáng chú ý trong số đó là đồ thị tái tổ hợp di truyền (Ancestral Recombination Graph - ARG), một công cụ quan trọng trong nghiên cứu di truyền quần thể và các bài toán liên quan đến tìm sự đa dạng trong hệ gen [1,58]. Với một tập các chuỗi nhiễm sắc thể, đồ thị ARG đầy đủ sẽ mô tả một cách chi tiết lịch sử di truyền, mối quan hệ của chúng với nhau và với một tổ tiên chung (common ancestor) thông qua ba sự kiện: đột biến (mutation), tái tổ hợp (recombination) và kết hợp (coalescence). Trong quá trình xây dựng đồ thị ARG, sự kiện tái tổ hợp và sự kiện đột biến là 2 sự kiện cốt lõi ảnh hưởng tới đồ thị kết quả, từ đó ảnh hưởng trực tiếp tới các ứng dụng liên quan như tìm vùng gen liên quan đến bệnh, đột biến gây bệnh, đặc trưng của quần thể quan sát, Tuy nhiên, số sự kiện tái tổ hợp và sự kiện đột biến cũng như vị trí thực sự xảy ra trong quá trình tiến hóa là không biết trước. Do đó, chúng ta không thể biết được ARG thực sự mà chúng ta chỉ có thể suy diễn chúng từ dữ liệu với các giả định tối ưu số sự kiện tái tổ hợp và sự kiện đột biến nhằm có được ARG với các sự kiện sát nhất với thực tế. 13
Nhiều phương pháp xây dựng đồ thị ARG đã được đề xuất [26], tập trung vào 2 cách tiếp cận chính: (1) xây dựng đồ thị ARG tối thiểu (minimal ARG), tức là đồ thị có chính xác số sự kiện tái tổ hợp nhỏ nhất; và (2) xây dựng đồ thị ARG hợp lý (plausible ARG), tức là đồ thị có số sự kiện tái tổ hợp tùy thuộc vào thuật toán xấp xỉ chúng. Tuy nhiên, các phương pháp xây dựng đồ thị ARG hiện tại vẫn gặp những hạn chế sau: - Đa số các phương pháp xây dựng đồ thị ARG mới chỉ giới hạn với những tập dữ liệu nhỏ và vừa hàng trăm trình tự [52,58,62]. - Các phương pháp xây dựng đồ thị ARG với hàm mục tiêu có chính xác số sự kiện tái tổ hợp ít nhất hiện thời còn tốn rất nhiều thời gian và chỉ khả thi với những tập dữ liệu rất nhỏ chứa vài chục trình tự [62,71]. Ngày nay, những thành tựu trong công nghệ giải trình tự gen thế hệ mới, sự phát triển và ngày càng hoàn thiện của các thư viện đặc tả biến dị di truyền trong quần thể người đã tạo tiền đề cho các nghiên cứu trên toàn hệ gen. Để có thể ứng dụng được vào các nghiên cứu về biến thể di truyền liên quan đến bệnh ở người một cách hiệu quả, các phương pháp phải có khả năng tính toán được trên dữ liệu liên quan đến hàng nghìn hệ gen. Từ đó, mục tiêu và kết quả của luận án đã đạt được là: 1. Nghiên cứu các phương pháp xây dựng đồ thị ARG hiện tại, từ đó đề xuất thuật toán gần đúng xây dựng đồ thị ARG cho hàng nghìn trình tự, thậm chí hàng nghìn hệ gen nhằm ứng dụng hiệu quả vào các bài toán thực tế trên các tập dữ liệu lớn. 2. Đề xuất thuật toán xây dựng đồ thị ARG với hàm mục tiêu tối thiểu hóa số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG bằng việc kết hợp thuật toán đề xuất trong (1) với một số đặc trưng của dữ liệu và kĩ thuật tối ưu được sử dụng trong các phương pháp tìm cận dưới tái tổ hợp và các phương pháp xây dựng đồ thị ARG tối thiểu đã có. 14
Các kết quả của luận án đã được công bố trong 1 bài tạp chí ISI (công trình khoa học số 1) và 2 báo cáo hội nghị quốc tế (công trình khoa học số 2 và 3). Ngoài phần kết luận, luận án được tổ chức như sau: Chương 1 đầu tiên giới thiệu khái quát về hệ gen người và các mạng phát sinh loài (phylogenetic networks). Sau đó là phần giới thiệu về bài toán xây dựng đồ thị ARG. Phần cuối của chương trình bày các cách tiếp cận giải bài toán xây dựng đồ thị ARG và ứng dụng của ARG trong nghiên cứu tương quan toàn hệ gen. Chương 2 đề xuất một thuật toán xây dựng đồ thị ARG cho dữ liệu lớn hàng nghìn trình tự độ dài hệ gen người. Để làm được điều đó, chúng tôi đưa ra các nhược điểm của các cách tiếp cận hiện có, đặc biệt là những hạn chế trong thuật toán Margarita xây dựng đồ thị ARG hợp lý được đề xuất bởi Minichiello và Durbin [52], từ đó đưa ra thuật toán đề xuất nhằm khắc phục các nhược điểm đó. Các kết quả thực nghiệm ở phần sau của chương đã chứng tỏ hiệu quả của thuật toán đề xuất. Phần cuối của chương giới thiệu kết quả ứng dụng thuật toán đề xuất vào bài toán tìm vùng gen liên quan đến bệnh sốt rét ở Châu Phi trên tập dữ liệu lớn gồm 5560 trình tự trên toàn nhiễm sắc thể 11. Các kết quả trong phần này đã khẳng định thêm hiệu quả, khả năng ứng dụng của thuật toán đề xuất trong các bài toán thực tế trên dữ liệu lớn. Chương 3 của luận án giới thiệu các phương pháp nhằm cực tiểu hóa số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Cụ thể, chúng tôi đề xuất hai phương pháp: (1) kết hợp một số đặc trưng của dữ liệu; (2) kết hợp kĩ thuật sử dụng trong các phương pháp xây dựng đồ thị ARG tối thiểu với chiến lược thực hiện sự kiện tái tổ hợp đề xuất trong chương 2 để tối ưu hóa số sự kiện tái tổ hợp. Các thực nghiệm trên các bộ dữ liệu khác nhau đã chứng tỏ hiệu quả của các phương pháp đề xuất. 15
Chương 1. GIỚI THIỆU 1.1. Giới thiệu chung Trong phần này luận án sẽ giới thiệu về hệ gen người, cụ thể là cấu trúc bộ gen người, các nguyên nhân dẫn tới các biến thể di truyền ở người, các loại biến thể di truyền phổ biến và một số loại dữ liệu hệ gen quan trọng. Luận án cũng giới thiệu sơ lược về các loại mạng phát sinh loài (phylogenetic networks), một công cụ quan trọng để biểu diễn các mối quan hệ tiến hóa trong nghiên cứu di truyền quần thể. 1.1.1. Hệ gen người Bộ gen người là tất cả vật liệu di truyền của một người được di truyền từ thế hệ này sang thế hệ khác. Bộ gen chứa các gen, mỗi gen là một đoạn DNA (deoxyribonucleic acid) mã hóa cho những sản phẩm riêng lẻ như các mRNA được sử dụng trực tiếp cho tổng hợp các enzim, các protein cấu trúc hay các chuỗi polypeptide để gắn lại tạo ra protein có hoạt tính sinh học. Các gen được đóng gói trong nhiễm sắc thể, nhiễm sắc thể nằm trong nhân tế bào, mỗi nhân tế bào có 23 cặp nhiễm sắc thể (Hình 1.1) [54]. Hình 1.1: Cấu trúc hệ gen người. Hệ gen người gồm 23 cặp nhiễm sắc thể, có khoảng 3 tỉ phân tử DNA, khoảng 20.000 đến 25.000 gen. Nguồn hình: 16
Chuỗi DNA người có cấu trúc xoắn kép, gồm 2 chuỗi DNA đơn cuộn xoắn vào nhau. Tại mỗi vị trí cụ thể trên một chuỗi DNA đơn là một trong 4 loại nucleotit: A, T, C, hoặc G. Hệ gen người có khoảng 3 tỉ phân tử DNA, trong đó có các vùng chứa thông tin mã hóa protein gọi là các gen. Con người có từ 20.000 đến 25.000 gen. Hầu hết các gen ở mọi người là như nhau, nhưng có khoảng 0.1% vị trí mà các nucleotit là khác nhau giữa 2 người gọi là các biến thể di truyền. Biến thể di truyền giúp cho mỗi người chúng ta là một cá thể duy nhất, không giống bất kì ai [54]. Đột biến và tái tổ hợp là 2 nguyên nhân chính của biến thể di truyền [22]. Đột biến là nguồn gốc của biến thể mới, xảy ra khi có lỗi trong quá trình sao chép DNA mà không được sửa chữa bởi các enzyme sửa chữa DNA. Trong khi tái tổ hợp di truyền là nguyên nhân chính của biến thể di truyền ở thế hệ con cái. Mỗi người có sự pha trộn các vật liệu di truyền từ cha mẹ. Tái tổ hợp góp phần vào biến đổi gen bằng cách xáo trộn DNA của cha mẹ và tạo ra các tổ hợp biến thể mới. Chi tiết về sự kiện tái tổ hợp được giới thiệu trong Mục 1.2.1. Hình 1.2: Các kiểu biến thể trình tự: (a) Thay thế một cặp bazơ đơn. Trong ví dụ, biến thể xuất hiện ở 2 vị trí so với trình tự tham chiếu, đó là thay thế nucleotit T↔A và G↔A. (b) Chuỗi GCA được chèn vào so với trình tự tham chiếu. (c) Chuỗi CG bị xóa so với trình tự tham chiếu. Biến thể di truyền có thể được phân loại thành biến thể trình tự và biến thể cấu trúc [20,57]. Các biến thể trình tự gồm dạng thay thế một cặp bazơ (base pair, viết tắt là bp) hay còn gọi là đa hình đơn nucleotit (Single Nucleotide Polymorphisms – SNP) và xóa hoặc thêm một đoạn DNA kích thước nhỏ hơn 1kb (1kb = 1000 bp) (Hình 1.2). Các trường hợp chèn và xóa phạm vi lớn hơn, cũng như các trường hợp đảo 17
ngược (inversion) hay lặp lại 2 lần (duplication) hoặc nhiều lần (copy-number variant) 1 đoạn DNA được gọi chung là các biến thể cấu trúc (Hình 1.3). Biến thể cấu trúc thường làm thay đổi cấu trúc của hệ gen, cấu trúc của protein tương ứng. Đoạn DNA biến thể có kích thước từ 1kb đến hơn 5Mb (1Mb = 106 bp). Hình 1.3: Các loại biến thể cấu trúc: xóa, thêm, lặp, đảo hay lặp nhiều lần 1 đoạn DNA. Đoạn đột biến cấu trúc có kích thước lớn hơn 1kb. Biến thể SNP là loại biến thể di truyền phổ biến nhất trong hệ gen người. Một biến đổi điểm có tần số xuất hiện trong quần thể lớn hơn 1% thì được gọi là SNP. Dữ liệu SNP Các dự án hệ gen người [12,13,40] đã chỉ ra có khoảng 10 triệu SNP trong hệ gen người và chúng đóng vai trò như là các dấu hiệu sinh học giúp phân biệt sự khác nhau giữa người với người. Chúng giải thích cho sự khác nhau về màu mắt, màu tóc, nhóm máu của con người. Một số SNP có thể ảnh hưởng tới nguy cơ phát triển một số bệnh hay rối loạn nào đó. Dữ liệu SNP đóng một vai trò đặc biệt quan trọng trong các nghiên cứu tương quan toàn hệ gen (Genome-Wide Association Study – GWAS) nhằm so sánh các vùng trong hệ gen người để định vị vùng gen và các biến thể di truyền có ảnh hưởng tới sức khỏe hay liên quan đến bệnh quan tâm, từ đó giúp cho quá trình chẩn đoán và điều trị [4,6,49,67]. 18
Hầu hết SNP ở người là 2 alen (biallelic SNP), tức là các vị trí SNP chỉ chứa alen tham chiếu và alen biến thể, chiếm đến hơn 99% tổng số SNP [8]. Ngoài ra còn có một số ít các SNP đa alen (multiallelic SNP), là các vị trí SNP chứa alen tham chiếu và 2 hoặc nhiều alen biến thể (Hình 1.4). Hình 1.4: Ví dụ dữ liệu SNP chứa biến thể 2 alen và nhiều alen. Có 8 vị trí SNP đều là 2 alen, gồm alen tham chiếu và 1 alen biến thể, ví dụ như A và G ở vị trí 1; T và C ở vị trí 2. Chỉ có vị trí 7 là 3 alen: alen tham chiếu (G) và 2 alen biến thể C, T. Dữ liệu haplotype Hình 1.5: Ví dụ 4 haplotype của 4 cá thể trên một vùng gen. Một haplotype được tạo thành từ sự kết hợp của các SNP được di truyền cùng nhau trong các đoạn DNA. 19
Haplotype là một nhóm các gen trong một sinh vật được di truyền cùng nhau từ bố hoặc mẹ của chúng. Nhóm gen này được di truyền cùng nhau do liên kết di truyền, hoặc hiện tượng các gen gần nhau trên cùng một nhiễm sắc thể thường được di truyền cùng nhau. Ngoài ra, thuật ngữ "haplotype" cũng còn được đề cập đến là nhóm các SNP được di truyền cùng nhau trong các đoạn DNA [13] (Hình 1.5). Dữ liệu SNP haplotype này là dữ liệu quan trọng trong các nghiên cứu di truyền quần thể và là dữ liệu đầu vào cho bài toán xây dựng đồ thị ARG. Trong hệ gen người (và các loài lưỡng bội nói chung), mỗi người có 2 haplotype trong một vùng xác định của hệ gen. Dữ liệu kiểu gen (genotype) Kiểu gen của một cá thể là tập hợp tất cả các alen – những dạng biến dị khác nhau của cùng một gen ở cá thể đó ( 234). Với các loài lưỡng bội, mỗi vị trí gen c sẽ có 2 alen. Nếu trạng thái alen tại vị trí c là P và Q, kí hiệu "P/Q" chỉ kiểu gen tại vị trí đó. Một vị trí được gọi là đồng hợp tử (homozygous) nếu kiểu gen tại vị trí đó mang 2 alen giống nhau, và được gọi là dị hợp tử (heterozygous) nếu kiểu gen tại vị trí đó mang 2 alen khác nhau. Ví dụ, ta có kiểu gen tại 5 vị trí tương ứng của 1 cá thể X là: A/A, G/A, C/C, T/T, A/T. Vị trí thứ 1, 3 và 4 được gọi là đồng hợp tử còn vị trí thứ 2 và vị trí thứ 5 được gọi là dị hợp tử. Nếu chỉ biết dữ liệu kiểu gen, ta không thể suy luận được 2 haplotype của cá thể X này vì sẽ có 2 cặp haplotype phù hợp với dữ liệu kiểu gen này do 2 vị trí dị hợp tử: Cặp 1: A G C T A Cặp 2: A A C T A A A C T T A G C T T 20
Bài toán tìm haplotype khi cho trước dữ liệu kiểu gen cũng như bài toán xác định kiểu gen và haplotype cho dữ liệu hệ gen thu được từ máy giải trình tự gen thế hệ mới là các bài toán đặc biệt quan trọng trong tin sinh [5,46,51]. 1.1.2. Mạng phát sinh loài Theo học thuyết tiến hóa của Darwin tất cả các loài sinh vật đều tiến hóa từ một tổ tiên chung. Mối quan hệ giữa các loài sinh vật được biểu diễn bởi một cây, gọi là cây phân loài (phylogenetic tree) với cấu trúc như Hình 1.6 [14,59]: Hình 1.6: Cây phân loài biểu diễn mối quan hệ tiến hóa của một số loài linh trưởng. Đười ươi và Khỉ đột rẽ nhánh sớm hơn các loài linh trưởng khác. Con người rẽ ra một nhánh riêng và nhánh còn lại cho ra Tinh tinh và vượn Bonobo. • Mỗi nút lá của cây biểu diễn cho một loài sinh vật hiện tại. • Mỗi nút bên trong của cây biểu diễn cho một loài sinh vật tổ tiên. Thông thường, chúng ta không có thông tin về các loài sinh vật tổ tiên này. • Một cạnh của cây nối hai nút của cây và biểu diễn mối quan hệ trực tiếp giữa hai loài sinh vật ở hai nút của cây. 21
• Độ dài của cạnh nối hai loài sinh vật trên cây cho biết khoảng cách tiến hóa giữa chúng. Khoảng cách này có thể được biểu diễn bằng thời gian, hay số lượng các biến đổi nucleotit giữa hai chuỗi DNA được sử dụng để so sánh hai loài. Cây phân loài là mô hình cơ bản nhất để biểu diễn quan hệ tiến hóa của các loài hoặc các gen. Tuy nhiên, mô hình dạng cây không thể biểu diễn các thông tin và hiện tượng sinh học khác như chuyển gen ngang (horizontal gene transfer), tái tổ hợp (recombination) hoặc lai ghép (hybridization). Trong những trường hợp đó, một số nhánh của cây kết hợp thành một nút mắt lưới (reticulation node) và cây trở thành mạng phát sinh loài (phylogentic network) [36]. Mạng phát sinh loài đang trở thành một công cụ quan trọng trong tiến hóa phân tử. Mạng phát sinh loài là đồ thị bất kì được sử dụng để biểu diễn các mối quan hệ tiến hóa (bằng các cạnh) giữa một tập hợp các nhãn (taxa) (bằng các nút lá) [37]. Với sự đa dạng dữ liệu sinh học hiện có, rất nhiều loại mạng phát sinh loài khác nhau đã ra đời. Có khoảng 20 loại mạng phát sinh loài khác nhau [36]. Một số mạng được đặt tên bởi các thuật toán tính toán chúng hoặc bởi các đặc tính toán học mà định nghĩa chúng, ví dụ như “neighbor-nets” hoặc “median networks”. Một số mạng khác được đặt tên theo các loại sự kiện tiến hóa mà họ mô hình hóa, ví dụ như “hybridization networks”, “recombination networks” hay “duplication-loss-transfer (DLT) networks”. Hình 1.7 minh họa một số mạng tổng quát hiện có. Mỗi mạng có vai trò khác nhau: cây phân loài mô tả mối quan hệ giữa các loài hoặc các gen; mạng phân tách mô tả sự khác nhau giữa các cây phát sinh loài; các sự kiện lai ghép hay tái tổ hợp được mô hình hóa trong các mạng lai ghép hay các mạng tái tổ hợp, Trong đó, sự kiện tái tổ hợp là sự kiện quan trọng thu hút được nhiều sự quan tâm của các nhà nghiên cứu, đặc biệt trong di truyền quần thể. Việc phân tích và xác định được các sự kiện tái tổ hợp giúp cho quá trình xác định đa dạng di truyền, tìm 22
hiểu các nguyên nhân dẫn đến các bệnh đa yếu tố như bệnh tiểu đường, ung thư, và là nền tảng nghiên cứu thuốc chữa bệnh [6]. Hình 1.7: Khái quát hóa các mạng phát sinh loài điển hình [36]. Trong luận án này, đối tượng nghiên cứu là đồ thị tái tổ hợp di truyền (đồ thị ARG), một loại mạng phát sinh loài mô hình hóa quan hệ di truyền giữa các trình tự của các cá thể được quan sát trong một quần thể khi có sự kiện tái tổ hợp xảy ra trong lịch sử tiến hóa của chúng. 1.2. Xây dựng đồ thị tái tổ hợp di truyền 1.2.1. Sự kiện tái tổ hợp Tái tổ hợp là một thành phần cơ bản trong quá trình truyền DNA từ trình tự này sang trình tự khác khi các nhiễm sắc thể được truyền từ thế hệ này sang thế hệ khác. Có 2 kiểu tái tổ hợp phổ biến là trao đổi chéo (crossing over) và chuyển đổi gen (gene conversion). Mỗi loài sinh vật có cơ chế tái tổ hợp khác nhau. Đối với loài người, trao đổi chéo là kiểu tái tổ hợp phổ biến nhất xảy ra trong quá trình giảm phân. Trao đổi chéo là 23
hiện tượng 2 trình tự DNA có độ dài bằng nhau có sự trao đổi lẫn nhau và sinh ra một trình tự tái tổ hợp thứ 3 có cùng độ dài, chứa phần đầu của một trình tự và theo sau bởi phần sau của trình tự còn lại (Hình 1.8a). Chuyển đổi gen là hiện tượng trình tự tái tổ hợp được tạo ra từ phần đầu của một trình tự, theo sau bởi phần giữa của trình tự thứ 2 và theo sau bởi phần cuối của trình tự đầu tiên (Hình 1.8b). Hình 1.8: Hai hiện tượng tái tổ hợp phổ biến của người: (a) trao đổi chéo và (b) chuyển đổi gen. Mặc dù đây là 2 hiện tượng khác nhau, nhưng trong nhiều mô hình, sự kiện chuyển đổi gen có thể được xét đến là 2 lần trao đổi chéo liên tiếp nhau. Tái tổ hợp xảy ra trong quá trình phân bào là sự kiện sinh học quan trọng tạo ra đa dạng hệ gen. Do sự tái tổ hợp trong tất cả các thế hệ trước, bộ gen mà bất kỳ cá thể nào thừa hưởng là sự pha trộn và phản ánh DNA của nhiều cá thể khác nhau qua các thế hệ tổ tiên [21]. Khả năng tạo chuỗi trình tự nhiễm sắc thể mới này cho phép các loài phản ứng nhanh chóng với những thay đổi trong môi trường và loại bỏ các đột biến gây hại. Do đó, tái tổ hợp là một đặc tính thích nghi quan trọng xảy ra ở hầu hết các loài sinh vật nhân chuẩn. Sự tồn tại của bộ gen tổ hợp phong phú như vậy thúc đẩy các nghiên cứu về sự biến đổi gen trong các quần thể để khám phá mối quan hệ giữa nội dung bộ gen và các đặc điểm quan tâm có ảnh hưởng từ yếu tố di truyền. 24
1.2.2. Đồ thị tái tổ hợp di truyền Từ dữ liệu trình tự của một quần thể quan sát, có nhiều câu hỏi liên quan chúng ta muốn biết như: đặc điểm di truyền của quần thể như thế nào? Lịch sử dân số hay nguồn gốc địa lý của quần thể? Hay tìm nguyên nhân cho sự mở rộng hoặc suy giảm dân số của quần thể, mức độ di cư như thế nào? Và quan trọng hơn là tìm mối liên hệ giữa kiểu hình quan sát (ví dụ như bệnh) trong các cá thể thuộc quần thể và dữ liệu trình tự để tìm ra các gen gây bệnh và các cơ chế liên quan. Đồ thị tái tổ hợp di truyền đóng một vai trò quan trọng trong việc trả lời các câu hỏi liên quan đến nghiên cứu di truyền quần thể và các bài toán liên quan đến tìm sự đa dạng trong hệ gen [1]. Khi xây dựng được đồ thị ARG, chúng ta không những xác định được các vùng liên quan đến bệnh quan tâm mà đồ thị còn cho ta cái nhìn tổng quan về đặc điểm của quần thể quan sát, nền tảng của đột biến gây bệnh, và từ đó có thể dự đoán và thay thế dữ liệu bị khuyết (imputing missing data) [53]. Trong nghiên cứu về bệnh dựa trên tập dữ liệu người bệnh và người không bệnh (case-control study), việc xây dựng đồ thị ARG giúp tìm được vị trí nhánh phân biệt rõ nhất giữa người bệnh và người không bệnh, từ đó xác định được vùng gen liên quan đến bệnh [52,71]. Đồ thị ARG còn được ứng dụng hiệu quả trong bài toán tìm SNP, một bài toán quan trọng được tập trung giải quyết trong dự án bản đồ hệ gen người [44]. Trong nghiên cứu di truyền quần thể, đồ thị ARG có ứng dụng trong bài toán xác định các dấu hiệu của chọn lọc tự nhiên [30]; nghiên cứu dòng gen (gene-flow) và sự di trú (migration) liên quan đến tổ tiên của người hiện đại [32]; bài toán phân biệt chuyển đổi gen với tái tổ hợp trao đổi chéo [61]; phát hiện dòng gen trong nấm men [38]; phát hiện đồng tiến hóa (coevolution) trong nấm [10], Tổng hợp nhiều ứng dụng của đồ thị ARG được giới thiệu trong [1,26,30]. Đồ thị ARG được xây dựng xuất phát từ lý thuyết kết hợp (coalescent theory) của Kingman năm 1982 [39]. Kingman đưa ra một cách mô hình quan hệ họ hàng của các chuỗi DNA khi không có sự kiện tái tổ hợp. Các sự kiện kết hợp và đột biến 25
được xem xét và biểu diễn dưới dạng cây. Tuy nhiên, ngoài sự kiện kết hợp và đột biến, sự kiện tái tổ hợp là một thực tế không thể loại bỏ của quá trình tiến hóa và di truyền. Do đó, lý thuyết kết hợp truyền thống đã được mở rộng để tính đến sự kiện tái tổ hợp dưới dạng đồ thị ARG [34]. Khi tái tổ hợp được xét đến trong một mô hình kết hợp, một trình tự được mô hình là có một hoặc hai trình tự cha trong thế hệ trước (là hai nếu tái tổ hợp xảy ra và là một trong trường hợp còn lại), và do đó chúng ta xét đến một đồ thị thay vì là một cây. Mục 1.2.2.1 dưới đây sẽ mô tả giả định được sử dụng để định nghĩa một mạng phát sinh loài là một đồ thị ARG. Từ đó dẫn tới mô tả dữ liệu vào cho thuật toán xây dựng đồ thị ARG trong mục 1.2.2.2 và cấu trúc của một đồ thị ARG trong mục 1.2.2.3. 1.2.2.1. Mô hình các vị trí vô hạn Sự kiện tái tổ hợp và sự kiện đột biến là hai sự kiện quan trọng dẫn tới các biến đổi trên hệ gen từ thế hệ này sang thế hệ khác. Sự kiện đột biến liên quan đến biến đổi trên một vị trí trên chuỗi DNA còn sự kiện tái tổ hợp liên quan đến biến đổi trên các đoạn DNA, dẫn tới tái cấu trúc lại hệ gen làm cho hệ gen của chúng ta là sự pha trộn và kế thừa di truyền từ các thế hệ trước. Trong chiều dài lịch sử tiến hóa, tại một vị trí trên tập các trình tự quan sát, sự kiện đột biến có thể xảy ra một hoặc nhiều lần. Trường hợp đột biến xảy ra tại một vị trí sau đó đột biến ngược lại trạng thái trước đó gọi là đột biến ngược (back mutation); trường hợp đột biến xảy ra tại một vị trí, sau đó lại xuất hiện lại tại vị trí đó một hoặc nhiều lần ở các nhánh tiến hóa (lineage) khác nhau gọi là đột biến lặp lại (recurrent mutation). Trường hợp mạng phát sinh loài biểu diễn sự kiện tái tổ hợp và đột biến, cho phép đột biến ngược hoặc lặp lại được gọi là mạng tái tổ hợp. Quá trình xây dựng đồ thị ARG gắn với giả định có nhiều nhất một sự kiện đột biến xảy ra tại mỗi vị trí trong toàn bộ lịch sử tiến hóa, không cho phép đột biến tái phát. Mô hình đột biến này gọi 26
là mô hình các vị trí vô hạn (infinite-sites model), mô tả sự tiến hóa của các chuỗi DNA rất dài với tỷ lệ đột biến thấp ở mỗi vị trí. Mô hình các vị trí vô hạn xuất phát từ quan sát dữ liệu trình tự thực tế cho thấy, số lượng các vị trí đột biến thường nhỏ so với số lượng vị trí giống hệt nhau. Nghiên cứu đã chỉ ra, tỉ lệ đột biến ở mỗi thế hệ người khoảng 1.2 x10-8 [60], tức là, với độ dài hệ gen người ~3 tỉ cặp base, mỗi thế hệ người chỉ được phép đột biến khoảng 3x109x1.2x10-8 = 36 cặp nucleotit/thế hệ. Vì vậy, khả năng xảy ra đột biến 2 lần tại một vị trí là rất thấp [28]. Đồ thị ARG dưới giả định các vị trí vô hạn còn được gọi là mạng phát sinh loài hoàn hảo có sự kiện tái tổ hợp [69] hay đồ thị ARG hoàn hảo (perfect ARG) [42]. Trên góc độ sinh học, khả năng xảy ra đột biến tái phát là có trong thực tế. Tuy nhiên, từ góc độ toán học và mô hình hóa, một đột biến ngược hoặc lặp lại có thể được mô hình bằng sự kiện trao đổi chéo [26]. 1.2.2.2. Dữ liệu đầu vào cho đồ thị ARG dựa trên mô hình các vị trí vô hạn Dưới giả định mô hình các vị trí vô hạn - chỉ cho phép 1 đột biến xảy ra tại một vị trí trong suốt lịch sử tiến hóa, dữ liệu đầu vào của bài toán xây dựng đồ thị ARG là dữ liệu SNP haplotype, chỉ tính đến các vị trí SNP 2 alen. Do đó, chỉ có nhiều nhất 2 nucleotit khác nhau xuất hiện ở mỗi vị trí trong dữ liệu. Như vậy, dữ liệu đầu vào có thể chuyển đổi thành dạng nhị phân với 2 trạng thái 0 và 1. Dữ liệu gen người được chuyển đổi thành dạng nhị phân bằng việc quy định alen tham chiếu là 0 và alen biến thể là 1 (Hình 1.9). Đối với quần thể người (và các loài lưỡng bội nói chung), 2 trình tự SNP haplotype trong mỗi người được coi là độc lập nhau, tương ứng với 2 trình tự trong dữ liệu đầu vào của bài toán xây dựng đồ thị ARG. 27
Hình 1.9: Biến đổi dữ liệu SNP thành dạng nhị phân. Vị trí có giá trị giống với tham chiếu là 0, giá trị khác tham chiếu là 1. 1.2.2.3. Cấu trúc đồ thị ARG Có 4 thành phần cần thiết để xác định một đồ thị ARG tổng quát cho 1 tập trình tự nhị phân D cho trước: đồ thị cơ sở, các nhãn cạnh, các nhãn nút, và các trình tự quan sát (xem Hình 1.10) [26]. − Đồ thị cơ sở: Cho một tập D gồm n trình tự nhị phân độ dài m, một đồ thị ARG cho D được xây dựng trên một đồ thị có hướng không có chu trình (directed acyclic graph – DAG) chứa chính xác một nút gốc không có cạnh đến, một tập các nút bên trong có cả cạnh đến và cạnh đi và n nút lá có một cạnh đến và không có cạnh đi. Một nút trong có một hoặc hai cạnh đến: nút trong với một cạnh đến gọi là nút cây; nút trong với 2 cạnh đến gọi là nút tái tổ hợp. Cạnh đến một nút tái tổ hợp gọi là cạnh tái tổ hợp; cạnh đến một nút cây gọi là cạnh cây; và cạnh đến nút lá gọi là cạnh lá. Nút gốc và nút trong có thể có một hoặc 2 cạnh đi, đại diện cho quá trình đột biến và sao chép. − Các nhãn cạnh: Mỗi cạnh có thể được gán nhãn bằng một tập các số nguyên từ 1 đến m, biểu thị vị trí trong D nơi một đột biến xảy ra. Với các cạnh không có đột biến xảy ra sẽ không có nhãn cạnh. 28
Hình 1.10: Đồ thị ARG cho tập dữ liệu M gồm 7 trình tự độ dài 5 [26]. Trình tự tổ tiên là “00000”; 5 sự kiện đột biến tại các vị trí tương ứng (1,2,3,4,5) được ghi trên các cạnh xảy ra đột biến của đồ thị; 2 sự kiện tái tổ hợp xảy ra tại vị trí 3 và 4. Dưới giả định có nhiều nhất một đột biến xảy ra tại một vị trí, một vị trí có đột biến sẽ chỉ xuất hiện trên một cạnh duy nhất trong suốt lịch sử tiến hóa. − Các nhãn nút: Mỗi nút trong đồ thị ARG được gán nhãn bởi một trình tự nhị phân độ dài m. Với một nút cây v, trình tự nút cây sv được lấy từ trình tự nút cha của v với sự thay đổi trạng thái từ 0 sang 1 hoặc từ 1 sang 0 tại các vị trí c là các nhãn cạnh hướng vào v. Việc tạo ra trình tự nút cây sv được lấy từ trình tự nút cha của v với sự thay đổi trạng thái từ 0 sang 1 hoặc từ 1 sang 0 tại một vị trí c gọi là sự kiện đột biến. Trong trường hợp sao chép, 2 trình tự con sinh ra sẽ giống hệt trình tự cha của chúng. Xét ngược chiều thời gian, sự kết hợp của 2 trình tự giống nhau thành 1 trình tự giống với 2 trình tự đó gọi là sự kiện kết hợp. 29
Với một nút tái tổ hợp x, đặt s và s’ là 2 trình tự cha độ dài m của nút x. Trình tự nút tái tổ hợp sx là một trình tự độ dài m với trạng thái tại mỗi vị trí c trong sx bằng với trạng thái tại vị trí c trong ít nhất một trong 2 trình tự s hoặc s’. Việc tạo ra trình tự sx từ s và s’ tại một nút tái tổ hợp được gọi là sự kiện tái tổ hợp. Điểm xảy ra sự kiện tái tổ hợp gọi là điểm cắt tái tổ hợp (breakpoint). Một điểm cắt tái tổ hợp liên quan đến một vị trí vật lý trên nhiễm sắc thể (Hình 1.11). Khi ta không biết vị trí chính xác của điểm cắt tái tổ hợp cho sự kiện tái tổ hợp liên quan đến nút tái tổ hợp x nhưng ta biết rằng nó phải xảy ra giữa vị trí c và vị trí d (d > c+1), ta nói rằng điểm cắt tái tổ hợp bx là trong khoảng (c,d], tức là bx lớn hơn c và nhỏ hơn hoặc bằng d và lựa chọn phần trình tự cha đóng góp vào trình tự tái tổ hợp được thay đổi tại vị trí bx. Cụ thể, khi một sự kiện tái tổ hợp xảy ra tại vị trí cắt tái tổ hợp bx, trình tự tái tổ hợp sẽ gồm phần đầu (prefix) từ vị trí 1 đến bx - 1 của trình tự s (s’) và theo sau là phần sau (suffix) từ vị trí bx đến m của trình tự s’ (s). Hình 1.11: Điểm cắt tái tổ hợp. − Trình tự quan sát: tức là các trình tự nhị phân trong tập D, là các trình tự ghi nhãn nút lá trong đồ thị ARG và nó được xác định duy nhất trong đồ thị ARG. Với các loài lưỡng bội, 2 trình tự SNP haplotype tương ứng với 2 nút lá trong đồ thị ARG. Nếu có một sự kiện tái tổ hợp thì sẽ có một chu trình tái tổ hợp (recombination cycle) tương ứng. Ví dụ, bắt đầu ở một nút tái tổ hợp và theo dõi lại dọc theo hai đường dẫn, do tất cả các đường dẫn cuối cùng đều kết hợp để hướng tới 1 nút gốc, hai đường dẫn cuối cùng phải kết hợp lại tại một nút, tại đó hai đường dẫn xác định 30
một chu trình tái tổ hợp. Như vậy, có bao nhiêu sự kiện tái tổ hợp thì sẽ có bấy nhiêu chu trình tái tổ hợp. Đặc điểm này của đồ thị ARG dẫn tới rất nhiều ý tưởng đề xuất xây dựng đồ thị ARG sau này. Một cách biểu diễn khác đồ thị ARG đó là chỉ mô tả các đoạn mang thông tin di truyền từ các trình tự quan sát đến một tổ tiên chung trong quá trình xây dựng đồ thị ARG. Cách biểu diễn này được minh họa như trong Hình 1.12, được lấy từ bài báo của Griffiths và Marjoram [23]. Khi đó, ngoài các đoạn mang thông tin di truyền, các nút trong đồ thị chứa các đoạn không có thông tin (các vị trí mang giá trị không xác định), đó là phần thông tin di truyền từ các tổ tiên khác do sự kiện tái tổ hợp. Sự kiện kết hợp Đột biến tại vị Sự kiện trí thứ 3 tái tổ hợp Hình 1.12: Một ví dụ đồ thị ARG cho 4 trình tự với các ký hiệu: ■: trạng thái di truyền, ◘: trạng thái di truyền đột biến, □: trạng thái không xác định. Xét ngược chiều thời gian, cho trước trình tự tái tổ hợp S và một sự kiện tái tổ hợp phân tách trình tự tái tổ hợp làm 2 trình tự S1 và S2, khi đó ta chỉ biết được phần đầu của S1 và phần sau của S2 mang phần thông tin di truyền từ trình tự tái tổ hợp, các phần còn lại của 2 trình tự này là các phần không có thông tin (phần không xác định). 31
Đồ thị ARG được biểu diễn và xây dựng theo cách tiếp cận như vậy có thể cho phép trình tự đầu vào trong D có một số giá trị khuyết [47,52]. Với một đồ thị ARG đầy đủ được mô tả như trong Hình 1.10 và Hình 1.12, mỗi vị trí c sẽ có một cây thành phần (cây biên - marginal tree) T(c) mô tả lịch sử của các cá thể cho vị trí đó. Từ tập trình tự ban đầu, với mỗi trình tự ta lần theo các cạnh của đồ thị tái tổ hợp di truyền cho vị trí c; khi một sự kiện tái tổ hợp xuất hiện, ta đi theo đường bên trái nếu vị trí tái tổ hợp xảy ra sau c và đi theo đường bên phải trong trường hợp ngược lại. Tập tất cả các cạnh đó sẽ định nghĩa T(c). Hình 1.13 minh họa các cây thành phần cho đồ thị ARG trong Hình 1.12. Bên cạnh các thuật toán xây dựng đồ thị ARG đầy đủ, rất nhiều thuật toán, đặc biệt theo cách tiếp cận thống kê thường xây dựng đồ thị ARG không đầy đủ, tức là đồ thị ARG được biểu diễn bằng tập các cây thành phần và các sự kiện tái tổ hợp. 1.2.3. Bài toán xây dựng đồ thị ARG Quá trình tái tổ hợp làm cho hệ gen của các cá thể trong quần thể bị thay đổi rất nhiều qua các thế hệ. Do đó, với một tập ngẫu nhiên các trình tự của các cá thể trong một quần thể quan sát ở hiện tại, ta không thể xác định được gia phả, lịch sử tiến hóa của chúng. Vì vậy, trong quá trình tái cấu trúc lại tổ tiên của tập trình tự quan sát, tức là xây dựng đồ thị ARG, số sự kiện tái tổ hợp và sự kiện đột biến cũng như vị trí thực sự xảy ra của chúng trong quá trình tiến hóa là không thể xác định được. Các hướng để giải quyết vấn đề này là thiết kế các mô hình gần đúng với quá trình tiến hóa với giả định tối ưu số sự kiện tái tổ hợp và sự kiện đột biến. 32
1. Cây thành phần cho vị trí 1 (2) Cây thành phần cho vị trí 3 (3) Cây thành phần cho vị trí 2 (4) Cây thành phần cho vị trí 4 Hình 1.13: Các cây thành phần (đường đậm nét) của đồ thị ARG trong Hình 1.12. Nguồn hình [43]. Bài toán xây dựng đồ thị ARG được chứng minh là một bài toán NP-khó [69]. Dưới giả định các vị trí vô hạn, bài toán xây dựng đồ thị ARG được phát biểu như sau: Cho một tập D gồm n trình tự nhị phân, mỗi trình tự có độ dài m, tìm một ARG hiển thị D với số sự kiện tái tổ hợp ít nhất. 33
Cụ thể: Dữ liệu vào: Dữ liệu đầu vào là một tập các trình tự nhị phân độ dài m. Các trình tự có độ dài bằng nhau. Tập các trình tự được ký hiệu là D = {S1, , SN}, trong đó N là số lượng trình tự, Sx là một trình tự trong tập D, 1 ≤ x ≤ N. Sx có độ dài m, Sx[i] biểu thị giá trị trạng thái của Sx tại vị trí i, Sx[i] có giá trị bằng 0 hoặc 1, 1 ≤ i ≤ m. Bài toán: Tìm đồ thị ARG mô tả mối quan hệ của các trình tự trong tập dữ liệu vào thông qua 3 sự kiện: đột biến, kết hợp và tái tổ hợp, với giả định chỉ có nhiều nhất một đột biến xảy ra tại mỗi vị trí. Do có nhiều phương pháp khác nhau cho kết quả với độ hợp lý cũng như thời gian thực hiện khác nhau, chúng ta cần đề xuất các phương pháp cho kết quả tốt dựa trên các tiêu chí về số sự kiện tái tổ hợp ít nhất, hình thái cây gần với cây thật nhất, khả thi với dữ liệu lớn hàng trăm đến hàng nghìn trình tự độ dài hệ gen, đồ thị có ứng dụng tốt trong các bài toán thực tế và có thời gian thực hiện khả thi. Dữ liệu đầu ra: Đồ thị ARG chứa các thông tin quan hệ dưới dạng 3 sự kiện cơ bản: đột biến, kết hợp và tái tổ hợp giữa các trình tự đầu vào (nút lá) với các trình tự trung gian được sinh ra trong quá trình xây dựng đồ thị (nút cây) và với một trình tự tổ tiên chung duy nhất (nút gốc). Rất nhiều nghiên cứu xây dựng đồ thị ARG khác nhau đã được đề xuất với các mô hình tái tổ hợp khác nhau phù hợp với quần thể quan sát và mục đích nghiên cứu khác nhau. Trong bài toán xây dựng đồ thị ARG cho các quần thể vi khuẩn, các sự kiện tái tổ hợp được xem xét và mô hình hóa là các sự kiện chuyển đổi gen [17,66]. Trong nghiên cứu di truyền quần thể người, sự kiện tái tổ hợp được mô hình hóa trong quá trình xây dựng đồ thị ARG hầu hết là sự kiện trao đổi chéo. Trong nhiều thuật toán, đặc biệt là các thuật toán tổ hợp tập trung vào đặc điểm cấu trúc của đồ thị, sự kiện chuyển đổi gen có thể được biểu diễn qua 2 sự kiện trao đổi chéo liên tiếp nhau [26]. 34
Luận án tập trung vào các thuật toán tổ hợp xây dựng đồ thị ARG đầy đủ có số sự kiện tái tổ hợp ít nhất dưới giả định mô hình các vị trí vô hạn. Sự kiện tái tổ hợp trong đồ thị ARG được đề cập đến chỉ sự kiện trao đổi chéo và được sử dụng như vậy trong suốt các phần tiếp theo của luận án. Dữ liệu trình tự được xét đến trong bài toán là dữ liệu SNP haplotype được biểu diễn ở dạng nhị phân. 1.3. Các phương pháp xây dựng đồ thị ARG Có 2 hướng nghiên cứu xây dựng đồ thị ARG: (1) Xây dựng đồ thị ARG tối thiểu (minimal ARG), tức là đồ thị có chính xác số sự kiện tái tổ hợp nhỏ nhất, và (2) xây dựng đồ thị ARG “hợp lý” (plausible ARG), tức là các thuật toán không cố gắng xây dựng ARG có chính xác số sự kiện tái tổ hợp ít nhất mà hướng đến việc xây dựng đồ thị ARG với số sự kiện tái tổ hợp được sinh ra phụ thuộc vào các phương pháp mô hình hóa sự kiện tái tổ hợp khác nhau. 1.3.1. Các phương pháp xây dựng đồ thị ARG tối thiểu Các cách tiếp cận theo hướng nghiên cứu này hầu hết đều dựa trên các phương pháp tìm kiếm vét cạn trên đồ thị để cực tiểu hóa số sự kiện tái tổ hợp nhằm đạt tới ARG tối thiểu. Trong đó, khái niệm cặp vị trí không tương thích được sử dụng trong hầu hết các thuật toán để xác định sự kiện tái tổ hợp: Cho một tập D gồm 4 hoặc nhiều hơn 4 trình tự, một cặp vị trí bất kì gọi là không tương thích nếu tồn tại 4 trình tự trong D lần lượt chứa 4 loại giao tử (0,0), (0,1), (1,0), (1,1) cho cặp vị trí đó. Dưới giả định các vị trí vô hạn (có nhiều nhất một đột biến xảy ra tại một vị trí), từ 1 tổ tiên chung, cách duy nhất để có cặp vị trí không tương thích này là ít nhất một sự kiện tái tổ hợp đã xảy ra trong lịch sử giữa 2 vị trí đó (Hình 1.14b). Trường hợp giữa 2 vị trí có ít hơn 4 loại giao tử trên được gọi là cặp vị trí tương thích, khi đó, từ 1 tổ tiên chung, dữ liệu quan sát có thể có được thông qua các sự kiện đột biến (Hình 1.14a). 35
(a) (b) Hình 1.14: (a) Ví dụ cặp vị trí tương thích: cặp vị trí này chỉ chứa 3 loại giao tử và có thể có được từ 1 tổ tiên chung thông qua 2 sự kiện đột biến. (b) Cặp vị trí không tương thích: cặp vị trí chứa 4 loại giao tử và trong trường hợp này phải có ít nhất 1 sự kiện tái tổ hợp xảy ra dưới giả định các vị trí vô hạn (kí hiệu * biểu thị vị trí không có thông tin). Khái niệm cặp vị trí không tương thích này là yếu tố cơ bản dẫn tới rất nhiều thuật toán tìm cận dưới tái tổ hợp và thuật toán xây dựng đồ thị ARG. Các phương pháp vét cạn hướng tới việc tìm ra các điểm cắt tái tổ hợp tối ưu, tức là, số sự kiện tái tổ hợp ít nhất để phá vỡ tất cả các vị trí không tương thích này. Song và cộng sự [62,63] tìm một chuỗi các cây, với mỗi cây cho mỗi vị trí và các sự kiện tái tổ hợp được yêu cầu để dịch chuyển cây tại một vị trí sang cây tại vị trí tiếp theo. Để làm điều này, tác giả xây dựng tất cả các cây có thể cho mỗi vị trí. Sau đó, các sự kiện tái tổ hợp cần thiết để chuyển từ tất cả các cây tại một vị trí sang tất cả các cây tại vị trí tiếp theo được tính toán. Các đồ thị ARG tối thiểu sau đó được xây dựng bằng cách lần theo các vị trí mà có số sự kiện tái tổ hợp ít nhất. Phương pháp này chỉ áp dụng được với tối đa 9 trình. Tuy nhiên, ý tưởng duyệt qua các cây là ý tưởng then chốt được sử dụng cho rất nhiều thuật toán xây dựng đồ thị ARG dựa trên thống kê sau này. Lyngsø và cộng sự [47] sử dụng phương pháp nhánh cận (branch and bound 36
approach) và đưa ra một cải tiến về tốc độ và bộ nhớ sử dụng. Thay vì tính toán từ trái qua phải dọc theo chuỗi trình tự, phương pháp làm việc ngược chiều thời gian, thực hiện các sự kiện đột biến, kết hợp và tái tổ hợp cho đến khi đến một tổ tiên chung tối ưu. Tìm kiếm phân nhánh được thực thi để khám phá tất cả các chuỗi sự kiện có thể, cố gắng tìm một chuỗi sự kiện với một số sự kiện tái tổ hợp cho trước. Nếu không tìm được, số sự kiện tái tổ hợp cho phép được tăng thêm một và cứ như vậy cho đến khi một đồ thị ARG được tìm thấy. Thuật toán sử dụng một số quy tắc để giảm kích thước của dữ liệu như: thu gọn các trình tự (các hàng) giống hệt nhau thành một, thu gọn các cột chứa cùng một giá trị trạng thái, Tuy nhiên, việc xác định cận dưới (lower bound) tái tổ hợp là một việc khó vì nếu chọn nhỏ quá thì quá trình xây dựng đồ thị sẽ phải lặp lại tốn nhiều thời gian. Hơn nữa, việc xây đồ thị với một lượng sự kiện tái tổ hợp cho trước đòi hỏi tốn nhiều bộ nhớ để lưu giữ các trường hợp đã thử và được đánh dấu để tiếp tục tìm các cạnh tiếp theo trong đồ thị. Phương pháp cũng chỉ chạy được với 10 mẫu trình tự. Gusfield và cộng sự [27] đề xuất thuật toán xây dựng một trường hợp đặc biệt của đồ thị ARG nếu có - đồ thị ARG có nút gốc cho trước là trình tự toàn trạng thái 0 với ràng buộc tất cả các chu trình tái tổ hợp không chung nút với nhau (node- disjoint). Khi đó, đồ thị ARG là một cây có nốt sùi (galled-tree) trong đó mọi chu trình tái tổ hợp là các nốt sùi (gall) thỏa mãn không nốt sùi nào chung nút với nốt sùi nào (Hình 1.15). Tác giả đã phát triển thuật toán với thời gian O(nm + n3) để xây dựng cây có số nốt sùi ít nhất (nếu tồn tại), tức là có ít số sự kiện tái tổ hợp nhất cho tập dữ liệu D với nút gốc sr cho trước. Thuật toán xác định các cặp vị trí không tương thích trên tập dữ liệu D' = D  sr, từ đó xác định các thành phần liên thông để xây dựng các gall cho các vị trí không tương thích và kết hợp các gall vào một galled-tree. Sau đó, tác giả đã mở rộng bài toán tìm galled-tree với nút gốc chưa biết [25]. 37
Hình 1.15: Một cây có nốt sùi cho tập trình tự giống với tập trong Hình 1.10 với 2 nốt sùi tương ứng với 2 chu trình tái tổ hợp không chung nút với nhau [27]. Chương trình SHRUB [64] xây dựng thuật toán tính cận trên tái tổ hợp Rub và đồ thị ARG cho tập dữ liệu D sử dụng chính xác Rub sự kiện tái tổ hợp bằng cách xây dựng đồ thị ARG lần lượt từ các nút lá. Các phép biến đổi kết hợp/thay thế các trình tự đầu vào được tiến hành song song tương ứng với các bước xây dựng đồ thị ARG cho đến khi đạt tới 1 nút chung duy nhất (chỉ còn lại một trình tự duy nhất qua các phép biến đổi). Lưu ý rằng trong trường hợp cận trên tái tổ hợp Rub tìm được bằng với cận dưới tái tổ hợp (có thể được tính sử dụng chương trình HapBound trong cùng nghiên cứu) thì SHRUB cho ra đồ thị ARG tối thiểu, trong trường hợp còn lại thì SHRUB sẽ cho ra đồ thị ARG hợp lý. SHRUB mở rộng để xử lý sự kiện chuyển đổi gen với tên gọi SHRUB-GC, được mô tả trong [61]. Trong đó, thuật toán sử dụng 2 sự kiện tái tổ hợp trao đổi chéo để biểu diễn chuyển đổi gen. Wu và cộng sự [71,72] đưa bài toán xây dựng đồ thị ARG về bài toán tìm số trình tự trung gian tối thiểu cần để xây dựng ARG. Thuật toán chạy được với dữ liệu hơn một trăm trình tự độ dài khoảng 100 SNP. Trong khoảng thời gian tiếp theo, nhiều nghiên cứu phát triển tập trung vào hướng xây dựng đồ thị ARG hợp lý dựa trên 38
thống kê. Gần đây, Cámara và cộng sự [7] đã đề xuất một kiểu đồ thị tổng hợp mới gọi là topological ARG dựa trên khái niệm barcode trong persistent homology [9]. Tác dụng chính của topological ARG là khả năng chứa các thông tin trong thời gian đa thức và cho phép ứng dụng được với hàng trăm trình tự. Tuy nhiên, bên cạnh việc không hoàn toàn xây dựng đồ thị ARG đầy đủ, thuật toán cũng mới chỉ áp dụng được với những trình tự ngắn, chưa khả thi với dữ liệu hệ gen người. 1.3.2. Các phương pháp xây dựng đồ thị ARG hợp lý Các phương pháp tìm ARG tối thiểu chỉ áp dụng được cho các bộ dữ liệu nhỏ và độ phức tạp tính toán lớn. Để tương tác được với dữ liệu lớn hơn, các phương pháp xây dựng đồ thị ARG hợp lý đã được đề xuất. Theo hướng nghiên cứu này, các phương pháp xây dựng đồ thị ARG thường theo 2 cách tiếp cận chính là dựa trên heuristic và dựa trên thống kê. Dựa trên ý tưởng từ thuật toán tìm ARG tối thiểu của Lyngsø và cộng sự [47], Minichiello và Durbin đã đề xuất chiến lược mới để xác định sự kiện tái tổ hợp, đó là sự kiện tái tổ hợp được thực hiện trên cặp trình tự có đoạn chung dài nhất [52]. Thuật toán chạy được với tập dữ liệu tối đa một nghìn trình tự có độ dài hàng trăm SNP. Ý tưởng độ dài đoạn chung giữa 2 cá thể cũng được khai thác trong thuật toán xây dựng đồ thị ARG hợp lý của Parida và cộng sự [55]. Một cách tiếp cận khác là lấy mẫu các ARG từ xác suất hậu nghiệm của các mô hình xấp xỉ quá trình kết hợp với tái tổ hợp (coalescent-with-recombination – CwR) [50]. Các thuật toán này cố gắng tích hợp quá trình kết hợp và tái tổ hợp vào các mô hình học máy để xây dựng các cây phả hệ. Trong mô hình này, cạnh của ARG có độ dài biểu thị thời gian trôi đi. CwR thường được mô tả như một quá trình thống kê theo thời gian, nhưng Wiuf và Hein [70] đã chỉ ra rằng nó có thể được thay đổi thành một quá trình toán học tương đương dọc theo trình tự hệ gen. Không giống quá trình theo thời gian, quá trình tuần tự (sequential process) này không phải quá trình Markov, đó là các nòi giống cho 39
các vị trí di truyền 1, , i-1 không thể bị “quên” khi sinh ra nòi giống ở vị trí i+1 dựa trên nòi giống ở vị trí i. Tức là, nòi giống tiếp theo phụ thuộc không chỉ vào nòi giống hiện tại mà còn phụ thuộc vào tất cả các nòi giống trước đó. Tuy nhiên, McVean and Cardin [50] chỉ ra rằng quá trình mô phỏng mà bỏ qua một cách đơn giản các đặc trưng không Markov của quá trình tuần tự cho ra các tập các trình tự mà rất giống với hầu hết các khía cạnh của các trình tự được sinh ra bởi CwR đầy đủ. Tức là, CwR bản chất hầu như là quá trình Markov theo nghĩa là các tương quan tầm xa là yếu và có ảnh hưởng nhỏ lên dữ liệu. Nghiên cứu cũng ước lượng xấp xỉ tới quá trình toàn cục gọi là quá trình kết hợp Markov tuần tự (sequentially Markov coalescent - SMC). SMC đã trở thành một hướng đầy tiềm năng cho phương pháp xấp xỉ suy luận thống kê các ARG [31,45,48]. Chìa khóa bên trong các phương pháp này là nếu không gian trạng thái liên tục cho chuỗi Markov (tức là tập tất cả các phả hệ có thể) được ước lượng bởi một tập hữu hạn, không quá lớn thường bằng việc đếm các cấu trúc liên kết cây và/hoặc thời gian rời rạc, sau đó suy luận có thể được thực hiện một cách hiệu quả sử dụng các thuật toán đã biết cho các mô hình Markov ẩn (HMM). Ví dụ đơn giản nhất của phương pháp này là kết hợp Markov tuần tự từng đôi của Li và Durbin [45]. Rassamusen và cộng sự [58] đã đề xuất thuật toán ARGweaver suy luận đồ thị ARG cho dữ liệu hàng trăm trình tự trên toàn nhiễm sắc thể bằng việc chia nhỏ và chạy song song trên các đoạn 2Mb dữ liệu. Ý tưởng chính của phương pháp là lấy mẫu một ARG của n trình tự với điều kiện trên một ARG của n- 1 trình tự, thao tác này gọi là “xâu chuỗi” (threading). Sử dụng phương pháp dựa trên HMM, người ta lấy mẫu các xâu chuỗi mới một cách hiệu quả từ chính phân phối có điều kiện mà đang quan tâm. Bằng việc lặp lại thao tác loại bỏ và xâu chuỗi lại các trình tự cá thể, thuật toán thu được một bộ lấy mẫu Gibbs hiệu quả cho ARG. Tài liệu kĩ thuật chi tiết hơn về thuật toán và cách dùng ARGweaver được mô tả trong [33]. Phiên bản mở rộng của ARGweaver là ARGweaver-D đã được đề xuất 40
trong [32] để tính đến cả sự kiện phân tách (split) và di trú (migration) trong quần thể. Heine và cộng sự [29] mới đây đã đề xuất một thuật toán MCMC mới trong quá trình suy luận các ARG gọi là bộ lấy mẫu MCMC bắc cầu cây (tree-bridging MCMC sampler). Đối với mỗi đoạn gen được xác định trước, các cây kết hợp tại các vị trí đầu và vị trí cuối được cố định trong khi một chuỗi cây mới (cây bắc cầu) được đề xuất tại các vị trí ở giữa tương thích với dữ liệu D và các cây đầu cuối. Thủ tục bắc cầu được thực hiện trên các đoạn gen như vậy giúp chia nhỏ bài toán suy luận trên toàn bộ dữ liệu D, giảm độ phức tạp tính toán và phù hợp với các hệ thống máy tính song song. Các phương pháp theo cách tiếp cận thống kê là một hướng tiếp cận được nhiều nhà nghiên cứu phát triển gần đây. Tuy nhiên, các phương pháp này không suy luận được các ARG đầy đủ mà chỉ là tập các cây biên với tập các sự kiện tái tổ hợp tương ứng. Các phương pháp này thường được dùng trong việc mô phỏng dữ liệu. Hơn nữa, cách tiếp cận này rất phức tạp, đòi hỏi chi phí tính toán lớn. Bên cạnh đó, việc thiết lập các tham số phù hợp cho các bộ dữ liệu khác nhau cũng là một rào cản lớn cho người dùng nên cách tiếp cận này vẫn chưa có được những ứng dụng thực tế trên những tập dữ liệu lớn. 1.3.3. Tổng hợp các phần mềm xây dựng đồ thị ARG Dưới đây là bảng tổng hợp các phần mềm xây dựng đồ thị ARG tiêu biểu được công khai cho cộng đồng nghiên cứu. Bảng 1.1: Các phần mềm xây dựng đồ thị ARG tiêu biểu. Tên phần mềm Địa chỉ truy cập Mô tả Các công cụ xây dựng đồ thị ARG tối thiểu Galledtree [27] - Mã nguồn mở. ~gusfield/galledtree.tar - Ngôn ngữ lập trình: Python - Yêu cầu: Python phiên bản 2. 41
SHRUB [64] - Chỉ cung cấp chương trình chạy. du/~yss/lu.html - Hỗ trợ chạy trên hệ điều hành MS Windows, Linux và Mac OS X. SHRUB-GC [61] - Phiên bản mở rộng của SHRUB du/~yss/lugc.html khi tính đến cả sự kiện chuyển đổi gen. - Chỉ cung cấp chương trình chạy. - Hỗ trợ chạy trên hệ điều hành Linux và Mac OS X. TARGet [7] - Mã nguồn mở b/TARGet - Ngôn ngữ lập trình: Python - Yêu cầu: Python phiên bản 2.7. Các công cụ xây dựng đồ thị ARG hợp lý Margarita [52] - Mã nguồn mở. urces/software/margarita/ - Ngôn ngữ lập trình: Java. ARGweaver [58] - Mã nguồn mở. argweaver - Yêu cầu cài đặt: Trình biên dịch C++ và Python. Arbores [29] - Mã nguồn mở. Arbores - Ngôn ngữ lập trình: C Bacter [66] - Chương trình xây dựng đồ thị cter/ ARG cho dữ liệu gen của vi khuẩn. - Mã nguồn mở. - Ngôn ngữ lập trình: Java. 1.4. Ứng dụng ARG trong nghiên cứu tương quan toàn hệ gen Trong nghiên cứu tương quan người bệnh - người không bệnh (case-control association study), tần số alen tại các vị trí quan tâm được so sánh trong các quần thể gồm các cá thể bị bệnh và các cá thể không bị bệnh. Tần số trong người bị bệnh mà cao hơn là minh chứng cho alen đó liên quan đến nguy cơ gây bệnh tăng lên. Bằng việc phân tích trên dữ liệu SNP haplotype, phân tích sự phân biệt của các alen SNP giữa các quần thể người bệnh và người không bệnh ta có thể xác định được các vị trí có liên quan một cách thống kê tới bệnh. 42
Có rất nhiều phương pháp khác nhau cho bài toán nghiên cứu vùng gen liên quan đến bệnh quan tâm như phương pháp kiểm thử sự liên đới áp dụng cho từng vị trí [16,56]; phương pháp dựa trên sự kết hợp thông qua bảng phả hệ [73]; phương pháp phân cụm haplotype [18]. Trong phần này, chúng tôi giới thiệu cách tiếp cận sử dụng bảng phả hệ dưới dạng đồ thị ARG cho bài toán tìm vùng gen liên quan đến bệnh quan tâm. Theo đúng lịch sử, đồ thị ARG đúng tạo ra các trình tự quan sát sẽ cho thấy các vị trí (cả về vị trí trong bộ gen và thời điểm xảy ra) của tất cả các đột biến và tái tổ hợp, và do đó, nó cho thấy rõ ràng các đoạn có trong người bệnh nhưng không có trong người khỏe mạnh. Do đó, xây dựng được ARG đúng theo lịch sử di truyền góp phần quan trọng trong các nghiên cứu tương quan toàn hệ gen. Di chuyển dọc theo nhiễm sắc thể, các hình thái của các cây biên liên tiếp nhau dịch chuyển theo tác động của các sự kiện tái tổ hợp mang tính lịch sử. Các sự kiện tái tổ hợp định nghĩa vùng nhiễm sắc thể mà mỗi cây biên mở rộng. Với một vị trí cho trước, cây biên có thể được trích xuất từ ARG bằng cách lần vết phả hệ của vị trí đó ngược chiều thời gian từ các nút lá. Khi một tái tổ hợp xuất hiện, phả hệ sẽ theo đường của thế hệ cha mẹ phía bên trái nếu điểm cắt tái tổ hợp nằm phía phải của vị trí được xét, và ngược lại thì đi theo thế hệ cha mẹ phía bên phải. Nếu có một đột biến nguy cơ gây bệnh tại một vị trí cụ thể trên nhiễm sắc thể (giả sử đột biến chỉ xuất hiện nhiều nhất một lần tại mỗi vị trí trong suốt lịch sử tiến hóa), nó sẽ xảy ra trên một số nhánh bên trong cây biên tại vị trí đó. Vì vậy, một cách để tìm các tương quan đến bệnh là kiểm tra các cây biên để tìm những cây có nhánh phân biệt rõ nhất giữa người bệnh và người không bệnh, tức là các nhánh mà ở đấy nhiều người bệnh và chỉ số rất ít hoặc không có người không bệnh thuộc nhánh đó. Cụm các người bệnh tập trung vào một nhánh như vậy sẽ gợi ý rằng có một đột biến gây bệnh xuất hiện trên nhánh đó (Hình 1.16). 43
Hình 1.16: (a) Đồ thị ARG cho tập 4 trình tự, trong đó trình tự s1, s2 là từ 2 cá thể khỏe mạnh, trình từ s3, s4 là từ 2 cá thể bị bệnh. (b) Đột biến 3 (vùng khoanh tròn) trên cây biên tại vị trí 3 của đồ thị ARG trong (a) cho ra sự phân biệt rõ nhất giữa các trình tự bệnh và trình tự không bệnh. Một lưu ý rằng cây biên T có thể được xác định tại một vị trí SNP hoặc là cho tất cả các vị trí trong một khoảng IT của các vị trí SNP liên tiếp nhau. Và điểm bắt đầu và kết thúc chính xác của khoảng vị trí vật lý mà một đột biến xuất hiện thì ta không thể biết vì ta chỉ biết cây biên tại vị trí SNP và điểm cắt tái tổ hợp có thể xuất hiện giữa 2 vị trí SNP liên tiếp nhau. Cách làm ánh xạ tương quan sử dụng đồ thị ARG được tóm tắt như sau: 1. Với một tập D các haplotype từ các cá thể người bệnh và người không bệnh, xây dựng đồ thị ARG G cho D sử dụng thuật toán xây dựng đồ thị ARG. 2. Tìm tập các cây biên T của G tại các vị trí của D. 3. Các cạnh e trong cây T T được tính điểm độ tốt trong việc phân biệt các lá 44
gán nhãn bệnh với các lá gán nhãn không bệnh. Sau đó, cạnh e với điểm cao nhất được thiết lập cho T. 4. Đặt T là cây biên trong T có độ tương quan lớn nhất, tức là chứa cạnh e với độ phân biệt lớn nhất giữa nhãn bệnh và nhãn không bệnh trong tất cả các cạnh trong tất cả các cây biên trong T. Nếu T đủ tốt (đạt một ngưỡng cho trước), kết luận là đột biến gây bệnh có khả năng xảy ra quanh vị trí cây biên T và đột biến xuất hiện có thể xảy ra trong thời gian được biểu thị bởi cạnh e được tìm thấy trong T. 1.5. Kết luận chương Đồ thị tái tổ hợp di truyền là một loại mạng phân loài mô hình hóa quan hệ tiến hóa của một tập trình tự quan sát trong một quần thể thông qua 3 sự kiện: kết hợp, tái tổ hợp và đột biến. Trong đó, sự kiện tái tổ hợp là sự kiện sinh học quan trọng tạo ra đa dạng hệ gen và là trọng tâm nghiên cứu trong đồ thị ARG. Xây dựng được đồ thị ARG đầy đủ giúp ta có cái nhìn tổng quan về quần thể quan sát, các đặc điểm và quan hệ di truyền giữa các cá thể trong quần thể, từ đó có thể trả lời các câu hỏi liên quan đến lịch sử tiến hóa, dấu hiệu chọn lọc tự nhiên, qua các thế hệ đến quần thể hiện tại. Đặc biệt, đồ thị ARG là công cụ hỗ trợ đắc lực trong các nghiên cứu tương quan (GWAS) nhằm xác định các vùng gen liên quan đến bệnh. Trong quá trình tái cấu trúc lại tổ tiên của tập trình tự quan sát, tức là xây dựng đồ thị ARG, số sự kiện tái tổ hợp và sự kiện đột biến cũng như vị trí thực sự xảy ra của chúng trong quá trình tiến hóa là không thể xác định được. Quan nghiên cứu dữ liệu thực tế cho thấy, khả năng xảy ra đột biến 2 lần tại một vị trí là rất thấp. Vì vậy, các phương pháp xây dựng đồ thị ARG hướng tới hàm mục tiêu tối ưu số sự kiện tái tổ hợp với giả định chỉ có nhiều nhất một sự kiện đột biến có thể xảy ra tại mỗi vị trí trong suốt lịch sử tiến hóa. 45
Có hai hướng nghiên cứu là xây dựng đồ thị ARG tối thiểu, tức là, đồ thị có chính xác số sự kiện tái tổ hợp ít nhất, và xây dựng đồ thị ARG hợp lý, có số sự kiện tái tổ hợp phụ thuộc vào cách mô hình hóa sự kiện tái tổ hợp của các thuật toán. Trong đó, hướng nghiên cứu xây dựng đồ thị ARG tối thiểu chỉ giới hạn ở các tập dữ liệu nhỏ. Hướng nghiên cứu xây dựng đồ thị ARG hợp lý có thể áp dụng được với các tập dữ liệu lớn hơn nhưng vẫn chưa có thuật toán nào thực sự hiệu quả với các bộ dữ liệu hàng nghìn hệ gen người. Ngày nay, những thành tựu trong công nghệ giải trình tự gen thế hệ mới, sự phát triển và ngày càng hoàn thiện của các thư viện đặc tả biến dị di truyền trong quần thể người đã tạo tiền đề cho các nghiên cứu trên toàn hệ gen. Để có thể ứng dụng được vào các nghiên cứu về biến thể di truyền liên quan đến bệnh ở người một cách hiệu quả, các phương pháp phải có khả năng tính toán được trên dữ liệu liên quan đến hàng nghìn hệ gen. Do đó, luận án hướng tới mục tiêu phát triển thuật toán xây dựng đồ thị ARG có khả năng chạy được với hàng nghìn hệ gen người. 46
Chương 2. THUẬT TOÁN ARG4WG XÂY DỰNG ĐỒ THỊ TÁI TỔ HỢP DI TRUYỀN HỢP LÝ CHO DỮ LIỆU HỆ GEN 2.1. Giới thiệu Khảo sát trên các phương pháp tìm ARG hợp lý cho thấy cách tiếp cận dựa trên heuristic của Minichiello và Durbin khả thi với tập dữ liệu một nghìn trình tự có độ dài hàng trăm SNP. Đồng thời, thuật toán Margarita của Minichiello và Durbin đã có những ứng dụng vào một số bài toán thực tế [44,52]. Do đó, luận án tập trung nghiên cứu theo hướng tiếp cận này, phân tích và tìm ra nguyên nhân dẫn tới hạn chế cho dữ liệu lớn hơn của thuật toán, từ đó đưa ra thuật toán đề xuất. Phần 2.1.1 sẽ giới thiệu về các định nghĩa được sử dụng trong thuật toán Margarita và cũng được sử dụng trong thuật toán đề xuất. Phần 2.1.2 trình bày về thuật toán Margarita đồng thời đưa ra hạn chế của Margarita trong việc xử lý dữ liệu lớn. 2.1.1. Các định nghĩa Đặt “*” là trạng thái không xác định (vị trí không được thừa kế giá trị từ bất kì trình tự nào) trong các nút trong của đồ thị ARG. Định nghĩa S1[i] khớp với S2[i] nếu hoặc cả 2 trình tự có cùng trạng thái hoặc trạng thái của ít nhất một trong 2 trình tự là *, tức là, S1[i] khớp với S2[i] khi và chỉ khi S1[i] = S2[i] hoặc S1[i] = * hoặc S2[i] = *. Một toán tử bù, ¬, được định nghĩa để nếu S[i] = 0 thì ¬S[i] = 1 và ngược lại, và * là phần bù của chính nó. Một trình tự S1 được coi là dài hơn trình tự S2 (L(S1) > L(S2)) nếu S1 chứa nhiều vật liệu di truyền hơn S2. Ta cũng định nghĩa (L(S1) > L(S2))[a,b] nếu S1 dài hơn S2 trong đoạn [a,b]. 47
2.1.2. Thuật toán Margarita xây dựng đồ thị ARG Ý tưởng then chốt của thuật toán Margarita là định nghĩa đoạn chung giữa các cặp trình tự và sử dụng đoạn chung dài nhất cho bước tái tổ hợp trong quá trình xây dựng đồ thị ARG. Thuật toán định nghĩa có một đoạn chung giữa trình tự S1 và S2 trên một đoạn liên tục từ vị trí a đến vị trí b nếu: 1. S1[i] khớp với S2[i] với mọi a ≤ i ≤ b; 2. Tồn tại ít nhất 1 vị trí i mà S1[i] = S2[i] ≠ *; 3. Nếu a > 1 thì S1[a-1] ≠ S2[a-1]; 4. Nếu b < m thì S1[b+1] ≠ S2[b+1]. Điều kiện thứ (1) yêu cầu 2 trình tự có cùng trạng thái hoặc trạng thái của ít nhất một trong 2 trình tự là * trên đoạn chung; điều kiện thứ (2) yêu cầu rằng có ít nhất một vị trí trong đoạn chung mà cả 2 trình tự đều có giá trị xác định; điều kiện thứ (3) và (4) yêu cầu rằng đoạn chung đó là cực đại. Để xây dựng đồ thị ARG, thuật toán Margarita tính ngược chiều thời gian, bắt đầu từ một tập các trình tự quan sát cho đến khi đến được một tổ tiên chung duy nhất. Lưu đồ thuật toán được biểu diễn như Hình 2.1. Các bước của thuật toán được diễn giải như sau: Bước 1: Margarita tìm các cặp trình tự S1 và S2 thỏa mãn S1[i] khớp với S2[i] với mọi 1 ≤ i ≤ m, chọn ngẫu nhiên ra một cặp trong số đó và thực hiện kết hợp 2 trình tự đó thành một trình tự cha. Quá trình này được lặp lại đến khi không thể thực hiện sự kiện kết hợp nào được nữa, thuật toán chuyển sang Bước 2. Bước 2: Các sự kiện đột biến sẽ được xét đến. Thuật toán tìm các vị trí đột biến, là các vị trí tại đó tồn tại duy nhất một trình tự có giá trị bằng 1 trong khi các trình tự khác có giá trị bằng 0 hoặc ngược lại. Thuật toán sẽ thay đổi giá trị tại các vị trí đột 48
biến này, và quay lại bước 1. Nếu không có sự kiện đột biến thì thuật toán chuyển sang Bước 3. Hình 2.1: Lưu đồ thuật toán Margarita. Bước 3: Thực hiện tái tổ hợp với chiến lược được mô tả như bên dưới. Sau đó quay lại bước 1. Nếu không thực hiện được tái tổ hợp thì thuật toán kết thúc và đạt đến một tổ tiên chung duy nhất SMCA. Trong bước tái tổ hợp, Margarita thực hiện tái tổ hợp trên cặp trình tự chứa đoạn chung liên tục dài nhất (đoạn chung dài nhất - longest shared tract) với xác suất 0.9 và thực hiện tái tổ hợp trên một cặp trình tự chứa đoạn chung bất kì với xác suất 0.1. Thuật toán thực hiện tái tổ hợp trên một trong 2 trình tự, phân tách trình tự đó ra để có được trình tự con chứa đoạn chung để thực hiện kết hợp với trình tự còn lại. Bài toán tìm xâu con chung dài nhất giữa 2 xâu cùng độ dài m đã được chứng minh có độ phức tạp O(m2) nếu sử dụng quy hoạch động và có độ phức tạp O(2m) nếu sử 49
dụng cây hậu tố tổng quát (generalized suffix tree) [24]. Trong bài báo, Margarita không công bố cụ thể thuật toán tìm đoạn chung dài nhất. Vì vậy, Margarita sẽ cần duyệt qua 2m vị trí để tìm đoạn con chung dài nhất giữa 2 trình tự độ dài m. Đáng chú ý, nếu đoạn chung được tìm thấy nằm bên trong trình tự, Margarita sẽ phải thực hiện 2 sự kiện tái tổ hợp, sinh ra 3 trình tự con từ 1 trình tự để có được trình tự chỉ chứa đoạn chung để thực hiện kết hợp với trình tự còn lại (Hình 2.2). Chiến lược này gây ra sự bùng nổ về số nút trong quá trình xây dựng đồ thị khi đoạn chung dài nhất luôn được tìm thấy bên trong trình tự. Điều này dẫn đến thuật toán không thể chạy được với số lượng lớn trình tự trên cả nhiễm sắc thể. Hình 2.2: Vấn đề trong việc thực hiện sự kiện tái tổ hợp của Margarita. Hai trình tự S1 và S2 với đoạn chung dài nhất giữa hai trình tự được biểu diễn bằng đoạn màu đen. Thuật toán thực hiện lần lượt 2 sự kiện tái tổ hợp R1 và R2 trên trình tự S1 để sinh ra 3 trình tự con S11, S12 và S13. Sau đó, trình tự con chứa đoạn chung dài nhất S13 sẽ được kết hợp với S2. Vì vậy, khi đoạn chung dài nhất được tìm thấy bên trong trình tự, thuật toán phải thực hiện 2 sự kiện tái tổ hợp trên một trình tự và từ 2 trình tự ban đầu (S1 và S2) sẽ thành 3 trình tự ở thế hệ tiếp theo (S11, S12 và S' (S' = S2)). 50
2.2. Thuật toán ARG4WG 2.2.1. Chiến lược tìm đoạn đầu chung dài nhất Cho trước một tập trình tự D và một trình tự s, ta sẽ chứng minh mệnh đề rằng số cực tiểu sự kiện tái tổ hợp để tách s thành các trình tự con mà kết hợp được với các trình tự trong D có thể đạt được bằng cách lặp lại việc lấy đoạn chung dài nhất từ một phía của s. Mệnh đề 1 dưới đây được phát biểu và chứng minh với trường hợp lặp lại việc lấy đoạn chung dài nhất từ phía trái của s. Ta có thể chứng minh tương tự với phía bên phải. Phần tiếp theo sẽ đưa ra một ví dụ chỉ ra rằng chiến lược lấy đoạn chung dài nhất bên trong trình tự không phải luôn cho ta số sự kiện tái tổ hợp ít nhất. Mệnh đề 1: Cho một tập các trình tự trong D, và 1 trình tự s có cùng độ dài m. Số cực tiểu sự kiện tái tổ hợp, f (s, D) , để tách s thành các trình tự con mà có thể kết hợp với các trình tự trong D có thể đạt được bằng cách lặp lại việc lấy các đoạn chung dài nhất từ phía trái của s. Chứng minh Mệnh đề 1: f (s, D) = 0 nếu s có thể kết hợp với một trình tự trong D. f (s, D) = 1 nếu s có thể được biểu diễn thành (s1,s2 ) trong đó s1 và s2 có thể kết hợp với các trình tự trong D. Ngoài ra, đặt s = s1 + s2 + s3 trong đó s1 và s2 có thể kết hợp với các trình tự trong D. f (s, D) có thể được viết dưới dạng: f (s, D) = min { f (s + s + s , D)} = s1 ,s2 ,s3 1 2 3 (1) min { f (s , D) + 2} s1 ,s2 ,s3 3 51
Đặt s = sl + sr trong đó sl là đoạn dài nhất phía bên trái của s mà có thể kết hợp với một trình tự trong D. Như vậy, tất cả các trình tự con từ phía bên trái của s mà có thể kết hợp với 1 trình tự trong D là một tập con của sl . Chúng ta sẽ chứng minh rằng: f (s, D) = min { f (s , D) + 2} = f (s , D) +1 (s1,s2 ,s3 ) 3 r (s*,s* ,s* ) s s Đặt 1 2 3 là giải pháp tối ưu của phương trình (1) và x là phần bù của l * * * * trong s1 + s2 , sx = s1 + s2 − sl . Hình 2.3: Tất cả các trình tự con từ phía bên trái của s mà có thể kết hợp với một trình tự trong D là một tập con của đoạn bên trái dài nhất của s ( sl ). * * Rõ ràng là sx là một chuỗi con của s2 vì s1 là một chuỗi con của sl . Do đó, nếu * s2 có thể kết hợp với một trình tự h trong D thì có thể kết hợp với h (xem Hình 2.3). Điều này dẫn tới: 52
* f (s, D) = f (sl + sx + s3 , D) = f (s3 *, D) + 2 = f (sr , D) +1 Nói cách khác, chúng ta có thể có được cực tiểu số sự kiện tái tổ hợp để tách s thành các trình tự con mà kết hợp được với các trình tự trong D bằng cách lặp lại việc lấy ra các đoạn chung dài nhất từ phía bên trái của s. Chứng minh tương tự với trường hợp lấy từ phía bên phải. Tuy nhiên, điều này là không đúng nếu chúng ta không chọn các đoạn chung dài nhất từ hai phía của s. Hình 2.4 mô tả giải pháp tối ưu mà chỉ cần một sự kiện tái tổ hợp (xem Kịch bản A). Tuy nhiên, nếu ta chọn đoạn chung dài nhất không phải từ 2 phía của s (ở đây là chọn đoạn chung dài nhất trong s) thì ta phải cần đến 2 sự kiện tái tổ hợp (Kịch bản B). Hình 2.4: Phân tách s bằng cách chọn các đoạn chung dài nhất trong s để kết hợp với các trình tự trong D có thể không dẫn tới số cực tiểu sự kiện tái tổ hợp. Từ mệnh đề trên, luận án định nghĩa đoạn đầu chung dài nhất (longest shared end) là đoạn chứa thông tin di truyền giống nhau liên tục dài nhất tính từ 2 đầu của các 53
trình tự. Từ đó, luận án đề xuất thuật toán ARG4WG sử dụng chiến lược tìm cặp trình tự có đoạn đầu chung dài nhất cho bước tái tổ hợp. Thuật toán cần duyệt (lR+lL) vị trí (với lR là độ dài đoạn đầu chung dài nhất từ phía bên phải, lL là độ dài đoạn đầu chung dài nhất từ phía bên trái, 0 ≤ lL,lR ≤ m) để tìm ra đoạn đầu chung dài nhất giữa 1 cặp trình tự. Chiến lược thực hiện tái tổ hợp sử dụng đoạn đầu chung dài nhất này giúp thuật toán ARG4WG luôn đảm bảo số nút được ổn định sau mỗi bước tái tổ hợp (minh họa trong Hình 2.5) và thuật toán ARG4WG có thể chạy được với dữ liệu lớn hàng nghìn trình tự độ dài hệ gen trong thời gian hợp lý. 2.2.2. Thuật toán ARG4WG Tương tự thuật toán Margarita, ARG4WG được xây dựng ngược chiều thời gian, từ một tập các trình tự cho tới khi đạt tới một tổ tiên chung duy nhất. ARG4WG gồm 3 bước chính: bước kết hợp, bước đột biến và bước tái tổ hợp. Lưu đồ thuật toán giống với Margarita như trong Hình 2.1. Bước kết hợp và bước đột biến được thực hiện tương tự như trong thuật toán Margarita. Trong bước tái tổ hợp, thuật toán ARG4WG định nghĩa đoạn đầu chung dài nhất là đoạn chứa chung phần vật liệu di truyền liên tục nhiều nhất tính từ phía bên trái hoặc phía bên phải của 2 trình tự. Thuật toán sẽ chọn một cặp trình tự có đoạn đầu chung dài nhất để thực hiện bước tái tổ hợp. Trong 2 trình tự được chọn, thuật toán thực hiện một sự kiện tái tổ hợp bằng việc tách một trình tự thành 2 trình tự con mới, trình tự con chứa đoạn chung sẽ được kết hợp với trình tự còn lại ngay sau đó (xem Hình 2.5). Như vậy, sự kiện tái tổ hợp được thực hiện theo chiến lược này chỉ thay thế một trình tự bằng một trình tự mới có ít vật liệu di truyền hơn. Điều này luôn đảm bảo số nút trong đồ thị được ổn định (luôn có xu hướng giảm) trong suốt quá trình xây dựng đồ thị ARG. Chi tiết thuật toán ARG4WG được trình bày trong phần dưới đây. 54
Hình 2.5: Sự kiện tái tổ hợp được biểu thị trong thuật toán ARG4WG. (a) Xét 2 trình tự S1 và S2, các đoạn đầu chung của 2 trình tự từ phía bên trái (hình lượn sóng) và từ phía bên phải (màu đen) được xác định. (b) Với 1 tập 3 trình tự S1, S2 và S3, các đoạn đầu chung của mỗi cặp được tính toán (hình lượn sóng) và đoạn đầu chung dài nhất được xác định được mô tả bằng đoạn màu đen giữa trình tự S1 và S2. (c) Một sự kiện tái tổ hợp được thực hiện trên trình tự S1 để sinh ra 2 trình tự con S11 và S12. S12 chứa đoạn đầu chung dài nhất sau đó sẽ được kết hợp với S2. Như vậy, ARG4WG luôn chỉ thực hiện 1 tái tổ hợp trên 1 trình tự và từ 2 trình tự ban đầu (S1, S2) sẽ thành 2 trình tự ở thế hệ tiếp theo (S11, S’), trong đó S’ = S2 và S11 có ít vật liệu di truyền hơn S1. 55
Với một cặp (S1, S2), đặt (S1, S2){d} là đoạn đầu chung của chúng. Cụ thể, (S1, S2){d=left} là phần chung của đầu bên trái của (S1, S2); (S1, S2){d=right} là phần chung của đầu bên phải của (S1, S2). Định nghĩa cặp có đoạn đầu chung cực đại (S1,S2){d,l} với độ dài đoạn chung l (0 < l ≤ m) của S1 và S2 nếu nó thỏa mãn các điều kiện sau: 1. Nếu d = left thì S1[i] khớp với S2[i] với mọi 1 i l và hoặc l = m hoặc S1[l+1] không giống S2[l+1]. 2. Nếu d = right thì S1[i] khớp với S2[i] với mọi m-l i m và hoặc l = m hoặc S1[m-l] không giống S2[m-l]. 3. Vùng giống nhau phải có ít nhất một vị trí i mà S1[i] = S2[i] *. Điều kiện đầu tiên và thứ 2 xác định rằng đoạn đầu chung liên tục của 2 trình tự là cực đại. Điều kiện thứ 3 nhấn mạnh rằng đoạn đầu chung giữa 2 trình tự có chung ít nhất một vị trí mang giá trị xác định. Điều này làm giảm số các nhánh dư thừa trong quá trình xây dựng ARG [50,52]. Xét cặp trình tự (S1,S2) có các đoạn đầu chung cực đại từ phía bên trái và từ phía bên phải tương ứng là lL và lR. Nếu đoạn đầu chung cực đại từ phía bên phải chứa nhiều phần vật liệu di truyền chung hơn đoạn đầu chung cực đại từ phía bên trái thì lR được xác định là đoạn đầu chung dài nhất của cặp trình tự (S1,S2). Xét ví dụ như hình trên đây, cặp trình tự (S1, S2) có lR = 4, lL = 5, tuy nhiên lR chứa nhiều phần vật liệu di truyền chung hơn (3 vị trí có giá trị xác định) so với lL (2 vị trí có giá trị xác định) nên cặp trình tự (S1, S2) được xác định có đoạn đầu chung dài nhất là từ phía bên phải. 56
Với 1 cặp trình tự được chọn cho bước tái tổ hợp, trình tự chứa ít vật liệu di truyền trong phần đoạn đầu chung hơn sẽ được chọn để thực hiện tái tổ hợp. Điều kiện này nhằm đảm bảo sự kiện tái tổ hợp chỉ thay thế một trình tự bằng một trình tự con của trình tự ở thế hệ trước đó, không phát sinh trình tự khác, giúp đảm bảo quá trình tính toán ổn định của thuật toán. Xét cặp trình tự (S1, S2) như trên, đoạn đầu chung dài nhất được xác định là từ phía bên phải lr = 4. Trong đó, S2 có ít vật liệu di truyền trong đoạn chung hơn (3 giá trị xác định) so với trình tự S1 (4 giá trị xác định). Do đó, S2 sẽ được chọn để thực hiện tái tổ hợp. Sau bước tái tổ hợp, S2 được phân tách thành 2 trình tự S21 = 100000001 và S22 = 011*. S22 chứa phần đoạn chung sẽ được kết hợp với S1 và sinh ra trình tự S' = S1. Như vậy, sau bước tái tổ hợp và thực hiện 1 kết hợp, từ 2 trình tự S1 và S2 sẽ thu được 2 trình tự S' = S1 và S21 chứa ít vật liệu di truyền hơn S2. Trong khi nếu ta chọn thực hiện tái tổ hợp trên S1 thì ta sẽ được 2 trình tự con S11 = 001000 và S12 = 0110. S12 chứa phần đoạn chung sẽ được kết hợp với S2 và sinh ra một trình tự mới S' = 1000000010110. Như vậy, sau bước tái tổ hợp và thực hiện 1 kết hợp sẽ sinh ra trình tự S11 chứa ít vật liệu di truyền hơn S1 và một trình tự mới S' là sự kết hợp của trình tự S2 và trình tự S12. Điều này làm tăng độ phức tạp tính toán của thuật toán. Thuật toán bắt đầu từ thời gian t = 1. Tập các trình tự tại thời gian t được ký hiệu là Dt (D1 = D). Với mỗi Dt, 3 danh sách ứng cử viên cho các sự kiện kết hợp, đột biến và tái tổ hợp được xây dựng như sau: • Danh sách kết hợp C: Với một cặp có đoạn đầu chung dài nhất (S1,S2){d,l}, nếu l = m, ta cho cặp này vào danh sách kết hợp. • Danh sách đột biến M: Với một vị trí i (1 ≤ i ≤ m), nếu tồn tại duy nhất một trình tự S1, và với mọi trình tự S2 trong Dt╲{S1} ta có S2[i] = ¬S1[i] thì S1[i] được cho vào danh sách đột biến. • Danh sách tái tổ hợp R: Với một cặp có đoạn đầu chung dài nhất (S1,S2){d,l}, nếu 0 < l < m, (S1,S2){d,l} được cho vào danh sách tái tổ hợp. 57
Khi một trong ba sự kiện xuất hiện, tập trình tự tiếp theo Dt+1 được tạo ra từ tập trình tự hiện thời Dt và 3 danh sách ứng cử viên được cập nhật. • Nếu một sự kiện kết hợp xảy ra giữa 2 trình tự S1 và S2 thì 2 trình tự S1 và S2 được hợp lại thành một tổ tiên chung S’. Loại bỏ S1 và S2 từ tập trình tự Dt; và thêm trình tự S’ vào Dt+1. • Nếu một sự kiện đột biến xảy ra trên trình tự S thì một trình tự mới S’ được tạo ra từ trình tự S với điểm đột biến. Thay thế trình tự S bằng trình tự mới S’ trong Dt+1. • Nếu một sự kiện tái tổ hợp xảy ra trên trình tự S1, trình tự S1 sẽ được tách ra thành 2 trình tự con mới S11, S12. Thay thế S1 bằng 2 trình tự con mới đó trong Dt+1. Thuật toán được mô tả đầy đủ trong Thuật toán 2.1 với các quy tắc sau: • Quy tắc 1: Thứ tự ưu tiên thực hiện các sự kiện là kết hợp, đột biến và cuối cùng là tái tổ hợp. Các thao tác giúp giảm số trình tự có thứ tự ưu tiên lớn hơn nhằm giúp cho quá trình xây dựng đồ thị nhanh hội tụ. 58
• Quy tắc 2: Nếu có nhiều hơn một ứng cử viên trong danh sách kết hợp, lấy ngẫu nhiên một ứng cử viên trong danh sách đó để thực hiện kết hợp. • Quy tắc 3: Nếu có nhiều hơn 1 ứng cử viên trong danh sách đột biến, một ứng cử viên được lấy ngẫu nhiên để thực hiện đột biến. Lưu ý rằng thuật toán chỉ thực hiện sự kiện đột biến duy nhất một lần tại mỗi vị trí trong suốt quá trình xây dựng đồ thị. • Quy tắc 4: Nếu một sự kiện tái tổ hợp cần được thực hiện, một cặp trình tự có đoạn đầu chung dài nhất sẽ được chọn để thực hiện tái tổ hợp. Nếu có hơn một ứng cử viên có cùng độ dài đoạn đầu chung dài nhất, một trong số đó sẽ được chọn một cách ngẫu nhiên. • Quy tắc 5: Với 1 cặp trình tự được chọn cho bước tái tổ hợp, trình tự chứa ít vật liệu di truyền trong phần đoạn đầu chung hơn sẽ được chọn để thực hiện tái tổ hợp. Sau bước tái tổ hợp, trình tự con chứa phần đoạn đầu chung dài nhất sẽ được thực hiện kết hợp với trình tự còn lại ngay sau đó. Điều này đảm bảo cho sự kiện tái tổ hợp chỉ thay thế một trình tự bằng một trình tự mới có ít vật liệu di truyền hơn trong thế hệ tiếp theo. Thuật toán ARG4WG ĐẦU VÀO: Tập dữ liệu D = {S1, , SN} trình tự (các haplotype) độ dài m, Sx[i] có giá trị bằng 0 hoặc 1, 1 ≤ x ≤ N, 1 ≤ i ≤ m. ĐẦU RA: một đồ thị ARG mô tả các mối quan hệ (các sự kiện kết hợp, đột biến, tái tổ hợp) giữa các nút (các trình tự) trong đồ thị đến một tổ tiên chung gần nhất SSCA. BEGIN t = 1; Dt = D; Gán danh sách kết hợp C = {tất cả các cặp (Sx,Sy){d,l} (1 ≤ x, y ≤ N) có l = m}; Gán danh sách đột biến M = {tất cả các trình tự chứa các vị trí đột biến đơn}; Gán danh sách tái tổ hợp R = {tất cả các cặp (Sx,Sy){d,l} (1 ≤ x, y ≤ N) có 0 < l < m }; while chưa đạt tới một tổ tiên chung duy nhất do if (danh sách kết hợp C không rỗng) then Lấy ngẫu nhiên một cặp trình tự có đoạn đầu chung (S1,S2); 59
Thực hiện kết hợp như sau: Gán S’ = S1 nếu L(S1) > L(S2); ngược lại S’ = S2; Dt+1 = (Dt\{S1,S2}) ∪ {S’}; Cập nhật 3 danh sách C, M, R; else if (danh sách đột biến M không rỗng) then Lấy ngẫu nhiên một trình tự S với một đột biến tại vị trí i; Thực hiện sự kiện đột biến như sau: Dt+1 = (Dt\{S})∪{S’} với S’[i] = S[i] và S’[j] = S[j] với mọi j ≠ i và 1 ≤ i,j ≤ m Cập nhật 3 danh sách C, M, R; else Lấy cặp trình tự có đoạn đầu chung dài nhất (S1,S2){d,l} từ danh sách tái tổ hợp; Thực hiện tái tổ hợp như sau: if d = left then // đoạn đầu chung dài nhất của (S1, S2) là từ đầu phía bên trái Gán SR = S1 nếu (L(S1) < L(S2))[1,l]; ngược lại SR = S2; SR1[i] = SR[i] với mọi 1 ≤ i ≤ l; SR1[j] = * với mọi l < j ≤ m SR2[i] = * với mọi 1 ≤ i ≤ l; SR2[j] = SR[j] với mọi l < j ≤ m else //d = right Gán SR = S1 nếu (L(S1) < L(S2))[m-l+1,m]; ngược lại SR = S2; SR1[i] = * với mọi 1 ≤ i ≤ m-l; SR1[j] = SR[j] với mọi m-l < j ≤ m SR2[i] = SR[i] với mọi 1 ≤ i ≤ m-l; SR2[j] = * với mọi m-l < j ≤ m endif Dt+1 = (Dt\{SR}) ∪ {SR1,SR2}; Cập nhật 3 danh sách C, M, R; endif endif endwhile END; Thuật toán 2.1: Thuật toán ARG4WG xây dựng một ARG từ một tập trình tự D cho trước. Như vậy, một sự kiện kết hợp làm giảm số lượng trình tự đi một. Một sự kiện đột biến chỉ làm thay đổi giá trị tại 1 vị trí của 1 trình tự, không làm tăng số lượng trình 60
tự. Một sự kiện tái tổ hợp chỉ thay thế một trình tự bằng một trình tự mới có ít vật liệu di truyền hơn con của nó. Chính vì vậy, ARG4WG luôn đạt tới được một tổ tiên chung duy nhất. Những lựa chọn ngẫu nhiên trong các bước của thuật toán dẫn tới việc sinh ra các đồ thị ARG khác nhau cho mỗi lần chạy. 2.3. Kết quả thực nghiệm Thuật toán ARG4WG được cài đặt sử dụng ngôn ngữ C++, có khả năng xử lý đa luồng. Mã nguồn của thuật toán có thể được truy cập tại địa chỉ: Do các phương pháp khác không khả thi cho các tập dữ liệu lớn nên ARG4WG được tiến hành thực nghiệm và so sánh với thuật toán Margarita. Các thực nghiệm được tiến hành trên hệ thống Rocks Cluster của Trung tâm tin học và tính toán, Viện Hàn lâm khoa học và công nghệ Việt Nam với cấu hình CPU 2600 MHz x 24 lõi, RAM 32GB, HDD 1TB. Lưu ý rằng các thực nghiệm so sánh ARG4WG với Margarita được tiến hành trên cùng hệ thống máy tính và ARG4WG không chạy đa luồng. Đầu tiên, các so sánh về thời gian chạy và số sự kiện tái tổ hợp của các ARG được xây dựng sử dụng Margarita và ARG4WG được tiến hành trên các tập dữ liệu thật, là các tập dữ liệu SNP haplotype gồm hàng nghìn trình tự và hàng nghìn đến chục nghìn SNP được trích xuất từ dự án 1000 hệ gen người (1kGP). Tiếp theo, hình thái cây được sinh ra từ 2 thuật toán được so sánh trên các tập dữ liệu mô phỏng (sử dụng các chương trình máy tính để tạo dữ liệu). 2.3.1. Các kết quả trên dữ liệu thật Thời gian chạy và số sự kiện tái tổ hợp của 2 thuật toán ARG4WG và Margarita được so sánh trên 3 tập dữ liệu trích xuất từ dự án 1000 hệ gen người (1kGP) [12] gồm 500, 1000, và 2000 trình tự với 1000, 2000, 5000, và 10000 SNP. Với mỗi cấu hình, chúng tôi thực hiện 3 thực nghiệm khác nhau, mỗi thực nghiệm tương ứng với 61
một ARG trên 1 vùng khác nhau của nhiễm sắc thể 1. Tổng hợp dữ liệu thử nghiệm được mô tả trong Bảng 2.1. Bảng 2.1: Tập dữ liệu trích xuất từ dự án 1000 hệ gen người. Tập dữ Số thử Số trình Số SNP liệu nghiệm tự 1 12 500 {1000, 2000, 5000, 10000} 2 12 1000 {1000, 2000, 5000, 10000} 3 12 2000 {1000, 2000, 5000, 10000} Trung bình thời gian chạy và số sự kiện tái tổ hợp được tính toán cho mỗi cấu hình. Do thuật toán Margarita sử dụng xác suất 0.9 cho việc chọn đoạn chung dài nhất và 0.1 cho việc chọn đoạn chung có độ dài bất kì cho bước tái tổ hợp, chúng tôi đã tiến hành thử nghiệm Margarita trên cả 2 trường hợp: cấu hình mặc định chọn đoạn chung có độ dài bất kì với xác suất 0.1 (Margarita) và cấu hình luôn chọn đoạn chung dài nhất (Margarita1.0). Hình 2.6 minh họa kết quả so sánh trung bình thời gian chạy giữa ARG4WG và Margarita, Margarita1.0. ARG4WG nhanh hơn hàng nghìn lần so với Margarita và Margarita1.0. Ví dụ, Margarita phải mất >105 giây để xây dựng 1 ARG cho 1000 trình tự và 10,000 SNP trong khi ARG4WG chỉ cần 140 giây (Hình 2.6b). Và Margarita không thể xây dựng được ARG cho 2000 trình tự và 10,000 SNP trong khi ARG4WG chỉ cần trung bình là 466 giây (Hình 2.6c). Ngoài ra, ta thấy Margarita1.0 hầu hết có thời gian chạy nhiều hơn so với Margarita . Điều này có thể được giải thích là do khi đó thuật toán mất nhiều chi phí thời gian cho việc tìm đoạn chung dài nhất giữa 2 trình tự để thực hiện sự kiện tái tổ hợp. 62
100000 10000 1000 100 (a) 10 Thời gian chạy (giây) chạy gian Thời 1 1000 2000 5000 10000 Số vị trí SNP 1000000 100000 10000 1000 100 (b) 10 Thời gian chạy (giây) chạy gian Thời 1 1000 2000 5000 10000 Số vị trí SNP 1000000 100000 10000 1000 100 (c) 10 1 Thời gian chạy (giây) chạy gian Thời 1000 2000 5000 10000 Số vị trí SNP Hình 2.6: Trung bình thời gian chạy của Margarita, Margarita1.0 và ARG4WG cho: (a) 500 haplotype; (b) 1000 haplotype; và (c) 2000 haplotype. 63
Hình 2.7 minh họa kết quả so sánh trung bình số sự kiện tái tổ hợp giữa ARG4WG và Margarita, Margarita1.0. Kết quả cho thấy thuật toán ARG4WG sinh ra các ARG có số sự kiện tái tổ hợp ít hơn rất nhiều so với Margarita. Số sự kiện tái tổ hợp của Margarita nhiều hơn trung bình là 1.4 lần so với ARG4WG trong tất cả các thực nghiệm. Ví dụ, Margarita sinh ra ARG trên 1000 trình tự và 10,000 SNP chứa trung bình 72,678 tái tổ hợp, trong khi ARG xây dựng trên tập dữ liệu đó sử dụng ARG4WG chỉ chứa trung bình 50,567 tái tổ hợp. Số sự kiện tái tổ hợp của Margarita1.0 ít hơn so với Margarita nhưng vẫn nhiều hơn so với ARG4WG từ 1.2 đến 1.3 lần. Các kết quả này một lần nữa khẳng định rằng chiến lược đoạn chung dài nhất không phải luôn cho ra đồ thị ARG có số sự kiện tái tổ hợp nhỏ nhất. 45000 40000 35000 30000 25000 20000 15000 (a) 10000 5000 Số sự kiện tái tổ hợp tổ tái kiện sự Số 0 1000 2000 5000 10000 Số vị trí SNP 80000 70000 60000 50000 40000 30000 (b) 20000 Số sự kiện tái tổ hợp tổ tái kiện sự Số 10000 0 1000 2000 5000 10000 Số vị trí SNP 64
100000 90000 80000 70000 60000 50000 40000 (c) 30000 20000 Số sự kiện tái tổ hợp tổ tái kiện sự Số 10000 0 1000 2000 5000 10000 Số vị trí SNP Hình 2.7: Trung bình số sự kiện tái tổ hợp của Margarita, Margarita1.0 và ARG4WG cho: (a) 500 haplotype; (b) 1000 haplotype; và (c) 2000 haplotype. Chúng tôi cũng đã kiểm tra khả năng thực thi của ARG4WG trên dữ liệu toàn nhiễm sắc thể. ARG4WG được chạy trên 4246 trình tự (2123 mẫu gen người) độ dài toàn nhiễm sắc thể 1 (174,234 SNP – nhiễm sắc thể dài nhất trong bộ gen người) từ dự án 1kGP. ARG4WG có thể sinh ra 1 ARG với thời gian ~ 4.5 giờ trong một lần chạy sử dụng 1 máy tính 16-thread. Kết quả này đã nhấn mạnh khả năng của ARG4WG trong việc thực thi trên dữ liệu hàng nghìn trình tự độ dài hệ gen người. 2.3.2. Các kết quả trên dữ liệu mô phỏng Trong thực nghiệm với dữ liệu thật, ta không biết được đồ thị ARG đúng (đồ thị ARG phản ánh đúng lịch sử tiến hóa của các trình tự quan sát). Do đó, để đánh giá độ hợp lý của đồ thị ARG sinh ra sử dụng thuật toán ARG4WG đề xuất, chúng tôi đã tiến hành các thực nghiệm trên dữ liệu mô phỏng. Để làm điều này, chúng tôi xây dựng các cây đúng (là cây được xác định ban đầu) và các bộ dữ liệu trình tự mô phỏng tương ứng sử dụng các chương trình máy tính. Sau đó, chúng tôi xây dựng đồ thị ARG trên các tập dữ liệu mô phỏng này sử dụng 2 thuật toán ARG4WG và Margarita. Tiếp theo, từ đồ thị ARG sinh ra từ 2 thuật 65
toán này chúng tôi sẽ trích xuất ra các cây biên và so sánh hình thái cây với cây đúng. MacS [11] được sử dụng để tạo ra 500 trình tự trên một vùng độ dài 1Mbp với kích thước quần thể hiệu quả là 15000. Tỉ lệ đột biến được thiết lập là 1.2x10-8 trên một vị trí trên một thế hệ. Giống với cách làm trong [58], 4 tập dữ liệu khác nhau được tạo ra với tỉ lệ của đột biến trên tái tổ hợp lần lượt là 1, 2, 4 và 6. Với mỗi tập dữ liệu, khoảng 800 SNP có vị trí gần với vị trí của SNP trên dữ liệu thật của 1Mbp dữ liệu của nhiễm sắc thể 1 (167,000,000-168,000,000) trong 1kGP được lựa chọn. Việc làm này nhằm mục đích lấy dữ liệu mô phỏng gần giống với dữ liệu thật nhất và có độ lớn phù hợp cho việc thử nghiệm. Mỗi thuật toán Margarita và ARG4WG được chạy để sinh ra 20 ARG. Sau đó, cây biên tại các vị trí SNP từ mỗi đồ thị ARG được trích rút ra. Các hình thái cây được sinh ra từ 2 thuật toán này được so sánh với cây đúng được sinh ra qua mô phỏng tại các vị trí tương ứng sử dụng khoảng cách Robinson-Fould (RF). Phần mềm PhyloNet [65] được sử dụng để tính khoảng cách RF này. Khoảng cách RF giữa cấu trúc của hai cây là tỷ lệ giữa số phân vùng chỉ có ở một trong hai cây trên tổng số phân vùng của cả hai cây. Khoảng cách RF có khoảng giá trị từ 0% đến 100%. Giá trị RF giữa hai cây càng nhỏ thì cấu trúc của hai cây càng giống nhau. Các kết quả khoảng cách RF trung bình với các tỉ lệ đột biến/tái tổ hợp được thể hiện trong Hình 2.8. Kết quả cho thấy Margarita có khoảng cách RF nhỏ hơn so với ARG4WG. Tuy nhiên, sự khác nhau giữa 2 thuật toán giảm dần với sự tăng lên của tỉ lệ đột biến/tái tổ hợp. Ví dụ, với tỉ lệ đột biến/tái tổ hợp bằng 1 thì độ chênh lệch giữa 2 thuật toán là ~2% nhưng với tỉ lệ bằng 6, khoảng cách RF của cả 2 thuật toán là gần như nhau (~0.14%). 66
64 62 60 58 56 ARG4WG 54 Margarita 52 Khoảng cách RF RF (%) cách Khoảng 50 1 2 4 6 Tỉ lệ đột biến/tái tổ hợp Hình 2.8: Khoảng cách RF của các cây được tạo ra bởi thuật toán Margarita và ARG4WG so với các cây đúng tương ứng trên các khoảng tỉ lệ đột biến và tái tổ hợp khác nhau. 2.4. Kết quả ứng dụng ARG4WG vào bài toán tìm vùng gen liên quan đến bệnh sốt rét ở Châu Phi Các kết quả thực nghiệm trong mục 2.3 đã chỉ ra khả năng nổi trội của thuật toán ARG4WG khi làm việc với dữ liệu lớn hàng nghìn trình tự trên toàn bộ nhiễm sắc thể. Để làm rõ hơn hiệu quả của thuật toán đề xuất cũng như chứng minh chất lượng của các đồ thị ARG sinh ra từ thuật toán, trong phần này, luận án trình bày ứng dụng ARG4WG vào bài toán ánh xạ tương quan giữa vùng gen liên quan tới bệnh quan tâm trên tập dữ liệu GWAS lớn. Cụ thể, luận án thực nghiệm ứng dụng ARG4WG vào bài toán tìm vùng gen liên quan đến bệnh sốt rét ở Châu Phi sử dụng tập dữ liệu tập Gambia [3]. Đây là tập dữ liệu case-control lớn gồm 2780 mẫu, trong đó có 1533 mẫu người khỏe mạnh và 1247 mẫu người bị bệnh sốt rét trên toàn nhiễm sắc thể 11 (22,028 SNP). Qua các khảo sát chúng tôi nhận thấy, chương trình Margarita vừa xây dựng đồ thị ARG hợp lý lại vừa xây dựng thuật toán tìm ánh xạ tương quan dựa trên đồ thị ARG hợp lý đó, nên đây là lựa chọn phù hợp cho việc thử nghiệm các đồ thị ARG kết quả từ thuật toán ARG4WG vào bài toán tìm ánh xạ tương quan toàn hệ gen. 67
Thuật toán tìm ánh xạ tương quan trong chương trình Margarita được thực hiện như sau [52]: Đầu vào: Tập dữ liệu D gồm các trình tự nhị phân (SNP haplotype). Các trình tự này có nhãn người bệnh và người không bệnh. Đầu ra: Điểm số về độ phân tách giữa người bệnh và người không bệnh tại từng vị trí SNP. Độ phân tách càng lớn thì gần vị trí SNP đó có nhiều khả năng liên quan đến bệnh. Do các lựa chọn ngẫu nhiên trong thuật toán xây dựng đồ thị ARG nên với cùng một tập dữ liệu đầu vào D, mỗi lần chạy thuật toán xây dựng đồ thị ARG trong Margarita sẽ cho ra một đồ thị ARG khác nhau. Để thực hiện ánh xạ tương quan, Margarita chạy n lần để có được n ARG khác nhau. Tại mỗi vị trí SNP c, n cây biên được trích rút ra từ n ARG. Và trong mỗi cây biên T, mỗi cạnh được tính điểm số về độ phân biệt giữa các cá thể bệnh với các cá thể khỏe mạnh. Sau đó, điểm số cao nhất trong số các cạnh đó được gán cho T. Cuối cùng, các điểm số liên quan đến n cây biên tại vị trí c được tính trung bình cho ra điểm số tương quan toàn cục cho vị trí c. Điểm số tương quan là thước đo hỗ trợ cho giả thuyết rằng có một đột biến gây bệnh xuất hiện gần vị trí c. Cụ thể, điểm số cho 1 cạnh e trong cây biên T cho vị trí c là một kiểm định tính độc lập sử dụng phân bố chi-bình-phương, sử dụng một bảng 2x2. Một trục của bảng là thông tin các cá thể là người bệnh hay người không bệnh; và một trục là thông tin các cá thể có lá nằm dưới e và các cá thể có lá không nằm dưới e trong T. Thống kê chi-bình-phương được tính cho e sẽ kiểm định xem việc bị bệnh và đặc tính nằm dưới cạnh e trong T có độc lập hay không. Từ đó, nó sẽ cung cấp một một điểm số cho giả thuyết rằng có một đột biến gây bệnh xuất hiện trên cạnh e trong T và một điểm số cho giả thuyết rằng đột biến xuất hiện gần vị trí c trong nhiễm sắc thể. Giá trị thống kê càng cao thì khả năng giả thuyết đúng càng lớn. Để kết hợp các điểm số 68
cho vị trí c, n kiểm định chi-bình-phương (tương ứng với n cây biên) sẽ được tính trung bình và cho ra điểm số tương quan cho vị trí c. Tiếp theo, một kiểm định hoán vị được sử dụng để đánh giá ý nghĩa thống kê của điểm số tương quan cho c. Nhiều lần (trong Margarita đặt mặc định là 10,000 lần), các nhãn người bệnh và người không bệnh được hoán vị ngẫu nhiên mà không thay đổi số lượng của từng loại. Một điểm số cho c sau đó được tính toán như trên cho mỗi hoán vị. P-value cho c được định nghĩa là tỉ lệ hoán vị trong đó điểm số tương quan được tính cho c cao hơn điểm số tương quan thu được sử dụng phân vùng đúng các nhãn. P-value cho vị trí c nhỏ dưới một ngưỡng nào đó sẽ ủng hộ cho giả thuyết rằng một đột biến gây bệnh xuất hiện gần vị trí c. Để tiến hành thực nghiệm, chúng tôi tiền xử lý dữ liệu Gambia, sử dụng công cụ SHAPEIT [15] để chuyển đổi dữ liệu kiểu gen sang dữ liệu haplotype. Dữ liệu haplotype (5560 haplotype cho 2780 mẫu) sẽ là đầu vào cho thuật toán ARG4WG để xây dựng các ARG. Với mỗi ARG, chúng tôi trích xuất ra các cây biên cho từng vị trí SNP và các cây biên này là đầu vào cho Margarita để thực hiện ánh xạ tương quan, tìm vùng gen liên quan đến bệnh sốt rét. ARG4WG mất khoảng 3.1 giờ để xây dựng một ARG sử dụng máy tính 16-thread cho toàn bộ tập dữ liệu Gambia. Chúng tôi thực hiện kiểm thử ánh xạ với 106 hoán vị với 3 kịch bản với số lượng SNP khác nhau như sau: • 10 ARG xây dựng từ toàn bộ tập dữ liệu Gambia bao gồm 5560 haplotype trên toàn nhiễm sắc thể 11 (22,028 SNP); • 30 ARG xây dựng từ 5560 haplotype độ dài 5000 SNP xung quanh gen HBB; • 50 ARG xây dựng từ 5560 haplotype và 1000 SNP xung quanh gen HBB; Các kết quả thực nghiệm được biểu thị trong Hình 2.9. 69
Hình 2.9: Sự tương quan đến bệnh từ 106 kiểm định hoán vị trên: (A) 10 ARG xây dựng trên toàn bộ NST 11; (B) 30 ARG xây dựng trên vùng 5000 SNP quanh gen HBB; và (C) Tổng hợp kết quả cho các thực nghiệm trên vùng 1000 SNP quanh gen HBB. Kết quả cho thấy cả 3 thực nghiệm đều phát hiện ra vùng có tín hiệu mạnh liên quan 70
đến bệnh là từ 4.43Mb tới 6.28Mb trên nhiễm sắc thể 11 với p-value ≤ 10-7. Kết quả này phù hợp với phân tích của nhóm tác giả Garvin Band [3] là vùng gen HBB (4.5Mb-5.5Mb) có nhiều khả năng liên quan đến bệnh sốt rét. Trong nghiên cứu của họ, họ đã chỉ ra giá trị P-value thấp nhất trong vùng này là 5.7x10-13 sử dụng phương pháp phân tích SNPTEST meta-analysis. Tuy nhiên, với dữ liệu vào như trên thì thuật toán Margarita không thực hiện được phân tích sâu đến >107. Sự thống nhất trong kết quả của 3 thực nghiệm trên cũng chỉ ra rằng ARG4WG rất ổn định trong việc xây dựng ARG với số lượng SNP khác nhau. Chúng tôi cũng đã kiểm tra thực thi của Margarita trong việc xác định vùng gen liên quan đến bệnh cho dữ liệu Gambia. Vì Margarita không thể thực thi trên cả nhiễm sắc thể nên ba thực nghiệm với 106 hoán vị trong vùng 4M đến 6M quanh gen HBB gồm 600 SNP đã được thực hiện: - Test1: 10 ARG sinh ra từ 2000 haplotype (1000 cá thể bệnh, 1000 cá thể khỏe mạnh). - Test2: 10 ARG sinh ra từ 3000 haplotype (1500 cá thể bệnh, 1500 cá thể khỏe mạnh). Test3: 10 ARG sinh ra từ 5560 haplotype (toàn bộ cá thể trong tập Gambia). Để xây dựng 10 ARG, Margarita mất 30 giờ cho Test1, 90 giờ cho Test2, và 420 giờ cho Test3. Các kết quả thực nghiệm trong Hình 2.10 chỉ ra rằng càng dùng nhiều dữ liệu thì Margarita càng cho kết quả tốt hơn. Với tập dữ liệu nhỏ, ít mẫu cá thể, Margarita chỉ cho ra các tín hiệu yếu trong vùng gen HBB. Tín hiệu trở nên mạnh hơn khi nhiều mẫu dữ liệu được sử dụng. Với toàn bộ mẫu dữ liệu như trong Test3, Margarita đã phát hiện được tín hiệu mạnh liên quan đến bệnh trong vùng từ 4.46Mb đến 5.98Mb. 71