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

Версия 17:29, 22 августа 2011

Имеем алгоритм поиска приближённых повторов в ДНК-последовательностях:

• Скользим по последовательности окном размера W и каждую букву заменяем на количество пуринов (A, G, U — аденин, гуанин, урацил) / пиримидинов (C, T — цитозин, тимин) в этом окне (W-окрестности буквы). Полученную числовую последовательность называем «профилем последовательности».
• Дальше, начиная с каждой из (N-A)/S позиций, где N — длина последовательности, A — второе окно, «окно аппроксимации», S — шаг аппроксимации (расстояние между соседними начальными позициями), выцепляем A чисел. Из этих A чисел методом ближайшего соседа выбираем K (K сильно меньше A, например, K=16, A=500), равномерно или в соответствии с сеткой разложения выбранного базиса. Эти K чисел раскладываем по выбранному базису, получаем в каждой из (N-A)/S позиций K коэффициентов, характеризующих часть последовательности длиной A.
• Получаем вместо длинной строки кучу векторов по K коэффициентов. Проделываем это с обеими сравниваемыми последовательностями (или с одной, если сравниваем её саму с собой).
• Сравниваем все вектора попарно и если они ближе, чем на ε, считаем, что в этой области последовательности, похоже, повтор.

Обоснование профиля:

Теперь представляем, что хотим это ускорить. Для простоты, за расстояние между векторами берём максимум модуля разности, а базис используем тригонометрический (Фурье). Хотим, по сути, уйти от сложности O(N²).

И сразу ловим 2 облома:

Первый: природа коэффициентов разложения — распределение похоже на нормальное, и что ещё хуже, первые коэффициенты сильно больше последних по модулю. Сильно — это примерно на 1-2 порядка. Причины две — во-первых, в реальных данных низкочастотные гармоники всегда больше, во-вторых, в нашем случае высокочастотные гармоники приглушены нами же — исходный вектор длины K => разница между двумя соседними значениями в нём по модулю не больше A/K, так как разница между двумя соседними значениями в профиле последовательности либо 0, либо 1, либо −1. Следовательно, максимум i-го коэффициента примерно (K+1-i)*A/K.

Второй: при линейном поиске и сравнении мы же не сравниваем вектора целиком. Мы их сравниваем последовательно от 1-ой координаты к K-ой, то есть в порядке убывания абсолютных значений коэффициентов, а значит, в среднем и разности. То есть большая часть неподходящих нам векторов отсекается вообще по 1-ой координате. Мало — по второй, ещё меньше — по третьей, и т. п. В реале получается, что по 8-16-ым не отсекается вообще почти ничего, а по где-то 7-и из них — совсем ничего. В то же время для полной проверки подходящего вектора нужно сделать точно K сравнений. То есть, 90 % неподходящих векторов дают часть сложности, сравнимую с полной проверкой оставшихся 10 %. Соответственно, даже если мы изобретём индекс, который позволит нам сразу отсечь почти все 90 % — мы убьём только примерно половину сложности. Потому что индекс никогда не выдаст сразу точно 10 %, а выдаст, скажем, 15 %, которые ещё нужно проверить и ещё 5 % отсечь. Соответственно и максимальное ускорение с помощью любого индекса зависит в основном от процента подходящих векторов и размерности. Чем больше потенциально подходит — тем хуже, чем больше размерность — тем хуже.

Мысли:

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