Luận án Tiến sĩ Công nghệ thông tin: 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

Tập các trình tự D

Số lượng trình tự trong một tập các trình tự N

độ dài của trình tự m

Trình tự thứ x trong một tập các trình tự Sx

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:

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

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 ( )....................... 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

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/.

cặp nhiễm sắc thể (Hình 1.1) [54].

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

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.

tái tổ hợp được giới thiệu trong Mục 1.2.1.

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.

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.

Đoạn DNA biến thể có kích thước từ 1kb đến hơn 5Mb (1Mb = 106 bp).

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

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.

và 2 hoặc nhiều alen biến thể (Hình 1.4).

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.

Dữ liệu haplotype

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ể đó (https://www.nature.com/scitable/definition/genotype-

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à

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.

cây phân loài (phylogenetic tree) với cấu trúc như Hình 1.6 [14,59]:

• 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ư, …

Hình 1.7: Khái quát hóa các mạng phát sinh loài điển hình [36].

và là nền tảng nghiên cứu thuốc chữa bệnh [6].

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

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.

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).

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ự

Hình 1.11: Điểm cắt tái tổ hợp.

s (s’) và theo sau là phần sau (suffix) từ vị trí bx đến m của trình tự s’ (s).

− 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ị

trí thứ 3

Sự kiện 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

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).

(a) (b)

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

Bảng 1.1: Các phần mềm xây dựng đồ thị ARG tiêu biểu.

công khai cho cộng đồng nghiên cứ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.

41

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.

- Mã nguồn mở. - Yêu cầu cài đặt: Trình biên dịch ARGweaver [58] https://github.com/mdrasmus/ argweaver 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.

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

Hình 2.1: Lưu đồ thuật toán Margarita.

sang Bước 3.

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

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)).

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ể.

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, , để 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:

nếu s có thể kết hợp với một trình tự trong D.

nếu s có thể được biểu diễn thành trong đó và có thể kết

hợp với các trình tự trong D.

Ngoài ra, đặt trong đó và có thể kết hợp với các trình tự

trong D. có thể được viết dưới dạng:

(1)

51

Đặt trong đó 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 .

Chúng ta sẽ chứng minh rằng:

Đặt là phần bù của là giải pháp tối ưu của phương trình (1) và

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 (

).

trong , .

Rõ ràng là là một chuỗi con của vì là một chuỗi con của . Do đó, nếu

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

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

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ái tổ hợp (Kịch bản B).

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ỉ:

https://github.com/thaontp711/arg4wg.

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

Bảng 2.1: Tập dữ liệu trích xuất từ dự án 1000 hệ gen người.

được mô tả trong Bảng 2.1.

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

) y â

10000

1000

i g ( y ạ h c n a

100

(a)

i g i

10

ờ h T

1

1000

2000

5000

10000

Số vị trí SNP

1000000

100000

) y â

10000

1000

i g ( y ạ h c n a

100

(b)

i g i

10

ờ h T

1

1000

2000

5000

10000

Số vị trí SNP

1000000

100000

10000

) y â

1000

100

(c)

