Об ускорении поиска повторов ОСАМ-методом
Имеем алгоритм поиска приближённых повторов в ДНК-последовательностях:
- Скользим по последовательности окном размера 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 коэффициентов. Проделываем это с обеими сравниваемыми последовательностями (или с одной, если сравниваем её саму с собой).
- Сравниваем все вектора попарно и если они ближе, чем на ε, считаем, что в этой области последовательности, похоже, повтор.
Во-первых, целесообразно оказывается использовать 3, максимум 5 коэффициентов разложения, так как по высшим гармоникам фильтрации почти не происходит. В простейшем случае 3 коэффициента - среднее, проекция на синус, проекция на косинус.
Теперь представляем, что хотим это ускорить. Для простоты, за расстояние между векторами берём максимум модуля разности, а базис используем тригонометрический (Фурье). Хотим, по сути, уйти от сложности O(N²).
И сразу ловим 2 облома:
Первый: природа коэффициентов разложения — распределение похоже на нормальное, и что ещё хуже, первые коэффициенты сильно больше последних по модулю. Сильно — это примерно на 1-2 порядка. Причины две — во-первых, в реальных данных низкочастотные гармоники всегда больше, во-вторых, в нашем случае высокочастотные гармоники приглушены нами же — исходный вектор длины K => разница между двумя соседними значениями в нём по модулю не больше A/K, так как разница между двумя соседними значениями в профиле последовательности либо 0, либо 1, либо −1. Следовательно, максимум i-го коэффициента примерно (K+1-i)*A/K. И всё это - распределено околонормально. Таким образом, разбиение пространства на какие-то области (построение индекса) уже становится геморройно - "дельных" измерений с большим относительно размера окрестности поиска диапазоном мало (1-2-3 максимум), бОльшая часть векторов в середине, середина маленькая, то есть при поиске придётся просматривать довольно много лишних векторов.
Второй: при линейном поиске и сравнении мы же не сравниваем вектора целиком. Мы их сравниваем последовательно от 1-ой координаты к K-ой, то есть в порядке убывания абсолютных значений коэффициентов, а значит, в среднем и разности. То есть большая часть неподходящих нам векторов отсекается вообще по 1-ой координате. Мало — по второй, ещё меньше — по третьей, и т. п. В реале получается, что по 8-16-ым не отсекается вообще почти ничего, а по где-то 7-и из них — совсем ничего. В то же время для полной проверки подходящего вектора нужно сделать точно K сравнений. То есть, 90 % неподходящих векторов дают часть сложности, сравнимую с полной проверкой оставшихся 10 %. Соответственно, даже если мы изобретём индекс, который позволит нам сразу отсечь почти все 90 % — мы убьём только примерно половину сложности. Потому что индекс никогда не выдаст сразу точно 10 %, а выдаст, скажем, 15 %, которые ещё нужно проверить и ещё 5 % отсечь. Соответственно если даже мы сделаем оптимальный неточный индекс, максимум ускорения с его помощью зависит только от процента подходящих векторов и размерности. Чем больше потенциально подходит — тем хуже, чем больше размерность — тем хуже.
Мысли:
- Разбить последовательность на слова и проиндексировать их с учётом возможных сдвигов.
- Использовать вместо классических базисов целочисленный вэйвлетоподобный базис, порождаемый из замены двух чисел на их сумму и разность.
- Но использовать сами вэйвлеты, то есть когда разности не раскладываются повторно, не совсем верно.
