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

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 2 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. 3 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 4 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 5 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 6 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 7 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 8 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: https://genomainternational.com/introduction-to-genomics/. ................................... 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 9 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 10 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 ( )....................... 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 ls 11 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 12 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 13 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ế. 14 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ó. 15 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. 16 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: https://genomainternational.com/introduction-to-genomics/. 17 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 18 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]. 19 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. 20 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 đề...há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 40 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 41 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] https://csiflabs.cs.ucdavis.edu/ ~gusfield/galledtree.tar - Mã nguồn mở. - Ngôn ngữ lập trình: Python - Yêu cầu: Python phiên bản 2. 42 SHRUB [64] https://people.eecs.berkeley.e du/~yss/lu.html - Chỉ cung cấp chương trình chạy. - Hỗ trợ chạy trên hệ điều hành MS Windows, Linux và Mac OS X. SHRUB-GC [61] https://people.eecs.berkeley.e du/~yss/lugc.html - Phiên bản mở rộng của SHRUB 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] https://github.com/RabadanLa b/TARGet - Mã nguồn mở - 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] ftp://ftp.sanger.ac.uk/pub/reso urces/software/margarita/ - Mã nguồn mở. - Ngôn ngữ lập trình: Java. ARGweaver [58] https://github.com/mdrasmus/ argweaver - Mã nguồn mở. - Yêu cầu cài đặt: Trình biên dịch C++ và Python. Arbores [29] https://github.com/heinekmp/ Arbores - Mã nguồn mở. - Ngôn ngữ lập trình: C Bacter [66] https://tgvaughan.github.io/ba cter/ - Chương trình xây dựng đồ thị 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. 43 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). 44 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á 45 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. 46 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. 47 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]. 48 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 49 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ử 50 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)). 51 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, ),( Dsf , để 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: 0),( =Dsf nếu s có thể kết hợp với một trình tự trong D. 1),( =Dsf nếu s có thể được biểu diễn thành ),( 21 ss trong đó 1s và 2s có thể kết hợp với các trình tự trong D. Ngoài ra, đặt 321 ssss ++= trong đó 1s và 2s có thể kết hợp với các trình tự trong D. ),( Dsf có thể được viết dưới dạng: }2),({min )},({min),( 3,, 321,, 321 321 + =++= Dsf DsssfDsf sss sss (1) 52 Đặt rl sss += trong đó ls 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 ls . Chúng ta sẽ chứng minh rằng: 1),(}2),({min),( 3),,( 321 +=+= DsfDsfDsf rsss Đặt ),,( * 3 * 2 * 1 sss là giải pháp tối ưu của phương trình (1) và x s là phần bù của trong * 2 * 1 ss + , lx ssss −+= * 2 * 1 . 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 ( ). Rõ ràng là xs là một chuỗi con của * 2s vì * 1s là một chuỗi con của ls . Do đó, nếu * 2s có thể kết hợp với một trình tự h trong D thì xs có thể kết hợp với h (xem Hình 2.3). Điều này dẫn tới: ls ls 53 1),( 2)*,( ),(),( 3 * 3 += += ++= Dsf Dsf DsssfDsf r xl 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 54 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. 55 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. 56 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. 57 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. 58 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ụ. 59 • 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 chun...12] với số trình tự và số SNP được mô tả trong Bảng 3.1. Chúng tôi tiến hành so sánh các thuật toán GAMARG, Margarita, ARG4WG, REARG về số sự kiện tái tổ hợp so với các thuật toán vét cạn trên các bộ dữ liệu nhỏ. Với bộ dữ liệu lớn hơn như bộ dữ liệu trích xuất từ dự án 1kGP, các thuật toán 89 vét cạn không có khả năng thực thi được, chúng tôi tiến hành so sánh, đánh giá 4 thuật toán trên về số sự kiện tái tổ hợp và thời gian chạy của các thuật toán. Bảng 3.1: Tập dữ liệu từ dự án 1kGP. Tập dữ liệu Số trình tự Số SNP {DS1, DS2, DS3} 100 2000 {DS1, DS2, DS3} 100 5000 {DS1, DS2, DS3} 100 10000 {DS1, DS2, DS3} 200 2000 {DS1, DS2, DS3} 200 5000 {DS1, DS2, DS3} 200 10000 REARG có 3 phiên bản là REARG_SIM, REARG_LEN, REARG_COM. Đối với các thực nghiệm mà kết quả của cả 3 phiên bản là như nhau thì chúng tôi để là kết quả của thuật toán REARG nói chung. 3.6.1. Kết quả trên các tập dữ liệu nhỏ Với mỗi thuật toán, 10000 lần chạy được thực hiện để sinh ra 10000 ARG và ARG được xây dựng với ít số sự kiện tái tổ hợp nhất được ghi lại. Đối với thuật toán GAMARG, nhiều giá trị khác nhau cho tham số ઠ đã được chúng tôi thực nghiệm và kết quả tốt nhất thu được là 18 ≤ 𝛿 < 43 cho tập dữ liệu Kreitman; 29 ≤ δ < 54 cho tập SDS1 và 11 ≤ δ < 45 cho tập SDS2. Các kết quả của các thuật toán được mô tả trong Bảng 3.2. 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ỏ. Kreitman SDS1 SDS2 MinARG 7 10 12 Margarita 8 14 18 ARG4WG 10 17 18 REARG 10 17 20 GAMARG 7 10 13 90 Kết quả cho thấy, chiến lược đoạn chung dài nhất của Margarita tỏ ra hiệu quả hơn so với chiến lược đoạn đầu chung dài nhất trong ARG4WG và REARG khi cho ra kết quả số sự kiện tái tổ hợp ít hơn trong trường hợp dữ liệu nhỏ này. Thuật toán REARG sử dụng các ràng buộc hạn chế hơn so với ARG4WG, đó là, thay vì lựa chọn ngẫu nhiên trong số các cặp trình tự có cùng độ dài đoạn chung dài nhất như trong thuật toán ARG4WG, REARG sử dụng thêm các ràng buộc về đặc trưng dữ liệu để lựa chọn ra ứng cử viên cho bước tái tổ hợp. Cách làm này dẫn đến REARG sinh ra các đồ thị ARG với số sự kiện tái tổ hợp luôn lớn hơn hoặc bằng số sự kiện tái tổ hợp trong ARG4WG. Các kết quả thử nghiệm trên cho thấy chiến lược sử dụng trong thuật toán REARG và ARG4WG không phù hợp với các tập dữ liệu nhỏ. Tuy nhiên, các kết quả của Margarita, ARG4WG, và REARG đều còn rất xa các kết quả tối ưu. Thuật toán GAMARG cho kết quả tốt nhất khi có thể đạt tới ARG tối thiểu với tập Kreitman và tập SDS1 và chỉ hơn 1 sự kiện tái tổ hợp so với phương pháp tối ưu ở tập SDS2. Kết quả này chỉ ra tính hiệu quả của phương pháp kết hợp kiểm thử 4 giao tử và đoạn đầu chung dài nhất trong thuật toán GAMARG khi chạy với tập dữ liệu nhỏ. 3.6.2. Các kết quả trên các tập dữ liệu từ dự án 1kGP Với mỗi thuật toán, 1000 lần chạy được thực hiện để sinh ra 1000 ARG và ARG được xây dựng với ít số sự kiện tái tổ hợp nhất được ghi lại. Do thuật toán Margarita mất rất nhiều thời gian để xây dựng 1 đồ thị ARG trên các tập dữ liệu này, chính vì vậy, ARG có ít số sự kiện tái tổ hợp nhất được tìm thấy sau 9 ngày chạy thuật toán được ghi lại là kết quả cho thuật toán Margarita. Đối với thuật toán GAMARG, các giá trị khác nhau (tức là, 5, 10, 15, 20, 25, và 30) cho tham số 𝛿 đã được thử nghiệm trên các tập dữ liệu có các kích thước khác nhau. 5000 ARG được xây dựng và ARG với số sự kiện tái tổ hợp ít nhất được ghi lại trên mỗi tập dữ liệu. Các kết quả cho thấy GAMARG cho ra các kết quả tương tự nhau 91 khi 𝛿 bằng một trong các giá trị 5, 10, 15 với độ dài 500 SNP. Tuy nhiên, với các trình tự dài hơn, tức là, 1000 và 2000 SNP thì thuật toán cho ra số sự kiện tái tổ hợp ít nhất với tham số 𝛿 = 5. Do đó, trong các thực nghiệm này, chúng tôi chạy GAMARG sử dụng 𝛿 = 5. Các kết quả thực nghiệm (xem Bảng 3.3 và Bảng 3.4) cho thấy các thuật toán REARG và GAMARG đề xuất đều cho ta các đồ thị ARG với số sự kiện tái tổ hợp ít hơn so với thuật toán Margarita và ARG4WG. Đặc biệt, thuật toán GAMARG đều có kết quả vượt trội trong tất cả các thực nghiệm. 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. Số SNP 2000 5000 10000 Margarita 30749 77574 13234 ARG4WG 2596 5837 10490 (a) REARG_SIM 2579 5786 10324 REARG_LEN 2560 5741 10336 REARG_COM 2560 5766 10307 GAMARG 2441 5610 10052 Số SNP 2000 5000 10000 Margarita 33820 86041 10874 ARG4WG 1924 4335 8741 (b) REARG_SIM 1949 4284 8623 REARG_LEN 1924 4288 8627 REARG_COM 1930 4279 8613 GAMARG 1824 4093 8298 92 Số SNP 2000 5000 10000 Margarita 32807 81829 11714 ARG4WG 1555 4120 9100 (c) REARG_SIM 1545 4087 8987 REARG_LEN 1527 4077 8960 REARG_COM 1545 4063 8934 GAMARG 1502 3994 8817 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. Số SNP 2000 5000 10000 Margarita 55788 139136 22692 ARG4WG 4218 9583 17083 (a) REARG_SIM 4176 9480 16860 REARG_LEN 4180 9437 16834 REARG_COM 4182 9451 16829 GAMARG 4099 9315 16715 Số SNP 2000 5000 10000 Margarita 60337 153145 18467 ARG4WG 3072 6998 14158 (b) REARG_SIM 3027 6909 13900 REARG_LEN 3020 6877 13887 REARG_COM 3021 6861 13925 GAMARG 2955 6679 13606 93 Số SNP 2000 5000 10000 Margarita 58313 144777 20047 ARG4WG 2620 6652 14813 (c) REARG_SIM 2586 6607 14620 REARG_LEN 2595 6564 14634 REARG_COM 2584 6600 14642 GAMARG 2521 6504 14583 Ta thấy, với các trình tự dài (10,000 SNP), Margarita sinh ra đồ thị ARG có số sự kiện tái tổ hợp nhiều hơn gấp khoảng 1.3 lần so với thuật toán REARG. Tuy nhiên, Margarita cần một lượng lớn số sự kiện tái tổ hợp, gấp khoảng 10 đến 20 lần so với REARG để xây dựng ARG cho các tập dữ liệu có trình tự ngắn hơn (2000 và 5000 SNP). Như vậy, Margarita cần số sự kiện tái tổ hợp để xây dựng 1 ARG cho dữ liệu độ dài 10,000 SNP thậm chí còn ít hơn rất nhiều so với các tập dữ liệu độ dài 2000 và 5000 SNP. Kết quả cho thấy, thuật toán Margarita không ổn định. Điều này có thể được lý giải là do chiến lược đoạn chung dài nhất của thuật toán. Margarita luôn cần đến gấp đôi số sự kiện tái tổ hợp để thực hiện tái tổ hợp khi đoạn chung dài nhất được tìm thấy bên trong trình tự. Vì vậy, khi đoạn chung dài nhất luôn được tìm phía bên trong trình tự thì thuật toán sẽ bị bùng nổ về số nút cũng như số sự kiện tái tổ hợp. Mặc dù cả 3 phiên bản của thuật toán REARG đều tốt hơn so với ARG4WG nhưng không có phiên bản nào tốt hơn hoàn toàn trong tất cả các thực nghiệm. REARG_LEN và REARG_COM tốt hơn REARG_SIM trên hầu hết các tập dữ liệu; tuy nhiên, REARG_SIM lại là thuật toán tốt nhất trong 3 phiên bản của REARG với 2 tập dữ liệu: tập DS1 với 200 trình tự và 2000 SNP và tập DS3 với 200 trình tự và 10,000 SNP. 94 Hình 3.6 minh họa các kết quả của 3 thuật toán GAMARG, REARG và ARG4WG. Trong đó, kết quả của thuật toán REARG là kết quả tốt nhất trong 3 phiên bản REARG_LEN, REARG_COM, và REARG_SIM trên mỗi thử nghiệm. 100 trình tự 2000 SNP 5000 SNP 10000 SNP 200 trình tự 2000 SNP 5000 SNP 10000 SNP 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. Thuật toán GAMARG cho kết quả tốt hơn hẳn các thuật toán còn lại trong tất cả các thực nghiệm. Sự vượt trội của GAMARG so với ARG4WG và REARG được thể hiện rõ trong các thực nghiệm với 100 trình tự. Với các tập dữ liệu lớn hơn với 1480 1780 2080 2380 2680 DS1 DS2 DS3 3900 4200 4500 4800 5100 5400 5700 DS1 DS2 DS3 8200 8500 8800 9100 9400 9700 10000 10300 10600 DS1 DS2 DS3 2500 2800 3100 3400 3700 4000 4300 DS1 DS2 DS3 6400 6700 7000 7300 7600 7900 8200 8500 8800 9100 9400 DS1 DS2 DS3 13500 13800 14100 14400 14700 15000 15300 15600 15900 16200 16500 16800 17100 DS1 DS2 DS3 95 nhiều trình tự hơn, độ đa dạng của dữ liệu tăng lên, do đó có nhiều cặp vị trí không tương thích hơn. Tuy nhiên, chỉ số ít trong số đó có thể thỏa mãn ràng buộc là ít nhất có một loại giao tử có tần số bằng 1. Trong trường hợp này, lợi thế của GAMARG so với ARG4WG và REARG là không đáng kể. Trung bình thời gian chạy để xây dựng một đồ thị ARG được tính toán cho mỗi thuật toán trên mỗi tập dữ liệu, các kết quả được tổng hợp trong Bảng 3.5 và Bảng 3.6. Do sự bùng nổ về số sự kiện tái tổ hợp trong thuật toán Margarita khiến cho thuật toán này có thời gian tính toán rất cao, gấp từ hàng trăm đến hàng nghìn lần so với các thuật toán đề xuất. Thời gian chạy của thuật toán REARG_LEN gần bằng với thời gian chạy của thuật toán ARG4WG. Do việc tính toán độ tương đồng mất nhiều thời gian hơn nên thuật toán REARG_SIM và REARG_COM tốn nhiều thời gian so với thuật toán ARG4WG và REARG_LEN. Thuật toán GAMARG có thời gian chạy chậm hơn so với ARG4WG và REARG_LEN nhưng nhanh hơn REARG_SIM và REARG_COM. 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. Số SNP 2000 5000 10000 Margarita 442.00 5690.00 201.50 ARG4WG 0.57 1.26 2.08 (a) REARG_SIM 3.27 17.66 60.64 REARG_LEN 0.62 1.35 2.46 REARG_COM 3.39 18.34 62.80 GAMARG 4.15 10.18 30.83 96 Số SNP 2000 5000 10000 Margarita 754.00 7709.00 140.00 ARG4WG 0.41 0.99 1.90 (b) REARG_SIM 2.41 12.98 53.50 REARG_LEN 0.44 1.06 2.15 REARG_COM 2.56 13.77 55.88 GAMARG 2.41 6.83 21.09 Số SNP 2000 5000 10000 Margarita 743.50 7627.50 166.00 ARG4WG 0.36 0.96 1.96 (c) REARG_SIM 2.04 13.01 53.63 REARG_LEN 0.39 1.02 2.08 REARG_COM 2.22 13.34 55.72 GAMARG 1.65 6.30 23.99 Các thuật toán này tuy không có khả năng chạy được với dữ liệu hàng nghìn trình tự độ dài hàng trăm nghìn SNP như ARG4WG nhưng có khả năng chạy được với dữ liệu hàng nghìn trình tự độ dài hàng chục nghìn SNP. 97 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. Số SNP 2000 5000 10000 Margarita 2826.51 31459.54 1157.00 ARG4WG 1.33 3.26 5.50 (a) REARG_SIM 8.84 49.85 168.09 REARG_LEN 1.43 3.51 5.76 REARG_COM 9.71 49.20 169.77 GAMARG 14.91 34.42 105.01 Số SNP 2000 5000 10000 Margarita 3953.000 42623.510 839.503 ARG4WG 0.998 2.442 4.739 (b) REARG_SIM 6.181 35.391 145.266 REARG_LEN 0.999 2.558 4.879 REARG_COM 7.001 35.246 146.843 GAMARG 5.886 20.043 59.926 98 Số SNP 2000 5000 10000 Margarita 3710.50 39230.55 1074.52 ARG4WG 0.96 2.42 4.94 (c) REARG_SIM 5.55 34.77 149.67 REARG_LEN 0.95 2.55 5.45 REARG_COM 6.28 34.46 153.76 GAMARG 3.95 18.62 60.35 3.7. Kết luận chương Chiến lược đoạn chung dài nhất của Margarita và chiến lược đoạn đầu chung dài nhất của ARG4WG đều không tối ưu số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Các thực nghiệm đã cho thấy chiến lược đoạn chung dài nhất trong Margarita phù hợp hơn so với chiến lược đoạn đầu chung dài nhất trong ARG4WG trên các bộ dữ liệu nhỏ. Tuy nhiên, Margarita trở nên không ổn định và cho ra các đồ thị ARG có số sự kiện tái tổ hợp nhiều hơn so ARG4WG trên các bộ dữ liệu vừa và lớn. Trong chương này, luận án đã giới thiệu hai thuật toán REARG và GAMARG cải tiến từ ARG4WG nhằm tối ưu thêm số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Thông qua việc kết hợp thêm các đặc trưng về độ tương đồng và độ dài của trình tự được chọn thực hiện tái tổ hợp bên cạnh ràng buộc đoạn đầu chung dài nhất, thuật toán REARG có thể giúp tìm ra các ARG có số sự kiện tái tổ hợp ít hơn so với ARG4WG với những tập dữ liệu vừa và lớn. Tức là, 2 đặc trưng này là phù hợp, có thể giúp quá trình tìm kiếm khu trú được vào các ARG có số sự kiện tái tổ hợp nhỏ nhanh hơn trong hữu hạn số lần chạy thuật toán. Tuy nhiên, không đặc trưng nào là 99 tốt nhất và việc tìm kiếm cần tiến hành song song giữa các phiên bản khác nhau của thuật toán REARG. Kiểm thử 4 giao tử là một kĩ thuật đơn giản nhưng hiệu quả trong việc tính toán đồ thị ARG tối thiểu cho các bộ dữ liệu nhỏ. Chiến lược đoạn đầu chung dài nhất trong ARG4WG rất hiệu quả đối với các bộ dữ liệu lớn. Sự kết hợp của chúng trong GAMARG làm cho thuật toán là tốt nhất so với ARG4WG và REARG về số sự kiện tái tổ hợp. GAMARG không những chạy được với dữ liệu hàng nghìn trình tự độ dài hàng nghìn đến chục nghìn SNP mà còn có khả năng so sánh được với các thuật toán xây dựng ARG tối thiểu trên các bộ dữ liệu nhỏ. Các kết quả nghiên cứu của chương này đã được công bố trong một bài kỉ yếu Hội nghị quốc tế NAFOSTED năm 2017 (công trình khoa học số 2) và một bài kỉ yếu Hội nghị quốc tế ICBBB năm 2019 (công trình khoa học số 3). 100 KẾT LUẬN Xác định nguồn gốc di truyền của bệnh bằng việc xác định các gen và alen nhạy cảm với bệnh là mục tiêu then chốt của nghiên cứu di truyền học con người. Đồ thị tái tổ hợp di truyền đóng một vai trò quan trọng trong nghiên cứu di truyền quần thể, đa dạng hệ gen và đa hình di truyền SNP. Tuy nhiên, bài toán xây dựng đồ thị ARG là một bài toán NP-khó và đòi hỏi tính toán khối lượng lớn nên ứng dụng vào thực tế còn hạn chế. Thông qua việc nghiên cứu các phương pháp xây dựng đồ thị ARG, tập trung theo hướng tiếp cận xây dựng đồ thị ARG có ít số sự kiện tái tổ hợp nhất và thuật toán Margarita, luận án đã đề xuất thuật toán ARG4WG xây dựng đồ thị ARG hợp lý cho dữ liệu lớn hàng nghìn trình tự trên toàn hệ gen. Bằng cách tiếp cận vấn đề theo cách của Margarita, cải tiến sử dụng đoạn đầu chung dài nhất cho bước tính toán sự kiện tái tổ hợp, thuật toán ARG4WG đề xuất không những giúp làm giảm đáng kể thời gian tìm kiếm các đoạn chung dài nhất sau mỗi lần thực hiện bước tái tổ hợp mà còn đảm bảo số nút trong đồ thị luôn được ổn định trong quá trình xây dựng đồ thị ARG. Kết quả thực nghiệm cho thấy thuật toán ARG4WG nhanh hơn hàng trăm đến hàng nghìn lần thuật toán Margarita. Đặc biệt, ARG4WG có thể chạy được với hàng nghìn trình tự trên toàn nhiễm sắc thể trong 1 lần chạy trong khoảng thời gian hợp lý thông qua xử lý đa luồng. Thuật toán ARG4WG cũng đã được thử nghiệm ứng dụng vào một bài toán thực tế xác định tương quan toàn bộ nhiễm sắc thể trên tập dữ liệu lớn, cụ thể là trong bài toán tìm vùng gen liên quan đến bệnh sốt rét ở Châu Phi trên 5560 trình tự độ dài toàn nhiễm sắc thể 11. Kết quả vùng tín hiệu bệnh sốt rét tìm được trùng với các kết quả phân tích đã có. Các kết quả này đã cho thấy khả năng ứng dụng của thuật toán ARG4WG vào các bài toán thực tế trên dữ liệu lớn. Thuật toán được cài đặt và để 101 dưới dạng mã nguồn mở cho cộng đồng nghiên cứu tại địa chỉ: https://github.com/thaontp711/arg4wg. Luận án cũng đã đề xuất 2 thuật toán cải tiến từ thuật toán ARG4WG là REARG và GAMARG nhằm tối ưu thêm số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Thuật toán REARG giúp quá trình xây dựng ARG khu trú được vào các ARG có số sự kiện tái tổ hợp nhỏ nhanh hơn ARG4WG trong hữu hạn số lần chạy thuật toán đối với các tập dữ liệu vừa và lớn. Tuy nhiên, thuật toán GAMARG tổng quát hơn. GAMARG có khả năng xây dựng được những ARG có chính xác hoặc gần chính xác số sự kiện tái tổ hợp nhỏ nhất. Trong thời gian tới, chúng tôi tiếp tục nghiên cứu các cách để xây dựng đồ thị ARG tối thiểu cho dữ liệu hệ gen. Hệ gen người có cấu trúc là các khối haplotype (haplotype blocks), mà trong các khối này không có sự kiện tái tổ hợp xảy ra và các khối được xác định bởi các sự kiện tái tổ hợp. Do đó, một hướng đi tiềm năng có thể xét đến là ứng dụng các thuật toán trong bài toán tìm các khối haplotype, từ đó tìm được các điểm nóng phân biệt các khối (hot spots) chính là các vị trí tiềm năng xảy ra các sự kiện tái tổ hợp, từ đó xây dựng đồ thị ARG dựa trên các vị trí tái tổ hợp đã được xác định đó. Ngoài ra, chúng tôi sẽ nghiên cứu kết hợp các phương pháp tối ưu tổ hợp vào thuật toán GAMARG để tối ưu số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Bên cạnh đó, chúng tôi sẽ tiếp tục nghiên cứu và triển khai các ứng dụng thuật toán ARG4WG, GAMARG vào các bài toán thực tế khác trên dữ liệu lớn như bài toán tìm đa hình di truyền đơn nucleotide, xử lý dữ liệu bị khuyết, 102 DANH MỤC CÁC CÔNG TRÌNH KHOA HỌC CỦA TÁC GIẢ LIÊN QUAN ĐẾN LUẬN ÁN 1. Nguyen, T. T. P., Le, V. S., Ho, H. B., & Si Le, Q. (2016), “Building ancestral recombination graphs for whole genomes”, IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 14(2), 478-483. 2. Nguyen, T. T. P., Le, V. S. (2017), "Building minimum recombination ancestral recombination graphs for whole genomes", The 4th NAFOSTED Conference on Information and Computer Science 2017, pp. 248-253. (IEEE conference). 3. Nguyen, T. T. P., Le, V. S. (2019), “A Hybrid Approach to Optimize the Number of Recombination in Ancestral Recombination Graphs”, In Proceedings of the 2019 9th International Conference on Bioscience, Biochemistry and Bioinformatics, pp. 36-42. (ACM conference). 103 TÀI LIỆU THAM KHẢO 1. Arenas M (2013), “The importance and application of the ancestral recombination graph,” Frontiers in genetics, Vol. 4, p.206. 2. Bafna V, Bansal V (2006), “Inference about recombination from haplotype data: lower bounds and recombination hotspots,” Journal of Computational Biology, Vol. 13(2), pp.501–521. 3. Band G, Le QS, Jostins L, Pirinen M, Kivinen K, Jallow M, et al. (2013), “Imputation-based meta-analysis of severe malaria in three African populations,” PLoS genetics, Vol. 9(5), p.e1003509. 4. Barsh GS, Copenhaver GP, Gibson G, Williams SM (2012), “Guidelines for genome-wide association studies,” PLoS genetics, Vol. 8(7), p.e1002812. 5. Browning SR, Browning BL (2011), “Haplotype phasing: existing methods and new developments,” Nature Reviews Genetics, Vol. 12(10), p.703. 6. Bush WS, Moore JH (2012), “Genome-wide association studies,” PLoS computational biology, Vol. 8(12), p.e1002822. 7. Cámara PG, Levine AJ, Rabadán R (2016), “Inference of ancestral recombination graphs through topological data analysis,” PLoS computational biology, Vol. 12(8), p.e1005071. 8. Cao M, Shi J, Wang J, Hong J, Cui B, Ning G (2015), “Analysis of Human Triallelic SNPs by Next-Generation Sequencing,” Annals of human genetics, Vol. 79(4), pp.275–281. 9. Carlsson G (2009), “Topology and data,” Bulletin of the American Mathematical Society, Vol. 46(2), pp.255–308. 10. Charlton ND, Carbone I, Tavantzis SM, Cubeta MA (2008), “Phylogenetic relatedness of the M2 double-stranded RNA in Rhizoctonia fungi,” 104 Mycologia, Vol. 100(4), pp.555–564. 11. Chen GK, Marjoram P, Wall JD (2009), “Fast and flexible simulation of DNA sequence data,” Genome research, Vol. 19(1), pp.136–142. 12. Consortium 1000 Genomes Project, others (2010), “A map of human genome variation from population-scale sequencing,” Nature, Vol. 467(7319), p.1061. 13. Consortium IH, others (2003), “The international HapMap project,” Nature, Vol. 426(6968), p.789. 14. Darwin C (2004), On the origin of species, 1859, Routledge. 15. Delaneau O, Zagury J-F, Marchini J (2012), “Improved whole-chromosome phasing for disease and population genetic studies,” Nature methods, Vol. 10(1), p.5. 16. Devlin B, Risch N (1995), “A comparison of linkage disequilibrium measures for fine-scale mapping,” Genomics, Vol. 29(2), pp.311–322. 17. Didelot X, Lawson D, Darling A, Falush D (2010), “Inference of homologous recombination in bacteria using whole-genome sequences,” Genetics, Vol. 186(4), pp.1435–1449. 18. Durrant C, Zondervan KT, Cardon LR, Hunt S, Deloukas P, Morris AP (2004), “Linkage disequilibrium mapping via cladistic analysis of single- nucleotide polymorphism haplotypes,” The American Journal of Human Genetics, Vol. 75(1), pp.35–43. 19. Fearnhead P, Harding RM, Schneider JA, Myers S, Donnelly P (2004), “Application of coalescent methods to reveal fine-scale rate variation and recombination hotspots,” Genetics, Vol. 167(4), pp.2067–2081. 20. Frazer KA, Murray SS, Schork NJ, Topol EJ (2009), “Human genetic variation and its contribution to complex traits,” Nature Reviews Genetics, Vol. 10(4), p.241. 105 21. Grelon M (2016), “Meiotic recombination mechanisms,” Comptes rendus biologies, Vol. 339(7–8), pp.247–251. 22. Griffiths AJF, Miller JH, Suzuki DT, others (2000), “Sources of variation,” An Introduction to Genetic Analysis 7th edition Available from: https://www.ncbi.nlm.nih.gov/books/NBK22012/,. 23. Griffiths RC, Marjoram P (1996), “Ancestral inference from samples of DNA sequences with recombination,” Journal of Computational Biology, Vol. 3(4), pp.479–502. 24. Gusfield D (1997), “Algorithms on stings, trees, and sequences: Computer science and computational biology,” Acm Sigact News, Vol. 28(4), pp.41–60. 25. Gusfield D (2005), “Optimal, efficient reconstruction of root-unknown phylogenetic networks with constrained and structured recombination,” Journal of Computer and System Sciences, Vol. 70(3), pp.381–398. 26. Gusfield D (2014), ReCombinatorics: the algorithmics of ancestral recombination graphs and explicit phylogenetic networks, MIT Press. 27. Gusfield D, Eddhu S, Langley C (2004), “Optimal, efficient reconstruction of phylogenetic networks with constrained recombination,” Journal of bioinformatics and computational biology, Vol. 2(01), pp.173–213. 28. Hein J, Schierup M, Wiuf C (2004), Gene genealogies, variation and evolution: a primer in coalescent theory, Oxford University Press, USA. 29. Heine K, Beskos A, Jasra A, Balding D, De Iorio M (2018), “Bridging trees for posterior inference on ancestral recombination graphs,” Proceedings of the Royal Society A, Vol. 474(2220), p.20180568. 30. Hejase HA, Dukler N, Siepel A (2020), “From Summary Statistics to Gene Trees: Methods for Inferring Positive Selection,” Trends in Genetics,. 31. Hobolth A, Christensen OF, Mailund T, Schierup MH (2007), “Genomic 106 relationships and speciation times of human, chimpanzee, and gorilla inferred from a coalescent hidden Markov model,” PLoS genetics, Vol. 3(2), p.e7. 32. Hubisz MJ, Williams AL, Siepel A (2019), “Mapping gene flow between ancient hominins through demography-aware inference of the ancestral recombination graph,” bioRxiv, p.687368. 33. Hubisz M, Siepel A (2020), “Inference of ancestral recombination graphs using ARGweaver,” In: Statistical Population Genomics, Humana, New York, NY, pp.231–266. 34. Hudson RR (1983), “Properties of a neutral allele model with intragenic recombination,” Theoretical population biology, Vol. 23(2), pp.183–201. 35. Hudson RR, Kaplan NL (1985), “Statistical properties of the number of recombination events in the history of a sample of DNA sequences,” Genetics, Vol. 111(1), pp.147–164. 36. Huson DH, Rupp R, Scornavacca C (2010), Phylogenetic networks: concepts, algorithms and applications, Cambridge University Press. 37. Huson DH, Scornavacca C (2011), “A survey of combinatorial methods for phylogenetic networks,” Genome biology and evolution, Vol. 3, pp.23–35. 38. Jenkins PA, Song YS, Brem RB (2012), “Genealogy-based methods for inference of historical recombination and gene flow and their application in Saccharomyces cerevisiae,” PloS one, Vol. 7(11), p.e46947. 39. Kingman JFC (1982), “The coalescent,” Stochastic processes and their applications, Vol. 13(3), pp.235–248. 40. Kitts A, Sherry S (2002), “The single nucleotide polymorphism database (dbSNP) of nucleotide sequence variation,” The NCBI Handbook McEntyre J, Ostell J, eds Bethesda, MD: US National Center for Biotechnology Information,. 107 41. Kreitman M (1983), “Nucleotide polymorphism at the alcohol dehydrogenase locus of Drosophila melanogaster,” Nature, Vol. 304(5925), p.412. 42. Lam F, Tarpine R, Istrail S (2010), “The imperfect ancestral recombination graph reconstruction problem: upper bounds for recombination and homoplasy,” Journal of Computational Biology, Vol. 17(6), pp.767–781. 43. Larribe F, Lessard S, Schork NJ (2002), “Gene mapping via the ancestral recombination graph,” Theoretical population biology, Vol. 62(2), pp.215– 229. 44. Le SQ, Durbin R (2011), “SNP detection and genotyping from low-coverage sequencing data on multiple diploid samples,” Genome research, Vol. 21(6), pp.952–960. 45. Li H, Durbin R (2011), “Inference of human population history from individual whole-genome sequences,” Nature, Vol. 475(7357), p.493. 46. Li Y, Chen W, Liu EY, Zhou Y-H (2013), “Single nucleotide polymorphism (SNP) detection and genotype calling from massively parallel sequencing (MPS) data,” Statistics in biosciences, Vol. 5(1), pp.3–25. 47. Lyngsø RB, Song YS, Hein J (2005), “Minimum recombination histories by branch and bound,” In: International Workshop on Algorithms in Bioinformatics, pp.239–250. 48. Mailund T, Dutheil JY, Hobolth A, Lunter G, Schierup MH (2011), “Estimating divergence time and ancestral effective population size of Bornean and Sumatran orangutan subspecies using a coalescent hidden Markov model,” PLoS genetics, Vol. 7(3), p.e1001319. 49. McCarthy JJ, Hilfiker R (2000), “The use of single-nucleotide polymorphism maps in pharmacogenomics,” Nature biotechnology, Vol. 18(5), p.505. 50. McVean GAT, Cardin NJ (2005), “Approximating the coalescent with 108 recombination,” Philosophical Transactions of the Royal Society B: Biological Sciences, Vol. 360(1459), pp.1387–1393. 51. Menelaou A, Marchini J (2012), “Genotype calling and phasing using next- generation sequencing reads and a haplotype scaffold,” Bioinformatics, Vol. 21(6), pp.84–91. 52. Minichiello MJ, Durbin R (2006), “Mapping trait loci by use of inferred ancestral recombination graphs,” The American Journal of Human Genetics, Vol. 79(5), pp.910–922. 53. Morris AP, Whittaker JC, Balding DJ (2002), “Fine-scale mapping of disease loci via shattered coalescent modeling of genealogies,” The American Journal of Human Genetics, Vol. 70(3), pp.686–707. 54. NIH (2007), “Understanding human genetic variation,” NIH Curriculum Supplement Series, ncbi nlm nih gov/books/NBK20363 (last accessed 13 October 2015),. 55. Parida L, Melé M, Calafell F, Bertranpetit J, Consortium G (2008), “Estimating the ancestral recombinations graph (ARG) as compatible networks of SNP patterns,” Journal of Computational Biology, Vol. 15(9), pp.1133–1153. 56. Pritchard JK, Przeworski M (2001), “Linkage disequilibrium in humans: models and data,” The American Journal of Human Genetics, Vol. 69(1), pp.1–14. 57. Rahim NG, Harismendy O, Topol EJ, Frazer KA (2008), “Genetic determinants of phenotypic diversity in humans,” Genome biology, Vol. 9(4), p.215. 58. Rasmussen MD, Hubisz MJ, Gronau I, Siepel A (2014), “Genome-wide inference of ancestral recombination graphs,” PLoS genetics, Vol. 10(5), 109 p.e1004342. 59. Salemi M, Vandamme A-M, Lemey P (2009), The phylogenetic handbook: a practical approach to phylogenetic analysis and hypothesis testing, Cambridge University Press. 60. Scally A, Durbin R (2012), “Revising the human mutation rate: implications for understanding human evolution,” Nature Reviews Genetics, Vol. 13(10), pp.745–753. 61. Song YS, Ding Z, Gusfield D, Langley CH, Wu Y (2007), “Algorithms to distinguish the role of gene-conversion from single-crossover recombination in the derivation of SNP sequences in populations,” Journal of Computational Biology, Vol. 14(10), pp.1273–1286. 62. Song YS, Hein J (2005), “Constructing minimal ancestral recombination graphs,” Journal of Computational Biology, Vol. 12(2), pp.147–169. 63. Song YS, Lyngso R, Hein J (2006), “Counting all possible ancestral configurations of sample sequences in population genetics,” IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), Vol. 3(3), pp.239–251. 64. Song YS, Wu Y, Gusfield D (2005), “Efficient computation of close lower and upper bounds on the minimum number of recombinations in biological sequence evolution,” Bioinformatics, Vol. 21(suppl_1), pp.i413--i422. 65. Than C, Ruths D, Nakhleh L (2008), “PhyloNet: a software package for analyzing and reconstructing reticulate evolutionary relationships,” BMC bioinformatics, Vol. 9(1), p.322. 66. Vaughan TG, Welch D, Drummond AJ, Biggs PJ, George T, French NP (2017), “Inferring ancestral recombination graphs from bacterial genomic data,” Genetics, Vol. 205(2), pp.857–870. 110 67. Visscher PM, Brown MA, McCarthy MI, Yang J (2012), “Five years of GWAS discovery,” The American Journal of Human Genetics, Vol. 90(1), pp.7–24. 68. Wall JD (2000), “A comparison of estimators of the population recombination rate,” Molecular Biology and Evolution, Vol. 17(1), pp.156– 163. 69. Wang L, Zhang K, Zhang L (2001), “Perfect phylogenetic networks with recombination,” Journal of Computational Biology, Vol. 8(1), pp.69–78. 70. Wiuf C, Hein J (1999), “Recombination as a point process along sequences,” Theoretical population biology, Vol. 55(3), pp.248–259. 71. Wu Y (2007), “Association mapping of complex diseases with ancestral recombination graphs: Models and efficient algorithms,” In: Annual International Conference on Research in Computational Molecular Biology, pp.488–502. 72. Wu Y, Gusfield D (2007), “Efficient computation of minimum recombination with genotypes (not haplotypes),” Journal of bioinformatics and computational biology, Vol. 5(02a), pp.181–200. 73. Zöllner S, Pritchard JK (2005), “Coalescent-based association mapping and fine mapping of complex trait loci,” Genetics, Vol. 169(2), pp.1071–1092.