i g ( y ạ h c n a

10

i g i

1

ờ h T

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.

p ợ h ổ t i

(a)

á t n ệ i k ự s ố S

45000 40000 35000 30000 25000 20000 15000 10000 5000 0

1000

10000

2000 5000 Số vị trí SNP

p ợ h ổ t i

(b)

á t n ệ i k ự s ố S

80000 70000 60000 50000 40000 30000 20000 10000 0

1000

2000

5000

10000

Số vị trí SNP

64

p ợ h ổ t i

(c)

á t n ệ i k ự s ố S

100000 90000 80000 70000 60000 50000 40000 30000 20000 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

ARG4WG

56

Margarita

54

52

F R h c á c g n ả o h K

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

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.

Tóm lại, Margarita cũng có thể tìm ra các tín hiệu liên đới nhưng mất rất nhiều thời

gian để làm việc đó. Ví dụ trong bài toán này, Margarita mất đến 420 giờ (17.5

ngày) để sinh ra 10 ARG cho 5560 haplotype trên vùng 2Mbp. Như vậy, trong

trường hợp vùng liên quan đến bệnh là chưa biết trước thì việc chia nhỏ hệ gen

thành các vùng dữ liệu và thực hiện ánh xạ dùng ARG từ Margarita là việc làm khó

khả thi trong thực tế. Điều này đã nhấn mạnh thêm sức mạnh, khả năng ứng dụng

của ARG4WG trong các bài toán thực tế với dữ liệu lớn.

2.5. Kết luận chương

Bước tái tổ hợp là nút thắt chính trong thời gian chạy của tất cả các thuật toán xây

dựng đồ thị ARG vì tái tổ hợp phân tách 1 trình tự thành 2 trình tự con sẽ làm tăng

số lượng nút trong đồ thị. Thuật toán Margarita đã thực hiện tái tổ hợp trên cặp trình

tự có đoạn chung dài nhất, phân tách 1 trong 2 trình tự đó nhằm có được trình tự

con chứa đoạn chung để kết hợp với trình tự còn lại. Mục đích của thuật toán nhằm

giảm lượng thông tin di truyền nhiều nhất qua bước tái tổ hợp, hướng tới mục tiêu

nhanh hội tụ về nút tổ tiên chung gần nhất. Tuy nhiên, chiến lược này làm bùng nổ

số nút trong đồ thị khi đoạn chung dài nhất luôn được tìm thấy bên trong trình tự.

Do đó, 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ó

72

đoạn đầu chung dài nhất (đoạn chung dài nhất từ phía bên phải hoặc phía bên trái

của trình tự) cho bước tái tổ hợp. Với chiến lược này, sự kiện tái tổ hợp chỉ thay thế

1 trình tự bằng một trình tự con chứa ít thông tin di truyền hơn. Điều này đảm bảo

số nút trong đồ thị luôn ổn định sau các bước tái tổ hợp.

Các thực nghiệm trên các bộ dữ liệu khác nhau đã đưa đến kết luận rằng, thuật toán

ARG4WG tuy có hình thái cây kém hơn một chút nhưng có thời gian tính toán

nhanh gấp hàng trăm tới hàng nghìn lần so với 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

một lần chạy trong khoảng thời gian hợp lý.

Ngoài ra, Mệnh đề 1 đã được chúng tôi đưa ra và chứng minh về tính tối ưu số sự

kiện tái tổ hợp của chiến lược tìm đoạn đầu chung dài nhất trên phạm vi cục bộ mỗi

bước thực hiện sự kiện tái tổ hợp. Điều này đã mang lại hiệu quả và cho ra các đồ

thị ARG có số sự kiện tái tổ hợp ít hơn Margarita trên các tập dữ liệu thử nghiệm

lớn.

ARG4WG đã được thử nghiệm ứng dụng 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 Gambia gồm 2780 người (1533 người khỏe

mạnh, 1247 người bị bệnh) trên toàn bộ nhiễm sắc thể 11. Các kết quả thực nghiệm

đã cho thấy khả năng ứng dụng của thuật toán đề xuất trong việc phát hiện ra vùng

gen liên quan đến bệnh cho các nghiên cứu tương quan toàn bộ nhiễm sắc thể trên

các tập dữ liệu lớn. Từ đó, ý tưởng sử dụng thuật toán ARG4WG để có nhanh được

kết quả bước đầu, khu trú được các vùng gen tiềm năng là một việc làm khả thi. Các

nhà phân tích dữ liệu sau đó có thể tiến hành các phân tích sâu trên các vùng gen

này. Điều này sẽ giúp giảm đáng kể chi phí và thời gian phân tích. Từ các kết quả

này, chúng tôi tin tưởng rằng ARG4WG có thể được ứng dụng hiệu quả cho nhiều

bài toán thực tế khác như bài toán tìm đa hình di truyền đơn nucleotit, bài toán xử lý

dữ liệu bị khuyết, … trên các tập dữ liệu lớn.

73

Các kết quả nghiên cứu của chương này đã được công bố trong một bài báo trên tạp

chí quốc tế IEEE/ACM Transactions on Computational Biology and Bioinformatics

năm 2017 (công trình khoa học số 1).

74

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

3.1. Giới thiệu

Hầu hết các lịch sử tiến hóa được cho là chứa một số lượng nhỏ sự kiện tái tổ hợp

có thể phát hiện được. Hơn nữa, trong tất cả các số liệu thống kê mà ta muốn xác

định liên quan đến lịch sử tái tổ hợp, số sự kiện tái tổ hợp nhỏ nhất là một số liệu có

thể được xác định cụ thể, và về nguyên tắc là có thể tính toán được.

Ngoài các ứng dụng của đồ thị ARG tối thiểu và gần tối thiểu, ngay cả cận dưới tái

tổ hợp cũng được chứng minh là có ích. Chúng có thể được dùng để trả lời các câu

hỏi về tái tổ hợp như: tìm các điểm nóng tái tổ hợp (recombination hotspots) tiềm

năng trong hệ gen [2,19]; ước lượng tỉ lệ tái tổ hợp trong các trình tự quan sát [68]

và xác định các thông số di truyền quần thể khác nhau [35].

Vì vậy, việc xây dựng đồ thị ARG với hàm mục tiêu cực tiểu hóa số sự kiện tái tổ

hợp luôn là mục tiêu hướng tới nhằm phát hiện ra các điểm tái tổ hợp được xem là

quan trọng phục vụ cho các phân tích. Quá nhiều sự kiện tái tổ hợp được phát hiện

sẽ gây khó khăn cho các tính toán dựa trên nó.

Thuật toán ARG4WG đề xuất được giới thiệu trong chương 2 (Công trình khoa học

số 01) có thể xây dựng được đồ thị ARG cho dữ liệu lớn hàng nghìn trình tự trên

toàn hệ gen người. Tuy nhiên, thuật toán không được thiết kế nhằm tối ưu hóa số sự

kiện tái tổ hợp. Việc sử dụng đoạn đầu chung dài nhất đã được chứng minh là tối ưu

về số sự kiện tái tổ hợp trên phạm vi cục bộ khi phân tách 1 trình tự để có thể kết

hợp với các trình tự còn lại. Tuy nhiên, điều này là không đúng trên phạm vi toàn

cục. Mục 3.3 sẽ trình bày về hạn chế này của ARG4WG. Từ đó, luận án trình bày 2

phương pháp: (1) kết hợp các đặc trưng của dữ liệu và (2) kết hợp các kĩ thuật tối

75

ưu vào thuật toán ARG4WG nhằm tối ưu hóa thêm số sự kiện tái tổ hợp trong quá

trình xây dựng đồ thị ARG.

3.2. Một số định nghĩa và khái niệm sử dụng trong các thuật

toán

Dưới giả định các vị trí vô hạn, ta gọi 2 vị trí i và j là không tương thích nếu chúng

chứa tất cả 4 loại giao tử 00, 01, 10, 11 [62]. Sẽ có ít nhất một sự kiện tái tổ hợp

giữa 2 vị trí không tương thích i và j. 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.

Xét một tập dữ liệu D(5) = {S1, S2, S3, S4, S5} của 5 trình tự với 5 vị trí như dưới

1 2 3 4 5

đây:

S1 = 1 1 0 0 1

S2 = 0 1 0 1 0

S3 = 0 0 0 1 0

S4 = 1 0 0 0 0

S5 = 1 0 0 0 1

Trong tập dữ liệu D(5), có 3 cặp vị trí không tương thích: vị trí 1 và vị trí 2; vị trí 2

và vị trí 4; vị trí 2 và vị trí 5. Trong trường hợp này, ít nhất 2 sự kiện tái tổ hợp phải

xảy ra trong lịch sử tiến hóa của các trình tự. Hình 3.1 là một ARG tối thiểu với 2

sự kiện tái tổ hợp (sự kiện thứ 1 và sự kiện thứ 4) biểu diễn lịch sử tiến hóa của tập

dữ liệu D(5).

76

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.

Đặt D = {S1, S2, …, SN} là tập N trình tự nhị phân có độ dài m, Sx[i] có giá trị bằng 0

hoặc 1, 1 ≤ x ≤ N, 1 ≤ i ≤ m; gọi * là trạng thái không di truyền, tức là không mang

thông tin di truyền từ dữ liệu quan sát. Một số định nghĩa được sử dụng giống như

trong thuật toán ARG4WG như sau:

• Xét 1 vị trí i, Sx[i] khớp với Sy[i] nếu Sx[i] = Sy[i] hoặc Sx[i] = * hoặc Sy[i] = *.

77

• (Sx,Sy){d,l} là một cặp trình tự Sx và Sy có đoạn đầu chung với độ dài tối đa l từ

phía bên trái (d = left) hoặc từ phía bên phải (d = right).

• (Sx,Sy){d,l} tồn tại nếu và chỉ nếu có ít nhất một vị trí i trong phần chung thỏa

mãn Sx[i] = Sy[i]  *.

• Cặp (Sx,Sy) gọi là cặp có đoạn đầu chung dài nhất nếu cặp đó chứa phần vật liệu

di truyền chung dài nhất trong đoạn đầu chung.

Với 1 cặp có đoạn đầu chung (Sx,Sy){d,l}, theo chiến lược đoạn đầu chung dài nhất

thì điểm cắt tái tổ hợp được xác định giữa:

• l và l + 1 khi d = left và Sx[i] khớp với Sy[i] với mọi 1  i  l và Sx[l+1] 

Sy[l+1].

• l -1 và l khi d = right và Sx[i] khớp với Sy[i] với mọi l  i  m và Sx[l-1] 

Sy[l-1].

Cũng giống với ARG4WG, các thuật toán đề xuất dưới đây đều hoạt động ngược

chiều thời gian và có cùng giả định có nhiều nhất một đột biến xảy ra tại mỗi vị trí

trong suốt quá trình xây dựng đồ thị ARG.

3.3. Hạn chế của thuật toán ARG4WG

ARG4WG trước tiên thực hiện tất cả các sự kiện kết hợp và đột biến có thể xảy ra.

Sau đó, thuật toán tìm một cặp trình tự có đoạn đầu chung dài nhất, tức là, đoạn

trong đó có chứa phần vật liệu di truyền giống nhau dài nhất từ phía bên trái hoặc

bên phải của 2 trình tự. Một sự kiện tái tổ hợp được thực hiện trên trình tự có ít vật

liệu di truyền trong đoạn chung hơn để tách trình tự đó thành 2 trình tự con. Trình

tự con chứa đoạn đầu chung sẽ được kết hợp với trình tự còn lại ngay sau bước tái

tổ hợp.

Chiến lược đoạn đầu chung dài nhất giúp cho ARG4WG chạy được với dữ liệu lớn

gồm hàng nghìn trình tự độ dài hệ gen. Tuy nhiên thuật toán hướng đến việc xây

78

dựng đồ thị ARG hợp lý chứ không thể xây dựng đồ thị ARG tối thiểu. Hình 3.2

minh họa ngắn gọn cách thuật toán ARG4WG thực hiện với tập dữ liệu D(5). Đoạn

đầu chung dài nhất được xác định là giữa trình tự S4 và S5 (trong khung chữ nhật),

ARG4WG sẽ luôn thực hiện tái tổ hợp trên trình tự S4 hoặc S5 trước. Lựa chọn này

không làm phá vỡ bất kì cặp vị trí không tương thích nào và thuật toán luôn cần đến

3 sự kiện tái tổ hợp để xây dựng đồ thị ARG. Trong khi số sự kiện tái tổ hợp ít nhất

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.

để xây dựng đồ thị ARG cho tập dữ liệu D(5) là 2.

3.4. Thuật toán REARG

3.4.1. Động cơ nghiên cứu

Xuất phát từ các quan sát trong quá trình làm thực nghiệm, chúng tôi nhận thấy rằng

việc lựa chọn cặp trình tự có độ dài đoạn đầu chung dài nhất cho việc thực hiện tái

tổ hợp trong thuật toán ARG4WG thường là không duy nhất. Nói cách khác,

ARG4WG thường phải chọn ngẫu nhiên 1 cặp trình tự cho việc thực hiện tái tổ hợp

từ rất nhiều cặp có cùng độ dài đoạn đầu chung dài nhất. Ví dụ, trong Hình 3.3, có 3

79

cặp ứng cử viên (S1, S2), (S2, S3) và (S4, S5) có cùng độ dài đoạn đầu chung dài nhất.

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.

Khi đó ARG4WG sẽ chọn ngẫu nhiên một trong 3 cặp để thực hiện tái tổ hợp.

Các phân tích thực nghiệm cho thấy, bên cạnh tiêu chí độ dài đoạn đầu chung dài

nhất, các yếu tố khác như độ tương đồng của cặp trình tự được chọn hay độ dài của

trình tự được chọn để thực hiện tái tổ hợp cũng có ảnh hưởng đáng kể đến số sự

kiện tái tổ hợp. Do đó, việc kết hợp các yếu tố này trong việc lựa chọn cặp trình tự

thích hợp nhất cho việc tái tổ hợp có thể giúp định hướng quá trình xây dựng đồ thị

ARG tới đồ thị ARG với số sự kiện tái tổ hợp tối ưu hơn khi chạy với dữ liệu lớn

với số lần chạy giới hạn.

3.4.2. Thuật toán REARG

Trong thuật toán này, chúng tôi định nghĩa:

Độ tương đồng của 2 trình tự S1 và S2:

Với

80

Độ dài của một trình tự S:

Với

Lưu đồ thuật toán REARG giống với lưu đồ thuật toán ARG4WG và Margarita như

trong Hình 2.1. Các bước của thuật toán được trình bày trong Thuật toán 3.1.

Thuật toán REARG

Đầu vào: Một tập N trình tự, mỗi trình tự độ dài m.

vào, các trình tự trung gian được sinh ra và trình tự tổ tiên chung duy nhất được tìm thấy.

• Bước 1: Thực hiện tất cả các khả năng kết hợp.

• Bước 2: Thực hiện một đột biến sau đó quay lại bước 1.

Nếu không thực hiện được đột biến nào nữa thì chuyển sang bước 3.

• Bước 3: Thực hiện 1 sự kiện tái tổ hợp với chiến lược được mô tả như bên dưới,

theo sau là 1 sự kiện kết hợp rồi quay lại bước 1.

Nếu không thực hiện được tái tổ hợp thì chuyển sang bước 4.

• Bước 4: Kết thúc thuật toán, đạt đến một trình tự tổ tiên chung duy nhất.

Thuật toán 3.1: Thuật toán REARG.

Đầu ra: Đồ thị ARG chứa các sự kiện kết hợp, đột biến và tái tổ hợp giữa các trình tự đầu

Khác với ARG4WG, để thực hiện tái tổ hợp, nếu có nhiều hơn 1 cặp trình tự có

cùng độ dài đoạn đầu chung dài nhất thì REARG sẽ sử dụng thêm tiêu chuẩn về độ

tương đồng giữa các cặp trình tự và độ dài của trình tự để lựa chọn ứng cử viên

cho bước tái tổ hợp. Phần dưới đây mô tả 3 phiên bản khác nhau của thuật toán

REARG: REARG_SIM, REARG_LEN và REARG_COM tương ứng với 3 tiêu

chuẩn lựa chọn ứng cử viên cho bước tái tổ hợp.

81

Bước tái tổ hợp của thuật toán REARG_SIM

• Bước 1: Tính độ dài của các đoạn đầu chung cho tất cả các cặp trình tự. Các

cặp trình tự có cùng độ dài đoạn đầu chung dài nhất được chọn là các cặp ứng

cử viên cho việc tái tổ hợp.

• Bước 2: Tính độ tương đồng của tất cả các cặp ứng cử viên. Chọn cặp ứng cử

viên có độ tương đồng cao nhất để thực hiện tái tổ hợp. Trong trường hợp có

nhiều ứng cử viên có cùng độ tương đồng cao nhất, một cặp trong số đó sẽ

được chọn ngẫu nhiên để thực hiện tái tổ hợp.

Bước tái tổ hợp của thuật toán REARG_LEN

• Bước 1: Tính độ dài của các đoạn đầu chung cho tất cả các cặp trình tự. Các

cặp trình tự có cùng độ dài đoạn đầu chung dài nhất được chọn là các cặp ứng

cử viên cho việc tái tổ hợp.

• Bước 2: Tính độ dài của trình tự ngắn hơn của tất cả các cặp ứng cử viên.

Chọn ứng cử viên có độ dài trình tự dài nhất để thực hiện tái tổ hợp. Trong

trường hợp có nhiều ứng cử viên có cùng độ dài trình tự dài nhất, một trong

số đó sẽ được chọn ngẫu nhiên để thực hiện tái tổ hợp.

Bước tái tổ hợp của thuật toán REARG_COM

• Bước 1: Tính độ dài của các đoạn đầu chung cho tất cả các cặp trình tự. Các

cặp trình tự có cùng độ dài đoạn đầu chung dài nhất được chọn là các cặp ứng

cử viên cho việc tái tổ hợp.

• Bước 2: Tính độ tương đồng của tất cả các cặp ứng cử viên và tính độ dài của

trình tự ngắn hơn của các cặp ứng cử viên.

• Bước 3: Chọn ngẫu nhiên một cặp ứng cử viên có độ tương đồng cao nhất

hoặc một ứng cử viên có độ dài trình tự dài nhất để thực hiện tái tổ hợp.

82

3.5. Thuật toán GAMARG

3.5.1. Động cơ nghiên cứu

Do chiến lược đoạn đầu chung dài nhất không dẫn đến số sự kiện tái tổ hợp cực

tiểu, nên ý tưởng đặt ra là kết hợp ARG4WG với các tiêu chí tối ưu khác để giảm số

sự kiện tái tổ hợp. Đáng chú ý, kiểm thử 4 giao tử (four-gamete test) là ý tưởng then

chốt dẫn đến nhiều thuật toán khác nhau trong bài toán tìm cận dưới số sự kiện tái

tổ hợp và trong bài toán xây dựng đồ thị ARG có chính xác số sự kiện tái tổ hợp

nhỏ nhất. Do đó, luận án đề xuất thuật toán GAMARG kết hợp giữa ràng buộc trong

kiểm thử 4 giao tử với chiến lược đoạn đầu chung dài nhất trong ARG4WG để tối

ưu hóa số sự kiện tái tổ hợp trong quá trình xây dựng đồ thị ARG. Các kết quả thực

nghiệm trên các tập dữ liệu khác nhau cho thấy GAMARG có thể chạy được với

hàng nghìn trình tự với hàng chục nghìn SNP và có thể so sánh được với các thuật

toán xây dựng đồ thị ARG tối thiểu trên các bộ dữ liệu nhỏ.

3.5.2. Thuật toán GAMARG

Trước khi đi vào chi tiết thuật toán đề xuất, chúng tôi đưa ra một số quan sát trong

quá trình xây dựng ARG sử dụng kiểm thử 4 giao tử, từ đó dẫn đến một số mở rộng

được đề xuất khi áp dụng kiểm thử 4 giao tử vào thuật toán.

Đặt FreqGametei,j = {freq00i,j, freq01i,j, freq10i,j, freq11i,j} là tần số của các loại

giao tử 00, 01, 10, 11 xuất hiện giữa vị trí i và vị trí j.

Tần số của 4 loại giao tử giữa các vị trí không tương thích trong tập dữ liệu D(5) lần

lượt là FreqGamete1,2 = {1,1,2,1}; FreqGamete2,4 = {2,1,1,1}; FreqGamete2,5 =

{2,1,1,1}.

Ta nhận thấy rằng có 3 loại giao tử có tần số bằng 1. Vì vậy, việc lựa chọn trình tự

chứa một trong số các loại giao tử đó giữa các vị trí không tương thích để thực hiện

tái tổ hợp sẽ cho ta kết quả tốt hơn. Ví dụ, việc thực hiện một sự kiện tái tổ hợp giữa

vị trí 1 và vị trí 2 trên một trong 3 trình tự S1, S2, S3 sẽ phá vỡ ngay cặp vị trí ko

83

tương thích này. Ngược lại, vì freq101,2 = 2 nên nếu ta thực hiện một sự kiện tái tổ

hợp giữa vị trí 1 và vị trí 2 trên một trong 2 trình tự S4, S5 thì ta chỉ làm giảm bớt 1

lần tần số xuất hiện của loại giao tử 10 và ta cần thêm 1 sự kiện tái tổ hợp nữa để

phá vỡ được cặp vị trí không tương thích này. Cụ thể, xét ví dụ như trong Hình 3.4,

việc chọn S4 để thực hiện tái tổ hợp sẽ cần đến 2 sự kiện tái tổ hợp giữa vị trí 1 và vị

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).

