Поиск повторов в ДНК на основе ОСАМ

Ссылки сюда (1) →

Или «Применение обобщенного спектрально-аналитического метода в задаче анализа биологических данных».

Ключевая задача анализа геномных последовательностей: поиск повторов. Прямых, обратных, симметричных. Что есть геномная последовательность? По сути, длинная строка в алфавите A, T, G, C (аденин, тимин, гуанин, цитозин — привет, биология за 10-й класс). T и C близки, это «пиримидиновые основания». G и A тоже близки, это «пуриновые основания». Методов куча, но есть Проблема: Последовательности Очень Длинные, анализ долгий. Если искать точные повторы, ещё более-менее, но как только переходим к поиску неточных повторов, сразу всё сильно замедляется. По поводу «обычных» методов — например, можно посмотреть программу UniPro DPView — творение неких Новосибирских коллег. Ещё и адовые проекты BioPerl, BioPython — большие сборники различных методов и библиотек решения биологических задач — в частности, и методов поиска повторов.

ОСАМ. Мысль проста: разложить сигнал по какому-нибудь классическому ортогональному базису, получить краткое описание, к тому же обладающее различными приятными свойствами (норма сохраняется; отсекая хвост, получаем приближения; есть мера точности) и обработать не сигнал, а описание. Применим в широком спектре задач распознавания.

Идея — применить ОСАМ к поиску повторов в ДНК, таким образом ускорив его. Как?! Во-первых, построить профиль последовательности, т. е. перевести её в длинный числовой вектор, выбрав w — окно профиля, и принимая за каждый элемент последовательности (количество пуринов в w-окрестности элемента) минус (количество пиримидинов в w-окрестности элемента). Далее, выбирая по N значений из полученной последовательности — (s — шаг аппроксимации) и раскладывая получаемые вектора из N чисел по k коэффициентам некоторого базиса, получить «индекс» последовательности. k << N, потому «индекс». Далее пробежаться по индексам обеих последовательностей (или одной и той же последовательности) и сравнить попарно все пары описаний на похожесть. А что такое похожесть? Критериев похожести можно выработать массу, среди них можно найти устойчивые к масштабу и т. п., однако у нас всё довольно просто: , где . Типа «нормированного L₂-расстояния». Здесь можно выиграть от т. н. «принципа дискриминантности», который гласит очевидную вещь: если уже при k < n, суммирование можно не продолжать, т. к. меньше сумма квадратов уже не станет.

Итак, от этого сравнения мы получим оценку «подобия» участков ДНК. Крупных или мелких, более или менее точное сравнение — это уже как захотим — для этого можно варьировать параметры. Задаём порог, можем пробежаться по результатам и сразу выявить участки, «подозрительные на повторы». То есть больше не нужно всё время искать повторы «везде»: сначала достаточно выявить крупные относительно похожие участки, а потом можно «увеличить масштаб» и выявить (или не выявить) точные координаты повторов. Единственное, для чего подход почти не подходит — для выявления «абсолютно точных» координат повторов. Это уже в «подозрительных» областях можно делать стандартными методами. Например, diff'оподобными алгоритмами.

Часть статьи

Для реализации программы поиска повторов с помощью ОСАМ был выбран язык C++. Такой выбор обусловлен сущностью процесса разложения функций, позволяющей с помощью объектно-ориентированного подхода разделить функционал на общий и зависящий от конкретного ортогонального базиса. Общий функционал — это функции подсчёта весовых коэффициентов, подсчёта интеграла на сетке Гаусса, подсчёта матрицы Грама заданного базиса, нормирования заданного базиса, интерполяции сигнала на заданную сетку, и воссоздания изначального сигнала по коэффициентам разложения. К базисо-зависимому функционалу относятся функции подсчёта сетки, весовых коэффициентов, и самих значений функции. Также такой подход, кроме всего прочего, даёт возможность оптимизировать части функционала отдельно друг от друга.

«Наивный» алгоритм

В целом основная задача программного обеспечения поиска повторов на основе ОСАМ — построение спектральной матрицы гомологии последовательности, в общем случае — двух последовательностей. При сравнении двух последовательностей каждый элемент спектральной матрицы гомологии отражает оценку подобия соответствующих участков последовательностей. Также последовательность можно сравнивать с самой собой.

Простейший «наивный» вариант алгоритма построения матрицы гомологии:

