Поиск одинаковых участков в нуклеотидных цепочках с помощью индексации

Данный пост написан в рамках конкурса Accelerate Your Code 2012.

Самый простой и самый не быстрый из методов поиска одинаковых участков в нуклеотидных цепочках - "наивный" перебор строк со смещением. Он хорошо подходит для коротких цепочек, так как не требует предварительной обработки и дополнительных объемов памяти.
Но требование O(nm) по времени в общем случае нас не устроит, потому будем подбирать другие варианты.

Подготовка данных
Для начала подумаем об ускорении используя подготовку данных.
Так как мы будем проводить поиск в нуклеотидных цепочках, мы можем ограничить алфавит или же понизить битность на один символ. Всего таких символов у нас 4, а именно АЦГТ (ACGT). Это дает возможность использовать всего 2 бита вместо 8 (ANSII). Каким же образом построить соответствие? Если у нас нет дополнительных наполеоновских планов, можно в любом, но если вспомнить что эти нуклеотиды составляют пары (А-Т, Ц-Г) лучше на будущее провести такие аналогии:

    • A = 002=0

    • T = 112=3

    • C = 012=1

    • G = 102=2



Таким незамысловатым способом мы бы могли решить одну маленькую проблему - проверку пары нуклеотидов на комплементарность. В общем случае на это приходилось бы 8 проверок, в нашем же - только одна, проверка суммы значений на равенство 3. Уже неплохое ускорение в 8 раз. Правда времени потратим О(N). Ну да ладно - не смертельно, все-равно надо от мусора (перевод строк или пробелы) избавляться при считывании из неподготовленных данных.

Теперь подумаем вот о чем: а почему бы нам каким-нибудь образом сравнивать не по одной паре, а сразу по несколько. Ну что же, на помощь к нам придут хеш-функции для задания соответствия некоторой цепочки нуклеотидов длиной W и целого натурального числа. Таким образом, если два значения хеш-функции не совпадают, значит и цепочки гарантированно разные.
С равенством возможны варианты. Если использовать хеш-функции без коллизий (например, представление цепочки как число по основанию 4), мы также сможем гарантировать и идентичность цепочек. Такая хеш-функция будет цикличной, что позволит использовать предыдущее рассчитанное значение при смещении на 1 нуклеотид. Но тут возникают некоторые трудности с хранением данных - в переменную размером К бит мы сможем сохранить только К/2 нуклеотида. Затраты на подготовку О(N), при использовании цикличных хеш-функций типа такой:

int32_t W_ = W - 1;

for (i = 0; i < W_; i++) {

	key = ((key << 2) + Seq[i]) & mask;

}

for (i = W_, p = 0; i < Seq_l; i++, p++) {

	key = ((key << 2) + Seq[i]) & mask;

	Table_Key[p] = key;

}


Для снятия ограничения на максимальную длину цепочки можем, например, использовать несколько различных хеш-функций с коллизиями для уменьшения вероятности совпадений коллизий. Или просто проверить совпадение на идентичность.

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

Первый способ
После превращения цепочки в массив хешей, каждый элемент записать как пару хеш-позиция. Дальше нам предстоит отсортировать этот массив по порядку значений хешей. Затраты по времени от O(N*Log(N)).
На этом можно и остановиться, а можно для ускорения поиска по нему преобразовать в двумерный массив. В нем элемент первого порядка будет соответствовать значению хеш-функции, а элементы второго порядка - позициям. При необходимости можно отсортировать и массив позиций - необходимо, если надо искать свертку одной строки на себя же (комплементарная связь, а не идентичность).
Такой способ подойдет в том случае если значения хеш-функции находятся в большом диапазоне величин, но используемых намного меньше от суммарного.
Поиск нужного значения хеш-функции возможен либо через бинарный поиск, либо через индексацию диапазонов с последующим бинарным поиском.

Второй способ
По-другому можем поступить в случае, если хеш-функция имеет ограниченный в разумных рамках диапазон величин.
Для примера используем хеш-функцию без коллизий, приведенную выше. В ней диапазон величин будет зависеть только от длины цепочки W, для которой мы строим соответствие. И этот диапазон составит от 0 до 4W-1. Таким образом мы сможем воспользоваться уже "линейной" сортировкой с временем О(N) для построения двумерного массива соответствий значений хеш-функции и позиции.
Смысл простой. Создаем два массива: линейный S размером 4W и двумерный T размером по первому индексу также 4W, а по второму как получится (см. ниже).
На первом этапе воспользуемся массивом S для построения гистограммы вхождений значений хеш-функции в массив за время O(N).
На втором этапе выделяем память для подмассивов Т на количество элементов, указанных в соответствующей ячейке S.
Для окончательного заполнения будем использовать массив S как счетчик заполненности. В Т первый индекс будет соответствовать значению хеш-функции, второй берем из S, значение - позиция хеша в начальном массиве. Время, затраченное на заполнение, также О(N).
При желании и необходимости этот массив можно "причесать" - для хеш-функций соответствующих однородным строкам (АА..АА и подобных) выкинуть часть позиций по некоторым правилам.
В итоге за линейное время получим двумерный массив-индекс Т для быстрого определения положений подстроки длиной W по части (так же длиной W) заданного шаблона.
Затраты по времени общей подготовки О(N). Затраты на дополнительную память зависят от выбранной длинны W: массив хешей будет больше начальной строки в 2-4-... раза в зависим от выбранного типа данных (16, 32 и больше бит). Массив S будет содержать 4W элементов размером от 8 до 32+ бит в зависимости от размеров строк. Массив Т будет занимать по объему примерно такой же объем как S (первый уровень индексации) и начальный массив суммарно.
Поиск нужного значения хеш-функции будет происходить за время О(1), но накладные затраты (массив S и первый уровень массива Т) могут быть больше чем полезная нагрузка (списки позиций).

