Multiple Global Alignment in C-Sibelia
C-Sibelia представляет собой инструмент для полногеномного выравнивания. На вход он принимает два генома, которые разбиваются на гомологичные блоки, и затем происходит глобальное выравнивание полученных блоков друг на друга. Было принято решение адаптировать существующий алгоритм для выравнивания множества последовательностей. При этом алгоритм разбиения на блоки успешно масштабируется на случай множества последовательностей, но дальнейшее выравнивание множества блоков является является сложной вычислительной задачей. Для решения задачи множественного глобального выравнивания существуют инструменты, однако время их работы оказалось неприемлемо долгим.
Таким образом, было необходимо:
- изучить существующие подходы к решению этой задачи, сравнить инструменты, предназначенные для глобального выравнивания;
- разработать более эффективный алгоритм для решения этой задачи;
- интегрировать полученный алгоритм в C-Sibelia
В рамках работы были изучены следующие инструменты: T-Coffie, Dialign-TX, MLagan, Clustal-Omega, MAFFT, MSARC. Был проведен замер производительности этих инструментов (указаны суммарный размер последовательностей, их количество, время в секундах; * означает, что замер неточен, поскольку программа выходила за рамки доступной физической памяти):
size(KB) |
SeqCnt |
mafft_def |
mlagan |
clustal-o |
mafft_ginsi |
t-coffee |
msarc |
172 |
29 |
1.092 |
50.171 |
126.246 |
123.928 |
624.200* |
1595* |
96 |
15 |
0.705 |
15.179 |
75.804 |
45.205 |
529.075* |
458* |
156 |
13 |
3.01 |
14.405 |
131.792 |
135.651 |
||
288 |
8 |
30.138 |
11.02 |
689* |
1381* |
||
372 |
8 |
36.314 |
13.058 |
1288* |
|||
148 |
8 |
1.573 |
7.516 |
154.678 |
|||
376 |
8 |
36.818 |
10.811 |
1375* |
|||
388 |
8 |
34.409 |
13.638 |
1410* |
|||
984 |
8 |
95.844 |
57.44 |
|
|||
1008 |
8 |
98.017 |
24.382 |
|
|||
1884 |
8 |
245.822 |
88.709 |
|
|||
236 |
8 |
4.897 |
9.384 |
387.911 |
|||
1260 |
8 |
149.253 |
126.594 |
|
|||
1244 |
8 |
150.605 |
86.158 |
|
|||
452 |
8 |
46.323 |
24.099 |
|
|||
100 |
8 |
0.631 |
5.911 |
69.049 |
|||
84 |
8 |
0.528 |
8.075 |
46.618 |
|||
176 |
8 |
1.366 |
13.355 |
216.709 |
|||
272 |
8 |
27.166 |
14.237 |
524* |
|||
1912 |
8 |
241.676 |
96.617 |
|
Как видно из таблицы, многие инструменты не могут обработать большую часть блоков, а обработка некоторых из них занимает много времени. Основная пробелема инструментов заключается в том, что они преднозначены для работы с произвольными геномами, однако в данном случае выравниваемые последовательности являются близкими, и для них могут быть эффективными совершенно иные подходы. И в этом предположении был разработан новый алгоритм.
Основная идея заключается в том, чтобы при выравнивании двух геномов находить для каждого k-мера из первого генома похожие на него во втором. Далее выравниваются похожие k-меры, промежутки между ними обрабатываются при помощи стандартного динамического программирования. Были сделаны следующие оптимизации:
- Каждый k-мер может выравниваться лишь на те, кто находится относительно близко от него. То есть выражение | p1 / len1 - p2 / len2| должно быть меньше определенного значения. Иными словами, относительные позиции двух k-меров в геномах должны быть близки.
- Посколько геномы могут содержать большие вставки и удаления, относительная позиция считалась лишь среди тех k-меров, для которых существуют похожие в другом геноме. То есть в таком примере: ABC AC. В последовательности ABC при расчете относительной позиции будет выбрасываться часть B.
- Для сравнения k-меров использовалось битовое сжатие (два бита на каждый нуклеотид), что позволяло быстро находить расстояние между двумя k-мерами.
- Был зафиксирован один нуклеотид (A), и для него для каждого генома был построен следующий индекс. Берется k-мер, в нем смотрится, на каких позициях находится A, и в зависимости от этого получается k-битная последовательностей из нулей и единиц. Это позволяет быстро находить возможные позиции для близких k-меров (с незначительным числом промахов), и затем определить их действительную степень близости.
Алгоритм был протестирован на приведенных выше данных, и были получены следующие результаты. На тесте, состоящем из 8 последовательностей длиной примерно 230000 нуклеотидов, MLagan работал 88 секунд, MAFFT - 245 секунд, а разработанный алгоритм - 1.87 секунд. Прочие алгоритмы были не способны обработать данный набор.
Качество работы было незначительно хуже: если у MLagan отношение колонок, в которых было хотя бы одно отличие, к общему числу колонок, было равно 55810 / 258355, то у разработанного алгоритма 56342 / 260357. При подсчете общего числа различных нуклеотидов, расположенных в одной колонке, получается 444268 и 465974 соответственно.
В дальнейшем планируется улучшить алгоритм прогрессивного выравнивания, поддержать обработку N в записи генома, интегрировать алгоритм в C-Sibelia и сравнить инструмент с прочими полногеномными выравнивателями.