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

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

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

В данной работе предлагается алгоритм поиска длинных разнесенных повторов. Лежащий в основе алгоритма обобщенный спектрально-аналитический метод, позволяет значительно ускорить процесс анализа последовательности за счет применения средств распаллеливания и векторизации. Также предлагается матрица спектральной гомологии генетических последовательностей. Близкая к точечной матрице гомологии, она предоставляет более быстрый инструмент для сравнительного анализа и визуализации внутренней структуры больших отрезков геномов (порядка 10e6 нуклеотидов), их тандемных и разнесенных повторов.

План презентации

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

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

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

Кстати, нужно использовать все современные возможности процессоров. Иначе будет обидно, если такую же программу написать на MATLAB’е и она — опа! — окажется быстрее в 5 раз. То есть нужно не забывать о многопоточности, не забывать об SIMD инструкциях, не забывать об аппаратном ускорении математических функций. Засчёт этого всего выигрываем в скорости ещё больше, реальная разница — в 10-20 раз (Core 2 Duo). Как?! Для многопоточности — голые нити (треды), никаких OpenMP! Так как это костылистая штуковина, приводит либо к сильному ухудшению структуры кода (причём фактическая логика получается аналогична голым тредам), либо к большим накладным расходам на распараллеливание — 5-15 %. Так что треды. Плюс библиотека Intel Integrated Performance Primitives для SIMD и аппаратного ускорения инструкций. А что это — IPP? А это такой векторный ассемблер, только на C. Библиотека, содержащая в себе оптимальные реализации большого спектра векторных операций (есть почти всё, что душе угодно — от сложений, умножений, корней и синусов, до узкоспециализированных функций ускорения декодирования аудио и видео, распознавания речи и т.д и т.п) для процессоров, имеющих различные расширения типа MMX / SSE1/2/3/4/5/ и т. д. Выражения над векторами там писать, к несчастью, нельзя, потому и получается код типа:

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

Вот где-то примерно это всё и было реализовано. Есть относительно простая программа, есть относительно хорошая библиотека для абстрагирования от деталей реализации конкретных базисов, есть сами базисы — Чебышева 1 и 2 рода, Якоби, Лежандра, Лагерра, Эрмита, Фурье, ДКП, ДСП. Она работает и рисует красивые картинки. [показать пару картинок и закончить]. Кстати, по поводу того, а какой же базис лучше? Вообще они все дают очень похожие результаты… Пока что «лучше» всех Чебышев 1-го рода. А что вообще такое «лучше»? «Лучше» — чисто умозрительно это «больше соотношение сигнал/шум» (в результатах). Как измерить? Ну, например, при одинаковых параметрах окон и глубине разложения подобрать eps такое, чтобы общее количество «похожих» участков было примерно равно, и посчитать, например, среднюю длину повторов. Можно и медиану тоже. Чем больше, тем лучше — мы ведь хотим найти как можно более длинные повторы. Начинали реализовывать с Чебышева 1-го рода, потом пробовали Лежандра, потом думали, что Чебышев 2-го рода произведёт революцию и всё будет гораздо лучше, т.к весовая функция выпуклая, центр отрезка учитывается сильнее, края меньше. Революции не произошло, результаты сильно похожие на 1-го рода, местами получше, местами похуже. Формально — похуже. Дальше есть табличка с «попугаями» по разным базисам. Тестовые данные — часть генома мыши (не спрашивайте какая, я не знаю) длиной 1.5 млн нуклеотидов. Сравнение приводилось при примерно одинаковых количествах найденных участков, «подозрительных» на повтор — в районе 5000. При выбранных настройках минимальная длина участка, подозрительного на повтор — 3500 нуклеотидов. Какие выводы? Лидирует Чебышев 1 рода. Базисы ДКП, ДСП и Фурье дают до жути похожие на него, практически идентичные, результаты. С небольшим отставанием за ними следует Лежандр, за ним — Чебышев 2 рода, а базисы Эрмита и Лагерра не подходят для поиска повторов вообще — что есть логичный факт, т. к. они оба работают на бесконечном интервале либо (0, +бесконечность), либо от — до + бесконечности. Вариантов значения медианной длины было всего 2: 3500 (минимально возможная) или 10000, она отражает, фактически, чистое количество шума — мелких отрезков, и гласит, что приемлемый уровень шума дают… Ясно кто.

-	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

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