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

Материал из YourcmcWiki
Перейти к: навигация, поиск
Строка 2: Строка 2:
  
 
В данной работе предлагается алгоритм поиска длинных разнесенных повторов. Лежащий в основе алгоритма обобщенный спектрально-аналитический метод, позволяет значительно ускорить процесс анализа последовательности за счет применения средств распаллеливания и векторизации. Также предлагается матрица спектральной гомологии генетических последовательностей. Близкая к точечной матрице гомологии, она предоставляет более быстрый инструмент для сравнительного анализа и визуализации внутренней структуры больших отрезков геномов (порядка 10e6 нуклеотидов), их тандемных и разнесенных повторов.
 
В данной работе предлагается алгоритм поиска длинных разнесенных повторов. Лежащий в основе алгоритма обобщенный спектрально-аналитический метод, позволяет значительно ускорить процесс анализа последовательности за счет применения средств распаллеливания и векторизации. Также предлагается матрица спектральной гомологии генетических последовательностей. Близкая к точечной матрице гомологии, она предоставляет более быстрый инструмент для сравнительного анализа и визуализации внутренней структуры больших отрезков геномов (порядка 10e6 нуклеотидов), их тандемных и разнесенных повторов.
 
== План презентации ==
 
  
 
Ключевая задача анализа геномных последовательностей: поиск повторов. Прямых, обратных, симметричных. Что есть геномная последовательность? По сути, длинная строка в алфавите A, T, G, C (аденин, тимин, гуанин, цитозин, привет, биология, 10-й класс). T и C близки, это «пиримидины». G и A тоже близки, это «пурины». Методов куча, но есть и Проблема: последовательности очень длинные, анализ долгий. Если искать точные повторы, ещё более-менее, но как только переходим к поиску неточных повторов, всё сразу сильно замедляется. По поводу «обычных» методов — например, можно посмотреть программу UniPro DPView — творение неких Новосибирских коллег. Ещё есть довольно адские проекты BioPerl, BioPython — большие сборники всяких методов и библиотек по поводу биологических задач, в частности, и методов поиска повторов, на скриптовых языках.
 
Ключевая задача анализа геномных последовательностей: поиск повторов. Прямых, обратных, симметричных. Что есть геномная последовательность? По сути, длинная строка в алфавите A, T, G, C (аденин, тимин, гуанин, цитозин, привет, биология, 10-й класс). T и C близки, это «пиримидины». G и A тоже близки, это «пурины». Методов куча, но есть и Проблема: последовательности очень длинные, анализ долгий. Если искать точные повторы, ещё более-менее, но как только переходим к поиску неточных повторов, всё сразу сильно замедляется. По поводу «обычных» методов — например, можно посмотреть программу UniPro DPView — творение неких Новосибирских коллег. Ещё есть довольно адские проекты BioPerl, BioPython — большие сборники всяких методов и библиотек по поводу биологических задач, в частности, и методов поиска повторов, на скриптовых языках.
Строка 11: Строка 9:
 
Идея: применить ОСАМ к поиску повторов в ДНК, таким образом ускорив его. Как?! Во-первых, построить профиль последовательности, т.&nbsp;е. перевести её в длинный числовой вектор, выбрав w — окно профиля, и принимая за каждый элемент последовательности (количество пуринов в w-окрестности элемента) минус (количество пиримидинов в w-окрестности элемента). Далее, выбирая по N значений из полученной последовательности — 0..N-1, s..N+s-1, 2s..N+2s-1, … (s — шаг аппроксимации) и раскладывая получаемые вектора из N чисел по k коэффициентам некоторого базиса, получить «индекс» последовательности. k << N, потому и «индекс». Далее пробежаться по всем полученным описаниям (по индексу) обеих последовательностей (или одной и той же последовательности) и сравнить попарно все пары описаний (на похожесть). А что такое похожесть? Критериев похожести можно выработать массу, среди них можно найти устойчивые к масштабу и т.&nbsp;п., однако у нас всё довольно просто: <m>\frac{|a-b|}{|a|+|b|}</m>, где <m>|x|=\sqrt{\sum {x}_{i}^{2}}</m>. Такое вот «нормированное L<sub>2</sub>-расстояние». Здесь, кстати, можно выиграть от т.&nbsp;н. «принципа дискриминантности», который гласит очевидную вещь: что если <m>\frac{\sqrt{{\sum }_{i=0}^{k}{({a}_{i}-{b}_{i})}^{2}}}{|a|+|b|}> \varepsilon</m> уже при k < n, то суммирование можно не продолжать, т.&nbsp;к. меньше сумма квадратов уже не станет. Итак, что мы получим от этого сравнения? Мы получим приближённые «близости» участков ДНК. Крупных или мелких, более или менее точное сравнение — это уже как захотим — для этого можно варьировать параметры. Задаём порог, можем пробежаться по результатам и сразу выявить «подозрительные на повторы» участки. Это есть важно, т.&nbsp;к. больше не нужно всё время искать повторы ВЕЗДЕ: сначала достаточно выявить крупные относительно похожие участки, а потом можно «увеличить масштаб» и выявить (или не выявить) точные координаты повторов. Кстати, единственное, для чего подход почти не подходит — для выявления «абсолютно точных» координат повторов. Это уже в «подозрительных» областях можно делать стандартными методами. Например, diffоподобным алгоритмом. :-)
 
