Тема создана как дискуссионная. В ней задано несколько вопросов, поэтому можно отвечать как на все, так и на один из них. На некоторые вопросы я дал свой ответ, не факт, что правильный. Если будет большой интерес, можно попытаться организовать чат для дискуссии.
Наиболее часто встречающаяся задача в биоинформатике: поиск подстроки в строке. В большинстве случаев это реализуется одним из стандартных способов, не требует много ресурсов и не заслуживает детального рассмотрения. Однако, есть класс задач, который мог бы быть, на мой взгляд, ускорен/упрощен применением нагруженных деревьев.
В данной теме рассмотрим только точный поиск.
Расшифровка терминологии:
референс: длинная ДНК-последовательность, как правило, содержащая расшифровку части хромосомы или хромосому/геном целиком. Представляет собой (как правило) строку, состоящую из A, G, T, C (обозначающих разные нуклеотиды) и N - символ означающий либо любой нуклеотид, либо отсутствие точной информации. Пример: AAGGCCTCGTCNTTNNGCCC...
паттерн (или рид): короткая ДНК-последовательность, строка, содержащая те же символы, что и выше, длиной от 1 (бесмыссленно, конечно, но бывает) и до разумных величин в несколько тысяч. Мы, в основном, работаем с длинами от 30 до 300 букв, но бывают и другие задачи.
палиндром - комплементарная строка, т.е. перевернутая по особым правилам последовательность (читаем строку с конца, заменяем буквы по правилу: A -> T, C -> G, G -> C, T -> A, N -> N)
Для разбора полетов в качестве референса будем рассматривать 8ю хромосому человека (A: 43333530, C: 29030173, G: 29103787, T: 43300646, N: 370500 штук, всего 145138636 букв). Риды, для простоты, будут длиной 100 букв.
Сразу скажу, что код я буду писать на Delphi. По большому счету, мне все равно - Delphi, c++, c# или java, но Delphi мне привычнее.
Вопрос 1. Создаем структуру дерева. Классы или связанный список.
2 варианта реализации дерева:
Класс
type
AGTCIndex = (A, C, G, T);
TTrie = class
fAmI: integer; // кто я - на всякий случай
fParent: TTrie; // ссылка на родительский узел
fChildren:array[AGTCIndex] of TTrie; // ссылка на узлы потомков
fCount: array of longint; // для адресации 145 миллионов хватит. Вот если бы работали с геномом секвойи, стоило бы использовать int64 :)
<...>
end;
или связанный список
type
AGTCIndex = (A, C, G, T);
pTrie = ^TTrie;
TTrie = record
fAmI: integer; // кто я - на всякий случай
fParent: pTrie; // ссылка на родительский узел
fChildren:array[AGTCIndex] of pTrie; // ссылка на узлы потомков
fCount: array of longint;
end;
Как мне видится, для нашей задачи связанный список подходит гораздо больше: меньше затрат памяти, да и работа будет быстрее (не будем терять время на лишних вызовах). Разумеется, для удобства обернём список в класс.
Вопрос 2. Затраты памяти
В стандартных методах потратится 145 с копейками мегабайтов памяти - столько будет занимать строка. Можно ухитриться и использовать в 2 раза меньше с помощью полубайтового представления (описание и способы работы).
В случае использования дерева затраты будут не менее, чем 100*145138636*16 или не менее 230 гигабайтов. Всё, овчинка не стоит выделки. Максимум, что можно потратить, это 32 гигабайта. Делаем вывод, что подобный подход применим при условии, что референс имеет длину не более 16 миллионов.
В принципе, полубайтовый подход и здесь реализуем. Он даёт даже большую, нежели 2 раза, выгоду. Но, боюсь, неудобств будет больше, чем выгоды. если есть альтернативное мнение - прошу озвучивать.
Вопрос 3. Быстродействие точного поиска
В стандартном подходе найти все включения паттерна в референс занимает примерно |length(референс)*ln(length(паттерн)|. Если распараллелить, то примерно в 1.6 раз быстрее на каждые удвоение ядер.
В случае использования дерева только на построение уйдет, по моим подсчетам, в лоб минимум O(100*length). Используя рамку считывания, можно сократить до O(30*length). Параллельный подход также ускорит процесс в 1.6 раз на каждое удвоение ядер.
Итого, имеем: на создание дерева уйдет где-то в 500 - 1000 раз больше времени, чем на поиск одного паттерна стандартными методами. Но, поскольку поиск по уже сформированному дереву гораздо быстрее, плюс, позволяет параллельно получать дополнительные данные, имеем: Выгода от создания дерева достигается в том случае, когда референс используется многократно, причем для разных паттернов. Оценка: не менее 500 раз. Область применения: точное картирование, потоковая обработка данных с прибора и т.п.
Правда, встает небольшая проблема комплементарного поиска. Фактически, нам придется идти от потомков к родителям, что будет несколько медленнее.
Если тема заинтересовала, можно перейти к нечеткому поиску, что, в целом, гораздо актуальнее и интереснее.
Update 1
Создал комнату в чате
Update 2
**Относительно побитовых представлений. Рассуждение вслух (к деревьям отношения почти не имеет) **
Почему не сэкономить память и, зачастую, время, и не перейти к полубайтовому, а то и четвертьбайтовому представлению?
- Точно не знаю. Не встречал пока ни одной GNU-софтины (насчет платных - см. ниже), которая бы использовала подобный подход. Все открещиваются совместимостью и нежеланием писать сложные рутины для побитовой работы. А в последнее время еще и отвечают - "память подешевела".
Моё мнение, почему все не используют побитовые подходы:
Совместимость не особо важна. В конце концов, преобразование без потерь: всегда можно преобразовать обратно. Но! Часто одни и те же рутины используются не только для поиска внутри ДНК, но и внутри белковых последовательностей, генных зависимостей и т.п. А там - и латинского алфавита будет мало.
AGTC - это очень частный случай. Есть еще метилированные нуклеотиды, дидезокси- нуклеотиды, "подозрительные" и т.п. Т.е. "алфавит" гораздо больше, нежели 4 буквы. Поэтому рутины создают "для всего".
И по первой ссылке выше я приводил пример работы с полубайтовым представлением, которое "разработал" в 2004 году, а полноценный геномный тритер (первый под Windows - чуток горжусь :) ) был создан в 2006 году. Тогда памяти было мало, а процессоры весьма тупили. Основная идея: A = 0001b, C = 0010b, G = 0100b, T = 1000b. Тогда очень удачно накладываются маски по поиску сразу нескольких вариантов: N = 1111b, R (A или G) =0101b, Y (C или T) = 1010b и т.п. А сдвиг на полубайт сразу у большого слова (64 или 128 бит) реализован в SSE2, что позволяет обрабатывать одновременно до 32 букв сразу.
Так вот касательно этой реализации у меня есть подозрение, что этот подход используется, как минимум, в одном коммерческом пакете :) Впрочем, мне не жалко, тем более, что подход лежал на поверхности.
Update 3
Опубликовал 2ю часть вопроса: Нагруженные деревья в биоинформатике. Часть 2. Нечеткий поиск