CC BY
13
Рубанов Л.И.1, Селиверстов А.В.2, Любецкий В.А.3
гИППИ РАН, г. Москва, к.т.н., в.н.с., rubanov@iitp .ru 2ИППИ РАН, г. Москва, к.ф.-м.н., в.н.с., slvstv@iitp.ru 3ИППИ РАН, г. Москва, д.ф.-м.н., зав.лаб., lyubetsk@ iitp . ru
ШИРОКОМАСШТАБНЫЙ ПОИСК УЛЬТРАКОНСЕРВАТИВНЫХ ЭЛЕМЕНТОВ
В ПОЛНЫХ ГЕНОМАХ
КЛЮЧЕВЫЕ СЛОВА
Большие данные, высокопроизводительные вычисления, полный геном, ультраконсервативный элемент, полулокальное выравнивание, плотный подграф, кластеризация.
АННОТАЦИЯ
Описывается метод совместного анализа большого числа полных геномов с целью нахождения в них ультраконсервативных элементов - похожих друг на друга (не обязательно тождественных) участков, встречающихся сразу во многих геномах. Задача связана с обработкой больших данных, для ее решения разработан ряд оригинальных параллельных алгоритмов с линейной или близкой к ней сложностью.
Введение и постановка задачи. Задача состоит в поиске ультраконсервативных элементов в данном наборе геномов, последовательностей порядка 3 миллионов - 6 миллиардов букв в четырёхбуквенном алфавите. Число геномов может достигать нескольких сотен. Нужно кластеризовать слова, подпоследовательности небольшой длины, взятые из данных последовательностей, таким образом, чтобы один кластер состоял из «почти одинаковых» слов. Последнее означает, что каждый кластер, называемый ещё ультраконсервативным элементом, допускает хорошее выравнивание входящих в него слов, т.е. редакционное расстояние между ними небольшое. Забегая вперёд, отметим, что типичная возникающая трудность состоит в том, чтобы не было слишком больших («гигантских») кластеров. Это зависит от правильного выбора параметров задачи, что само по себе является нетривиальной задачей.
Задача относится к числу типичных трудных задач обработки больших данных, о чем свидетельствуют характерные размеры полных геномов, а также число уже известных геномов и темпы расшифровки новых. На неформальном уровне каждый такой элемент представляет собой участок ДНК, который встречается одновременно во многих геномах, но не обязательно в виде точной копии, а с возможностью замен, вставок и удалений отдельных букв. Таким образом, задача соединяет в себе две важнейшие задачи - кластеризацию слов и построение их множественного выравнивания.
Более строго, рассматривается следующая задача. Даны М символьных последовательностей над конечным алфавитом 2 . Любую конечную последовательность
А = {а1а2...аы} будем называть строкой, а произвольный ее отрезок между позициями i и ]
включительно - подстрокой (или словом), которую обозначим А [ г.. j А . Длины строк могут различаться, но не превосходят заданной величины N■; без ограничения общности можно считать длину каждой строки равной N. Для пары слов определено редакционное расстояние [1, 2] 1, используя фиксированные цены элементарных операций, которые предполагаются симметричными и однородными, т.е. цена вставки равна цене удаления (8 1 = 8й = 8 ^) и цены всех замен равны между собой. Часто предполагают также, что цена вставки/удаления больше цены замены, 8й > 8Х, но для простоты изложения будем временно считать 8й = 8Х = 1 ■ при этом редакционное расстояние между двумя словами равно минимальному числу операций, преобразующих одно слово в другое.
Определение. Будем говорить, что слово W = ^W1W2•••W/^ встречается в данной строке
А = {а1а2 ... ам], если в ней существуют две позиции 1<г< ]< N , такие, что редакционное расстояние [г,]])<£ , где £ - заданный порог. Для таких случаев будем использовать
обозначения w е ^ А [г.. ^ .
Задача состоит в том, чтобы найти все слова длины не менее I = N, которые
встречаются, по меньшей мере, в т < М строках, с указанием их позиций в каждой строке.
Согласно определению, здесь под «встречей» понимается не только точное совпадение слов, т.е. равенство алфавитных символов на соответственных позициях, но и приблизительное совпадение, когда слова неодинаковы, но близки в смысле редакционного расстояния. При этом порог £ предполагается небольшим, например, не более 6 редакционных операций (замен, вставок, удалений) при длине I = 150 ■ Основная сложность вытекает из размерности задачи: длина строк N очень велика (общая длина полного генома составляет порядка 109), а число строк М может достигать сотен. Трудоемкость прямого перебора квадратичная и составляет порядка
0(М2 N'^12) операций сравнения, т.е. свыше 1025, а для параллельного перебора необходима
общая память размером по меньшей мере О(МЫ) - терабайтного объема. Это не позволяет решать такую задачу даже с использованием современных суперкомпьютеров. Требуется более эффективная технология, основанная на применении быстрых алгоритмов и лучше приспособленная к распараллеливанию. С теоретико-графовых позиций предлагаемый метод состоит из двух этапов.
На первом этапе анализируются две строки и в них выделяются все пары искомых слов-кандидатов (для определенности - слов второй строки, которые встречаются в первой строке). Каждая найденная пара слов задает ребро будущего графа, в котором вершины соответствуют словам, а ребрами соединены вершины, для которых редакционное расстояние между соответствующими словами не превосходит £ . Для поиска разработан быстрый алгоритм с линейной (или почти линейной) сложностью. При наличии в составе вычислительной системы М(М — 1) / 2 процессоров данный этап может выполняться параллельно для каждой неупорядоченной пары строк, независимо от прочих пар. Однако алгоритм наиболее эффективен, когда в каждой параллельной ветви выполняется сравнение одной строки со всеми оставшимися, для чего достаточно М — 1 процессоров. Возможны также и промежуточные по числу ветвей варианты, позволяющие лучше сбалансировать неоднородную вычислительную нагрузку между процессорами; этот вопрос исследован нами в [3]. Результатом первого этапа являются независимые ребра, которые в совокупности образуют М-дольный граф низкой плотности (в котором число ребер ближе по порядку к числу вершин, чем к квадрату числа вершин), обычно состоящий из большого числа связных компонент.
На втором этапе проводится уплотнение начального графа с последующим выделением в нем плотных подграфов (в идеале - являющихся кликами) с вершинами, принадлежащими т и более различным долям. Каждый такой т-плотный подграф в целом описывает одно найденное искомое слово; вершина этого подграфа, относящаяся к к-й доле начального графа, определяет позицию этого слова в к-й исходной строке. Совокупность всех найденных т-плотных подграфов дает решение задачи. Уплотнение графа выполняется путем склейки вершин, которым соответствуют одни и те же или сильно пересекающиеся (не менее чем на длине слова, с одновременным удалением возникающих кратных ребер. Для поиска т-плотных подграфов в большом графе разработан новый параллельный алгоритм, построенный по принципу клеточного автомата, в котором каждой вершине исходного графа соответствует свой многошаговый процесс, однако эти процессы исполняются с единой синхронизацией шагов. Как показали эксперименты на реальных данных, даже для очень больших графов (например, порядка 107 вершин и 109 ребер) алгоритм завершается после небольшого числа шагов (обычно менее сотни), и гибко масштабируется для любого числа доступных процессоров. В следующих разделах излагаются подробности указанных двух этапов метода и описываются лежащие в их основе алгоритмы.
Поиск слов, встречающихся в двух строках. Вначале уточним ряд моментов, опущенных выше для краткости. Во-первых, из введенного определения встречаемости следует, что для слова w е А любое его подслово также встречается в этой строке: ^с w ^и е А . Поэтому для уменьшения числа решений задачи естественно потребовать, чтобы искомые слова были
непродолжаемыми, т.е. никакое найденное слово не было подстрокой другого найденного слова. Однако такое требование может вступить в конфликт с основной целью задачи. В самом деле, легко представить ситуацию, когда некоторое слово w встречается в m исходных строках, а его подслово uQ w - в m ' >m строках. Оба эти решения могут представлять практический интерес, поэтому ограничим требование непродолжаемости слова пределами каждого выявленного множества из m и более строк, в которых встречается это слово.
Во-вторых, на практике часто представляют интерес не любые слова, широко встречающиеся в исходных строках, а только обладающие достаточной колмогоровской сложностью. Например, в полном геноме иногда присутствуют длинные подстроки, состоящие из многих тысяч повторений одного-двух алфавитных символов, например, АТ-повторы. Найденные в таких подстроках элементы обычно не представляют интереса и должны быть отброшены. Как показано ниже, это также помогает снизить затраты памяти и ускорить работу алгоритма. Для оценки сложности найденных кандидатов использован алгоритм сжатия данных Зива-Лемпеля [4] в реализации GNU Gzip [5] с установленным порогом для коэффициента компрессии r; величина порога определялась эмпирически.
Основная трудность этого этапа связана с квадратичной (от длины строки) сложностью непосредственного сравнения слов двух последовательностей. Известные быстрые алгоритмы поиска подстроки в строке в данном случае неприменимы, поскольку способны отыскивать только точно совпадающие слова. Требовался алгоритм более низкой сложности, который был разработан на базе трех принципов:
1) Индексация строк. Наивный квадратичный алгоритм проверяет для каждого слова строки B, встречается ли оно в одной и той же строке A. Ясно, что такую проверку можно в среднем существенно ускорить с помощью надлежащей индексации строки A, которая проводится один раз, а затем многократно используется. Эта идея хорошо известна и широко используется, например, поисковыми машинами Интернет. Подробности разработанных для данной задачи индексов обсуждаются ниже.
2) Предварительный поиск по ключу. В рассматриваемой постановке ищутся элементы с высокой консервативностью, так что порог допустимого редакционного расстояния между искомыми словами предполагается низким: различия между любыми двумя словами устраняются не более чем за £ операций вставки, удаления или замены одного символа. Отсюда следует, что искомые слова гарантированно содержат точно совпадающие подслова некоторой длины k, которые можно использовать в качестве обычных поисковых ключей. Например, для указанных
выше значений параметров £ = 6, I = 150 длина ключа может составлять k » 30 .
3) Полулокальное выравнивание. После поиска по точному ключу каждое подслово-кандидат подлежит расширению в обе стороны с соблюдением порога £ , пока не будет достигнута общая длина не менее l, иначе кандидат отвергается. Для нахождения редакционного расстояния между словами обычно применяется алгоритм Нидлмана-Вунша [6] со сложностью
0(12). Этот алгоритм глобального выравнивания предполагает концы слов зафиксированными,
тогда как в нашем случае известно только начало или конец слова, в зависимости от того, в какую сторону выполняется расширение. Использование общих алгоритмов локального выравнивания (Смита-Уотермена, BLAST и др.) в рассматриваемом случае явно избыточно, коль скоро один из концов каждого слова фиксирован (примыкает к ключу). На основе алгоритма [6] был разработан полулокальный алгоритм выравнивания слов с одним закрепленным концом, причем благодаря учету порога £ и ограничению числа подряд идущих делеций D, созданный алгоритм имеет
пространственно-временную сложность 0(DI), т.е. является линейным по l, в отличие от
исходного квадратичного алгоритма.
Пусть k - длина ключа, определенная в соответствии со вторым принципом.
Проиндексируем строку A следующим образом. Строится хеш-таблица HA, в которую поочередно
заносятся все присутствующие ключи, т.е. слова длины k в составе строки A, вместе с позицией их начала в строке:
Ha = [К(i) е (Ka, j) | Ka = A[j..(j + k -1)], i = Hash(KA)},
где Hash ( Ka) - хеш-функция, преобразующая значение очередного ключа в адрес внутри хеш-таблицы. При построении хеш-таблицы неизбежно возникают коллизии, которые приводят к
снижению производительности алгоритма и в данном случае имеют двоякую природу: (а) возникающие из-за несовершенства хеш-функции и (б) вызванные повторениями одного и того же слова длины k в исходной строке. Частоту появления коллизий первого рода можно пытаться снизить путем уменьшения коэффициента заполнения таблицы и/или выбора другой хеш-функции, но коллизии второго рода неустранимы в принципе, хотя с увеличением k их частота естественным образом снижается. Частоту коллизий второго рода для реальных больших данных
(полный геном Sarcocystis neurona, длина строки N = 124377056, размер алфавита |S| = 4)
иллюстрирует табл. 1. В этом примере при длине ключа k >24 коллизии возникают в среднем менее чем в 4% случаев, но при уменьшении длины ключа до 16 эта доля превосходит 20%. Многократно повторяющиеся ключи обычно принадлежат участкам генома с низкой колмогоровской сложностью, которые в дальнейшем будут отброшены, однако для ускорения мы используем порог t для числа встреч ключа в анализируемой последовательности. Поскольку коллизии неизбежны, хоть и достаточно редки, для построения индекса выбран вариант хеш-таблицы с цепочками. Размер хеш-таблицы определяется исключительно длиной строки A и желаемым коэффициентом заполнения. Поскольку строка B на этом этапе не участвует, индекс для каждой из M исходных строк может быть построен заранее (и параллельно), а затем многократно использоваться со всевозможными B.
Таблица 1. Число повторяющихся ключей в строке большой длины
Число вхождений ключа k=16 k=24 k=32 k=48
1 94919216 119034327 121193592 121823700
2 5284353 764077 401812 242260
3, 4 1438861 190900 65945 24494
5-8 486137 60275 15503 4875
9-16 183799 20869 4129 1116
17-32 72581 7429 1202 194
33-64 29861 2996 437 155
65-128 11196 1217 208 85
129-256 4986 498 86 41
257-512 1776 148 9 3
513-1024 739 67 8 6
Более 1024 447 90 3 1
Число различных ключей 102433952 120082826 121682934 122096930
Среднее число вхождений 1.21408 1.03559 1.02191 1.01833
С построенной хеш-таблицей НА, для каждой строки В выполняется следующий
алгоритм.
1. Поочередно от каждой позиции/ строки В берем слово К^ = В[ j..(] + k — 1)], которое будет использоваться в качестве ключа для обращения к хеш-таблице.
2. Проверяем, присутствует ли этот ключ в НА; в случае отсутствия переходим на шаге 1 к следующей позиции /.
3. Если существует такое 5, что hA (5) = (Kj, /), то найдена новая пара кандидатов (/, /): на
позиции / в строке А и на позиции / в строке В. Если элемент ^ (5) содержит цепочку
позиций] и имеет место коллизия второго рода, то перебираются все пары (/,/), пропуская элементы с другими значениями ключа (при коллизии первого рода).
4. Для найденной пары кандидатов проверяем возможность расширения ключей до искомых слов и в случае успеха запоминаем их. Алгоритм расширения подробно описывается ниже, здесь просто предположим, что его сложность равна некоторой константе С.
5. Переходим на шаге 1 к следующей позиции/ строки В вплоть до конца строки.
Если не учитывать коллизии, сложность алгоритма О(N)-С , т.е. в среднем линейная от длины строки. Используемая память 0( N) также линейная, но и это может создавать трудности при очень больших N. Однако имеется важный резерв экономии памяти, связанный с представлением ключей в хеш-таблице. В случае короткого алфавита хешируемые ключи можно перекодировать с использованием меньшего числа бит на символ и хранить в упакованном виде в
меньшем числе байтов. Поскольку решаемая задача относится к строкам из 4-буквенного алфавита нуклеотидов (А, С, G, Т), то этим способом затраты памяти на хранение ключей в хеш-таблице удается снизить в четыре раза. В этом случае удобно выбирать к кратным 4. Перед проверкой по хеш-таблице слова-кандидаты подвергаются такому же преобразованию.
В п. 4 алгоритма для найденной пары (/, ]) позиций ключа в строках А и В проверяется возможность расширения ключей до искомых слов, а именно, в каждой строке ищутся положения начала, г1<1,]1< } и конца, г2>г + к— 1, }2— } +к — 1, такие что 12—11 — 1—1, }2 — Л — 1—1 и
А|Х.. Ч] » В[*
1 • ^.7*2]. Согласно определению выше, последнее условие подразумевает редакционное расстояние не более £ между расширенными словами. Вариантов такого расширения может быть много; мы хотим выбрать тот, который доставляет максимум величины Ь = i2 — ^ + *2 — 71 (суммарная длин слов), что обеспечивает непродолжаемость. Для этого одновременно с расширением ключа в каждую сторону строится полулокальное выравнивание добавляемых подстрок и запоминаются те позиции, на которых редакционное расстояние между ними изменяется. Расширение выполняется параллельно в обе стороны и продолжается до тех пор, пока не будет превышен порог £ или достигнуто начало (конец) строки. В двух полученных списках позиций выбирается лучший по суммарной длине вариант расширения вправо и влево, удовлетворяющий всем вышеприведенным условиям; отсутствие такого варианта означает неуспех проверки. Учитывая малость £, перебор всех вариантов не создает трудностей.
Опишем алгоритм расширения ключа с полулокальным выравниванием в одну сторону, для определенности - от концов ключей в сторону конца строк. Как известно, классический алгоритм [6] выравнивания двух строк с длинами т и п опирается на вычисление двух матриц F, Р порядка (т + 1)(п +1), откуда и проистекает его квадратичная пространственная и временная сложность. Этот алгоритм неэффективен для нашего применения, потому что позиции строк А и В, до которых расширяется ключ, заранее неизвестны, а матрицы требуется вычислить полностью, несмотря на наличие известного порога £ для максимума редакционного расстояния. Для устранения этого недостатка внесены следующие изменения (рис. 2): в каждой матрице выделена область, прилегающая к главной диагонали и содержащая дополнительно по D соседних диагоналей сверху и снизу (рис. 2, а). Смысл параметра D - максимально допустимое число подряд идущих делеций. Назовем такую область D-поясом. Матрица редакционного предписания нужна для составления вышеупомянутого списка позиций расширения, на которых редакционное расстояние возрастает.
Заметим, что траектория обратного хода алгоритма оптимального выравнивания обязана лежать внутри D-пояса независимо от позиции, до которой расширяется ключ в данную сторону, иначе нарушается ограничение на число подряд идущих делеций, т.к. каждый сдвиг в сторону от диагонального направления в матрице соответствует вставке или удалению символа.
Разместим выделенные D-пояса в двух матрицах G и Q порядка (АD + 1 )х(У + 1) , где V -заранее неизвестная граница расширения вправо (рис. 2, б)). Переход от координат и, V в матрицах G, Q к координатам т, п в матрицах F, Р и строках А, В осуществляется по формулам
т = V + (и — D)+, П = V + (D — и)+ где (D — u)+ означает положительную часть. С учетом
этого преобразования алгоритм выравнивания строится как обычно, но очередность вычисления элементов матриц и смысл направлений указателей соответствуют указанным на рис. 2, в). Временная и пространственная сложность этого алгоритма равны О(КО) = О(1), поскольку D = I. Таким образом, для проверки и расширения слов-кандидатов построен линейный от длины слов алгоритм.
Суммируя изложенное, заключаем, что общая сложность первого этапа нашего метода составляет О( N1) для одной пары строк. Для всех пар строк эта величина умножается еще на
О(М2), но пары строк можно обрабатывать параллельно - по отдельности или группами, как было указано выше. Перейдем теперь ко второму этапу нашего метода.
0 0 1 2 3 4 5 6 7 3 9 10
0,0 1,0 ОД 1,1 0,2
1 1,2 1,3
2 2,0 2Д 2,2 2,3 2 А ШАй ир
3 3,1 3,2 3,3 3,4 3,5 1_ЕП М
4 4,2 4,3 4,4 4,5 4,6
5 5,3 5,4 5,5 5,6 5,7 (а)
6 6,4 6,5 6,6 6,7 6,8
7 7,5 7,6 7,7 7,8 7,9
3 8,6 3,7 3,8 8,9 8,10
9 9,7 9,8 9,9 9,10
10 10,3 10,9 10,10
11 11,9 11,10
Г и Р
0
0,2 1,3 2,4 3,5 4,6 5,7 6,3 7,5 8,10
1 од 1,2 2,3 3,4 4,5 5,6 6,7 7,8 3,9 9,10
2 0,0 1,1 2,2 3,3 4,4 5,5 6,6 7,7 8,8 9,9 10,10
3 1,0 2Д 3,2 4,3 5,4 6,5 7,6 3,7 9,В 10,9 11Д0 (б)
4 2,0 3,1 4,2 5,3 6,4 7,5 3,6 9,7 10,3 11,9
0 1 2 3 4 5 6 7 3 9 10
0
1
2 т 0 < -> 5 ^ ^10 ^715 < ^20 7>30< ^35^40< >45
3 % £ /3 < ■4"47Щ (в)
4 'и 1 Щ ^141 —'24 ^ ¿К*
с и а
(ЛЕП
1_ЕП" ш
[>1_ЕП" к. Т □ОММ
Рис. 2. Преобразование матриц алгоритма Нидлмана-Вунша: (а) вид исходной матрицы с выделенным D-поясом для т=11, п=10, D=2; в клетках записаны координаты элементов; (б) расположение D-пояса в преобразованной матрице; (в) очередность вычисления элементов преобразованной матрицы и направления указателей обратного хода в матрице редакционного предписания.
Нахождение т-плотных подграфов в М-дольном графе. После первого этапа имеем совокупность пар слов, которые принадлежат двум различным исходным строкам А и В и
совпадают с точностью до редакционного расстояния е : А\_11.. /2] » В[]1..у2]. Сопоставим каждой
такой паре ребро графа; этот граф содержит М долей, вершины которых соединены ребрами (внутри долей ребер нет). Рассматриваемая задача в целом эквивалентна нахождению в этом графе т-клик, т.е. подграфов с вершинами, которые принадлежат т долям, и каждая вершина соединена ребрами со всеми вершинами прочих т — 1 долей. Это известная ЫР-трудная задача, для которой не существует быстрого алгоритма. Поэтому в рамках нашего метода решается близкая задача нахождения т-плотных подграфов, прежде всего - с наибольшим суммарным весом ребер. В качестве веса ребра естественно использовать монотонную функцию от длины и/или степени консервативности слов. Предполагается, что среди найденных решений этой
задачи можно по очереди отобрать присутствующие среди них клики, либо уже перебором выделить внутри них подграфы, являющиеся т-кликами.
Поскольку расширение точно совпадающих слов, т.е. построение вершин на конце ребра выше проводилось независимо от других ребер, на конце нескольких ребер могут стоять вершины, относящиеся к одному и тому же месту некоторой строки, но с несовпадающими концами слов и потому формально различные. Это снижает плотность всего графа и ведет к ненужному увеличению числа компонент связности, поэтому вначале проводится операция склейки вершин. При этом требуется, чтобы все объединяемые в одну вершину участки последовательности не просто пересекались попарно, а для любой пары содержали один и тот же общий фрагмент длиной не менее й.
После того, как выполнены все возможные операции склейки вершин, выполняется собственно алгоритм нахождения т-плотных подграфов. Этот алгоритм характеризуется глубоким внутренним параллелизмом и построен аналогично клеточному автомату, у которого в роли клеток выступают вершины графа, а роль соседей каждой клетки играют инцидентные ей вершины графа. В каждой вершине графа запускается независимый процесс, который состоит из последовательности шагов двух видов и начинается с шага 1. После каждого шага выполняется синхронизация процессов во всех вершинах графа, так что очередной шаг каждого процесса начинается только после окончания предыдущего шага во всех вершинах. Содержание шагов:
• Если вершина соединена ребрами менее чем с (т — 1) долями графа, то сама эта вершина и все инцидентные ей ребра удаляются.
• Если вершина соединена с некоторой долей только одним ребром, то помечаем его; в дальнейшем помеченное ребро может быть удалено только вместе с одним из его концов.
• Если при выполнении шага 1 в графе произошли какие-либо изменения, то снова выполняется шаг 1 во всех вершинах графа. В противном случае в каждой вершине однократно выполняется шаг 2.
• Если инцидентное вершине ребро не помечено и его вес строго меньше весов всех других не помеченных инцидентных ребер (или оно единственное), то ребро удаляется.
• Если при выполнении шага 2 в графе произошли какие-либо изменения, то все вершины снова переходят к шагу 1, иначе конец алгоритма. Каждая компонента связности получившегося графа есть искомый т-плотный подграф.
Данный алгоритм предоставляет гибкие возможности масштабирования вплоть до числа процессоров, равного числу вершин графа, и не требует общей памяти. Если число доступных процессоров меньше числа вершин, они равномерно распределяются между имеющимися процессорами, а затем каждый из них выполняет очередной шаг процесса по очереди для каждой назначенной ему вершины. Для этого алгоритма трудно получить аналитические оценки сложности, но это не имеет большого значения: практические расчеты с реальными большими данными показали, что второй этап предложенной технологии выполняется значительно быстрее первого.
Однако серьезная трудность состоит в том, что при неудачном выборе параметров или способа построения исходного графа в нём возникает гигантская связная компонента, объединяющая значительную часть вершин графа. Критериями успеха служат отсутствие гигантской компоненты и наличие т-плотных подграфов для значений т, больших половины от числа долей графа, то есть числа геномов. Как показали эксперименты, риск появления и размеры гигантской компоненты уменьшаются, если предварительно удалить из графа все ребра с весом меньше некоторого порога, но вопрос требует дальнейшего изучения. В любом случае, существование гигантской компоненты в исходном графе не обязательно является препятствием к выделению т-плотных подграфов с малым числом вершин при достаточно большом т.
Пример использования. Анализировались полные геномы 22 видов из надтипа А^еоМа (максимальная длина генома составляла 1,24* 108). Счет проводился для значений параметров
I = 60, k = 16, X = 200, Ss = 1, 8Ш = 2.1, £ = 17.5, D = 2, г = 2.2, ё = 40, m = 3.
Первый этап изложенного метода выполнялся в МСЦ РАН на суперкомьютерах МВС-100К и МВС-10П [7] с использованием 256-512 процессоров и потребовал около 200 часов. В результате было найдено 7.75-108 подходящих пар участков (общее число участков 4.54-107). Второй этап выполнялся в ИППИ РАН на 64-ядерном сервере с 256 Гб оперативной памяти. Сборка графа потребовала 13 час. на 16 процессорах, в результате был построен граф с 4977108 вершинами (и
тем же числом ребер 7.75-108 , т.к. кратные ребра удалялись в процессе поиска плотных подграфов). Поиск плотных подграфов продолжался 13 час. на 22 процессорах, в результате чего было построено 845 кластеров, в которые вошли 3924 вершины (помимо гигантской компоненты, содержащей около 3 млн. вершин). Данные о числе кластеров в зависимости от числа долей, представленных в кластере: 2 с представителями из 13 долей, 1 - из 12 долей, 4 - из 11 долей, 8 -из 10 долей, 17 - из 9 долей, 24 - из 8 долей, 35 - из 7 долей, 75 - из 6 долей, 131 - из 5 долей, 290 -из 4 долей, 258 - из 3 долей. Подавляющая часть кластеров (787) содержит не более одной вершины из каждой доли, и только 9 кластеров содержат более двух вершин из какой-нибудь доли (максимум 5 вершин). Отметим, что значительная часть машинного времени была потрачена на чтение и запись огромных файлов. Анализ обнаруженных ультраконсервативных элементов и полученные биологические результаты представлены в отдельном докладе.
В заключение укажем, что описанный метод поиска ультраконсервативных элементов может быть применен и в других областях, где используется сопоставление любых текстов, в том числе на естественном языке, с целью определения авторства, стиля, заимствований и т.д. [8].
Работа выполнена за счёт гранта Российского научного фонда (проект № 14-50-00150).
Литература
1. Левенштейн В.И. Двоичные коды с исправлением выпадений, вставок и замещений символов // Доклады Академии Наук СССР. 1965, т. 163, № 4, стр. 845-848.
2. Gusfield D. Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology. Cambridge University Press, Cambridge, UK, 1997.
3. Рубанов Л.И. Параллельное моделирование Монте-Карло на системах с распределённой памятью // International Journal of Open Information Technologies, 2014, т. 2, № 2, стр. 12-20.
4. Ziv J., Lempel A. A universal algorithm for sequential data compression // IEEE Transactions on Information Theory. 1977, vol. 23, no. 3, pp. 337-343.
5. GNU Operating System. GNU Gzip. http://www.gnu.org/software/gzip/
6. Needleman S.B., Wunsch C.D. A general method applicable to the search for similarities in the amino acid sequence of two proteins / / Journal of Molecular Biology. 1970, vol. 48, no. 3, pp. 443-453.
7. Межведомственный суперкомпьютерный центр РАН. http://www.jscc.ru/scomputers.shtml
8. Herranz J., Nin J., Solé M. Optimal Symbol Alignment Distance: A New Distance for Sequences of Symbols // IEEE Transactions on Knowledge & Data Engineering, 2011, vol. 23, no. 10, pp. 1541-1554.