Идея: применить ОСАМ к поиску повторов в ДНК, таким образом ускорив его. Как?! Во-первых, построить профиль последовательности, т.&nbsp;е. перевести её в длинный числовой вектор, выбрав w — окно профиля, и принимая за каждый элемент последовательности (количество пуринов в w-окрестности элемента) минус (количество пиримидинов в w-окрестности элемента). Далее, выбирая по N значений из полученной последовательности — 0..N-1, s..N+s-1, 2s..N+2s-1, … (s — шаг аппроксимации) и раскладывая получаемые вектора из N чисел по k коэффициентам некоторого базиса, получить «индекс» последовательности. k << N, потому и «индекс». Далее пробежаться по всем полученным описаниям (по индексу) обеих последовательностей (или одной и той же последовательности) и сравнить попарно все пары описаний (на похожесть). А что такое похожесть? Критериев похожести можно выработать массу, среди них можно найти устойчивые к масштабу и т.&nbsp;п., однако у нас всё довольно просто: <m>\frac{|a-b|}{|a|+|b|}</m>, где <m>|x|=\sqrt{\sum {x}_{i}^{2}}</m>. Такое вот «нормированное L<sub>2</sub>-расстояние». Здесь, кстати, можно выиграть от т.&nbsp;н. «принципа дискриминантности», который гласит очевидную вещь: что если <m>\frac{\sqrt{{\sum }_{i=0}^{k}{({a}_{i}-{b}_{i})}^{2}}}{|a|+|b|}> \varepsilon</m> уже при k < n, то суммирование можно не продолжать, т.&nbsp;к. меньше сумма квадратов уже не станет. Итак, что мы получим от этого сравнения? Мы получим приближённые «близости» участков ДНК. Крупных или мелких, более или менее точное сравнение — это уже как захотим — для этого можно варьировать параметры. Задаём порог, можем пробежаться по результатам и сразу выявить «подозрительные на повторы» участки. Это есть важно, т.&nbsp;к. больше не нужно всё время искать повторы ВЕЗДЕ: сначала достаточно выявить крупные относительно похожие участки, а потом можно «увеличить масштаб» и выявить (или не выявить) точные координаты повторов. Кстати, единственное, для чего подход почти не подходит — для выявления «абсолютно точных» координат повторов. Это уже в «подозрительных» областях можно делать стандартными методами. Например, diffоподобным алгоритмом. :-)
  
