15

Тема создана как дискуссионная. В ней задано несколько вопросов, поэтому можно отвечать как на все, так и на один из них. На некоторые вопросы я дал свой ответ, не факт, что правильный. Если будет большой интерес, можно попытаться организовать чат для дискуссии.

Наиболее часто встречающаяся задача в биоинформатике: поиск подстроки в строке. В большинстве случаев это реализуется одним из стандартных способов, не требует много ресурсов и не заслуживает детального рассмотрения. Однако, есть класс задач, который мог бы быть, на мой взгляд, ускорен/упрощен применением нагруженных деревьев.
В данной теме рассмотрим только точный поиск.

Расшифровка терминологии:

референс: длинная ДНК-последовательность, как правило, содержащая расшифровку части хромосомы или хромосому/геном целиком. Представляет собой (как правило) строку, состоящую из 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-софтины (насчет платных - см. ниже), которая бы использовала подобный подход. Все открещиваются совместимостью и нежеланием писать сложные рутины для побитовой работы. А в последнее время еще и отвечают - "память подешевела".

  • Моё мнение, почему все не используют побитовые подходы:

    1. Совместимость не особо важна. В конце концов, преобразование без потерь: всегда можно преобразовать обратно. Но! Часто одни и те же рутины используются не только для поиска внутри ДНК, но и внутри белковых последовательностей, генных зависимостей и т.п. А там - и латинского алфавита будет мало.

    2. 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. Нечеткий поиск

3
  • Все же, добавьте чат, пожалуйста. Не знаю, как у остальных, а у меня точно есть много вопросов, но мало ответов. 4 дек 2017 в 7:07
  • @Viktor Tomilov коментарии к этой статье видели?
    – Vasek
    5 дек 2017 в 12:01
  • @Vasek Да, видел. 5 дек 2017 в 12:56

4 ответа 4

7

Деревья хорошо себя зарекомендовали тем, что у них достаточно экономное расходование памяти, в случаях, когда существует только часть возможных ключей. Скажем нам надо сделать дерево на 3 уровня. В случае если у нас будут только варианты AAA, ACA, GCC у нас на первом уровне окажется только два узла 'A' и 'G', на втором уровне 3 узла, 'A','C'(для A),'C'(для G) и на третьем уровне 3 узла. При 64 возможных вариантах ключа из 3 нуклеотидов, всего 8 узлов - отличный результат.

А теперь проведем статистический анализ текста 8й хромосомы... Нас интересует в первую очередь соотношение всех возможных значений ключей, к реально существующим. До 11 буквенных последовательностей в референсе существуют практически все возможные варианты расстановок нуклеотидов. В 11-буквенных из 4х миллионов возможных значений отсутствуют только 240 тысяч (порядка 6%). Т.е. нам потребуется создать порядка 4 миллионов узлов для индексации цепочек из 11 нуклеотидов.

При этом, делая типизированные массивы (of TTrie) вы видимо предполагаете, что там будут храниться ссылки. А ссылки - это указатели на адрес в памяти процесса. Если вы создадите индекс и тут же будете его использовать, то все будет отлично. Но при больших затратах на индексацию хорошо бы его уметь хранить на диске. А на диск адреса в физической памяти не положишь, вам надо будет как то сериализовать объекты, а при чтении с диска опять воссоздавать с физическими адресами ... Кроме того, если память распределяется отдельно на каждый объект то появляются накладные расходы менеджера памяти. В итоге один узел вам обойдется не в 16, как вы предположили, а в лучшем случае байт в 30...

Таким образом на первый вопрос "какой из приведенных способов хранения лучше" ответ - ни один из вышеперечисленных. На второй вопрос вы ответили сами - затраты памяти безумны. В 3м вопросе, собственно вопрос отсутствует, это констатация факта, что индексация потребует много времени (кстати оценка через "O" уж слишком приблизительна. С учетом накладных расходов на поддержания таких структур средствами языка)

