Быстрые алгоритмы поиска уникальных наборов в двух очень длинных последовательностях текста
Мне нужно сравнить последовательности ДНК X и Y хромосом и найти паттерны (составленные из примерно 50-75 пар оснований), которые являются уникальными для Y-хромосомы. Обратите внимание, что эти части последовательности могут повторяться в хромосоме. Это должно быть сделано быстро (BLAST занимает 47 дней, нужно несколько часов или меньше). Существуют ли какие-либо алгоритмы или программы, особенно подходящие для такого сравнения? Опять же, здесь ключевым является скорость.
Одна из причин, по которой я поставил это на SO, заключалась в том, чтобы получить представление от людей за пределами конкретной предметной области, которые могут выдвигать алгоритмы, которые они используют при сравнении строк при повседневном использовании, которые могут к нашему использованию. Так что не стесняйтесь!
4 ответа
- Построить суффиксное дерево S в последовательности X.
- Для каждой начальной позиции i в последовательности Y найти строку Y[i..i+75] в S. Если совпадение не может быть найдено, начиная с позиции i (т. е. если поиск не удается после того, как j < 75 совпадающих нуклеотидов), то вы нашли строку длины j, уникальную для Y.
- Наименьшая такая строка среди всех начальных позиций i является самой короткой уникальной строкой (или просто остановитесь после того, как найдете любую такую строку, если вы не заинтересованы в минимизации длины).
Общее время: O(|X| + m|Y|), где m — максимальная длина строки (например, m = 75).
Вероятно, существуют еще более эффективные алгоритмы, основанные на обобщенных деревьях суффиксов.
Я предполагаю, что у вас есть один X и один Y для сравнения. Объедините их, разделив символом маркера, который не появляется ни в одном из них, чтобы сформировать, например. XoY. Теперь сформируйте http://en.wikipedia.org/wiki/Suffix_array за линейное время.
Вы получаете массив указателей на позиции в объединенной строке, где указатели расположены так, что подстроки, на которые они указывают, появляются в массиве в алфавитном порядке. Вы также получаете массив LCP, задающий длину самого длинного общего префикса, разделяемого между суффиксом и суффиксом непосредственно перед ним в массиве, который является суффиксом, который сортируется чуть меньше, чем он. На самом деле это самый длинный общий префикс, разделяемый между этой позицией и ЛЮБОЙ подстрокой меньше ее, потому что все, что имеет более длинный общий префикс и меньше строки[i], будет сортироваться между ней и текущей строкой[i - 1].
Вы можете сказать, на какую исходную строку указывает указатель, по его позиции, потому что X предшествует Y. Вы можете разрезать массив на чередующиеся части указателей X и Y. Длина общего префикса, разделяемого между pos[i] и pos[i - 1], равна lcp[i]. Длина префикса, разделяемого между pos[i] и pos[i-2], равна min(lcp[i], lcp[i-1]). Поэтому, если вы начинаете со значения Y непосредственно перед диапазоном X, вы можете вычислить количество символов префикса между этим Y и всеми X по очереди, уменьшая раздел, выполняя минимальную операцию на каждом шаге. Точно так же вы можете вычислить количество символов префикса, разделяемых между всеми этими X и Y, которые появляются следующими в массиве суффиксов, по цене одной минуты за X.То же самое, конечно, для диапазонов Y. Теперь сделайте максимальное значение для каждой записи, чтобы определить самый длинный префикс, общий для каждой позиции в X (или Y) и любой позиции в Y (или X).
Я думаю, вам нужны подстроки внутри X или Y, которые имеют небольшие префиксы, общие для них и любой другой подстроки другого пола, потому что строки на один символ длиннее этого, начиная с той же позиции, не появляются в другом поле. вообще. Я думаю, что после того, как вы выполнили приведенные выше вычисления min(), вы можете извлечь N наименьших подстрок префикса, используя кучу, чтобы отслеживать N наименьших записей. Я думаю, что все здесь требует времени, линейного по |X| + |Y| (если только N не сравнимо с |X| или |Y|).
У меня есть интересный ответ, он будет технологический. Основная идея заключается в том, что сравнение подпоследовательностей должно выполняться на GPU, поскольку GPU современных видеокарт представляет собой высокопараллельную вычислительную среду (подобно маленькому суперкомпьютеру). Таким образом, мы можем закодировать пару оснований как один пиксель, учитывая, что Х-хромосома составляет 154 миллиона пар - мы получаем изображение для Х-хромосомы, состоящее из 154 миллионов пикселей - размер этого изображения будет около 500 МБ. Для Y-хромосомы получаем изображение размером 160 МБ. Таким образом, эти два изображения (500+160) МБ могут быть очень эффективно обработаны спусковой видеокартой. (Просто возьмите видеокарту с >= 1 ГБ видеопамяти).
Следующим шагом является написание программы для графического процессора, возможно, с помощью Pixel Shader, Cuda или OpenCL
Программа для графического процессора будет простой – в основном она будет сравнивать 50-75 соседних пикселей изображения Y-хромосомы ко всем пикселям изображения X-хромосомы. Таким образом, эта программа для графического процессора должна иметь максимум 75 * 154 миллиона операций, которые будут вычисляться на современном графическом процессоре примерно за ЧАС. Потому что все подпоследовательности Y будут проверяться параллельно.
надеюсь поможет
В этом документе могут быть некоторые альтернативы адаптации BLAST для повышения его производительности (подробно разделение проблемного пространства AFAIKS).