• Загрузить входные файлы последовательностей.
• По всем подпоследовательностям 1-ой последовательности:
  • Подсчитать коэффициенты разложения подпоследовательности по выбранному ОНБ.
  • Вычислить норму вектора коэффициентов.
  • По всем подпоследовательностям 2-ой последовательности:
    • Подсчитать коэффициенты разложения подпоследовательности 2-ой последовательности по выбранному ОНБ.
    • Вычислить норму вектора коэффициентов разложения подпоследовательности 2-ой последовательности.
    • Подсчитать L₂-расстояние между векторами коэффициентов разложения подпоследовательностей.
    • Поделить подсчитанное расстояние на сумму норм векторов коэффициентов.
    • Сохранить подсчитанное значение как (i, j)-ый элемент матрицы гомологии.
• Записать матрицу гомологии в выходной файл.

Подготовительный этап:

• Подсчитать сетку Гаусса (то есть, корни n+1-ой функции базиса).
• Подсчитать весовые и нормировочные коэффициенты.

«Наивный» вариант алгоритма разложения:

• Интерполировать выбранную подпоследовательность длины N > n на подсчитанную сетку алгоритмом «ближайшего соседа».

То есть, по сути, не интерполировать её никак. Практика показала, что любая предварительная интерполяция никак не улучшает разложение по причине большой плотности точек в исходном сигнале и маленькой — в раскладываемом массиве.

• Подсчитать в цикле , где:

 — j-ый коэффициент разложения сигнала .

 — i-ый весовой коэффициент.

 — значение j-ой функции базиса в i-ой точке сетки.

 — j-ый нормировочный коэффициент.

Оптимизированный для рекуррентных соотношений алгоритм разложения:

• Интерполировать выбранную подпоследовательность длины N > n на подсчитанную сетку алгоритмом «ближайшего соседа».
• В цикле по i = 1..n:
• В цикле по i = 1..n:
  • Вычислить и сохранить в памяти все значения с помощью рекуррентных соотношений.
  • В цикле по j = 1..n:

Оптимизация

При реализации системы поиска повторов в виде программы учитывалась необходимость использования всех современных возможностей процессоров — ведь нужно понимать, что в наше время процессоры уже давно не i386, все суперскалярные, поддерживающие многопоточность, SIMD-инструкции (Single Instruction, Multiple Data) — инструкции, позволяющие за один такт выполнить несколько одинаковых операций сразу, аппаратно ускоренные математические функции и другие возможности поднятия производительности. Также не следует забывать, что большинство из этих возможностей успешно используется математическими пакетами вроде Matlab и Maple, популярными при тестировании и исследованиях математических методов. Поэтому, если забыть об этих возможностях в программе, можно испытать разочарование от скорости работы по сравнению с той же программой, реализованной с помощью математического пакета. К счастью, общий алгоритм разложения дискретизированных сигналов по классическим ортогональным базисам, являющийся просто алгоритмом вычисления соответствующего интеграла Гаусса, весьма прост и допускает оптимизацию также с помощью простых методов.

Кроме того, ОСАМ позволяет и производить практически идеальное распараллеливание алгоритма по причине небольшого объёма необходимой памяти в случае, если не используется т. н. «индексация последовательности» — такой подход может быть полезен при вычислениях с массовым параллелизмом. Индексацией называется процесс предварительного разложения сравниваемой последовательности по выбранному ортогональному базису и сохранения в памяти всех векторов коэффициентов разложения для последующего использования. Достоинство индексации — отсутствие необходимости производить большой объём вычислений во вложенном цикле; её недостаток — существенное увеличение объёма используемой оперативной памяти и увеличение требований к пропускной способности памяти. Последнее особенно важно при массивно-параллельных вычислениях — отдельные процессоры, ядра или узлы кластера могут вообще не иметь общего доступа ко всей оперативной памяти системы, не говоря уже о существенном замедлении обмена данных между вычислителями и памятью в случае конкуретной работы с большой области памяти. Такая проблема присутствует даже на многоядерных стандартных настольных компьютерах и серверах нижнего класса — оперативная память обычно работает приблизительно со скоростью, равной четверти скорости процессоров и, начиная с определённого количества ядер/процессоров, индексация становится менее выгодной, чем могла бы быть, так как чипсет и оперативная память не могут обеспечить требуемую скорость обмена.

Тем не менее, на обычных ПК и серверах нижнего класса наличие индексации хотя бы одной последовательности всё равно выгодно, поэтому при реализации был выбран следующий подход: индексация одной последовательности и разложений второй на лету. Соответственно, в любом случае — как в случае сравнения последовательности с самой собой, так и в случае сравнения двух последовательностей — вычисления коэффициентов разложения последовательностей происходят только 1 раз: первой при индексации, а второй во внешнем цикле.