В связи с вышеизложенным я пришел к выводу, что единственно компактным способом хранения начальных уровней такого дерева будет один общий массив с указателями на 4 миллиона возможных ключей. пустовать в нем будет только 6% элементов. Индексом в таблице будет служить ключ, составленный из 11 букв, представленных в системе счисления с основанием в размер алфавита (в нашем случае 2 бита на букву). В случае встречи в тексте одиночных N мы будем расценивать их как все 4 возможных ключа. Участки с более чем одной N мы будем хранить в виде списка областей (на рассмотренном примере 8й хромосомы таких участков всего 12). Указатели в массиве могут указывать на узлы дерева, "продолжающие" те 11-буквенные последовательности. Или, например, сразу на список адресов, где встречены такие последовательности. Только по некоторым ключам списки очень велики (выделяются 11 A и 11 T с 106к вхождений). Всего есть около 740 ключей (из 4 млн) списки за которыми больше 3000 адресов.

Я реализовал индексацию и поиск по выше предложенному способу. За указателем храню непосредственно списки адресов и при поиске использую их быстрое сравнение (параллельный, синхронный сдвиг указателей). Работу с N-областями предусмотрел, но до конца реализовывать не стал. Требуемая память для основного индекса 32 Мб. Для хранения списков адресов используется по 4 байта на адрес, т.е. длина референса * 4. Суммарный объем файла индекса составил 581Мб для 141Мб входного файла. Время индексации (за 2 прохода по референсу) порядка 30 секунд (включая сброс на диск). Кстати, такой индекс можно использовать для быстрого построения деревьев. Алгоритм поиска не оптимизирован и обрабатывает только порядка 10-15 тысяч запросов в секунду на одном ядре. Комплиментарный поиск так же не реализован, вижу такой вариант, как создание 16 Mb таблицы перевода комплиментарных ключей и поиск по тому же индексу с пересчетом смещений в паттерне на противоположный порядок. Исходный код: Индексация, Поиск (Собиралось gcc под linux)

6
  • Забавно, один из коммерческих продуктов при импорте референсной последовательности превращает её в файл размером как раз примерно чуть больше, чем в 4 раза :) 6 дек 2017 в 10:11
  • @ViktorTomilov Работая над этим ru.stackoverflow.com/a/768244/194569 я вдруг подумал, что предложенный в этой статье способ можно дополнить, расположив внутри нашего словаря смещения к областям данных, где они встречены не в порядке возрастания смещения (как было нужно для поиска слиянием списков), а в порядке последующих букв идущих за выбранным ключом. Тогда мы можем не сравнивать списки, а искать по одному списку, выбранному для начала паттерна, двоичным поиском, заглядывая в данные.
    – Mike
    9 янв 2018 в 13:26
  • Как оценивать затраты в этом случае? 9 янв 2018 в 13:39
  • @ViktorTomilov Все бы вам оценить :) вопрос как всегда в размере списка для конкретного ключа. Как там считается скорость двоичного поиска от размера словаря ... вроде O=log2(N) где N размер списка. НО это число гипотетически. А фактически для осуществления поиска нам придется прыгать в непредсказуемом порядке по файлу данных, что бы посмотреть больше там значение чем искомое или меньше, что может привести к существенному вводу-выводу из довольно далеких областей диска (хотя на ssd это не особо критично)
    – Mike
    9 янв 2018 в 13:50
  • 1
    @ViktorTomilov В этом случае в изначальном алгоритме индексации надо будет убрать проверку на "слишком большой список", который там не индексировался. таких списков очень мало и их объем при 11 буквенном ключе, если не изменяет память, был порядка 150к элементов (для ключей вроде 11 A). Это около 16 операций сравнения. причем сравнивать можно начиная сразу с 12 буквы (потому что первые 11 гарантированно равны для одного списка).
    – Mike
    9 янв 2018 в 16:02
4

