CC BY
23
№ 5 (41) 2012
Ю. М. Суворова, аспирант Центра «Биоинженерия» РАН, г. Москва
М. А. Короткова, канд. техн. наук, доцент Национального ядерного исследовательского университета «МИФИ», г. Москва
Е. В. Коротков, докт. биол. наук, профессор Центра «Биоинженерия» РАН, г. Москва
Изучение точек разладки триплетной периодичности в нуклеотидных последовательностях генов
Задача разбиения неоднородной последовательности или временного ряда точками, разделяющими однородные участки, впервые возникла в контексте контроля качества на производстве и получила название задачи о разладке (change point problem). Она оказалась актуальной для самых различных областей, таких как контроль качества, оценка надежности, эконометрика и медицина, так как точки разладки отражают внутренние изменения исследуемых процессов.
Введение
Позиции, разделяющие два однородных сегмента, называются точками разладки (change points) [1]. Последовательности ДНК также не являются абсолютно однородными по составу, и неоднородность нельзя объяснить лишь случайными флуктуациями [2]. В результате возникла задача разбиения последовательности ДНК на однородные участки, для чего были применены некоторые методы из теории поиска точек разладки (ТР) [3]. Большинство методов создали для выявления неоднород-ностей на уровне геномов [4], однако они могут присутствовать и на уровне отдельных генов.
Основным объектом данной работы являются последовательности ДНК, которые кодируют белки. В ДНК встречается четыре вида азотистых оснований: аденин («а»), гуанин («д»), тимин («t»), цитозин («с»). Ген кодирует аминокислотную последовательность белка тандемными триплетами, называемыми кодонами, всего существует 64 различных кодона. Каждой из 20 аминокислот соответствует, по крайней мере, один
кодон. Последовательность кодонов транслируется в аминокислоты в соответствии с действующей рамкой считывания, всего есть три возможных рамки [5]. Мутации в генах осуществляются посредством замен оснований ДНК, а также путем вставок и делений, как отдельных оснований ДНК, так и фрагментов ДНК различной длины. В настоящее время популярна теория о том, что на определенном этапе эволюции последовательности не создаются заново, а собираются из частей уже существующих последовательностей [6]. Перестановка более мелких блоков, так называемых доменов, может создавать белки с новыми свойствами. Анализ показал, что число таких независимых эволюционных блоков намного меньше общего числа белков организма. Разнообразие белков может быть результатом рекомбинации доменов, осуществляющейся путем так называемых склеек (часть одного гена присоединяется к другому, образуя новую последовательность) [7] или вставок (когда часть «чужого» гена вставляется внутрь родительской последовательности).
Установлено, что кодирующие последовательности всех живых организмов обла-
№ 5 (41) 2012
дают свойством триплетной периодичности (ТП) [8-9]. Классификационный анализ показал, что ТП большинства известных кодирующих последовательностей принадлежит одному из 2500 различных классов [8]. Это намного меньше общего числа последовательностей. Можно предположить, что, если в процессе эволюции склеиваемые участки имеют различное эволюционное происхождение, они могут иметь и разный тип ТП. В таком случае точка изменения типа ТП (точка разладки ТП) может быть сигналом того, что данный ген был образован в результате склейки [10] или вставки. Для исследования таких генов необходимо разработать метод, позволяющий выделить участки гена, внутри которых ТП однородна и между которыми она различна. Это позволит создать инструмент эволюционного анали-§ за последовательностей и выявления мута-^ ций типа вставки и склейки, не основанный Е на выравниваниях и гомологии. | Основное внимание в настоящем иссле-¡3 довании уделено поиску генов, содержа-§ щих парные точки разладки (ПТР) триплет-<| ной периодичности. Эти гены можно разде-| лить на три последовательных участка, таких й что первый и последний участок обладают Ц схожим типом ТП, а ТП среднего участка ото личается от них. Таким образом, первая ТР || разделяет первый и второй участок, дру-* гая — второй и третий. Последовательно-£ сти с ПТР эволюционно могли быть получены | либо в процессе вставки среднего участка, г| либо путем двойной склейки, причем первый § и последний участки изначально имели по-^ хожую периодичность. Для того чтобы найти | такие точки, нужно иметь возможность опре-| делять сходство и различие ТП на двух участках гена. Для этого введена мера сходства и различия ТП участков гена. Меры основа-<| ны на сравнении частотных матриц. Кроме <| того, следует учесть тот факт, что при встав-^ ках оснований длины, некратной трем, про-е исходит так называемый сдвиг рамки считы-| вания, который приводит к сдвигу фазы ТП. ^ Задачей данного исследования являлась § разработка алгоритма и соответствующего
программного обеспечения для поиска точек разладки ТП (в первую очередь парных) в последовательностях ДНК, на основании мер подобия и различия ТП и с учетом возможного сдвига рамки считывания.
Обзор программного инструментария
В настоящее время благодаря новейшим технологиям секвенирования банки данных биологических последовательностей (ДНК, РНК и аминокислот), такие как GenBank, UniProt, KEGG, достигли очень больших размеров, а объемы ежегодно добавляемых данных растут. К сожалению, пока скорость получения новых последовательностей намного превышает скорость понимания смысла этих последовательностей. Проблема стала настолько очевидна, что ее обсуждение появляется даже в обычной прессе. При таком потоке данных их обработка без использования компьютерных программ уже давно стала невозможной. За последние 30-40 лет для анализа последовательностей создано большое число методов и программ, направленных на решение различных задач. Их можно условно разделить на следующие категории:
1. Сравнение (выравнивание) двух или более последовательностей для установления сходства между ними.
2. Сборка геномов и поиск последовательностей, кодирующих белки.
3. Предсказание структуры и функции белка на основании последовательности (аннотация).
4. Поиск регуляторных последовательностей.
5. Эволюционные исследования (филогения).
Гомология — одно из центральных понятий анализа последовательностей, означает сходство, обусловленное происхождением от общего предка, подразумевает, что символьное подобие последовательностей отражает функциональное, структурное и эволюционное сходство между ними, а существующие различия были приобре-
№ 5 (41) 2012
тены ими в процессе эволюции. Таким образом, если есть некоторая последовательность, о которой ничего не известно, и есть другая последовательность, похожая на нее с известными свойствами, то на основании этого сходства можно предсказать структуру, функции или эволюционное происхождение неизвестной последовательности. Поэтому одной из первых задач стало создание алгоритмов и программ, позволяющих выявлять уровень гомологии между двумя (и более) последовательностями (задача выравнивания).
В упрощенном виде задачу выравнивания можно представить следующим образом: две последовательности записываются одна под другой так, чтобы совпадало максимальное число символов в обеих последовательностях (точнее, чтобы достигался максимум весовой функции), при этом допускаются пропуски символов и сдвиги частей обеих последовательностей. Существуют методы глобального выравнивания, когда последовательности выравниваются от начала до конца, и локального выравнивания, которое подразумевает поиск наиболее похожих участков. Оба алгоритма основаны на динамическом программировании с использованием весовой матрицы, представляющей собой симметричную квадратную матрицу, ячейки которой — это веса, задающие уровень сходства между отдельными символами алфавита. Вид выравнивания зависит от этой матрицы. Задача состоит в том, чтобы использовать такую весовую матрицу, которая наиболее точно описывает эволюционные изменения в последовательности. В настоящее время существует большое количество различных модификаций и оптимизаций методов парного выравнивания и реализующих их программ, которые входят в состав практически любого пакета, предназначенного для анализа последовательностей.
Динамическое программирование позволяет построить оптимальное выравнивание двух последовательностей, однако для этого требуется время порядка O(nm) (где m и п — длины последовательностей), что
не позволяет использовать эти программы § для поиска подобий в базах данных. Для по- ^ иска в них были созданы специальные про- ^ граммы, основанные на эвристических алго- и ритмах. Наиболее известной из них являет- § ся программа BLAST (Basic Local Alignment S Search Tool), которая позволяет осуществ- ¿э лять поиск подобий в базе данных в течение ^ нескольких секунд на компьютере пользо- ^ вателя или онлайн на сайте Национального | института здоровья (http://blast.ncbi.nlm.nih. & gov/Blast. cgi). Другая популярная программа поиска подобий в базах данных, суще- Ц ствующая даже дольше, чем BLAST, — про- ^ грамма FASTA.
Существуют также программы, реализующие множественные выравнивания (выравниваются три и более последовательностей). Построение множественного выравнивания алгоритмически намного более сложная задача. Применение методов динамического программирования вычислительно затруднено из-за возрастания размерности пространства. Существующие программы для построения множественных выравниваний основываются на парных выравниваниях и эвристических алгоритмах. Наиболее известные из них: ClustalW, T-Coffee (Tree-based Consistency Objective Function For alignment Evaluation) и MUSCLE (multiple sequence comparison by log-expectation). Предсказания, сделанные на основании множественного выравнивания, считаются более надежными, чем те, которые сделаны на основании парного. По результатам множественных выравниваний строится профиль семейства белков для поиска более удаленных гомологов в программе PSI-BLAST. К сожалению, эвристические методы, используемые современными программами, не гарантируют глобальный оптимум, так как для построения множественного выравнивания используется парное выравнивание и глобальный оптимум может быть пропущен. Был даже создан проект Phylo, авторы которого считают, что человек может построить множественное выравнивание лучше, чем компьютер-
№ 5 (41) 2012
ные программы, и на своем сайте предлагают пользователю поправить построенные программами множественные выравнивания белковых семейств так, чтобы увеличилась итоговая весовая функция.
С появлением большого числа современных секвенаторов (систем для определения генетических последовательностей) очень распространенной задачей стала сборка геномов в единую последовательность из множества отдельных кусочков, выдаваемых этими машинами, а также поиск в последовательностях участков, кодирующих белки. Далее не будем рассматривать такие программы. Они обычно входят в состав оборудования самих секвенаторов. После того как найдены кодирующие последовательности, возникает вопрос — какой именно белок кодирует каждая из них? На него пы-§ таются ответить программы и методы, пред-^ сказывающие третичную структуру и функ-Е цию белка по его последовательностям. Эта | задача (задача аннотации) до сих пор оста-¡3 ется одной из самых сложных. Для ее ре-§ шения часто используется понятие гомо-<5 логии. На основании множественных вы-| равниваний сходные последовательности й объединяются в семейства. Составляются Ц ресурсы, которые объединяют по функциям о в семейства последовательности целиком, || такие как PIRSF, или выделяя в них отдель-* ные мотивы (функционально значимые уча-£ стки аминокислотной последовательности) | на основе множественных выравниваний, г| PROSITE или скрытых марковских моделей Pfam. В настоящее время создано большое ^ количество методов предсказания структу-| ры белков, использующих самые разные ! подходы.
Ц После завершения полного секвениро-вания первых геномов было обнаружено, <| что гены, кодирующие белки, составляют <| лишь малую часть эукариотичеких гено-^ мов, а их число в геноме человека (примерна но 21 тыс.) незначительно отличается от ге-| номов более простых организмов. Тогда для ^ объяснения сложности организмов исследо-§ ватели сосредоточили внимание на так на-
зываемых регуляторных последовательностях. Эти некодирующие участки ДНК часть аппарата регуляции транскрипции — механизма увеличения или подавления уровня транскрипции определенного гена или группы генов, для регуляции уровня соответствующих белков в клетке. Программы поиска таких участков делятся по типу используемой информации на программы априорного поиска, которые применяют только саму исследуемую последовательность, и обучающиеся, которые дополнительно имеют библиотеку известных мотивов, а также по типу регуляторных последовательностей, для поиска которых они созданы (например, программа TESS (Transcription Element Search System) для поиска сайтов связывания транскрипции или программа Promoter для поиска промоторов в геномах позвоночных). Для полного изучения некодирующей части генома в 2003 г. был создан проект под названием Encyclopedia of DNA Elements (ENCODE), в ходе которого различные методы (вычислительные и экспериментальные) объединяют для определения того, что помимо генов обусловливает сложность генома человека.
Еще одной задачей анализа последовательностей является филогенетический анализ. Его цель — выявление родства и взаимосвязей между популяциями, организмами, геномами или генами, изучение происхождения отдельных генов и белков. Результаты таких исследований обычно представляются в виде дерева. Существует большое число программ для построения филогенетических деревьев на основании биологических последовательностей. Они различаются математическими алгоритмами построения деревьев, мерами расстояния и типами используемых последовательностей. Обзор около 300 программ филогенетического анализа представлен на сайте http://evolution.genetics. washington.edu/phylip/software.html. Часто используемым является свободно распространяемый пакет программ PHYLIP (Phylogeny Inference Package). Другая из-
ПРИКЛАДНАЯ ИНФОРМАТИКА /-
' № 5 (41) 2012
вестная программа — CLUSTRAL — позволяет строить филогенетические деревья на основании множественных выравниваний.
Программы выравниваний не только исторически являются одними из первых программ анализа биологических последовательностей, но также лежат в основе многих других программ в этой области. Так, программа BLAST, пожалуй, самый популярный инструмент современной биоинформатики — согласно Google на исходную статью 1990 г. насчитывается 41 876 ссылок. Однако, несмотря на это, возможности выравниваний ограничены. Программы аннотации на основе анализа последовательностей могут покрыть до 70% белков, представленных в банке данных аминокислотных последовательностей UniProt. Это ограничение связано с тем, что в процессе эволюции последовательностей пред-ковые формы не сохранились или исходное сходство было утеряно в результате мутаций. Поэтому в последнее время все чаще предпринимаются попытки создания алгоритмов анализа последовательностей, не использующих выравнивания так называемых alignment-free-методов. Создатели таких программ используют для решения задач только статистические свойства самой символьной последовательности. Этот подход особенно часто используется при поиске регуляторных последовательностей, так как для многих из них известны специфические мотивы. Существуют такие программы и для филогенетических исследований. Однако по сравнению с методами, использующими выравнивания, их доля очень невелика. Разработке метода анализа последовательностей без использования выравниваний посвящена данная работа.
Методы
Мера различия ТП двух участков
Рассмотрим нуклеотидную последовательность S длины L (L кратна трем).
Выберем в в две координаты х1 и х2 §
(0 < х1 < х2 < L -I) и выделим участки длиной ^
I [х1, х1 + I) и [х2, х2 + I). Для каждого из них по- ^
строим частотные матрицы и
М = М( х1, I) = [тД/, У)]4хз ||
и М2 = М(х2, I) = [т2(/, У)]4хз. ¿|
Элемент матрицы т(/, ]) равен числу элемен- ^ тов типа / (/ = 1 для 'а', / = 2; для Т, / = 3, для | 'д' и / = 4 для 'с'), находящихся в позиции & ] кодона (] = 1,2,3), для соответствующего участка. Например, элемент т1 (1,2) равен Ц числу символов Т, стоящих на вторых пози- ® циях кодонов на участке [х1, х1 + I).
В качестве меры различия двух частотных матриц будем использовать величину
1(М1,М2) = 11(М1, М2) + 12(М1,М2) + 13(М1, М2), (1)
где I], (] = 1,2,3) — информационная мера различия [11] для соответствующих столбцов рассматриваемых матриц
] (ММ) = Х тД/, У)1п(т1</, ])) +
/=1
4
+Х т2(/, ])1п(т2(/, ])) -
/=1
4 (2)
-X т (/, ]) + тг (/, У))!п(т1 (/, ]) + тг (/, ])) +
/=1
+ ]) + sг( ] ))!п(в1( ]) + S2( ])) -- S1(У)!п(в1(]))- sг(У)1п(Б2(])),
4
где sk (]) = Х тк (/, ]).
/=1
Величина 21] имеет асимптотическое распределение хи-квадрат с тремя степенями свободы [11]. Следовательно 21(М1,М2) ~%2^), где df = 3 + 3 + 3 = 9. Тогда, используя аппроксимацию нормального распределения
7(М1,М2) = .7 4 • 1(М1, М2)-V2 • df -1, (3)
получим величину7(М1,М2) ~ N(0,1).
-ч ПРИКЛАДНАЯ ИНФОРМАТИКА
№ 5 (41) 2012 ' -
со о
ё
13
u §
sa is
0
1
0 с й
1 г
0
t
со
is
u §
is §
t
1
i
il
SI *
is is
I
I
§
Для того чтобы учесть возможный сдвиг рамки считывания после точки х2, введем две дополнительные матрицы на этом участке:
Мг' = М(Х2 +1,/) = [< (Ш4хз и М" = М(х2 + 2,/) = [т2"(Ш4хз.
Заметим, что матрицы могут быть получены циклическим поворотом матрицы М2. Используя формулы (1-3), вычислим различия между матрицей М1 и матрицами М'2 и М2". Тогда в качестве меры различия ТП двух участков гена длины / с координатами начала х1 и х2 будем использовать величину
D( х1,х2) = тт[7(М1,М2)7(М1,М2')7(М1,М2")]. (4)
Мера подобия ТП двух участков
Аналогично предыдущему пункту рассмотрим матрицы М1 и М2, соответствующие участкам длины /, начинающимся в точках х1 и х2, соответственно. Рассмотрим нулевую гипотезу Н0 о том, что они являются случайными, нескоррелированными матрицами. Введем меру подобия двух матриц ТП. Для этого приведем матрицы к нормальному виду, используя следующее поэлементное преобразование:
nk ^ i) =
mk (i, i) - lpk (i, i)
фрк (i, i)(1 - Pk (i, i)) (£ mk (i, i)) ■ (£ mk (i, i))
(5)
Pk (', i) =
i=1
/2
где к = 1,2 пк(/,у) ~ N(0,1). Полученные в результате такого преобразования матрицы обозначим N1 и N2. Затем построим новую матрицу 7 = ^ (/, у)]4 . 3, перемножая соответствующие элементы матриц ^и N2
z(i, i) = n(i, i) ■ n2(i, i).
(6)
Так как случайная величина, являющаяся произведением двух стандартных нормально распределенных случайных вели-
чин, имеет плотность распределения [12] f(b) = ят1К0(1 Ь I), где К0 — модифицированная функция Бесселя второго рода, Ь — переменная закона распределения, то для любого z(/, у) можно рассчитать вероятность Р(Ь > z(i,у)).
Затем, используя обратную функцию нормального распределения, можно найти соответствующую величину у(/, у), удовлетворяющую условию Р(д > у(/,у)) = Р(Ь > z(/,у)), где д — переменная закона распределения. Затем, суммируя значения, рассчитанные для отдельных ячеек, получим
s( Х1.Х2) = ££ у (i, i).
(7)
i=1 i=1
Далее в рамках нулевой гипотезы полагаем S(x1,x2) ~ N(0,6), а величина P(N(0,6) > S(x1,x2)) есть вероятность того, что сходство двух матриц обусловлено случайными совпадениями. При больших значениях S(x1, x2) эта вероятность становится достаточно мала, что позволяет отвергнуть нулевую гипотезу.
Метод поиска парных ТР
Последовательность S разобьем точками xk = step ■ (k-1) +1; k = 1,2...K. Для каждой позиции xk заполним матрицы M(xk, l). Всего K = |_(L -1) / stepJ +1 матриц должно быть рассчитано (размер рассматриваемых участков l был выбран равным 60 символам, величина шага step = 9). Затем, сравнивая K частотных матриц по принципу «каждая с каждой», заполняем две большие матрицы Sim = [sim(i, i)] K ■ K и Dif = [dif(i, i)] K ■ K
sim(i, i) = S( xt, xi) dif (i, i) = D( x,, xi)
(8)
Элементы матрицы Sim рассчитываются по формуле (7) и отражают подобие, а элементы Dif, рассчитанные по формуле (4), — различие соответствующих участков.
Далее для всех k1 и k2 (1 < k1 < k2 < K) вычислялись значения
№ 5 (41) 2012
Вд,к2) = £ £ s/m(/,j>Г £ £ б№(Ц)+
1</<к К / <к2 1</<к к < ¡<к2
+ £ £ s/m(/,j)+ £ £ s/m(/,j)+ (9)
1</<к, к2 < ¡<К к,</<к2 к,< ¡<к2
+г £ £ ш(/,1)+ £ £ »т(/Д
к,</<к2 к2 < ¡<К к2 < ¡<Кк2 < ¡<К
Для того чтобы проиллюстрировать смысл выражения (9), рассмотрим случай, когда рассматриваемая последовательность содержит вставку, кратную трем, между позициями, соответствующими к, и к2 (случай не кратной трем вставки будет рассмотрен ниже). Тогда первое, четвертое и седьмое слагаемые уравнения (9) отражают подобие внутри интервалов (\х„),( х^М Хк2, хк). Второе и пятое слагаемые отражают различие ТП между участками (1, хк) - (хк,хк ) и (хк,хк ) - (хк,хК), соответственно, а третье слагаемое — подобие ТП между участками (1,хк )и (хк ,хК). Коэффициент г был выбран эмпирически, чтобы уравновесить вклад величин подобия и различия в итоговое значение (г = 7).
Для учета общей однородности исследуемой последовательности в была введена поправка
М = £ £ эта,¡),
1</<К 1< ¡ < К
(10)
которая отражает случай отсутствия вставки ТП в последовательности Э. С учетом этой поправки будем далее для поиска ПТР использовать величину
Определение статистической значимости Шшу.
С целью определения статистической значимости найденного значения нужно иметь в распоряжении контрольную выборку. Ее должны составлять случайные последовательности такой же длины и обладающие тем же уровнем ТП, как реальные последовательности анализируемых генов. Для этого каждый рассматриваемый ген Э был разделен на три последовательности, таким образом, что в первую последовательность (С1) вошли символы, стоящие на первых позициях кодонов Э( х1 : / = 1 + 3п; п = 0,1,2...), во вторую (С2) — символы, находящиеся на вторых позициях (х1 : / = 2 + 3п; п = 0,1,2...), а в третью (С3) — на третьих (х1 : / = 3 + 3п; п = 0,1,2...). Затем из каждой последовательности С( путем случайного перемешивания была получена последовательность Я.. Полученные последовательности Я1, Я2 и Я3 вновь объединили в одну (Я), в соответствии с позицией в кодо-не. Новая последовательность Я имеет ту же длину, частоты символов и тот же уровень ТП, что и исходная последовательность Э, но при этом перемешивание должно выровнять все имеющиеся неоднородности ТП.
Описанным выше способом для каждой последовательности Э было получено 500 случайных последовательностей Я, для каждой из которых по формуле (11) рассчитана величина 1Мтах. Затем рассчитано среднее (1Жтах) и дисперсия (0(Мтах)). Итоговая статистика для последовательности Э выражается формулой
и
I
I
0
1
со
о &
£
м к кг) = ед, кг) - м
(11)
7 = Мтах - М тах
(12)
Найдем такие позиции к1 и к2, которые дают максимум величины
Мтах = тт ах(М (к1, к,)),
— это потенциальные точки разладки ТП. После того как позиции максимума найдены, нужно проверить статистическую значимость такого максимума.
Также возможен сдвиг рамки считывания на один или два символа после второй точки разладки (что соответствует вставке, длина которой не кратна трем). Для того чтобы учесть такую возможность, в формуле (10) в третьем слагаемом нужно заменить матрицу первой рамки на втором участке, на матрицу, соответствующую второй или третьей рамке. Таким образом, переменную 7, рас-
81
№ 5 (41) 2012
о ё
13
и §
е §
0
1
о с
Й £
г
0
1
со
¡5
и §
£ ¡5
I
!
I
I
£ е
и
! §
считанную без учета сдвига, обозначим с учетом сдвига на один символ — 72, и для сдвига на два символа — 7Ъ.
Ввиду триплетной структуры исследуемых последовательностей величина, рассчитанная по формуле (12), имеет распределение, отличное от центрального нормального. Поэтому требуется провести дополнительное моделирование, чтобы выбрать уровни значимости для величин 72 и 7Ъ.
Поиск одинарных ТР
Важно заметить, что гены, содержащие одинарные точки разладки, подробно описанные в работе [10], также могут приводить к превышению порога одной из величин 710, 720 или 730 для данной последовательности. Следовательно, каждый ген, в котором было найдено значение 710, 720 или 730, превышающее соответствующий порог, должен быть дополнительно исследован на предмет одинарной точки разладки, прежде чем можно будет сделать вывод о наличии ПТР в этой последовательности. Процесс поиска одинарных склеек похож на поиск парных, но здесь рассматривается только одна координата к1 и вычисление величины W1(k1, к2) заменяется величиной W1(k1):
ВД) = X X ту)+
1</<к11< у <к1
+ г XX d/f(/,у) + X X ту). (13)
1< / < к1 к1 < у <К к1 </<К к1 < у <К
Далее производятся аналогичные вычисления с использованием формул (10-12). Итоговую величину обозначим V. Таким образом, если одно из значений Z■¡ (/ = 1,2,3) превосходит значение V, это событие рассматривается как двойная точка разладки, иначе как одинарная.
Определение пороговых значений для поиска ТР
В качестве данных для исследования были использованы гены 17 бактериальных геномов, загруженных с сайта базы данных KEGG/Genes [13]. Они содержат 69 936 ко-
дирующих последовательностей. На основе этих данных, путем перемешивания, описанным в пункте определения статистической значимости, была создана случайная выборка. Как уже указывалось, такой способ перемешивания позволяет создавать искусственные последовательности, максимально соответствующие интересующим нас свойствам реальных генов. Появление ПТР в таких сгенерированных последовательностях может быть вызвано только случайными факторами. Основываясь на этой выборке, пороговые значения меры поиска парных склеек 710, 720 и 730) были определены таким образом, чтобы число случаев ПТР, найденных на случайной выборке составляло не более 18% от числа случаев, найденных в реальных кодирующих последовательностях. Пороговые значения составили 710 = 3,1, 720 = 2,4 и 730 = 2,5. Аналогичным образом было выбрано пороговое значение для поиска одинарных точек разладки V, оно составило V, = 3,8. Следовательно, общий порог составляет также 3,8. То есть в случае если максимальным окажется величина V (V > 7 > V0,i = 1,2,3), будем считать, что этот ген содержит одинарную точку разладки. Если же одно из значений 7 превосходит V (7 > V > V0), то последовательность содержит ПТР, соответственно, без сдвига рамки или со сдвигом на один или два символа.
Контурные диаграммы различия ТП в генах
Для иллюстрации различия ТП в последовательностях генов использовались контурные графики меры различия D(x1,x2) (формула (3)), участков последовательности в длиной /. Варьируя независимо координаты х1 и х2 (х1 < х2) с шагом, равным трем, перебираем все возможные сочетания двух участков длины / для данной последовательности, на каждом шаге вычисляя различия ТП этих участков. Контурные графики, созданные таким образом, симметричны относительно главной диагонали. Чем темнее цвет некоторой области, тем сильнее ТП участка отличается от остальной последовательности.
82
№ 5 (41) 2012
Результаты и обсуждение
В предыдущей части был описан алгоритм поиска парных и одинарных точек разладки в кодирующих последовательностях. Этим методом проанализированы как реальные гены, так и искусственные склейки и вставки.
Анализ искусственно созданных последовательностей оснований ДНК
Для начала был проведен контрольный поиск в ПТР в искусственных периодических последовательностях. Сгенерирована искусственная периодическая последовательность в = ад)160. При обработке искусственной последовательности нашей программой, как и ожидалось, никаких точек разладки не было найдено. Затем проверили последовательность реального гена (KEGG идентификатор BSU26890). Данный ген также не содержит ТР триплетной периодичности. В его последовательность после 240-й позиции была произведена вставка из 180 нуклеотидов, содержащая другой тип ТП. Контурный график различий ТП на разных участках гена BSU26890 со вставкой из 180 нуклеотидов после 240-го символа представлен на рис. 1.
100 200 300 400 500 600 700 300 X
Рис. 1. Контурная диаграмма различия ТП
На диаграмме можно видеть, что ТП на уча- § стке примерно от 240 до 400 заметно отли- ^ чается от периодичности остальной после- ^ довательности. и
Для того чтобы проверить, как точно § данный метод определяет события вставки ¡5 и склейки генов, были созданы два набора ¿э модельных последовательностей. Первый ч; из них содержал искусственно склеенные ^ последовательности, второй — последова- | тельности с искусственной вставкой. При & создании этих выборок из генов 17 бактериальных геномов случайным образом вы- Ц бирались два гена. Затем последователь- ^ ности выбранных генов случайным образом делились на две части, и часть одного гена склеивалась с частью другого. В случае вставки часть одного гена была вставлена в случайно выбранную позицию другого. Эта процедура повторялась 10 000 раз в обеих выборках. Полученные последовательности обработали нашей программой. В результате определили 8018 случаев, из них 6306 случаев были позже определены как ОТР, оставшиеся 1712 — как значимые случаи ПТР. При обработке выборки, содержащей случайные склейки, 7566 случаев определены как ОТР и только 127 как ПТР. Результаты моделирования указывают на то, что в работе была определена только нижняя граница возможного числа ПТР, так как значительная часть вставок определяется как ОТР. В то же время довольно точно определены ОТР в случае склейки (только около 2% склеек были определены как ПТР).
Анализ последовательностей оснований ДНК различных бактериальных генов
Разработанным алгоритмом были проанализированы гены 17 бактериальных геномов (список геномов представлен в табл. 1). Всего было найдено 9159 генов, в которых одна из величин 7( превосходит пороговое значение (см. «Методы»), подробная статистика случаев также приведена в табл. 1.
№ 5 (41) 2012
Таблица 1
Статистика результатов поиска парных и одиночных точек разладки в 17 геномах
Геном Число генов с одиночной ТР Число генов с ПТР, без сдвига Число генов с ПТР, сдвиг равен 1 Число генов с ПТР, сдвиг равен 2 Всего
A. butzleri 227 50 23 20 320
A. vinelandii_Ent 477 92 71 64 704
B. avium 232 72 32 14 350
B. mallei 847 150 94 87 1 178
B. subtilis 444 114 49 20 627
E. coli 357 70 35 32 494
L. fermentum 170 41 15 19 245
M. capsulatus 281 78 25 26 410
P. aeruginosa 635 142 96 98 971
S. aureus_COL 221 51 17 18 307
S. enterica_ Choleraesuis 417 88 60 33 598
S. pneumoniae 150 29 13 8 200
S. sonnei 396 71 35 30 532
S. typhimurium 392 95 50 43 580
V. cholerae 246 48 31 17 342
X. campestris 604 91 63 33 791
Y. pseudotuberculo-sis_YPIII 363 80 48 19 508
Всего 6 459 1362 757 581 9 159
I
I
§
После дополнительной проверки было Гены с ПТР подробно изучены в рабо-
установлено, что 6459 генов содержат ОТР, те [10], поэтому произвели лишь сравнение
что составляет около 10% общего числа ге- результатов, полученных с ОТР, с прежними
нов, и 2700 генов содержат ПТР. данными. Внутри тех же 17 геномов в преды-
84
№ 5 (41) 2012
дущей работе было найдено 2100 случаев ОТР ТП, из них 1109 — это гены, обнаруженные и в настоящей работе. Разницу в результатах можно объяснить тем, что в предыдущей работе при поиске точек разладки использовалась только мера различия ТП смежных участков, в то время как в настоящей работе была дополнительно введена мера подобия участков гена (см. «Методы»), что позволило более аккуратно определять точки разладки.
Основной интерес для нас составляют оставшиеся 2700 случаев ПТР ТП, которые составляют около 4% размера исходной выборки. Примеры двух генов из генома бактерии B. subtilis, содержащих ПТР, представлены на рис. 2 (ген BSU02140) и 3 (ген BSU1632).
На рисунке 2 видно, что координаты ПТР составляют около ~600 и ~700 нуклеотидов соответственно.
Из рисунка 3 видно, что парные точки разладки приходятся на позиции ~650 и ~800 нуклеотидов. В обоих случаях легко увидеть участок, периодичность которого отличается от ТП остального гена.
Для проведения дальнейшего анализа ген, содержащий парные точки разладки ТП, удобно представить в следую-
щем виде: S(L) = S1 + S2 + S3, ^ = S[0,CP1]; S2 = [СР1,СР2]; S3 = [СР2^], где СР1 и СР2 — координаты точек разладки ТП. Рассмотрим две возможные гипотезы эволюционного происхождения таких генов. Первая состоит в том, что последовательность S была образована путем вставки последовательности S2 в предковую ^ + S3). Вторая гипотеза — это формирование гена S двумя последовательными склейками участков S1, S2 и S3, изначально принадлежавших разным генам (однако участки S1 и S3 обладали схожей периодичностью).
Для проверки обеих гипотез был проведен поиск возможных предковых последовательностей исследуемых генов с использованием поиска подобий в базе данных. В рамках первой гипотезы производился поиск последовательностей, имеющих подобие с участками S1 и S3, но без центральной части S2. Если же справедлива вторая гипотеза, то должны существовать три независимых гена, каждый из которых подобен только одному из участков S1, S2 или S3. Однако, для того чтобы можно было говорить о достоверности таких подобий, требовалось, чтобы как минимум для двух участков рассматриваемого гена нашлись независимые подобия. Для поиска подобия ис-
и
I
S
в
о t
S
О £
Рис. 2. Контурная диаграмма различия ТП для гена BSU02140
Рис. 3. Контурная диаграмма различия ТП для гена BSU16320
85
№ 5 (41) 2012
Таблица 2
Пример выравнивания для гена BMA0552
Выравнивание
Первая часть
Score = 241 bits (615), Expect = 2e-62 Identities = 118/282 (41%); Positives = 174/282 (61%); Gaps = 1/282 (0%); Frame = +1
Query: 34 RFEPKPILAQLPHLPGVYRYYDVQDAVLYVGKARDLKKRVSSYFTKTQLSPXXXXXXXXX 213
+F+ K L + PGVYR YD V+YVGKA+DLKKR+SSYF S
Sbjct: 4 QFDAKAFLKTVTSQPGVYRMYDAGGTVIYVGKAKDLKKRLSSYFRSNLASRKTEALVAQI 63
Query: 214 XXXETTVTRSXXXXXXXXXXXXXXXXPRYNILFRDDKSYPYLKLTGHRFPRMAYYRGAVD 393
+ TVT + PRYN+L RDDKSYP++ L+G PR+A +RGA
Sbjct: 64 QQIDVTVTHTETEALLLEHNYIKLYQPRYNVLLRDDKSYPFIFLSGDTHPRLAMHRGAKH 123
Query: 394 KKNQYFGPFPSAWAVRESIQILQRVFQLRTCEDSVFNNRTRPCLLHQIGRCSAPCV-GAI 570
K +YFGPFP+ +AVRE++ +LQ++F +R CE+SV+ NR+RPCL +QIGRC PCV G + Sbjct: 124 AKGEYFGPFPNGYAVRETLALLQKIFPIRQCENSVYRNRSRPCLQYQIGRCLGPCVEGLV 183
Query: 571 GEEDYARDVDNASRFLLGRQGEVMGELERKMHAFAAELKFEQAAAVRNQMSSLAKVLHQQ 750
EE+YA+ V+ FL G+ +V+ +L +M + L+FE+AA +R+Q+ ++ +V +Q Sbjct: 184 SEEEYAQQVEYVRLFLSGKDDQVLTQLISRMETASQNLEFEEAARIRDQIQAVRRVTEKQ 243
Query: 751 AIDVGGDSDVDILAVVAQGGRVCVNLAMVRGGRHLGDKAYFP 876
+ GD D+D++ V G CV++ +R G+ LG ++YFP Sbjct: 244 FVSNTGD-DLDVIGVAFDAGMACVHVLFIRQGKVLGSRSYFP 284
Вторая часть
Score = 247 bits (630); Expect = 3e-64; Identities = 155/346 (44%); Positives = 203/346 (58%); Gaps = 11/346 (3%); Frame = +1
Query: 1240 TEVLEAFIAQHYLGNR----VPPVLVVSHAPANRELI-DLLVEQAGHKVAVVRQPQGQKR 1404
+EV+E F+ Q YL +P +++ +++ L+ D L E AG K+ V +P+G +
Sbjct: 293 SEVVETFVGQFYLQGSQMRTLPGEILLDFNLSDKTLLADSLSELAGRKINVQTKPRGDRA 352
Query: 1405 AWLTMAEQNARLALARLLSEQGSQQARTRSLADVLGYESDDLAQLRIECFDISHTMGEAT 1584
+L +A NA AL LS+Q + R +LA VL R+ECFDISHTMGE T Sbjct: 353 RYLKLARTNAATALTSKLSQQSTVHQRLTALASVLKLPEVK----RMECFDISHTMGEQT 408
Query: 1585 QASCVVYHHHRMQSSEYRRYNIAGITPGDDYAAMRQVLTRRYEKMVXXXXXXXXXXXXXG 17 64
ASCVV+ + +EYRRYNI GITPGDDYAAM QVL RRY K Sbjct: 409 VASCVVFDANGPLRAEYRRYNITGITPGDDYAAMNQVLRRRYGK---------------- 452
Query: 17 65 IDGNAVHAAASAGRLPNVVLIDGGRGQVEIARQVFSELGLDIS------MLVGVAKGEGR 192 6
A ++P+V+LIDGG+GQ+ A+ VF+EL D+S +L+GVAKG R
Sbjct: 453 --------AIDDSKIPDVILIDGGKGQLAQAKNVFAEL--DVSWDKNHPLLLGVAKGADR 502
Query: 1927 KVGLETLIFADGRAPLELGKESAALMLVAQIRDEAHRFAITGMRAKRAKTRQTSRLEELE 2106
K GLETL F L +S AL ++ IRDE+H AI G R KRAK + TS LE +E
Sbjct: 503 KAGLETLFFEPEGEGFSLPPDSPALHVIQHIRDESHDHAIGGHRKKRAKVKNTSSLETIE 562
Query: 2107 GVGAKRRQRLLARFGGLRGVVAASVDELASVEGISRALAEQIYRQL 2244
GVG KRRQ LL GGL+G+ ASV+E+A V GIS+ LAE+I+ L Sbjct: 563 GVGPKRRQMLLKYMGGLQGLRNASVEEIAKVPGISQGLAEKIFWSL 608
ПРИКЛАДНАЯ ИНФОРМАТИКА /-
' № 5 (41) 2012
пользовали программу BLAST [14] с порогом E-value 0.001, и базу данных белковых последовательностей Swiss-Prot [15] (объем базы 532473 аминокислотных последовательностей). Используя каждую из 2700 последовательностей с ПТР в качестве запроса, искали подобия, соответствующие одной из описанных выше гипотез. Поскольку координаты ТР, найденные разработанным методом, нельзя рассматривать как абсолютно точные, для сравнения координат ТР и позиций выравниваний была введена погрешность, равная 5% от длины рассматриваемой последовательности.
В результате в 73 случаях были найдены выравнивания, соответствующие первой гипотезе (происхождение гена путем вставки).
Для второй гипотезы в 18 случаях были найдены независимые подобия, соответствующие всем трем участкам, и в 273 случаях соответствующие подобия нашлись для двух из трех участков. Эти результаты показывают, что и события вставки и двойной склейки могли являться причиной ПТР в последовательностях генов. Однако тот факт, что предковые последовательности могли быть утеряны в процессе эволюции и не представлены в современных базах данных, не позволяет установить ее более точно.
Для того чтобы дополнительно проверить гипотезу о том, что ПТР соответствуют вставкам, был проведен поиск ТР в известных последовательностях со вставками. Для этого использовалась база данных вставок в белковых структурах [16]. С сайта базы данных было загружено 2137 белковых структур (PDB). Для 1676 из них с использованием онлайн-инструмента [17] были найдены соответствующие кодирующие последовательности. После удаления избыточности из выборки (последовательностей, схожих между собой более, чем на 95%) остались 232 последовательности, которые обработали программой. Значимые случаи точек разладки ТП были обнаружены в 55 последовательностях (35 —
ОТР, 9 — ПТР без сдвига рамки, 6 — ПТР g
со сдвигом рамки на один, и 5 — со сдви- ^
гом на два символа). Процент найденных ^
случаев (20 из 232) сравним с процентом, и
полученным ранее при моделировании со- §
бытий вставки (1712 из 10 000). Расхожде- S
ние может быть обусловлено различной Ja
представленостью классов ТП в модель- <
ной выборке и в реальных последователь- ^
ностях со вставкой. §
О
&
Программная реализация
Программная реализация алгоритмов по- Ц иска парных и одиночных точек разладки, ® а также моделирования для определения уровня статистической значимости найденных случаев, была произведена с использованием языка программирования Fortran. Этот язык выбрали для обеспечения эффективности алгоритма, который содержит большое количество вычислительных операций. Дополнительно для ускорения вычислений на больших объемах реальных данных и при моделировании были разработаны параллельные версии программ с использованием технологии MPI. Вычисления производились на кластере Центра «Биоинженерия» РАН. Программы для дальнейшего анализа результатов, сравнение с другими банками данных реализованы на языке Python. Все программы написаны для операционной системы Linux.
Заключение
В результате исследования был разработан алгоритм и реализующий его комплекс программ для поиска парных и одинарных точек разладки триплетной периодичности в последовательностях ДНК, кодирующих белки. Также создано программное обеспечение для визуализации наличия точек разладки в последовательностях генов.
Проведенное моделирование позволяет сделать вывод о том, что найденные случаи точек разладки в генах не могут являться результатами случайных событий. Кроме того, точки разладки не могут отражать перехо-
-N ПРИКЛАДНАЯ ИНФОРМАТИКА
№ 5 (41) 2012 ' -
ды в структурах кодируемых белков, так как в таком случае они встречались бы практически в каждой рассматриваемой последовательности. С другой стороны, процент парных точек разладки, которые были найдены в исследуемой выборке, выше, чем ожидалось на основании процента одинарных ТР, найденных в той же выборке. Это показывает, что независимое происхождение парных ТР является статистически менее вероятным событием, чем одновременное происхождение, которое возникает при вставках в ген фрагментов ДНК с другим типом периодичности.
Разработанная программа позволяет находить следы эволюционных событий, таких как вставки и склейки фрагментов ДНК, только на основании свойств последовательностей, без использования допол-§ нительных данных. Важно заметить, что ^ программа может определять эти события Е только в тех случаях, когда исходные гены | обладали различными типами периодично-§ сти. Это позволяет сделать предположение § о том, что реальное число последователь-<| ностей, образованных путем вставки или | склейки, может превосходить значение, най-й денное в настоящей работе, в 5-10 раз.
4 Дальнейшее улучшение разрешающей о способности разработанного алгоритма по-|| иска точек разладки в последовательностях * ДНК может быть достигнуто двумя способа-£ ми: во-первых, независимым определением | порогового уровня методом статистических г| испытаний Монте-Карло для каждой рас-^ сматриваемой последовательности; во-вто-^ рых, дополнительным исследованием одно-g родности периодичности для других длин
^ периодов, отличных от трех. |
| Список литературы
il
1. Bhattacharya P. K. Some aspects of change-point
^ analysis. In the book: Change-Point Problems. IMC
5 Lecture notes. Hayward CA.: IMC. 1994.
js 2. Li W. The study of correlation structures of DNA
[S sequences — a critical review // Computers &
§ Chemistry. 1997. № 21.
3. Braun J. V, Muller H. G. Statistical methods for DNA segmentation // Statistical Science. 1998. № 13.
4. Li W., Bernaola-Galvan P. A., Haghighi F, Grosse I. Applications of recursive segmentation to the analysis of DNA sequences // Computers & Chemistry. 2002. № 26.
5. Graur D., Li W-H. Fundamentals of Molecular Evolution, 2-nd Edition. Sinauer Associates, Sunderland. 2000.
6. Moore A. D, Bjorklund A. K, Ekman D, BornbergBauer E, Elofsson A. Arrangements in the modular evolution of proteins // Trends in Biochemical Sciences. 2008. № 33.
7. Li W. H, Gu Z, Wang H, Nekrutenko A. Evolutionary analyses of the human genome // Nature 2001. № 409.
8. Френкель Ф. Е., Коротков Е. В. Классификация триплетной периодичности последовательностей, собранных в банке данных KEGG // Молекулярная биология. 2008. № 42.
9. Frenkel F. E, KorotkovE. V. Using triplet periodicity of nucleotide sequences for finding potential reading frame shifts in genes // DNA Research. 2009. № 16.
10. Suvorova Y. M, Rudenko V. M, Korotkov E. V. Detection change points of triplet periodicity of gene // Gene. 2012. № 491.
11. Кульбак С. Теория информации и статистика. M.: Наука, 1967.
12. Craig C. C. On the frequency function of xy // The Annals of Mathematical Statistics. 1963. № 7.
13. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes // Nucleic Acids Research. 2008. № 28.
14. Altschul S. F, Gish W, Miller W, Myers E. W, Lip-man D. J. Basic local alignment search tool // Journal of Molecular Biology. 1990. № 215 (3).
15. Boeckmann B., Bairoch A., Apweiler R, et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL // Nucleic Acids Research. 2003. № 31.
16. Aroul Selvam R., Hubbard T., Sasidharan R. Domain Insertions in Protein Structures // Journal of Molecular Biology. 2004. № 338.
17. Hovmoller S., Zhou T. 2004, Protein shape strings and DNA sequences [Online]. Available: http://www.fos.su.se/~pdbdna/.
