13 521
правка
Изменения
м
Нет описания правки
Или '''«Применение обобщенного спектрально-аналитического метода в задаче анализа биологических данных»'''.
Ключевая задача анализа геномных последовательностей: поиск повторов. Прямых, обратных, симметричных. Что есть геномная последовательность? По сути, длинная строка в алфавите A, T, G, C (аденин, тимин, гуанин, цитозин — привет, биология за 10-й класс). T и C близки, это «[[rupedia:Пиримидин|пиримидиныпиримидиновые]]»основания». G и A тоже близки, это «[[rupedia:Пурин|пуриныпуриновые]]»основания». Методов куча, но есть '''Проблема: Последовательности Очень Длинные''', анализ долгий. Если искать точные повторы, ещё более-менее, но как только переходим к поиску неточных повторов, сразу всё сильно замедляется. По поводу «обычных» методов — например, можно посмотреть программу UniPro DPView — творение неких Новосибирских коллег. Ещё и адовые проекты [http://www.bioperl.org/ BioPerl], [http://www.biopython.org/ BioPython] — большие сборники различных методов и библиотек решения биологических задач — в частности, и методов поиска повторов.
'''ОСАМ.''' Мысль проста: разложить сигнал по какому-нибудь классическому ортогональному базису, получить краткое описание, к тому же обладающее различными приятными свойствами (норма сохраняется; отсекая хвост, получаем приближения; есть мера точности) и обработать не сигнал, а описание. Применим в широком спектре задач распознавания.
Идея — применить ОСАМ к поиску повторов в ДНК, таким образом ускорив его. Как?! Во-первых, построить профиль последовательности, т. е. перевести её в длинный числовой вектор, выбрав w — окно профиля, и принимая за каждый элемент последовательности ''(количество пуринов в w-окрестности элемента) минус (количество пиримидинов в w-окрестности элемента)''. Далее, выбирая по N значений из полученной последовательности — <m>(0 \ldots N-1), (s \ldots N+s-1), (2s \ldots N+2s-1), \ldots</m> (s — шаг аппроксимации) и раскладывая получаемые вектора из N чисел по k коэффициентам некоторого базиса, получить «индекс» последовательности. k << N, потому «индекс». Далее пробежаться по индексам обеих последовательностей (или одной и той же последовательности) и сравнить попарно все пары описаний на похожесть. А что такое похожесть? Критериев похожести можно выработать массу, среди них можно найти устойчивые к масштабу и т. п., однако у нас всё довольно просто: <m>\frac{|a-b|}{|a|+|b|}</m>, где <m>|x|=\sqrt{\sum {x}_{i}^{2}}</m>. Типа «нормированного L<sub>2</sub>-расстояния». Здесь можно выиграть от т. н. «принципа дискриминантности», который гласит очевидную вещь: если <m>\frac{\sqrt{{\sum }_{i=0}^{k}{({a}_{i}-{b}_{i})}^{2}}}{|a|+|b|}> \varepsilon</m> уже при k < n, суммирование можно не продолжать, т. к. ''меньше'' сумма квадратов уже не станет.
Итак, от этого сравнения мы получим оценку «подобия» участков ДНК. Крупных или мелких, более или менее точное сравнение — это уже как захотим — для этого можно варьировать параметры. Задаём порог, можем пробежаться по результатам и сразу выявить участки, «подозрительные на повторы». То есть больше не нужно всё время искать повторы "«''везде''"»: сначала достаточно выявить крупные относительно похожие участки, а потом можно «увеличить масштаб» и выявить (или не выявить) точные координаты повторов. Единственное, для чего подход почти не подходит — для выявления «абсолютно точных» координат повторов. Это уже в «подозрительных» областях можно делать стандартными методами. Например, diff'оподобными алгоритмами.
== Часть статьи ==
Для реализации программы поиска повторов с помощью ОСАМ был выбран язык C++. Такой выбор обусловлен сущностью процесса разложения функций, позволяющей с помощью объектно-ориентированного подхода разделить функционал на общий и зависящий от конкретного ортогонального базиса. Общий функционал — это функции подсчёта весовых коэффициентов, подсчёта интеграла на сетке Гаусса, подсчёта матрицы Грама заданного базиса, нормирования заданного базиса, интерполяции сигнала на заданную сетку, и воссоздания изначального сигнала по коэффициентам разложения. К базисо-зависимому функционалу относятся функции подсчёта сетки, весовых коэффициентов, и самих значений функции. Также такой подход, кроме всего прочего, даёт возможность оптимизировать части функционала отдельнодруг от друга.
=== «Наивный» алгоритм ===
*** Сохранить подсчитанное значение как (i, j)-ый элемент матрицы гомологии.
* Записать матрицу гомологии в выходной файл.
Подготовительный этап:
* Подсчитать сетку Гаусса (то есть, корни n+1-ой функции базиса).
* Подсчитать весовые и нормировочные коэффициенты.
=== Алгоритм разложения ===
«Наивный» вариант алгоритма разложения:
* Интерполировать выбранную подпоследовательность длины N > n на подсчитанную сетку алгоритмом «ближайшего соседа».
: То есть, по сути, не интерполировать её никак. Практика показала, что любая предварительная интерполяция никак не улучшает разложение по причине большой плотности точек в исходном сигнале и маленькой — в раскладываемом массиве.
* Подсчитать в цикле <m>c_j = \sum_{i=1}^{n} y_i \cdot w_i \cdot f_j(x_i) \cdot r_j, j=1 \ldots n</m>, где:
: <m>c_j</m> — j-ый коэффициент разложения сигнала <m>y_i</m>.
: <m>w_i</m> — i-ый весовой коэффициент.
: <m>f_j(x_i)</m> — значение j-ой функции базиса в i-ой точке сетки.
: <m>r_j</m> — j-ый нормировочный коэффициент.
Оптимизированный для рекуррентных соотношений алгоритм разложения:
* Интерполировать выбранную подпоследовательность длины N > n на подсчитанную сетку алгоритмом «ближайшего соседа».
* В цикле по ''i = 1..n'':
** <m>c_i = 0</m>
* В цикле по ''i = 1..n'':
** Вычислить и сохранить в памяти все значения <m>f_j(x_i), j = 1 \ldots n</m> с помощью рекуррентных соотношений.
** В цикле по ''j = 1..n'':
*** <m>c_j = c_j + y_i \cdot f_j(x_i) \cdot r_j \cdot w_i</m>
Псевдокод оптимизированного с учётом векторных операций алгоритма разложения здесь не приведён по причине его объёма. Кратко можно описать два момента: во-первых, циклы сменены местами — внешний цикл идёт по коэффициентам разложения, а не по функциям базиса, и во-вторых, на всех этапах используются векторные операции — сложения, умножения, возведения в квадрат и т. п.
=== Оптимизация ===
Кроме того, ОСАМ позволяет и производить практически идеальное распараллеливание алгоритма по причине небольшого объёма необходимой памяти в случае, если не используется т. н. «индексация последовательности» — такой подход может быть полезен при вычислениях с массовым параллелизмом. ''Индексацией'' называется процесс предварительного разложения сравниваемой последовательности по выбранному ортогональному базису и сохранения в памяти всех векторов коэффициентов разложения для последующего использования. Достоинство индексации — отсутствие необходимости производить большой объём вычислений во вложенном цикле; её недостаток — существенное увеличение объёма используемой оперативной памяти и увеличение требований к пропускной способности памяти. Последнее особенно важно при массивно-параллельных вычислениях — отдельные процессоры, ядра или узлы кластера могут вообще не иметь общего доступа ко всей оперативной памяти системы, не говоря уже о существенном замедлении обмена данных между вычислителями и памятью в случае конкуретной работы с большой области памяти. Такая проблема присутствует даже на многоядерных стандартных настольных компьютерах и серверах нижнего класса — оперативная память обычно работает приблизительно со скоростью, равной четверти скорости процессоров и, начиная с определённого количества ядер/процессоров, индексация становится менее выгодной, чем могла бы быть, так как чипсет и оперативная память не могут обеспечить требуемую скорость обмена.
Тем не менее, на обычных ПК и серверах нижнего класса наличие индексации хотя бы одной последовательности всё равно выгодно, поэтому при реализации был выбран следующий подход: индексация одной последовательности и разложений разложение второй на лету. Соответственно, в любом случае — как в случае сравнения последовательности с самой собой, так и в случае сравнения двух последовательностей — вычисления коэффициентов разложения последовательностей происходят только 1 раз: первой при индексации, а второй во внешнем цикле.
Реальный выигрыш в производительности засчёт чисто программной оптимизации достигает 10-20 раз на стандартных двухъядерных процессорах архитектуры Core 2.
ippsMulC_64f(wn, -1, tn, n);
ippsSqrt_64f_I(tn, n);
И последний важный момент — принцип «дискриминантности». Напомним, что расстояние между двумя векторами коэффициентов разложения определяется как <m>\frac{|a-b|}{|a|+|b|}</m>, где <m>|x|=\sqrt{\sum {x}_{i}^{2}}</m> Принцип «дискриминантности» же гласит очевидную вещь: если <m>\frac{\sqrt{{\sum }_{i=0}^{k}{({a}_{i}-{b}_{i})}^{2}}}{|a|+|b|}> \varepsilon</m> уже при k < n, суммирование можно не продолжать, т.к. ''меньше'' ε сумма квадратов уже не станет. Эта идея также использовалась при оптимизации алгоритма. Однако здесь возникает определённое препятствие: суммирование с постоянными условными проверками не векторизуется, т.е., при подсчёте нормы с учётом принципа "дискриминантности" IPP использовать мы уже не можем. Но так как IPP даёт весьма неплохой прирост производительности, можно применить следующий нетривиальный ход: сначала суммировать до ''k = d'', где d - делитель n, больший 1, с использованием векторных операций, потом проверять, не превышен ли порог, потом до ''k = 2d'', потом до ''k = 3d'', и т.д.
=== Алгоритм с учётом индексации ===
С учётом выбранного подхода — индексации одной последовательности и разложения другой «на лету» — алгоритм принимает следующий вид:
* Загрузить входные файлы последовательностей.
* ''Подсчитать и сохранить в памяти коэффициенты разложения всех подпоследовательностей 1-ой последовательности по выбранному ОНБ.''
* ''Подсчитать и сохранить в памяти нормы всех векторов коэффициентов разложения этих подпоследовательностей.''
* По всем ''сохранённым коэффициентам разложения подпоследовательностей'' 1-ой последовательности:
** По всем подпоследовательностям 2-ой последовательности:
*** Подсчитать коэффициенты разложения подпоследовательности 2-ой последовательности по выбранному ОНБ.
*** Вычислить норму вектора коэффициентов разложения подпоследовательности 2-ой последовательности.
*** Подсчитать L<sub>2</sub>-расстояние между векторами коэффициентов разложения подпоследовательностей.
*** Поделить подсчитанное расстояние на сумму норм векторов коэффициентов.
*** Сохранить подсчитанное значение как (i, j)-ый элемент матрицы гомологии.
* Записать матрицу гомологии в выходной файл.
=== Алгоритм с учётом параллелизма ===
Изменения с учётом параллелизма тривиальны: наиболее внешние циклы разделяются на ''M'' частей и для обработки каждой части работы создаётся собственный поток. Далее главный поток приложения ожидает завершения всех созданных, т.е., ожидает окончания очередного этапа работы.
* Загрузить входные файлы последовательностей.
* ''Создать требуемое число M вычислительных потоков, далее, для каждого из них:''
** ''Подсчитать и сохранить в памяти коэффициенты разложения своей <m>\frac{1}{M}</m>-ой части подпоследовательностей 1-ой последовательности по выбранному ОНБ.''
** ''Подсчитать и сохранить в памяти нормы своей <m>\frac{1}{M}</m>-ой части векторов коэффициентов разложения этих подпоследовательностей.''
* ''Создать требуемое число M вычислительных потоков, далее, для каждого из них:''
** ''По своей <m>\frac{1}{M}</m>-ой части сохранённых коэффициентов разложения подпоследовательностей 1-ой последовательности'':
*** По всем подпоследовательностям 2-ой последовательности:
**** Подсчитать коэффициенты разложения подпоследовательности 2-ой последовательности по выбранному ОНБ.
**** Вычислить норму вектора коэффициентов разложения подпоследовательности 2-ой последовательности.
**** Подсчитать L<sub>2</sub>-расстояние между векторами коэффициентов разложения подпоследовательностей.
**** Поделить подсчитанное расстояние на сумму норм векторов коэффициентов.
**** Сохранить подсчитанное значение как (i, j)-ый элемент матрицы гомологии.
* Записать матрицу гомологии в выходной файл.
=== Сравнение векторов с учётом векторных операций и дискриминантности ===
* Вычислить относительный порог <m>l = (\varepsilon \cdot (s_1 + s_2))²</m>, где s<sub>1</sub> и s<sub>2</sub> — нормы векторов.
* Начальное значение ''f = 0''.
* В цикле:
** С помощью функции IPP <code>ippsNormDiffL2_64f</code> (или 32f, в зависимости от требуемой точности) вычислить норму разности очередных участков длины ''d'' сравниваемых векторов.
** Добавить к ''f'' квадрат полученного значения.
** Если ''f > l'', принять, что вектора «не подобны».
* Если цикл завершился без принятия того, что вектора «не подобны», принять, что вектора подобны.
=== Сравнение ОНБ ===
Ниже приводится табличка с замерами средней длины найденных повторов на различных базисах и части генома мыши длиной 1.5 млн нуклеотидов в качестве тестовых данных. Сравнение производилось при приблизительно равных количествах найденных «подобных» участков — 5000. При выбранных настройках минимально возможная найденная длина подобного участка — 3500 нуклеотидов.
<tab sep="tab" border="1" class="simpletable" head="topleft">
- Eps Среднее Медиана
Чебышева 1 рода .025 '''3978''' '''10000'''
Каковы выводы? По средней длине повтора лидирует базис Чебышев 1 рода, а базисы ДКП, ДСП и Фурье дают чрезвычайно похожие на него, практически идентичные, результаты. С небольшим отставанием следует базис Лежандра, далее — базис Чебышева 2 рода, а базисы Эрмита и Лагерра для поиска подобных участков не подходят вообще, чему есть простое математическое обоснование — оба они действуют на бесконечной полупрямой — либо <m>(0, +\inf)</m>, либо <m>(-\inf, +\inf)</m>. Вариантов значения медианной длины при этом было всего 2: 3500 (минимально возможная) или 10000. Медианная длина в данном случае отражает, фактически, «чистое» количество шума — мелких отрезков, и гласит, что приемлемый уровень шума дают базисы Чебышева 1 рода, ДКП, ДСП, Фурье и Лежандра.
[[Категория:УчёбаСтатьи]][[Категория:Биоинформатика]]