CC BY
410
УДК 004.05
А.В. Григорьев ВЫРАВНИВАНИЕ ДВУХ ДНК ПОСЛЕДОВАТЕЛЬНОСТЕЙ, ОПТИМАЛЬНОЕ ДЛЯ РЕАЛИЗАЦИИ НА ГРАФИЧЕСКОМ ПРОЦЕССОРЕ
В работе рассмотрены подходы к решению задачи о нахождении наибольшей общей подпоследовательности структуры ДНК. Приведён обзор известных алгоритмов выравнивания последовательностей, выделаются основные достоинства и недостатки разных подходов. Сформулированы требования к алгоритму, наиболее полно использующему существующее современное оборудование. Предложен алгоритм, удовлетворяющий этим требованиям.
Генетика, ДНК, выравнивание последовательностей, наибольшая общая подпоследовательность, многопоточность, CUDA
A.V. Grigoryev ALIGNMENT OF TWO DNA SEQUENCES OPTIMIZED FOR GPU CALCULATIONS
Methods for solve DNA longest common subsequence problem are present.
Known methods of alignments are discussed and demonstrated. Dynamic programming techniques (using memoization through tabular computation) are discussed. Requirements for algorithm, witch full quality using a present hardware, are present. New method of alignments is presented and applied to the implementations.
Genetics, DNA sequence alignment, multithreading, longest common subsequence, CUDA
Введение
Г енная инженерия в настоящее время весьма востребована, одной из важных ее задач является определение степени схожести последовательностей ДНК. Определение этого соответствия в генетике называется выравниванием ДНК последовательностей.
Задача нахождения наибольшей общей подпоследовательности - поиск последовательности, которая является подпоследовательностью нескольких последовательностей. Подпоследовательность можно получить из некоторой конечной последовательности, если удалить из последней некоторое множество её элементов (возможно, пустое). Спектр задач, в которых требуется выполнять выравнивание, весьма обширен: реконструкция эволюции организмов, поиск активного центра, моделирование 3D структуры белка, выявление паттерна функциональных семейств и сигналов в ДНК. В последнее время задача выравнивания находит широкое применение в проблеме секвенирования ДНК.
Особое распространение получило множественное выравнивание - определение наибольшей общей подпоследовательности среди многих последовательностей. Графически это можно представить как несколько гомологических последовательностей, записанных друг под другом оптимальным способом (рис. 1, [1]). Оптимальным считается расположение в одном
285
столбце аминокислот, которые находятся в одинаковом пространственном положении, имеют одинаковую функциональную нагрузку или похожи. В некоторых случаях возможны ситуации, когда выравниваемые последовательности имеют разную длину. В таком случае меньшая последовательность дополняется так называемыми разрывами (на рис. 1 обозначены как «-») -элементами, не несущими функциональной нагрузки в последовательности. Указанные разрывы вводятся условно, для удобства представления результатов выравнивания.
<1> 28-68 <4> 333-373 <8> 154-193 <б> 615-647 <4> 502-533 <1> 844-872 <8> 194-220 <5> 451-480
AAC AAGCА-А-ACTTTTATССATGGTCGTGGTACAGAGGGGTС AACAAGCA-A-ACTTTT AT CCATGGTCGTGGT ACAGAGGGGTC AACAAGCА-А-ACTTTTАТССATGGTCGTGGTACAGAGGGGT-
AACAAGCAGA-ACTTTTАТССATGGTCGTGGTAC----------
AACAAGCA-ACCCTTTTATCCATGGTCGTGGTA-----------
AACAAGCA-A-ACTTTTATCCATGGTCGTGG-------------
---------A-ACTTTTATCCATGGTCGTGGTACAGA-------
------------CTTTСA-ACGTGGTCGTGGTACAAAGGGGTС
Рис. 1. Пример представления результатов множественного выравнивания.
В строках располагаются ДНК последовательности. Слева в таблице указаны условные названия организмов, справа сама последовательность, где латинскими буквами обозначаются белковые основания, символ «-» обозначает разрыв
Для решения задачи выравнивания предложены алгоритмы Ниделмана-Вунша, Мил-лера-Маерса и их модификации. Рассмотрим их более подробно.
Алгоритм Ниделмана-Вунша
В алгоритме Ниделмана-Вунша [2] используется подход динамического программирования. Требуется найти максимальное вхождение одной последовательности ДНК в другую, длины последовательностей составляют п1 и п2, соответственно. Пусть существуют решения для всех подзадач (т1, m2), меньших заданной, где m1 и m2 - длины последовательностей такие, что 0 < m1 < п1 и 0 < m2 < п2 . Тогда задача ( п1 , п2 ) сводится к меньшим подзадачам следующим образом:
substitute(n1, n2) =
0, n1 = 0 v n2 = 0,
substitute (n1 -1, n2 -1) +1, s1[n1] = s2[n2], (1)
max(substitute(n1 -1, n2),substitute(n1, n2 -1)), s1[n1] Ф s2[n2].
Графически это можно изобразить в виде матрицы М (рис. 2), включающей в виде элементов i, j решения подзадач (i, j). Стрелками на схеме указаны переходы между решениями подзадач. Фактурными стрелками обозначены переходы, при которых увеличивается длина общей подпоследовательности, белыми стрелками - при которых образуются или удлиняются разрывы одной из выравниваемых последовательностей. В каждой ячейке di,j матрицы М хранится длина наибольшей подпоследовательности, которую можно получить из первых i элементов первой и j элементов второй последовательностей. Метод Ниделмана
- Вунша позволяет учитывать разрывы и сводить их к минимуму через введение функции штрафов gap_penalty(gj), которая задает размер штрафа за вставленный разрыв, где g - массив всех разрывов.
Эта функция вычисляется для разрывов в обеих строках выравнивания следующим образом:
gap _penalty(gj) = costbegin + length(gj)*costextention. (2)
Здесь costbegin - штраф за инициализацию разрыва, costextention - штраф за удлинение на один символ, length(gj) - длина разрыва gj.
Рис. 2. Фрагмент матрицы М
Рис. 3. Матрица М. Отмечена полоса шириной W, внутри которой происходят вычисления
С учетом штрафа за разрыв (1) принимает вид:
score(dt-1 j-1) + substitute(di j), score(dt j) = max<! score(di_1 j) - gap _ penalty(g), (3)
score(dt j-1) - gap _ penalty(g),
Далее выполняется прямой ход алгоритма - по формуле (3) рекуррентно заполняется матрица score с первых элементов до последних. Таким образом, в каждой ячейке i, j матрицы М записывается самое лучшее решение подзадачи (mi, mj). Затем производят обратный ход алгоритма: матрица М проходится от элемента n1, n2 к элементу 0, 0. При обратном проходе восстанавливается искомое выравнивание.
Асимптотическая оценка времени работы алгоритма будет О(пх*п2) [3], где n1 и n2 -длины последовательностей ДНК. Так как информация записывается в матрицу, то данный алгоритм будет занимать в памяти объем, равный ni*n2.
Оптимизация памяти. Алгоритм Миллера-Маерса
При полном выравнивании ДНК, рассматриваются последовательности, длина которых исчисляется в миллионах аминокислот. Для работы алгоритма Ниделмана - Вунша потребуются колоссальные объемы памяти. Наиболее простой способ уменьшить требуюмую память - хранить в памяти не всю матрицу М, а только часть, в которой возможно найти наибольшую подпоследовательность (рис. 3). Наибольшая из всех возможных подпоследовательностей будет лежать на главной диагонали. Наиболее длинные подпоследовательности (А, на рис. 3 графически они выражаются наиболее короткими линиями) будут примыкать к главной диагонали. Подпоследовательности D1 и D2 (на рис. 3) не представляют интереса, т.к. содержат, в основном, разрывы последовательности, известные в генетике как делеции [4]. Поэтому матрицу М приводят к ленточной матрице и в памяти ЭВМ хранятся только элементы, примыкающие к главной диагонали. Следует помнить, что данный метод не гарантирует нахождение наилучшего решения, т.к. он учитывают не все подзадачи.
Наиболее экономный по потребляемой памяти (будет занимать в памяти объем, равный ni) - алгоритм Миллера-Маерса [5]. Принцип его работы состоит в разбиении одной из последовательностей на две равные части (рис. 4 слева). Для каждой точки x линии раздела находим веса по формуле (3) оптимальных выравниваний из начала (0, 0) в x и из конца (n1, n2) в x: W(x), W(x). Вес оптимального выравнивания, проходящего через точку x, равен: W(x) = W+(x) + W(x). Вес оптимального выравнивания равен W = maxx (W(x)). Таким образом, нашли точку, через которую проходит оптимальное выравнивание.
Рис. 4. Первый (слева) шаг для обсчета квадранта матрицы М по методу Миллера-Маерса; справа приведен конечный результат применения алгоритма Миллера-Маерса
Найденная точка x разбивает матрицу М на четыре квадранта (подматрицы матрицы М), два из которых заведомо не содержат оптимального выравнивания. Для двух квадрантов, содержащих оптимальный путь, можно применить аналогичный прием, запомнить точки x' и x". Продолжая процедуру деления квадрантов пополам, найдем все точки, через которые проходит оптимальный путь (рис. 4, справа).
У этого метода есть три существенных достоинства:
- все полученные квадранты можно обсчитывать независимо;
- каждый квадрант может обсчитываться двумя независимыми потоками;
- метод требует минимум памяти - щ*2 вместо пх*п2.
Эти преимущества делают его идеальным алгоритмом для выполнения в многопоточной среде.
Формулировка альтернативного алгоритма и его реализация
Приведенная выше оптимизация памяти требует многократного расчета одних и тех же элементов матрицы М. В связи с этим разработан метод, позволяющий объединить преимущества однократного прохода по матрице (в отличие от алгоритма Миллера-Маерса) и экономное использование памяти ЭВМ. В отличие от метода вычисления по ленточной матрице, предлагаемый обеспечивает большую достоверность решения. Средний расход памяти составит и*^(и) или меньше, в худшем случае расход памяти и2. Изложим суть метода. Рассмотрим произвольный столбец матрицы М, в котором записываются наилучшие длины подстрок. Имеем неубывающую последовательность, состоящую преимущественно из повторяющихся значений. Не обязательно хранить в памяти все значения, достаточно запомнить позицию, на которой улучшен результат. Таким образом, ряд {1, 1, 1, 1, 4, 4, 4, 6, 6, 8, 8, 8, 8} будет записан как {[0, 1], [4, 4], [7, 6], [9, 8]}, где первый элемент - позиция, второй - значение. Из (3) следует, что длина наибольшей подпоследовательности будет увеличиваться в случае, когда аминокислоты в позициях I первой последовательности и ] второй последовательности равны. Для оптимизации поиска нужных аминокислот, запишем вторую входную последовательность в виде множества списков. Каждый список будет содержать позиции определенной аминокислоты. Пример фрагмента программы, написанной на языке С++, решающий задачу выравнивания, приведен ниже:
//Формирование первой последовательности: vector<byte> sq1;
std::ifstream in(FileName.c_str()); char ch; in >> ch;
while (!in. eof ()){ switch (ch){
case
case
case
case
'A'
'T'
'G'
'C'
case
case
case
case
sq1.push_back(0); break; sq1.push_back(3); break; sq1.push_back(2); break; sq1.push_back(1); break;
}
in
>> ch;
}
in.close () ; //Процедура чтения unsigned int i=0; while (!in.eof())
второй последовательности выглядит следующим образом:
//Основной алгоритм:
vector<unsigned int> *vPos = new vector<unsigned int> [csq1],
*vVal = new vector<unsigned int> [csq1]; maxInCol = new unsigned int [n]; maxInRect = new unsigned int [n]; for (int i=0; i<n; i++){ maxInCol[i] = 0; maxInRect[i] = 0;
if (px!=0)
lev.push(px); if (py!=0)
lev.push(-py);
Теперь переменная lev содержит информацию о конечном выравнивании двух входных последовательностей:
- 0 - если необходимо записать в результат нуклеотиды и первой и второй последовательности,
- положительное значение, если необходимо записать нуклеотиды второй последовательности и писать пропуски вместо первой,
- отрицательное значение, если необходимо записать пропуски вместо второй последовательности и записать нуклеотиды из первой.
Для ненулевых элементов модуль этого значения показывает количество символов последовательности, которые надо записать.
Алгоритм для CUDA процессора
Описанный выше метод хорошо подходит для выполнения на многопроцессорном суперкомпьютере, однако момент для параллельных вычислений чаще применяют технологии вычисления на графических процессорах. Современные видеокарты почти не уступают суперкомпьютерам по скорости и количеству ядер, у них имеется один существенный недостаток: в графическом процессоре арифметические устройства группируются в блоки и все устройства, находящиеся в одном блоке, выполняют одну операцию [6]. Этот подход идеально применим для компьютерной графики, где требуется обработать массив точек или выполнить операции над матрицами.
Метод Миллера-Маерса будет неоптимальным для выполнения на видеокарте. Можно запустить в 2 потока расчет одного квадранта: первый будет выполнять вычисления с начала, другой с конца матрицы М и все операции для расчетов будут идентичными. При параллельном вычислении нескольких квадрантов, сталкиваемся со следующей проблемой: имеются квадранты обычно разной размерности, и когда один квадрант вычислен, то потоки будут
вынуждены ждать обработки оставшихся квадрантов. Это понижает производительность всего алгоритма в целом. Необходимо распараллелить вычисление внутри одного квадранта. Это можно сделать, изменив порядок вычисления каждого элемента. Обычно предполагается последовательный алгоритм вычисления матрицы М [7]. Последовательно обходим матрицу по строчкам сверху вниз и по столбцам слева направо. В предлагаемом алгоритме будем осуществлять в матрице М независимое вычисление элементов, стоящих на линии, параллельной побочной диагонали матрицы. Указанный порядок представлен на рис. 5 справа. Из него видно, что классический метод (представленный слева) выполняется за 16 этапов, а предложенный в альтернативном алгоритме - за 7.
1 2 3 4 1 2 3 4
I 5 6 7 8 2 3 4 5
9 10 11 12 3 I 4 5 6
13 14 15 16 4 5 6 7
Рис. б. Метод Нидельмана-Вунша представлен слева, альтернативный - справа
Заключение
Благодаря дешевизне и быстрому качественному росту видеокарт с их помощью начинают решать все больше задач. Структура процессора и схема управления в видеокарте отличается от обычного центрального процессора компьютера и это накладывает особые требования на алгоритм. Предложенный в работе подход позволяет наиболее полно использовать ресурсы видеокарты, чтобы повысить эффективность реализации известных и зарекомендовавших алгоритмов выравнивания ДНК последовательностей имеющих большое значение в современной молекулярной генетике.
ЛИТЕРАТУРА
1. URL: http://www.pnas.Org/content/102/5/1285/F1.large.jpg
2. Needleman S.B., Wunsch C.D., A general method applicable to the search for similarities in the amino acid sequence of two proteins // Journal of Molecular Biology. 48. 1970. P.443-453.
3. Грин Д., Кнут Д. Математические методы анализа алгоритмов. М.: Мир, 1987. 120 с.
4. URL: http://www.biology-online.org/dictionary/Deletion_mutation
5. Mayers E.W., Miller W. Optimal alignments in linear space // Computer Applications in the Biosciences. 4(1). 1988. P. 11-17.
6. Sanders J., Kandrot E. CUDA by Example: An Introduction to General-Purpose GPU Programming. Addison-Wesley Professional, 2010. 312 p.
7. URL: http://www.ibm.com/developerworks/java/library/j-seqalign/index.html?ca=dgr-jw17dynamicj ava&S_TACT= 105AGX59&S_CMP=GR
Григорьев Антон Владимирович -
аспирант кафедры «Информационная безопасность автоматизированных систем» Саратовского государственного технического университета имени Г агарина Ю. А.
Статья поступила в редакцию 3.02.12, принята к опубликованию 12.03.12