Учитывая что в данном конкретном случае ключевым моментом (на мой взгляд) является размер данных, предлагаю оптимизировать именно этот фактор; (Для разбора полетов в качестве референса будем рассматривать 8ю хромосому человека (A: 43333530, C: 29030173, G: 29103787, T: 43300646, N: 370500 штук, всего 145138636 букв)) Создадим пять массивов A,C,G,T,N размерностью 145138636 бит ... бит 1 в соответствующем массиве будет означать наличие нуклеотида в этой позиции;

1
  • 1
    Интересно. В таком ключе я не рассматривал проблему. Возможно, достаточно 4х массивов - в случае N каждый из массивов A, G, T, C получает единицу. Но это, получается, вообще третий подход: и не прямой поиск, и не деревья. 4 дек 2017 в 7:29
1

Я бы предложил следующий способ. Поскольку имеется всего 4 нуклеотида: A, C, G, T и + пятый неопределенный нуклеотид N, то логично напрашивается 5-теричная система счисления базирующаяся на 0=A, 1=C, 2=G, 3=T, 4=N, соответственно 8-ю хромосому содержащую 145 млн. последовательностей можно закодировать в 5-теричное число, которое занимает в памяти 145/8*5 млн.= 90 мегабайт.

В этом случае поиск подстроки вырождается в поиск 5-теричной битовой маски

На Java это относительно просто реализуется на основе BigInteger

2
  • 2
    боюсь маска не 5-ти ричная. N фактически означает "любая буква". под запрос ATNC должны попасть 4 варианта ATAC, ATGC, ATTC, ATCC. И все усугубляется тем, что N бывают и в референсе и соответственно последовательность с одной N в референсе дложна быть найдена по 4 разным ключам ...
    – Mike
    4 дек 2017 в 10:27
  • @Mike прав. В комментариях неудобно писать, - внес добавление (Update 2) по этому поводу непосредственно в вопрос 4 дек 2017 в 13:18
1

Я дам интересней направление: нечеткий поиск. Возможно вы не видели этот https://habrahabr.ru/post/115147/ метод ("MinHash — выявляем похожие множества"). Если немного модифицировать и считать геном человека как много-много перекрывающихся окон со словами, то после индексирования мы сможем относительно быстро искать нужное нам окно (в котором будет не обязательно точное совпадение). Этот метод ложится и на SQL, правда у меня скорость получилась не очень быстрая. (Т.е. лучше заранее посчитать вычислительную сложность.)

3
  • :) Цитата: "В данной теме рассмотрим только точный поиск." Нечеткий поиск был заранее запланирован в Части 2. 5 дек 2017 в 9:56
  • 2
    В приведенной статье нет ни слова о поиске, наличие хеша само по себе ничего не даст. какое перекрытие окон вы предполагаете, если в 1 байт, то мы получаем список кешей длиной в весь исходный текст и теперь нам надо сравнить хеш от полученной строки со всем списком хешей окон, как это сделать за время значительно меньшее времени перебора входного массива (объем массива хешей в 4 раза больше входного текста)
    – Mike
    5 дек 2017 в 14:10
  • "Риды, для простоты, будут длиной 100 букв." Это, вроде, маловато. Но можно и 100. Окно делаем того же порядка. Сдвиг окна (видимо) половина от длины (50). Для каждого окна делаем (к примеру) 50 хэшей (хэш берется от "слов"). За слово будем принимать (к примеру) 5-9 идущих подряд букв, при этом слова можно (нужно?) брать с перекрытием. Размер нашей структуры будет немаленьким, скорость поиска умеренно быстрая. Поиск нечеткий. В результате поиска будут выданы координаты окон, которые предположительно имеют хорошее совпадение.
    – Vlad Chapl
    6 дек 2017 в 15:07

Ваш ответ

Нажимая на кнопку «Отправить ответ», вы соглашаетесь с нашими пользовательским соглашением, политикой конфиденциальности и политикой о куки

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.