Вот где-то примерно это всё и было реализовано. Есть относительно простая программа, есть относительно хорошая библиотека для абстрагирования от деталей реализации конкретных базисов, есть сами базисы — Чебышева 1 и 2 рода, Якоби, Лежандра, Лагерра, Эрмита, Фурье, ДКП, ДСП. Она работает и рисует красивые картинки. [показать пару картинок и закончить]. Кстати, по поводу того, а какой же базис лучше? Вообще они все дают очень похожие результаты… Пока что «лучше» всех Чебышев 1-го рода. А что вообще такое «лучше»? «Лучше» — чисто умозрительно это «больше соотношение сигнал/шум» (в результатах). Как измерить? Ну, например, при одинаковых параметрах окон и глубине разложения подобрать eps такое, чтобы общее количество «похожих» участков было примерно равно, и посчитать, например, среднюю длину повторов. Можно и медиану тоже. Чем больше, тем лучше — мы ведь хотим найти как можно более длинные повторы. Начинали реализовывать с Чебышева 1-го рода, потом пробовали Лежандра, потом думали, что Чебышев 2-го рода произведёт революцию и всё будет гораздо лучше, т.к весовая функция выпуклая, центр отрезка учитывается сильнее, края меньше. Революции не произошло, результаты сильно похожие на 1-го рода, местами получше, местами похуже. Формально — похуже. Дальше есть табличка с «попугаями» по разным базисам. Тестовые данные — часть генома мыши (не спрашивайте какая, я не знаю) длиной 1.5 млн нуклеотидов. Сравнение приводилось при примерно одинаковых количествах найденных участков, «подозрительных» на повтор — в районе 5000. При выбранных настройках минимальная длина участка, подозрительного на повтор — 3500 нуклеотидов. Какие выводы? Лидирует Чебышев 1 рода. Базисы ДКП, ДСП и Фурье дают до жути похожие на него, практически идентичные, результаты. С небольшим отставанием за ними следует Лежандр, за ним — Чебышев 2 рода, а базисы Эрмита и Лагерра не подходят для поиска повторов ''вообще — ''что есть логичный факт, т.&nbsp;к. они оба работают на бесконечном интервале либо (0, +бесконечность), либо от — до + бесконечности. Вариантов значения медианной длины было всего 2: 3500 (минимально возможная) или 10000, она отражает, фактически, чистое количество шума — мелких отрезков, и гласит, что приемлемый уровень шума дают… Ясно кто.
+
Вот где-то примерно это всё и было реализовано. Есть относительно простая программа, есть относительно хорошая библиотека для абстрагирования от деталей реализации конкретных базисов, есть сами базисы — Чебышева 1 и 2 рода, Якоби, Лежандра, Лагерра, Эрмита, Фурье, ДКП, ДСП. Она работает и рисует красивые картинки.
  
<tab sep="tab" border="1">
+
== Часть статьи ==
- 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'''
+
</tab>
+
  
 
Для реализации программы поиска повторов с помощью ОСАМ был выбран язык C++. Такой выбор обусловлен сущностью процесса разложения функций, позволяющей с помощью объектно-ориентированного подхода разделить функционал на общий и зависящий от конкретного ортогонального базиса. Общий функционал — это функции подсчёта весовых коэффициентов, подсчёта интеграла на сетке Гаусса, подсчёта матрицы Грама заданного базиса, нормирования заданного базиса, интерполяции сигнала на заданную сетку, и воссоздания изначального сигнала по коэффициентам разложения. К базисо-зависимому функционалу относятся функции подсчёта сетки, весовых коэффициентов, и самих значений функции. Также такой подход, кроме всего прочего, даёт возможность оптимизировать части функционала отдельно.
 
Для реализации программы поиска повторов с помощью ОСАМ был выбран язык C++. Такой выбор обусловлен сущностью процесса разложения функций, позволяющей с помощью объектно-ориентированного подхода разделить функционал на общий и зависящий от конкретного ортогонального базиса. Общий функционал — это функции подсчёта весовых коэффициентов, подсчёта интеграла на сетке Гаусса, подсчёта матрицы Грама заданного базиса, нормирования заданного базиса, интерполяции сигнала на заданную сетку, и воссоздания изначального сигнала по коэффициентам разложения. К базисо-зависимому функционалу относятся функции подсчёта сетки, весовых коэффициентов, и самих значений функции. Также такой подход, кроме всего прочего, даёт возможность оптимизировать части функционала отдельно.
Строка 44: Строка 32:
 
  ippsMulC_64f(wn, -1, tn, n);
 
  ippsMulC_64f(wn, -1, tn, n);
 
  ippsSqrt_64f_I(tn, n);
 
  ippsSqrt_64f_I(tn, n);
 +
 +
=== Сравнение ОНБ ===
 +
 +
Учитывая, что поиск повторов может осуществляться по выбору с использованием любого из ортогональных базисов, и что в библиотеке функций разложения их было реализовано 9 различных - базис Чебышева 1 рода, базис Чебышева 2 рода, дискретные косинусное и синусное преобразования, базис Фурье, базис Лежандра, базис Лагерра, базис Якоби и базис Эрмита - очевидным образом встаёт вопрос: а какой же из них "лучше" в задаче поиска повторов в последовательностях? А кроме того, каковы в целом критерии качества, по которым требуется производить сравнение базисов?
 +
 +
Очевидным подходом к данному вопросу является критерий «максимум соотношения сигнал/шум в найденных в итоге повторах».
 +
 +
Другой вариант - максимум средней длины найденных подобных участков, т.к. цель поиска повторов заключается в том, чтобы найти как можно более длинные подобные участки. Как можно оценить эту длину? Опишем простейший подход. Во-первых, нужно выбрать ширину скользящих окон и глубину разложения и выбрать некоторые тестовые данные, содержащие широкий спектр различных повторов - здесь хорошо подходит часть реальной ДНК-последовательности. Далее, используя различные базисы и подбирая порог сравнения (<m>\varepsilon</m>) такой, чтобы общее число найденных подобных участков было приблизительно равно, подсчитывать среднюю длину найденных подобных участков. Как вариант — можно вычислять медианное значение.
 +
 +
В процессе реализации программы вначале был выбран базис Чебышева 1-го рода; потом пробовали базис Лежандра. Потом было высказано предположение о том, что базис Чебышева 2-го рода произведёт «революцию» по той причине, что имеет выпуклую весовую функцию и сильнее учитывает центр сравниваемого отрезка, чем края, но революции не произошло, результаты базиса Чебышева 2-го рода сильно похожи на базис Чебышева 1-го рода, и даже немного хуже, в том числе и по средней длине найденных повторов.
 +
 +
Ниже приводится табличка с замерами средней длины найденных повторов на различных базисах и части генома мыши длиной 1.5 млн нуклеотидов в качестве тестовых данных. Сравнение производилось при приблизительно равных количествах найденных «подобных» участков — 5000. При выбранных настройках минимально возможная найденная длина подобного участка — 3500 нуклеотидов.
 +
 +
<tab sep="tab" border="1">
 +
- 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'''
 +