Реальный выигрыш в производительности засчёт чисто программной оптимизации достигает 10-20 раз на стандартных двухъядерных процессорах архитектуры Core 2.

Очевидными вариантами достижения параллелизма в алгоритме поиска повторов являются библиотека OpenMP и ручная реализация распараллеливания на основе потоков — в UNIX-среде pthreads (POSIX threads — потоки POSIX), а в Windows-среде функций WINAPI. Можно было бы предположить, что использование библиотеки OpenMP упростит переносимость программы, однако, при переопределении всего лишь двух функций — создания потока и ожидания завершения потока (т. н. «join») — ручной подход достигает в точности такой же идеальной переносимости программы. Собственно говоря, функции создания потока и ожидания завершения потока являются настолько базовыми в любой библиотеке работы с потоками на любой платформе, поддерживающей потоки, что при реализации можно не бояться их потенциального отсутствия, тем более, когда на дворе 2009-ый год. Вместе с тем как раз реализация OpenMP потенциально существует не для всех ОС.

Главным же минусом библиотеки OpenMP является то, что её работа построена на директивах компилятора, и в итоге транслируется обычно в код, постоянно создающий и завершающий вычислительные потоки, для каждой итерации распараллеливаемого цикла. Таким образом при использовании OpenMP либо приходится учитывать такое поведение, распараллеливая циклы с небольшими (по крайней мере, относительно) количествами итераций, ухудшая структуру кода и фактически сводя его логику к логике ручного распараллеливания, либо мириться с накладными расходами на распараллеливание, в нашем случае достигавшими 5-15 %.

Таким образом, для параллелизма использовалось ручное разделение задачи на подзадачи и ручное управление вычислительными потоками.

Для использования аппаратно-ускоренных и векторных (SIMD) инструкций использовалась библиотека Intel Integrated Performance Primitives (IPP). Ближайшая сравнение IPP — «векторный язык ассемблера», содержащий простые векторные «инструкции», а точнее оптимизированные функции-обёртки, для весьма широкого спектра задач — от сложений, умножений, корней и синусов, до узкоспециализированных функций ускорения декодирования аудио и видео, распознавания речи и т. п. Библиотека IPP даёт преимущества при использовании любых x86-процессоров, имеющих расширения наборов команд MMX, SSE, SSE2, SSE3 и т. п. Нужно отметить, что IPP сравнима в первую очередь действительно с языком ассемлера, так как не поддерживает трансляцию выражений над векторами, а только сами операции, реализованные в виде функций (аналог инструкций). Это, к сожалению, приводит к неочевидному «ассемблерному» коду следующего вида:

ippsCopy_64f(xn, wn, n);
ippsSqr_64f_I(wn, n);
ippsAddC_64f_I(-1, wn, n);
ippsMulC_64f(wn, -1, tn, n);
ippsSqrt_64f_I(tn, n);

И последний важный момент — принцип «дискриминантности». Напомним, что расстояние между двумя векторами коэффициентов разложения определяется как , где Принцип «дискриминантности» же гласит очевидную вещь: если уже при k < n, суммирование можно не продолжать, т.к. меньше ε сумма квадратов уже не станет. Эта идея также использовалась при оптимизации алгоритма. Однако здесь возникает определённое препятствие: суммирование с постоянными условными проверками не векторизуется, т.е., при подсчёте нормы с учётом принципа "дискриминантности" IPP использовать мы уже не можем. Но так как IPP даёт весьма неплохой прирост производительности, можно применить следующий нетривиальный ход: сначала суммировать до k = d, где d - делитель n, больший 1, с использованием векторных операций, потом проверять, не превышен ли порог, потом до k = 2d, потом до k = 3d, и т.д.

Алгоритм с учётом индексации

С учётом выбранного подхода - индексации одной последовательности и разложения другой "на лету" - алгоритм принимает следующий вид:

