Поиск повторов в ДНК на основе ОСАМ — различия между версиями

Версия 19:18, 29 августа 2009

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

Ключевая задача анализа геномных последовательностей: поиск повторов. Прямых, обратных, симметричных. Что есть геномная последовательность? По сути, длинная строка в алфавите 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)-ый элемент матрицы гомологии.
• Записать матрицу гомологии в выходной файл.

Алгоритм с учётом параллелизма

Изменения с учётом параллелизма тривиальны: наиболее внешние циклы разделяются на M частей и для обработки каждой части работы создаётся собственный поток. Далее главный поток приложения ожидает завершения всех созданных, т.е., ожидает окончания очередного этапа работы.

• Загрузить входные файлы последовательностей.
• Создать требуемое число M вычислительных потоков, далее, для каждого из них:
  • Подсчитать и сохранить в памяти коэффициенты разложения своей -ой части подпоследовательностей 1-ой последовательности по выбранному ОНБ.
  • Подсчитать и сохранить в памяти нормы своей -ой части векторов коэффициентов разложения этих подпоследовательностей.
• Создать требуемое число M вычислительных потоков, далее, для каждого из них:
  • По своей -ой части сохранённых коэффициентов разложения подпоследовательностей 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 нуклеотидов.

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