</tab>
 +
 +
Каковы выводы? По средней длине повтора лидирует базис Чебышев 1 рода, а базисы ДКП, ДСП и Фурье дают чрезвычайно похожие на него, практически идентичные, результаты. С небольшим отставанием следует базис Лежандра, далее — базис Чебышева 2 рода, а базисы Эрмита и Лагерра для поиска подобных участков не подходят вообще, чему есть простое математическое обоснование — оба они действуют на бесконечной полупрямой — либо <m>(0, +\inf)</m>, либо <m>(-\inf, +\inf)</m>. Вариантов значения медианной длины при этом было всего 2: 3500 (минимально возможная) или 10000. Медианная длина в данном случае отражает, фактически, «чистое» количество шума — мелких отрезков, и гласит, что приемлемый уровень шума дают базисы Чебышева 1 рода, ДКП, ДСП, Фурье и Лежандра.
  
 
[[Категория:Учёба]]
 
[[Категория:Учёба]]

Версия 14:33, 29 августа 2009

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

В данной работе предлагается алгоритм поиска длинных разнесенных повторов. Лежащий в основе алгоритма обобщенный спектрально-аналитический метод, позволяет значительно ускорить процесс анализа последовательности за счет применения средств распаллеливания и векторизации. Также предлагается матрица спектральной гомологии генетических последовательностей. Близкая к точечной матрице гомологии, она предоставляет более быстрый инструмент для сравнительного анализа и визуализации внутренней структуры больших отрезков геномов (порядка 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, потому и «индекс». Далее пробежаться по всем полученным описаниям (по индексу) обеих последовательностей (или одной и той же последовательности) и сравнить попарно все пары описаний (на похожесть). А что такое похожесть? Критериев похожести можно выработать массу, среди них можно найти устойчивые к масштабу и т. п., однако у нас всё довольно просто: , где . Такое вот «нормированное L2-расстояние». Здесь, кстати, можно выиграть от т. н. «принципа дискриминантности», который гласит очевидную вещь: что если уже при k < n, то суммирование можно не продолжать, т. к. меньше сумма квадратов уже не станет. Итак, что мы получим от этого сравнения? Мы получим приближённые «близости» участков ДНК. Крупных или мелких, более или менее точное сравнение — это уже как захотим — для этого можно варьировать параметры. Задаём порог, можем пробежаться по результатам и сразу выявить «подозрительные на повторы» участки. Это есть важно, т. к. больше не нужно всё время искать повторы ВЕЗДЕ: сначала достаточно выявить крупные относительно похожие участки, а потом можно «увеличить масштаб» и выявить (или не выявить) точные координаты повторов. Кстати, единственное, для чего подход почти не подходит — для выявления «абсолютно точных» координат повторов. Это уже в «подозрительных» областях можно делать стандартными методами. Например, diffоподобным алгоритмом. :-)

Вот где-то примерно это всё и было реализовано. Есть относительно простая программа, есть относительно хорошая библиотека для абстрагирования от деталей реализации конкретных базисов, есть сами базисы — Чебышева 1 и 2 рода, Якоби, Лежандра, Лагерра, Эрмита, Фурье, ДКП, ДСП. Она работает и рисует красивые картинки.

Часть статьи

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

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

Реальный выигрыш в производительности засчёт чисто программной оптимизации достигает 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);

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

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

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

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

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

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

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

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