• Загрузить входные файлы последовательностей.
• Подсчитать и сохранить в памяти коэффициенты разложения всех подпоследовательностей 1-ой последовательности по выбранному ОНБ.
• Подсчитать и сохранить в памяти нормы всех векторов коэффициентов разложения этих подпоследовательностей.
• По всем сохранённым коэффициентам разложения подпоследовательностей 1-ой последовательности:
  • По всем подпоследовательностям 2-ой последовательности:
    • Подсчитать коэффициенты разложения подпоследовательности 2-ой последовательности по выбранному ОНБ.
    • Вычислить норму вектора коэффициентов разложения подпоследовательности 2-ой последовательности.
    • Подсчитать L₂-расстояние между векторами коэффициентов разложения подпоследовательностей.
    • Поделить подсчитанное расстояние на сумму норм векторов коэффициентов.
    • Сохранить подсчитанное значение как (i, j)-ый элемент матрицы гомологии.
• Записать матрицу гомологии в выходной файл.

Сравнение векторов с учётом векторных операций и дискриминантности

• Вычислить относительный порог , где s₁ и s₂ — нормы векторов.
• Начальное значение f = 0.
• В цикле:
  • С помощью функции IPP ippsNormDiffL2_64f (или 32f, в зависимости от требуемой точности) вычислить норму разности очередных участков длины d сравниваемых векторов.
  • Добавить к f квадрат полученного значения.
  • Если f > l, принять, что вектора «не подобны».
• Если цикл завершился без принятия того, что вектора «не подобны», принять, что вектора подобны.

Сравнение ОНБ

Учитывая, что поиск повторов может осуществляться по выбору с использованием любого из ортогональных базисов, и что в библиотеке функций разложения их было реализовано 9 различных — базис Чебышева 1 рода, базис Чебышева 2 рода, дискретные косинусное и синусное преобразования, базис Фурье, базис Лежандра, базис Лагерра, базис Якоби и базис Эрмита — очевидным образом встаёт вопрос: а какой же из них «лучше» в задаче поиска повторов в последовательностях? А кроме того, каковы в целом критерии качества, по которым требуется производить сравнение базисов?

Очевидным подходом к данному вопросу является критерий «максимум соотношения сигнал/шум в найденных в итоге повторах».

Другой вариант — максимум средней длины найденных подобных участков, так как цель поиска повторов заключается в том, чтобы найти как можно более длинные подобные участки. Как можно оценить эту длину? Опишем простейший подход. Во-первых, нужно выбрать ширину скользящих окон и глубину разложения и выбрать некоторые тестовые данные, содержащие широкий спектр различных повторов — здесь хорошо подходит часть реальной ДНК-последовательности. Далее, используя различные базисы и подбирая порог сравнения () такой, чтобы общее число найденных подобных участков было приблизительно равно, подсчитывать среднюю длину найденных подобных участков. Как вариант — можно вычислять медианное значение.

В процессе реализации программы вначале был выбран базис Чебышева 1-го рода; потом пробовали базис Лежандра. Потом было высказано предположение о том, что базис Чебышева 2-го рода произведёт «революцию» по той причине, что имеет выпуклую весовую функцию и сильнее учитывает центр сравниваемого отрезка, чем края, но революции не произошло, результаты базиса Чебышева 2-го рода сильно похожи на базис Чебышева 1-го рода, и даже немного хуже, в том числе и по средней длине найденных повторов.

Ниже приводится табличка с замерами средней длины найденных повторов на различных базисах и части генома мыши длиной 1.5 млн нуклеотидов в качестве тестовых данных. Сравнение производилось при приблизительно равных количествах найденных «подобных» участков — 5000. При выбранных настройках минимально возможная найденная длина подобного участка — 3500 нуклеотидов.

-	Eps	Среднее	Медиана
Чебышева 1 рода	.025	3978	10000
Чебышева 2 рода	.0285	3882	3500
ДКП	.025	3978	10000
ДСП	.021	3975	10000
Фурье	.025	3978	10000
Эрмита	.0015	3502	3500
Лагерра	.0063	3505	3500
Лежандра	.0225	3966	10000

Каковы выводы? По средней длине повтора лидирует базис Чебышев 1 рода, а базисы ДКП, ДСП и Фурье дают чрезвычайно похожие на него, практически идентичные, результаты. С небольшим отставанием следует базис Лежандра, далее — базис Чебышева 2 рода, а базисы Эрмита и Лагерра для поиска подобных участков не подходят вообще, чему есть простое математическое обоснование — оба они действуют на бесконечной полупрямой — либо , либо . Вариантов значения медианной длины при этом было всего 2: 3500 (минимально возможная) или 10000. Медианная длина в данном случае отражает, фактически, «чистое» количество шума — мелких отрезков, и гласит, что приемлемый уровень шума дают базисы Чебышева 1 рода, ДКП, ДСП, Фурье и Лежандра.