trí 2 (R1,2(1) và R1,2(2)) để phá vỡ cặp vị trí không tương thích này.

Vì việc xây dựng ARG tối thiểu là một bài toán NP-khó nên không thể tìm được

ARG tối thiểu với các tập dữ liệu lớn. Để vượt qua gánh nặng tính toán, luận án đơn

giản hóa chiến lược kiểm tra 4 giao tử bằng cách chỉ xem xét các cặp vị trí không

tương thích có ít nhất một loại giao tử có tần số bằng 1. Giả định này đảm bảo rằng

thuật toán sẽ luôn phá vỡ ít nhất một cặp vị trí không tương thích khi thực hiện một

tái tổ hợp giữa một cặp vị trí không tương thích i và j.

Đặt ઠ là kích thước của cửa sổ trượt mà ta sẽ quét để tìm tất cả các cặp vị trí không

tương thích trong vùng đó. Cụ thể, tất cả các vị trí sẽ được kiểm tra. Với mỗi vị trí i

(0 ≤ i < m), chúng tôi sẽ quét để tìm tất cả các cặp vị trí không tương thích trong

phạm vi [i, i+ ઠ].

Đặt Sx(i,j) là một trình tự có loại giao tử có tần số bằng 1 giữa cặp vị trí không

𝑓𝑟𝑒𝑞00𝑖,𝑗 > 0 𝑎𝑛𝑑 𝑓𝑟𝑒𝑞01𝑖,𝑗 > 0 𝑎𝑛𝑑 𝑓𝑟𝑒𝑞10𝑖,𝑗 > 0 𝑎𝑛𝑑 𝑓𝑟𝑒𝑞11𝑖,𝑗 > 0 { 𝑓𝑟𝑒𝑞00𝑖,𝑗 = 1 𝑜𝑟 𝑓𝑟𝑒𝑞01𝑖,𝑗 = 1 𝑜𝑟 𝑓𝑟𝑒𝑞10𝑖,𝑗 = 1 𝑜𝑟 𝑓𝑟𝑒𝑞11𝑖,𝑗 = 1

tương thích i và j (0 ≤ i < m, j - i ≤ ઠ). Tức là, Sx(i,j) thỏa mãn điều kiện sau:

84

Cho một trình tự ứng cử viên Sx(i,j), ta cần tìm điểm cắt tái tổ hợp tốt nhất trong

phạm vi [i,j]. Vấn đề này lại được giải quyết sử dụng chiến lược đoạn đầu chung dài

nhất. Đoạn đầu chung dài nhất giữa trình tự này với tất cả các trình tự khác được

tính toán. Nếu tồn tại một trình tự Sz mà (Sx, Sz){d,l} thỏa mãn i ≤ l ≤ j, thì điểm cắt

tái tổ hợp được xác định là giữa l và l + 1 khi d = left hoặc giữa l - 1 và l khi d =

right. Nếu không tồn tại cặp có đoạn đầu chung nào trong phạm vi [i,j] thì điểm cắt

tái tổ hợp được chọn ngẫu nhiên giữa vị trí i và i+1 hoặc giữa vị trí j - 1 và j.

Hình 3.5: Lưu đồ thuật toán GAMARG.

Lưu đồ thuật toán được biểu diễn trong Hình 3.5.

Thuật toán GAMARG bắt đầu từ thời điểm t = 1. Tập các trình tự tại thời điểm t

được kí hiệu là Dt (D1=D). Với mỗi Dt, các danh sách 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:

85

• Danh sách kết hợp C: Đối với một cặp trình tự Sx và Sy có đoạn đầu chung

(Sx,Sy){d,l}, nếu l = m thì (Sx,Sy){d,l} được thêm 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 Sx[i] = 1 và ∀𝑆𝑦 ∈ 𝐷𝑡 ∖

{𝑆𝑥}: 𝑆𝑦[𝑖] ≠ 1 hoặc Sx[i] = 0 và ∀𝑆𝑦 ∈ 𝐷𝑡 ∖ {𝑆𝑥}: 𝑆𝑦[𝑖] ≠ 0, thì Sx[i] được thêm

vào danh sách đột biến.

• Danh sách giao tử G: Đối với một cặp vị trí không tương thích (i,j) (0 ≤ i < m,

j-i ≤ ઠ), nếu tồn tại một trình tự Sx chứa loại giao tử có tần số bằng 1 thì Sx(i,j)

được thêm vào danh sách giao tử.

• Danh sách đoạn đầu chung S: Với một cặp trình tự Sx và Sy có đoạn đầu chung

(Sx,Sy){d,l}, nếu 0 < l < m thì (Sx,Sy){d,l} được thêm vào danh sách đoạn đầu

chung.

Khi một trong 3 sự kiện xảy ra, tập trình tự tiếp theo Dt+1 được tạo ra từ tập trình tự

Dt hiện thời như mô tả dưới đây và 4 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ự Sx và Sy, hai trình tự Sx và Sy được

kết hợp thành 1 trình tự tổ tiên chung S’:

𝐷𝑡+1 = (𝐷𝑡\{𝑆𝑥, 𝑆𝑦}) ∪ {𝑆′}.

• Nếu một sự kiện đột biến xảy ra trên trình tự S, một trình tự mới S’ được tạo ra

từ trình tự S với đột biến:

𝐷𝑡+1 = (𝐷𝑡\{𝑆}) ∪ {𝑆′}.

• Nếu một sự kiện tái tổ hợp xảy ra trên một trình tự Sx(i,j), một điểm cắt tái tổ

hợp được thực hiện trong đoạn [i,j]. Hai trình tự con mới Sx1 và Sx2 được tạo ra

từ Sx:

𝐷𝑡+1 = (𝐷𝑡\{𝑆𝑥}) ∪ {𝑆𝑥1, 𝑆𝑥2}

• Nếu một tái tổ hợp xuất hiện trên một cặp trình tự có đoạn đầu chung (Sx,

Sy){d,l}, thực hiện tái tổ hợp trên trình tự có ít thông tin di truyền trong phần

86

đoạn chung hơn. Giả sử trình tự Sx có ít thông tin di truyền trong phần đoạn

chung hơn Sy, Sx sẽ được tách thành 2 trình tự con Sx1 và Sx2:

𝐷𝑡+1 = (𝐷𝑡\{𝑆𝑥}) ∪ {𝑆𝑥1, 𝑆𝑥2}

Các bước thực hiện thuật toán được trình bày trong Thuật toán 3.2.

Thuật toán GAMARG

Đầu vào: Một tập N trình tự nhị phân độ dài m.

Đầu ra: Một đồ thị ARG chứa các sự kiện kết hợp, đột biến, tái tổ hợp giữa các trình

tự, các trình tự trung gian được sinh ra và trình tự tổ tiên chung duy nhất được tìm thấy.

• Bước 1: Nếu danh sách kết hợp C không rỗng, thực hiện tất cả các kết hợp có thể.

• Bước 2: Nếu danh sách đột biến M không rỗng, thực hiện tất cả các đột biến có thể sau

đó chuyển sang Bước 1. Nếu không có đột biến nào, chuyển sang Bước 3.

• Bước 3: Nếu danh sách giao tử G không rỗng, thực hiện một tái tổ hợp sau đó chuyển

sang Bước 1. Nếu G rỗng, chuyển sang Bước 4.

• Bước 4: Nếu danh sách đoạn đầu chung S không rỗng, thực hiện một tái tổ hợp theo sau

là một sự kiện kết hợp. Chuyển đến Bước 1. Nếu S rỗng, chuyển sang Bước 5.

• Bước 5: Kết thúc thuật toán, đạt đến một trình tự tổ tiên chung duy nhất.

Thuật toán 3.2: Thuật toán GAMARG.

Các ứng cử viên từ 4 danh sách được chọn như sau:

• Ứng cử viên từ danh sách kết hợp và đột biến để thực hiện kết hợp và đột biến

được lấy ngẫu nhiên.

• Trong danh sách Giao tử G, nếu một trình tự ứng cử viên Sx(i, j) có khoảng cách

ngắn nhất từ vị trí i đến vị trí j, tức là, (j – i) có giá trị nhỏ nhất thì Sx có thứ tự

ưu tiên hàng đầu để thực hiện tái tổ hợp. Nếu có nhiều hơn một ứng cử viên có

cùng khoảng cách ngắn nhất thì một trong số ứng cử viên đó sẽ được chọn ngẫu

nhiên.

87

• Trong danh sách có đoạn đầu chung S, cặp trình tự có đoạn đầu chung dài nhất

(có đoạn thông tin di truyền chung nhiều nhất trong đoạn đầu chung) sẽ được

chọn thực hiện tái tổ hợp đầu tiên. Nếu có nhiều hơn một ứng cử viên có cùng

đoạn đầu chung dài nhất thì một trong số đó sẽ được chọn ngẫu nhiên.

Các lựa chọn ngẫu nhiên trong thuật toán GAMARG dẫn đến các ARG khác nhau

cho các lần chạy khác nhau.

3.6. Kết quả thực nghiệm

Để đánh giá hiệu suất của REARG và GAMARG, các thử nghiệm trên các bộ dữ

liệu khác nhau đã được tiến hành. Đầu tiên, chúng tôi so sánh các thuật toán

REARG, GAMARG, Margarita, và ARG4WG với kết quả từ các thuật toán xây

dựng đồ thị ARG tối thiểu (MinARG) trên các tập dữ liệu nhỏ gồm:

- Tập dữ liệu của Kreitman [41] bao gồm 11 trình tự độ dài 43. Tập dữ liệu

nhỏ này là một bộ dữ liệu chuẩn được sử dụng để đánh giá hiệu suất của

nhiều thuật toán tìm cận dưới tái tổ hợp cũng như các thuật toán xây dựng

ARG tối thiểu [47,62].

- Hai tập dữ liệu mô phỏng: SDS1 bao gồm 50 trình tự độ dài 54 và SDS2

gồm 75 trình tự độ dài 45 được công khai tại địa chỉ

https://people.eecs.berkeley.edu/~yss/lu.html.

Thứ 2, bốn thuật toán REARG, GAMARG, Margarita và ARG4WG được thử

nghiệm và so sánh trên 18 bộ dữ liệu trích xuất từ 3 vùng dữ liệu khác nhau của

Nhiễm sắc thể 1 từ dự án 1kGP [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

88

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

Bảng 3.1: Tập dữ liệu từ dự án 1kGP.

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.

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

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ỏ.

được mô tả trong Bảng 3.2.

Kreitman

SDS1

SDS2

MinARG

7

10

12

Margarita

8

14

18

ARG4WG

10

17

18

REARG

10

17

20

GAMARG

7

10

13

89

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

90

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

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.

đều có kết quả vượt trội trong tất cả các thực nghiệm.

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

2441

5610

10052

GAMARG

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

1824

4093

8298

GAMARG

91

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

1502

3994

8817

GAMARG

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

4099

9315

16715

GAMARG

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

2955

6679

13606

GAMARG

92

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

2521

6504

14583

GAMARG

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.

93

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

10600

2680

5700

10300

5400

10000

2380

9700

5100

2080

9400

4800

9100

REARG_LEN, REARG_COM, và REARG_SIM trên mỗi thử nghiệm.

100 trình tự

4500

1780

8800

4200

8500

1480

3900

8200

DS1

DS2

DS3

DS1

DS2

DS3

DS1

DS2

DS3

2000 SNP

5000 SNP

10000 SNP

4300

4000

3700

3400

200 trình tự

3100

2800

2500

9400 9100 8800 8500 8200 7900 7600 7300 7000 6700 6400

17100 16800 16500 16200 15900 15600 15300 15000 14700 14400 14100 13800 13500

DS1

DS2

DS3

DS1

DS2

DS3

DS1

DS2

DS3

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

94

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

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.

REARG_SIM và REARG_COM.

Số SNP

2000

5000

10000

442.00

5690.00

201.50

Margarita

0.57

1.26

2.08

ARG4WG

(a)

3.27

17.66

60.64

REARG_SIM

0.62

1.35

2.46

REARG_LEN

3.39

18.34

62.80

REARG_COM

4.15

10.18

30.83

GAMARG

95

Số SNP

2000

5000

10000

754.00

7709.00

140.00

Margarita

0.41

0.99

1.90

ARG4WG

(b)

2.41

12.98

53.50

REARG_SIM

0.44

1.06

2.15

REARG_LEN

2.56

13.77

55.88

REARG_COM

2.41

6.83

21.09

GAMARG

Số SNP

2000

5000

10000

743.50

7627.50

166.00

Margarita

0.36

0.96

1.96

ARG4WG

(c)

2.04

13.01

53.63

REARG_SIM

0.39

1.02

2.08

REARG_LEN

2.22

13.34

55.72

REARG_COM

1.65

6.30

23.99

GAMARG

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.

96

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

2826.51 31459.54

1157.00

Margarita

1.33

3.26

5.50

ARG4WG

(a)

8.84

49.85

168.09

REARG_SIM

1.43

3.51

5.76

REARG_LEN

9.71

49.20

169.77

REARG_COM

14.91

34.42

105.01

GAMARG

Số SNP

2000

5000

10000

3953.000 42623.510

839.503

Margarita

0.998

2.442

4.739

ARG4WG

(b)

6.181

35.391

145.266

REARG_SIM

0.999

2.558

4.879

REARG_LEN

7.001

35.246

146.843

REARG_COM

5.886

20.043

59.926

GAMARG

97

Số SNP

2000

5000

10000

3710.50 39230.55

1074.52

Margarita

0.96

2.42

4.94

ARG4WG

(c)

5.55

34.77

149.67

REARG_SIM

0.95

2.55

5.45

REARG_LEN

6.28

34.46

153.76

REARG_COM

3.95

18.62

60.35

GAMARG

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à

98

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).

99

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à để

100

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, …

101

DANH MỤC CÁC CÔNG TRÌNH KHOA HỌC CỦA TÁC

GIẢ LIÊN QUAN ĐẾN LUẬN ÁN

for whole genomes”, recombination graphs

1. Nguyen, T. T. P., Le, V. S., Ho, H. B., & Si Le, Q. (2016), “Building IEEE/ACM ancestral 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 in Ancestral Recombination Graphs”, In Number of Recombination Proceedings of the 2019 9th International Conference on Bioscience, Biochemistry and Bioinformatics, pp. 36-42. (ACM conference).

102

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,”

103

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.

104

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

105

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,.

106

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

107

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, http://www 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),

108

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.

109

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.

110