Поиск совпадений
Для коротких строк достаточно построение массива хеш-функций подстрок длиной W, с последующим расширением совпадающим участков по 1 или W нуклеотидов.
Для длинных строк возможно будет иметь смысл сравнивать два массива индексов - только те позиции подстрок, хеши которых будут присутствовать в обоих.
Разделение строк на длинные и короткие разное для разных случаев. И зависит от выигрыша в скорости работы.

Сравнение с другими методами
Казалось бы, что приведенный выше метод не очень хорош, так как дает только "точки соприкосновения" двух подстрок. И дальнейшее расширение до необходимой минимальной длины может "съесть" все ускорение. Но:

    • иногда минимальная необходимая длина достаточно маленькая (W = Wmin), что бы массив индексов был не столь уж и большим.

    • алгоритм поиска с помощью суффиксных деревьев/массивов точно потребует объем в много раз превышающий размер входных данных, при сравнимых скоростях поиска




Масштабируемость
В наше время когда корабли... О чем это мы... ах да... в наше время когда в каждом ПК минимум по 2 ядра, можно подумать и о параллельной обработке.
Сразу оговорюсь, что использовал OpenMP. Ибо до этого параллельным программированием не пользовался, а данная технология оказалось весьма и весьма... простой в понимании и в реализации...
Рассмотрим основные этапы алгоритма и мое решение по распараллеливанию.

1. Создание массива хешей подстрок
Ну тут все довольно просто - делим начальную строку на участки, каждый участок своему потоку. Учитываем при этом, что на построение N хешей необходимо N+W-1 нуклеотидов. То есть нарезка происходит с перехлестами. Запись готовых результатов тоже разделена в памяти. Как дополнительные вычисления на каждый поток будет просчет первых W-1 хешей недостаточной длины. Но это явно мелочи по сравнению с выгодой в скорости. Но и кидаться в крайности тоже не стоит - должно выполняться условие ограничение разбиения P*W << N, где P - количество участков.

2. Индексация подстрок
Эту проблему я решил так же разбиением начального массива уже готовых хешей на несколько участков. Для разделения участков записи в конечных индекс воспользовался массивами счетчиками для каждого потока, с последующим суммирование нужных и использованием значения как сдвиг начала при записи. Основной проигрыш - использование дополнительных массивов под гистограммы и доступ к нужной ячейке с помощью сдвига. Тут тоже стоит себя ограничивать в бездумном дроблении - каждый массив под гистограмму имеет 4W элементов от 8 до 64 бит. Что при большом количестве потоков может занимать значимый объем памяти.

3. Поиск совпадений
Тут еще проще - всю работу можно разделить на задания которые между собой не сильно связаны. Таким заданием, в случае сравнения массива и индекса, будет проверка совпадений для подстроки из шаблона и списка позиций возможных совпадений из индекса. Их будет приличное количество - динамическое распределение работы между потоками позволит адекватно балансировать нагрузку.
Для сравнения двух индексов неделимым заданием может быть как сравнения списка позиций для совпадающих хешей, так и более мелким - просчет результата для пары из возможного списка комбинаций позиций. Это решит проблему перекоса в случае если в обоих индексах некоторые значения хеш-функций будут содержать значительно большее количество позиций, чем остальные.
Еще вариант при котором будем отсекать некоторые варианты, если необходимо найти совпадение подстрок длиной намного больше чем 2*W. При этом возьмем две подстроки длиной W на таком расстоянии L, так что бы 2*W + L = Wmin. Дальше будем сравнивать массивы позиций из индекса для хешей этих двух строк. Если расстояние между определенными элементами будет также равно L - возможно совпадение больших подстрок, что подтвердится полной проверкой. Оценка по времени такой проверки на совпадения О(Log(n1)*Log(n2)), где ni - длины списков.
Чтобы не использовать критические секции, для каждого потока можно использовать отдельные массивы для хранения найденных результатов. С последующим слиянием и сортировкой/выборкой за время O(К + К*Log(К)).

Общее время алгоритма стремится к О(N + K + K*Log(K)), где N общий объем входных данных, К - объем не фильтрованных выходных данных.

Данный пост написан в рамках конкурса Accelerate Your Code 2012.
Спасибо за внимание. ЦУ буду рад прочитать в комментариях.

Для получения подробной информации о возможностях оптимизации компилятора обратитесь к нашему Уведомлению об оптимизации.