Оптимизация поиска гомологов нуклеотидных последовательностей
Автор работы: Пользователь скрыл имя, 01 Февраля 2013 в 13:00, дипломная работа
Краткое описание
Задача поиска гомологов нуклеотидных последовательностей в банке данных является одной из важнейших задач биоинформатики. В случае кодирующих последовательностей для проведения такого поиска успешно используется программа TBLASTN.
Содержание
ВВЕДЕНИЕ
4
1 ЛИТЕРАТУРНЫЙ ОБЗОР
6
1.1 Программы поиска гомологов нуклеотидных последовательностей . . . .
6
1.1.1 FASTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.1.2 BLASTN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1.1.3 YASS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.1.4 Nhunt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.2 Подходы к решению проблемы малой сложности . . . . . . . . . . . . . .
14
1.2.1 Маскировка участков малой сложности . . . . . . . . . . . . . . .
14
1.2.2 Корректировка матрицы замен аминокислотных остатков (метод
Ю – Альтшуля) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2 МАТЕРИАЛЫ И МЕТОДЫ
19
2.1 Средства программной реализации . . . . . . . . . . . . . . . . . . . . . .
19
2.2 Материалы тестирования методов поиска гомологов . . . . . . . . . . . .
19
2.2.1 Выборки последовательностей запроса . . . . . . . . . . . . . . . .
19
2.2.2 Банки данных . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3 РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
25
3.1 Сравнение программ поиска гомологов . . . . . . . . . . . . . . . . . . . .
25
3.1.1 Сценарии для сравнения качества программ поиска гомологов . .
25
3.1.2 Примеры результатов сравнения . . . . . . . . . . . . . . . . . . . .
27
3.2 Адаптация метода Ю – Альтшуля для матрицы замен нуклеотидов . . .
33
3.3 Новые методы корректировки матрицы замен . . . . . . . . . . . . . . . .
35
3.3.1 Сохранение ожидаемого счета . . . . . . . . . . . . . . . . . . . . .
36
3.3.2 Сохранение относительной энтропии . . . . . . . . . . . . . . . . .
3.4 Учет частот нуклеотидов в программе Nhunt . . . . . . . . . . . . . . . .
3.5 Сравнение разных методов корректировки матрицы замен . . . . . . . .
4 ЗАКЛЮЧЕНИЕ
ВЫВОДЫ
Список литературы
ПРИЛОЖЕНИЯ
56
3
Прикрепленные файлы: 1 файл
МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
имени М. В. ЛОМОНОСОВА
ФАКУЛЬТЕТ БИОИНЖЕНЕРИИ И БИОИНФОРМАТИКИ
Оптимизация поиска гомологов
нуклеотидных последовательностей
Дипломная работа
студента 5 курса
Пекова Ю. А.
Научный руководитель:
кандидат физико-математических наук,
старший научный сотрудник
Спирин С. А.
Москва 2012
СОДЕРЖАНИЕ
ВВЕДЕНИЕ
4
1 ЛИТЕРАТУРНЫЙ ОБЗОР
6
1.1 Программы поиска гомологов нуклеотидных последовательностей . . . .
6
1.1.1 FASTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.1.2 BLASTN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1.1.3 YASS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.1.4 Nhunt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.2 Подходы к решению проблемы малой сложности . . . . . . . . . . . . . .
14
1.2.1 Маскировка участков малой сложности . . . . . . . . . . . . . . .
14
1.2.2 Корректировка матрицы замен аминокислотных остатков (метод
Ю – Альтшуля) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2 МАТЕРИАЛЫ И МЕТОДЫ
19
2.1 Средства программной реализации . . . . . . . . . . . . . . . . . . . . . .
19
2.2 Материалы тестирования методов поиска гомологов . . . . . . . . . . . .
19
2.2.1 Выборки последовательностей запроса . . . . . . . . . . . . . . . .
19
2.2.2 Банки данных . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3 РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
25
3.1 Сравнение программ поиска гомологов . . . . . . . . . . . . . . . . . . . .
25
3.1.1 Сценарии для сравнения качества программ поиска гомологов . .
25
3.1.2 Примеры результатов сравнения . . . . . . . . . . . . . . . . . . . .
27
3.2 Адаптация метода Ю – Альтшуля для матрицы замен нуклеотидов . . .
33
3.3 Новые методы корректировки матрицы замен . . . . . . . . . . . . . . . .
35
3.3.1 Сохранение ожидаемого счета . . . . . . . . . . . . . . . . . . . . .
36
3.3.2 Сохранение относительной энтропии . . . . . . . . . . . . . . . . .
36
3.4 Учет частот нуклеотидов в программе Nhunt . . . . . . . . . . . . . . . .
37
3.5 Сравнение разных методов корректировки матрицы замен . . . . . . . .
39
2
4 ЗАКЛЮЧЕНИЕ
50
ВЫВОДЫ
52
Список литературы
52
ПРИЛОЖЕНИЯ
56
3
ВВЕДЕНИЕ
Задача поиска гомологов нуклеотидных последовательностей в банке данных яв-
ляется одной из важнейших задач биоинформатики. В случае кодирующих последова-
тельностей для проведения такого поиска успешно используется программа TBLASTN.
Но для поиска гомологов некодирующих последовательностей нет инструмента, кото-
рый был бы полностью удовлетворителен по качеству работы. Для этого используется
ряд программ. Самые распространенные из них — FASTA [4, 13], BLASTN [5, 6, 14] и
discontiguous MEGABLAST [16], однако каждая из этих программ обладает каким-либо
существенным недостатком. В 2004 году появилась программа YASS [2, 3], в которой
применяется ряд новых способов улучшения качества поиска гомологов. В ходе пред-
варительной работы нами была создана компьютерная программа Nhunt [17] для поиска
гомологов нуклеотидных последовательностей, которая по чувствительности превосхо-
дит программы FASTA, BLASTN и YASS.
Все известные программы поиска гомологов выдают выравнивание найденных
последовательностей с последовательностью запроса. Для выявления биологически осмыс-
ленных результатов обычно применяют оценку статистической значимости каждого по-
лученного выравнивания. В свою очередь, статистическая значимость рассчитывается
на основе веса выравнивания.
Вес выравнивания вычисляют как сумму трех слагаемых. Первое слагаемое —
количество позиций выравнивания, в которых буквы (нуклеотиды) в двух последова-
тельностях совпадают, умноженное на некоторое положительное число — цену совпа-
дения (match). Второе слагаемое — количество позиций выравнивания, в которых бук-
вы не совпадают, умноженное на некоторое отрицательное число — цену несовпадения
(mismatch). Третье слагаемое — штраф за делеции, неположительное число, которое
зависит от наличия и количества делеций в выравнивании.
Однако, при таком подходе к вычислению веса выравнивания возникает так на-
зываемая “проблема малой сложности”. Она состоит в том, что выравнивание двух
последовательностей малой сложности (то есть имеющих смещенный нуклеотидный
состав) будет иметь высокий вес в силу случайных совпадений часто встречающих-
ся нуклеотидов, а не в силу гомологии последовательностей. Аналогичная проблема
4
возникает также при поиске гомологов аминокислотных последовательностей. Перво-
начально для решения этой проблемы как для нуклеотидных, так и для аминокис-
лотных последовательностей использовалась маскировка участков малой сложности в
последовательностях банка данных, которая предотвращала появление этих участков
в результатах поиска. Затем для аминокислотных последовательностей был разрабо-
тан подход, основанный на корректировке различных матриц замен аминокислотных
остатков. Его основная идея состояла в уменьшении цены совпадения для часто встре-
чающихся остатков, что позволяло корректно оценивать значимость выравнивания в
целом. Такой алгоритм корректировки был встроен в пополярную программу поиска
гомологов аминокислотных последовательностей BLASTP.
Тем не менее, нам не известно ни одной программы поиска гомологов нуклео-
тидных последовательностей, которая бы использовала подобный метод корректировки
матрицы замен при расчете веса найденных выравниваний. Целью дипломной работы
является анализ и сравнение различных (в том числе разработанных нами) способов
корректировки матрицы замен нуклеотидов и применение этих способов в программе
Nhunt для решения проблемы малой сложности.
5
Глава
1
ЛИТЕРАТУРНЫЙ
ОБЗОР
1.1
Программы поиска гомологов нуклеотидных
последовательностей
Существует два основных подхода к построению выравнивания двух последова-
тельностей — глобальное и локальное выравнивание. Оба этих подхода предназначены
для построения выравнивания с максимальным весом. Различие состоит в том, что в
первом случае выравнивание включает в себя обе последовательности целиком, а во
втором для достижения максимального веса разрешается ограничиться лишь частями
последовательностей, которые выравниваются наилучшим образом.
Диагональю при сравнении двух последовательностей называется диагональ мат-
рицы — таблицы локального сходства. Диагональ также можно описать как сдвиг одной
последовательности относительно другой, или как полное выравнивание с делециями,
которые находятся только по краям. Все эти определения эквивалентны.
Для построения локального выравнивания обычно используется алгоритм Смита
– Ватермана [7] или его модификации. Основное же различие между существующими
программами поиска гомологов нуклеотидных последовательностей заключается в ме-
тодах предварительного поиска участков возможного локального сходства.
При построении выравнивания важной задачей является оценка значимости най-
денного выравнивания. Для случая локального выравнивания без разрывов предлагает-
ся ( [1]) использовать статистический подход, основанный на распределении максимума
весов выравнивания по независимым случайным последовательностям. Чем меньше ве-
роятность, что этот максимум будет больше, чем вес найденного выравнивания, тем
более значима данная находка.
Ожидаемое число выравниваний по случайным последовательностям с весом,
6
превышающим S, может быть рассчитано по формуле [1]
𝐸 = 𝐾𝑚𝑛𝑒
−𝜆𝑆
,
(1.1)
где 𝑚 и 𝑛 — длины сравниваемых последовательностей, а 𝐾 и 𝜆 — некоторые постоян-
ные. Полученное значение 𝐸 в литературе часто обозначается как “ожидаемая величи-
на” или “expected value” (далее по тексту — 𝐸𝑣𝑎𝑙𝑢𝑒). 𝜆 вычисляется как положительный
корень уравнения
∑︁
𝑖,𝑗
𝑝
𝑖
𝑝
𝑗
exp(𝜆𝑠
𝑖𝑗
) = 1,
где 𝑝
𝑖
, 𝑝
𝑗
— частоты встречаемости букв 𝑖 и 𝑗, а 𝑠
𝑖𝑗
— соответствующий элемент матрицы
замен. Постоянная 𝐾 также зависит только от значений 𝑝
𝑖
, 𝑝
𝑗
и 𝑠
𝑖𝑗
, но вычисляется по
более сложной формуле.
Необходимо отметить, что формула (1.1) применима лишь при отрицательном
математическом ожидании цены сопоставления двух случайно выбранных букв. Кроме
того, алгоритм Смита – Ватермана выдает осмысленный результат только при выпол-
нении того же условия (хотя формально может быть применён и без его выполнения).
Вероятность того, что существует вес выравнивания больший, чем 𝑆, рассчиты-
вается по формуле:
𝑃(𝑥 > 𝑆) = 1 − 𝑒
−𝐸
При значениях 𝐸 < 0,01 вероятность 𝑃 практически совпадает с 𝐸.
Ниже рассмотрены основные принципы работы нескольких распространенных
программ поиска гомологов нуклеотидных последовательностей. Программы перечис-
лены в хронологическом порядке, в соответствии с годом первой публикации, касаю-
щейся данной программы.
1.1.1
FASTA
Одной из самых ранних программ поиска гомологов стала программа FASTA [4],
предназначенная для эвристического поиска последовательностей (в том числе нуклео-
тидных). Для нахождения локального выравнивания с высоким весом она использует
процедуру, состоящую из четырех этапов. На первом этапе создается таблица, в которую
заносится расположение всех одинаковых слов длины 𝑘𝑡𝑢𝑝 в двух последовательностях.
Для нуклеотидных последовательностей 𝑘𝑡𝑢𝑝 может принимать значения от 1 до 6. За-
тем отбираются некоторое количество лучших (с наибольшим количеством совпадений
слов) диагоналей.
На втором этапе в отобранных диагоналях ищутся безразрывные участки с мак-
симальным весом, так называемые “затравки”. На третьем этапе проверяется, можно
7
ли соединить эти затравки, используя разрывы и учитывая штрафы за них. Наконец,
для последнего этапа отбираются диагонали, лежащие вокруг одного участка с самым
высоким весом, на которых вновь проводится локальное выравнивание с помощью мо-
дифицированных алгоритмов Нидельмана – Вунша и Смита – Ватермана.
По неясным причинам FASTA часто находит не все гомологи с высоким весом
выравнивания.
1.1.2
BLASTN
Самой распространенной на сегодняшний день программой поиска гомологов яв-
ляется программа BLASTN [5,6] и ее модификации, такие как Discontiguous MEGABLAST
[16]. Программа BLASTN перед началом поиска индексирует последовательности из
банка данных. При этом создается список, содержащий все слова фиксированной дли-
ны 𝑁, которые локально выравниваются с этими последовательностями с весом выше
порогового значения. Для каждого слова указывается расположение всех участков, с
которыми оно выровнялось. Длина слова 𝑁 для нуклеотидных последовательностей
больше или равна четырем (по умолчанию — 11).
Затем происходит собственно поиск, в ходе которого все слова той же длины из
последовательности запроса сравниваются со словами, найденными в базе данных при
индексации. Если найдена пара одинаковых слов, то начинается процесс расширения
совпадающего участка. Расширение происходит с учетом возможных разрывов, в обо-
их направлениях и до достижения максимального веса. Получающиеся в ходе такого
расширения выравнивания выдаются в качестве результата.
Такой подход, однако, приводит к потенциальному уменьшению чувствительно-
сти: если в гомологичной последовательности не встретится участок длиной 𝑁 букв,
полностью совпадающий с участком входной последовательности, то выравнивание
этих двух последовательностей не обнаружится.
Программа Discontiguous MEGABLAST по сравнению с BLASTN предназначена
для поиска более дальних гомологов. При индексации банка данных записываются не
целые слова, а паттерны длиной 𝑡 (𝑡 = 16, 18 или 21). Каждому типу паттернов соот-
ветствует число 𝑊: считается, что паттерн подобен слову, если в них совпадают хотя
бы 𝑊 букв в определенных местах (𝑊 = 11 или 12). После индексации банка данных
каждое слово из последовательности запроса сравнивается с записанными паттернами,
и при нахождении подобных паттернов соответствующий участок расширяется подобно
тому, как это происходит в программе BLASTN.
Для поиска используется несколько различных типов паттернов, предназначен-
ных для поиска гомологов кодирующих или некодирующих последовательностей раз-
личной степени сходства. Все перечисленные особенности программы призваны уве-
8
личить чувствительность, но на практике в большинстве случаев чувствительность
Discontiguous MEGABLAST в сравнении с BLASTN заметно уменьшается.
1.1.3
YASS
Алгоритм, лежащий в основе программы YASS [2, 3], во многом подобен алго-
ритму программ семейства BLASTN, однако содержит два существенных отличия.
Во-первых, для поиска лучших диагоналей используется некоторое количество
затравок — коротких участков высокого сходства между последовательностями запроса
и банка данных. Это количество не является фиксированным, как в других програм-
мах подобного рода, а может быть задано пользователем. При этом затравки могут
перекрываться друг с другом и лежать на разных, хотя и близких, диагоналях. Такие
затравки образуют так называемые “группы затравок”. Для дальнейшего поиска гомо-
логов выбираются такие группы затравок, для которых число совпадающих нуклеоти-
дов в последовательностях запроса и банка данных превышает заданный пользователем
порог.
Вторая особенность программы YASS основывается на известном свойстве боль-
шинства геномных последовательностей: транзиции (то есть мутационные замены одно-
го пурина на другой пурин или одного пиримидина на другой пиримидин) происходят
заметно чаще, чем трансверсии (мутационные замены пурина на пиримидин или наобо-
рот). Поэтому в качестве затравок используются паттерны, которые в определенных
позициях позволяют сопоставить как пару совпадающих нуклеотидов, так и транзи-
цию нуклеотида последовательности банка данных по отношению к последовательно-
сти запроса. Такие позиции по весу ценятся вдвое меньше, чем обычные совпадения
нуклеотидов.
Авторы программы утверждают, что такой подход заметно повышает как чув-
ствительность, так и селективность поиска гомологов нуклеотидных последовательно-
стей по сравнению с программами семейства BLASTN.
1.1.4
Nhunt
В ходе предварительной работы нами была создана программа Nhunt [17] для
поиска гомологов нуклеотидных последовательностей. Показано, что на большинстве
примеров она превосходит программы BLASTN, Discontiguous MEGABLAST и FASTA в
чувствительности, а также превосходит программу FASTA в скорости. Основные прин-
ципы работы программы Nhunt заключаются в следующем:
9
Поиск наилучших диагоналей
Вначале программа индексирует последовательность запроса. При этом для каж-
дого возможного слова длины 𝑤 + 1 вычисляются позиции в последовательности, в ко-
торых встретились оно само или его “соседи”. Слово A является соседом слова B, если
одновременно выполняются следующие условия:
— длины слов A и B совпадают
— первые буквы слов A и B совпадают
— хотя бы 𝑤 букв в одинаковых позициях в словах A и B совпадают
Точно так же индексируется последовательность, комплементарная к последо-
вательности запроса. После этого начинается проход по последовательностям банка
данных. При этом программа поочередно просматривает все позиции каждой последо-
вательности банка, сдвигаясь на один нуклеотид. Для каждой позиции находится слово
длиной 𝑤 + 1, которое начинается в данной позиции.
Затем для каждого слова, найденного в позиции 𝑦 последовательности банка
данных, подбирается идентичное ему слово из результата индексации запроса. Пусть
данное слово (сопоставленное с точно таким же словом или с соседом другого слова)
встречалось 𝑘 раз в запросе в позициях 𝑏
𝑖
. Тогда для всех 𝑖 от 1 до 𝑘 счет для каждой
диагонали, численно задаваемой выражением
𝑦 − 𝑏
𝑖
+ 𝑙𝑒𝑛𝑞𝑢𝑒𝑟𝑦 − 𝑤,
увеличивается на единицу. Здесь 𝑙𝑒𝑛𝑞𝑢𝑒𝑟𝑦 — длина последовательности запроса.
После того, как подобным образом обрабатывается очередная последователь-
ность банка данных, из всех возможных диагоналей отбираются те, для которых вы-
полняется условие:
𝐷
𝑘
≥ 𝐸
𝑞
+ 𝜎 · 𝑝,
(1.2)
где 𝐷
𝑘
— счет 𝑘-ой диагонали, 𝑝 — задаваемый пользователем параметр, 𝐸
𝑞
— ожи-
даемый счет и 𝜎 — среднеквадратичное отклонение. Эта оценка основана на предпо-
ложении, что в каждой позиции каждой диагонали 𝑃
𝑐𝑜𝑖𝑛𝑐
— вероятность случайного
совпадения слова из последовательности запроса и слова из банка данных — одина-
кова и зависит только от длины слова. Такая вероятность может быть вычислена по
формуле
𝑃
𝑐𝑜𝑖𝑛𝑐
=
1 + 3 · 𝑤
4
𝑤+1
,
(1.3)
где 𝑤 + 1 — длина слова. Тогда ожидаемый счет для каждой диагонали может быть
рассчитан как произведение длины запроса на вероятность совпадения слов в одной
10
позиции:
𝐸
𝑞
= 𝑃
𝑐𝑜𝑖𝑛𝑐
· 𝑙𝑒𝑛𝑞𝑢𝑒𝑟𝑦,
(1.4)
В рамках предложенной модели количество соседей для каждой диагонали имеет би-
номиальное распределение, для которого среднеквадратичное отклонение задается сле-
дующей формулой:
𝜎 =
√︀
𝑙𝑒𝑛𝑞𝑢𝑒𝑟𝑦 · 𝑃
𝑐𝑜𝑖𝑛𝑐
· (1 − 𝑃
𝑐𝑜𝑖𝑛𝑐
)
Таким образом, чем больше счет диагонали отличается от ожидаемого, тем более
вероятно, что она не является случайной и на ней следует провести локальное вырав-
нивание.
Локальное выравнивание для лучших диагоналей
Для отобранных вышеописанным способом лучших диагоналей в несколько ите-
раций проводится локальное выравнивание с помощью упрощённого алгоритма Смита
– Ватермана. Покажем, как это происходит, на примере одной короткой диагонали.
1. Обе последовательности в соответствии с номером диагонали записываются друг
под другом:
C A C C G C G A A
G C A A T C A G C C T A T A C G G
2. Удаляются краевые участки, в которых представлена только одна последователь-
ность и в которых выравнивание невозможно:
C A C C G C G A A
A T C A G C C T A
3. С правого края добавляется “лишняя” пара букв (F), а все позиции (пары букв)
нумеруются, начиная с левого края:
1 2 3 4
5 6 7
8 9 10
C A C C G C G A A F
A T C A G C C T A F
4. Каждой позиции такого выравнивания присваивается число в соответствии со
следующими правилами:
— последней по счету позиции (с парой букв F) присваивается число 0.
11
— если в позиции под номером 𝑛 буквы совпадают, то ей присваивается число
𝑆
𝑛+1
+𝑚𝑎𝑡𝑐ℎ, где 𝑆
𝑛+1
— число, присвоенное позиции 𝑛+1, 𝑚𝑎𝑡𝑐ℎ — натуральное
число, которое по смыслу является добавленным весом за совпадение букв.
Параметр 𝑚𝑎𝑡𝑐ℎ задается пользователем.
— если в позиции под номером 𝑛 буквы не совпадают, то ей присваивается число
𝑆
𝑛+1
+ 𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ, где 𝑆
𝑛+1
— число, присвоенное позиции 𝑛 + 1, 𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ — от-
рицательное целое число, которое по смыслу является штрафом за несовпадение
букв. Однако если 𝑆
𝑛+1
+ 𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ < 0, то позиции присваивается число 0.
Параметр 𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ задается пользователем.
Для примера в приведенном ниже выравнивании использованы следующие зна-
чения: 𝑚𝑎𝑡𝑐ℎ = 5,𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ = −4.
1 2
3
4
5
6 7
8 9 10
C A C C G C G A A F
A T C A G C C T A F
3 7 11 6 10 5 0
1 5
0
5. Выбирается позиция, которой присвоено наибольшее число (в данном случае это
позиция под номером 3). От этой позиции вправо начинается расширение вырав-
нивания, которое заканчивается первой встреченной позицией, которой присвоено
число 0 (в данном случае это позиция 7). Таким образом, в найденное в данном
случае локальное выравнивание попадают позиции от 3 до 6 включительно.
6. Каждому локальному выравниванию сопоставляется вес 𝑆, равный числу, при-
своенному первой позиции этого выравнивания (в нашем примере это число 11).
Среди найденных участков локального сходства отбираются те, для которых вес
𝑆 ≥ 𝑧 · 𝑚𝑎𝑡𝑐ℎ,
(1.5)
где 𝑧 — натуральное число, которое задается пользователем.
Если данное условие выполнено, то участок с найденным выравниванием “выре-
зается”, а участки, оставшиеся по краям, проверяются на длину. Если их длина
больше заданного пользователем параметра 𝐿 (по умолчанию 𝐿 = 10), то с ними
по отдельности проводится та же процедура, что и с исходным выравниванием,
и так далее, пока не останется ни одного непроверенного кусочка, большего или
равного 𝐿.
12
Слияние участков локального сходства
Отобранные со всех диагоналей участки локального сходства сортируются по
номеру позиции того нуклеотида, с которого начинаются эти участки (номер вычисля-
ется по последовательности банка данных). В том случае, если пара соседних участков
расположена достаточно близко, из запроса и из банка данных вырезаются последова-
тельности максимальной длины, лежащие между началом первого участка и концом
второго. Эти две последовательности выравниваются алгоритмом Смита — Ватермана
для локального выравнивания двух последовательностей с афинными штрафами за де-
леции. Афинные штрафы означают, что штраф за участок идущих подряд 𝑛 делеций
рассчитывается как сумма −𝐺𝑎𝑝𝑂𝑝𝑒𝑛 + (𝑛 · −𝐺𝑎𝑝𝐸𝑥𝑡𝑒𝑛𝑠𝑖𝑜𝑛), где 𝐺𝑎𝑝𝑂𝑝𝑒𝑛 — штраф
за открытие участка делеций, 𝐺𝑎𝑝𝐸𝑥𝑡𝑒𝑛𝑠𝑖𝑜𝑛 — штраф за каждую делецию в отдель-
ности. Оба этих параметра являются неотрицательными целыми числами и задаются
пользователем.
Если вес получившегося выравнивания больше, чем вес каждого из первоначаль-
ных участков, первый участок из пары удаляется, а второй заменяется на новое вырав-
нивание. Далее проверяется следующая пара, в которой второй участок из предыдущей
пары становится первым, и так далее до проверки последней пары.
Все оставшиеся участки, которые не были слиты с соседями, перевыравниваются
алгоритмом Смита — Ватермана.
Далее по формуле, предложенной в [1], для каждого выравнивания вычисляется
𝐸𝑣𝑎𝑙𝑢𝑒:
𝐸𝑣𝑎𝑙𝑢𝑒 = 𝐾𝑚𝑛𝑒
−𝜆𝑆
,
где 𝑚 и 𝑛 — длины последовательностей запроса и банка данных, а 𝐾 и 𝜆 — некоторые
постоянные, описанные в главе 2.
В выходной файл записываются все участки локального выравнивания с 𝐸𝑣𝑎𝑙𝑢𝑒,
значение которого не превышает пороговое значение, которое также задается пользо-
вательским параметром.
Программа Nhunt имеет некоторое количество параметров, которые позволяют
увеличить чувствительность программы, одновременно увеличив время ее работы. Од-
нако в большинстве случаев пользователю может понадобится только два таких пара-
метра. Первый — порог для отбора лучших диагоналей 𝑝 (см. формулу 1.2), который
имеет значение по умолчанию 10. Второй — минимальный вес для участка локального
сходства 𝑧 (см. формулу 1.5), который имеет значение по умолчанию 8.
Программа реализована на языке программирования C.
13
1.2
Подходы к решению проблемы малой сложности
Так называемая “проблема малой сложности” состоит в том, что выравнивание
двух последовательностей с сильно смещенным нуклеотидным составом будет иметь
высокий вес в силу случайных совпадений часто встречающихся нуклеотидов, а не в
силу гомологии последовательностей. Это может негативно отразиться на результатах
поиска, так как негомологичные последовательности могут иметь такой же или даже
более высокий вес по сравнению с гомологичными последовательностями.
Ниже описаны основные подходы, которые используются для решения проблемы
малой сложности как для нуклеотидных, так и для аминокислотных последовательно-
стей.
1.2.1
Маскировка участков малой сложности
Хронологически первым способом решения проблемы малой сложности была
маскировка участков малой сложности. Основная идея этого метода состоит в том,
что в последовательности банка данных заранее определяются участки малой слож-
ности, внутри которых поиск гомологов просто не проводится. Для аминокислотных
последовательностей данный способ считается устаревшим, однако для нуклеотидных
последовательностей он используется достаточно широко.
В частности, для решения проблемы малой сложности в программах семейства
BLASTN используется программа DUST [12]. Для любой последовательности из более
чем двух нуклеотидов вводится специальная мера сложности, которая основана на том,
сколько раз каждый из 64 возможных триплетов встретился в данной последователь-
ности. Затем при помощи динамического программирования во всех последовательно-
стях банка данных ищутся такие подпоследовательности, которые не содержат в себе
подпоследовательности с меньшей мерой сложности и которые одновременно обладают
мерой сложности, не превышающей определенный порог. Такие подпоследовательности
объявляются участками малой сложности, в которых не будет осуществляться поиск
гомологов.
1.2.2
Корректировка матрицы замен аминокислотных остатков
(метод Ю – Альтшуля)
Маскировка участков малой сложности, однако, не решает проблему малой слож-
ности в том случае, если интересующий нас гомолог находится как раз внутри тако-
го участка. Кроме того, на практике маскировка часто оказывается недостаточной и
поиск гомологов последовательности с сильно сдвинутым аминокислотным составом
14
выдаёт много ложных находок. Поэтому для аминокислотных последовательностей Ю
(Y.K.Yu), Альтшулем (S.K.Altschul) и их соавторами ( [8], [10], [11]) был разработан
алгоритм корректировки матрицы замен аминокислотных остатков.
Цель такой корректировки состоит в том, чтобы совпадение более частых букв
давало меньший вклад в вес выравнивания, чем совпадение более редких букв. Это
позволяет избежать нежелательного появления в выдаче программы поиска гомологов
выравниваний, которые имеют высокую значимость в основном за счет совпадения ча-
сто встречающихся букв. Кроме того, такой подход позволяет повысить значимость тех
выравниваний, в которых много совпадений редко встречающихся букв: такие совпа-
дения могут быть истолкованы как дополнительный признак гомологии.
В том случае, когда вероятности появления любого аминокислотного остатка
в каждой позиции последовательностей постоянны, произвольная исходная матрица
замен аминокислотных остатков может быть записана в виде:
𝑠
𝑖𝑗
=
1
𝜆
ln(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
),
(1.6)
где 𝑝
𝑖
, 𝑝
𝑗
— частоты встречаемости остатков 𝑖 и 𝑗, 𝜆 — параметр нормировки, а 𝑞
𝑖𝑗
(элемент матрицы частот 𝑞) — предполагаемая частота позиций правильного выравни-
вания, в которых остаток 𝑖 выровнен с остатком 𝑗 . Эта частота традиционно выводится
из частот сопоставленных пар остатков в наборе эталонных выравниваний; например,
для популярной серии матриц BLOSUM в качестве эталонных используются наборы
множественных выравниваний из базы данных BLOCKS.
В невырожденных случаях одна произвольная матрица замен аминокислотных
остатков соответствует только одному варианту распределения частот остатков в двух
выравниваемых последовательностях [8]. Это означает, что стандартная матрица замен
не подходит для расчета веса выравнивания последовательностей с заметно сдвинутым
аминокислотным составом.
Пусть 𝑃
𝑖
— частоты остатков в первой последовательности выравнивания, а 𝑃
′
𝑗
— во второй. Пусть также 𝑄
𝑖𝑗
— элемент новой, скорректированной матрицы частот 𝑄,
которую мы будем использовать вместо 𝑞. На 𝑄 накладываются естественные ограни-
чения, следующие из биологического смысла матрицы частот:
∑︁
𝑗
𝑄
𝑖𝑗
= 𝑃
𝑖
(1.7)
и
∑︁
𝑖
𝑄
𝑖𝑗
= 𝑃
′
𝑗
(1.8)
Это 2𝑁 − 1 независимое ограничение. В частности, из них следует, что ∑︀𝑄
𝑖𝑗
= 1.
15
Однако естественные ограничения 3.3 и 3.4 сами по себе еще не позволяют од-
нозначно получить матрицу 𝑄. Поэтому на 𝑄 также может быть наложено дополни-
тельное ограничение — минимизация величины относительной энтропии старой и новой
матриц частот аминокислотных остатков [9]. Это позволяет однозначно рассчитать но-
вые (скорректированные) значения матрицы замен аминокислотных остатков на основе
стандартной матрицы, используя вычисленные частоты остатков в двух выравнивае-
мых последовательностях.
Относительная энтропия старой и новой матриц частот вычисляется по формуле
𝐷(𝑄,𝑞) =
∑︁
𝑖,𝑗
𝑄
𝑖𝑗
ln(
𝑄
𝑖𝑗
𝑞
𝑖𝑗
)
Следовательно, искомая матрица 𝑄 должна удовлетворять условию минимума функции
𝐷(𝑄,𝑞):
𝑓(𝑄) =
∑︁
𝑖,𝑗
𝑄
𝑖𝑗
ln(
𝑄
𝑖𝑗
𝑞
𝑖𝑗
) → min
Кроме того, 𝑄 должна удовлетворять 2𝑁 − 1 естественным ограничениям, которые
могут быть записаны как функции
𝑔
𝑖
(𝑄) =
∑︁
𝑗
𝑄
𝑖𝑗
= 𝑃
𝑖
𝑔
𝑗
(𝑄) =
∑︁
𝑖
𝑄
𝑖𝑗
= 𝑃
′
𝑗
Для решения задачи нахождения матрицы 𝑄 предлагается применить метод мно-
жителей Лагранжа. Составим функцию Лагранжа в виде линейной комбинации функ-
ции 𝑓 и функций 𝑔
𝑖
(𝑄) и 𝑔
𝑗
(𝑄), взятыми с коэффициентами 𝛼
𝑖
и 𝛽
𝑗
соответственно,
называемыми множителями Лагранжа:
𝐿(𝑄,𝛼,𝛽) = 𝑓(𝑄) −
∑︁
𝑖
[𝛼
𝑖
(
∑︁
𝑗
𝑄
𝑖𝑗
− 𝑃
𝑖
)] −
∑︁
𝑗
[𝛽
𝑗
(
∑︁
𝑖
𝑄
𝑖𝑗
− 𝑃
′
𝑗
)]
Так как общее число ограничений равно 2𝑁 − 1, будем считать, что 𝛽
𝑁
= 0.
Затем мы приравниваем нулю частные производные функции Лагранжа по каждой из
переменных 𝑄
𝑖𝑗
и по каждому из множителей Лагранжа 𝛼
𝑖
и 𝛽
𝑗
:
𝛿𝐿
𝛿𝑄
𝑖𝑗
= 𝑙𝑛
𝑄
𝑖𝑗
𝑞
𝑖𝑗
+ 1 − 𝛼
𝑖
− 𝛽
𝑗
= 0
(1.9)
𝛿𝐿
𝛿𝛼
𝑖
= 𝑃
𝑖
−
∑︁
𝑗
𝑄
𝑖𝑗
= 0
(1.10)
16
𝛿𝐿
𝛿𝛽
𝑗
= 𝑃
′
𝑗
−
∑︁
𝑖
𝑄
𝑖𝑗
= 0
(1.11)
Сделаем замену переменных: 𝐴
𝑖
= 𝑒
𝛼
𝑖
−1
и 𝐵
𝑗
= 𝑒
𝛽
𝑗
. Тогда из 1.9 следует, что
𝑄
𝑖𝑗
= 𝐴
𝑖
𝐵
𝑗
𝑞
𝑖𝑗
(1.12)
Подставляя полученное выражение для 𝑄
𝑖𝑗
в 1.10 и 1.11, получаем:
⎧
⎪
⎪
⎨
⎪
⎪
⎩
𝑃
𝑖
= 𝐴
𝑖
∑︁
𝑗
𝐵
𝑗
𝑞
𝑖𝑗
𝑃
′
𝑗
= 𝐵
𝑗
∑︁
𝑖
𝐴
𝑖
𝑞
𝑖𝑗
(1.13)
Решение такой системы относительно переменных 𝐴
𝑖
и 𝐵
𝑗
может быть получено
в несколько итераций методом многомерного приближения Ньютона. После этого зна-
чения 𝑄
𝑖𝑗
выводятся в соответствии с 1.12. Затем элементы 𝑆
𝑖𝑗
новой матрицы замен 𝑆
по аналогии с формулой 1.6 вычисляется так:
𝑆
𝑖𝑗
=
1
𝜆
ln(
𝑄
𝑖𝑗
𝑃
𝑖
𝑃
′
𝑗
),
(1.14)
Существует также модификация данного метода корректировки матрицы замен.
Она состоит в том, чтобы дополнительно к требованию минимизации относительной
энтропии старой и новой матриц сохранять также энтропию самой матрицы частот [10].
Энтропия матрицы частот, или ее информационная наполненность, вычисляется
по формуле
𝐻 =
∑︁
𝑖,𝑗
𝑞
𝑖𝑗
𝑙𝑛(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
)
Тогда дополнительным требованием к скорректированной матрице частот 𝑄 станет со-
хранение информационной наполненнности 𝐻 исходной матрицы 𝑞:
∑︁
𝑖,𝑗
𝑄
𝑖𝑗
𝑙𝑛(
𝑄
𝑖𝑗
𝑃
𝑖
𝑃
′
𝑗
) = 𝐻
Новая матрица замен вычисляется аналогично тому, как это происходило в предыду-
щем случае, когда информационная наполненность матрицы частот не сохранялась.
Отличия состоят лишь в том, что в функции Лагранжа учитывается еще одно ограни-
чение для сохранения информационной наполненности, а многомерный метод Ньюто-
на применяется непосредственно к системе частных производных функции Лагранжа.
Этот метод, в частности, по умолчанию используется для корректировки матрицы за-
мен в популярной программе поиска гомологов аминокислотных последовательностей
BLASTP [9].
17
Однако показано [10], что при поиске гомологов аминокислотных последователь-
ностей применение дополнительного ограничения — сохранения информационной на-
полненности — не приводит к значимому улучшению результатов поиска.
Таким образом, для аминокислотных последовательностей существует метод, ко-
торый более эффективен для борьбы с малой сложностью, чем метод маскировки. Этот
метод позволяет скорректировать матрицу замен аминокислотных остатков и более
точно оценивать значимость найденных выравниваний. При этом методы, которые ис-
пользуются для корректировки матрицы замен аминокислотных остатков, формально
полностью применимы и к матрице замен нуклеотидов. Тем не менее, все известные
нам программы поиска гомологов нуклеотидных последовательностей либо никак не
решают проблему малой сложности, либо маскируют участки малой сложности таким
образом, чтобы эти участки не появлялись в результатах поиска.
18
Глава
2
МАТЕРИАЛЫ
И
МЕТОДЫ
2.1
Средства программной реализации
Программа поиска гомологов нуклеотидных последовательностей Nhunt реали-
зована на языке программирования C. Кроме стандартных кроссплатформенных биб-
лиотек, была использована также библиотека unistd.h, которая является стандартной
для Unix-платформ.
Пакет сценариев для сравнения различных программ поиска гомологов нуклео-
тидных последовательностей реализован на языке программирования Python 2.7.
2.2
Материалы тестирования методов поиска
гомологов
Для тестирования различных программ поиска гомологов нуклеотидных после-
довательностей применялось несколько наборов данных. Одни наборы данных были
использованы как последовательности запроса, а другие — как последовательности
банка данных. Для удобства каждый набор данных имеет короткий идентификатор
на латинице, который используется в главе 3 при описании процедуры и результатов
тестирования. Ниже для каждого из этих наборов данных указано, что за последова-
тельности в них содержатся, каков их размер и GC-состав (то есть суммарный процент
нуклеотидов G и С среди общего числа нуклеотидов последовательности).
2.2.1
Выборки последовательностей запроса
В качестве последовательностей запроса были взяты выборки аннотированных
РНК из геномов различных организмов.
19
1. Выборка тРНК из генома E. coli
В данную выборку входят все аннотированные транспортные РНК из генома E.
coli str. K-12 (запись EMBL AP009048). В случаях, если две или более РНК вы-
борки полностью совпадали по последовательности, в выборке оставлялась только
одна из этих РНК.
Общая длина последовательностей выборки — 3808 нуклеотидов. Количество по-
следовательностей — 48. Средний GC-состав по всем последовательностям — 58.9%,
минимальное значение — 49.3%, максимальное значение — 67.5%.
В описании результатов идентификатор этой выборки — ecoli_tRNA.
2. Выборка тРНК из генома митохондрии человека
В данную выборку входят все аннотированные транспортные РНК из митохон-
дриального генома Homo sapiens (запись EMBL FJ157845). В случаях, если две
или более РНК выборки полностью совпадали по последовательности, в выборке
оставлялась только одна из этих РНК.
Общая длина последовательностей выборки — 1512 нуклеотидов. Количество по-
следовательностей — 22. Средний GC-состав по всем последовательностям — 51.4%,
минимальное значение — 40.9%, максимальное значение — 61.9%.
В описании результатов идентификатор этой выборки — mit_tRNA.
3. Выборка тРНК из генома малярийного плазмодия Plasmodium knowlesi
В данную выборку входят все аннотированные транспортные РНК из генома
Plasmodium knowlesi. В случаях, если две или более РНК выборки полностью сов-
падали по последовательности, в выборке оставлялась только одна из этих РНК.
Общая длина последовательностей выборки — 3219 нуклеотидов. Количество по-
следовательностей — 43. Средний GC-состав по всем последовательностям — 56.5%,
минимальное значение — 49.3%, максимальное значение — 62.5%.
В описании результатов идентификатор этой выборки — plasm_tRNA.
4. Выборка фрагментов рРНК из генома малярийного плазмодия Plasmodium
falciparum
В данную выборку входят 18S и 25S рибосомальные РНК из генома Plasmodium
falciparum (запись EMBL X95275). При этом каждая из последовательностей де-
лилась на фрагменты длиной 300 нуклеотидов. Длина последовательностей не
кратна 300, поэтому их оставшиеся фрагменты, имеющие длину более 220 нук-
леотидов, также включались в выборку.
20
Общая длина последовательностей выборки — 4122 нуклеотидов. Количество по-
следовательностей — 14. Средний GC-состав по всем последовательностям — 21%,
минимальное значение — 13.3%, максимальное значение — 33.7%.
В описании результатов идентификатор этой выборки — plasm_rRNA.
2.2.2
Банки данных
В качестве последовательностей банков данных были взяты выборки полных
геномов различных прокариотических организмов. Каждая выборка состоит либо из
геномов бактерий, либо из геномов архей. Одновременно каждая выборка содержит
либо последовательности геномов с высоким GC-составом (GC-богатые геномы), либо
с низким GC-составом (GC-бедные), либо с GC-составом, близким к среднему (50%).
В таблицах, которые описывают банки данных, применяются следующие обозна-
чения: Организм — название организма, геном которого находится в выборке. EMBL
ID — идентификатор записи банка EMBL, которая содержит данный геном. Длина —
длина последовательности данного генома. GC-состав — GC-состав последовательно-
сти данного генома.
1. Выборка GC-богатых геномов бактерий
В данную выборку входят несколько полных геномов бактерий, последовательно-
сти которых имеют относительно высокий GC-состав. Полный список этих гено-
мов и их характеристики приведены в таблице (2.1).
Таблица 2.1
Геномы, входящие в выборку GC-богатых геномов бактерий
Организм
EMBL ID
Длина
GC-состав
Sphingobium sp.
AP012222 4199332
65.57
Stenotrophomonas maltophilia HE798556 4769156
66.75
Xanthomonas campestris
CP000050 5148708
64.96
Caulobacter crescentus
AE005673 4016947
67.21
Kineococcus radiotolerans
CP000750 4761183
74.44
Mycobacterium intracellulare
CP003324 5501090
68.07
Общая длина последовательностей выборки — 28 396 416 нуклеотидов.
В описании результатов идентификатор этой выборки — rich_bact.
2. Выборка GC-богатых геномов архей
В данную выборку входят несколько полных геномов архей, последовательности
21
которых имеют относительно высокий GC-состав. Полный список этих геномов и
их характеристики приведены в таблице (2.2).
Таблица 2.2
Геномы, входящие в выборку GC-богатых геномов архей
Организм
EMBL ID
Длина
GC-состав
Halobacterium salinarum AM774415 2000962
68.01
Halopiger xanaduensis
CP002839 3668009
65.98
Halobacterium sp.
AE004437 2014239
67.91
Haloferax volcanii
CP001956 2847757
66.64
Общая длина последовательностей выборки — 10 530 967 нуклеотидов.
В описании результатов идентификатор этой выборки — rich_ar.
3. Выборка геномов бактерий, имеющих GC-состав, близкий к среднему
В данную выборку входят несколько полных геномов бактерий, последовательно-
сти которых имеют GC-состав, близкий к среднему. Полный список этих геномов
и их характеристики приведены в таблице (2.3).
Таблица 2.3
Геномы выборки бактерий, имеющие GC-состав, близкий к среднему
Организм
EMBL ID
Длина
GC-состав
Xylella fastidiosa
AE003849 2679306
52.67
Yersinia pestis
AE009952 4600755
47.64
Shigella flexneri
AE005674 4607202
50.89
Pantoea ananatis
HE617160 4605545
53.45
Anaplasma marginale
CP001079 1202435
49.77
Acetobacter pasteurianus AP011121 2907495
53.04
Общая длина последовательностей выборки — 20 602 738 нуклеотидов.
В описании результатов идентификатор этой выборки — normal_bact.
4. Выборка геномов архей, имеющих GC-состав, близкий к среднему
В данную выборку входят несколько полных геномов архей, последовательности
которых имеют GC-состав, близкий к среднему. Полный список этих геномов и их
характеристики приведены в таблице (2.4).
22
Таблица 2.4
Геномы выборки архей, имеющие GC-состав, близкий к среднему
Организм
EMBL ID
Длина
GC-состав
Haloquadratum walsbyi
AM180088 3132494
47.86
Methanothermobacter thermautotrophicus AE000666 1751377
49.54
Methanocella paludicola
AP011532 2957635
54.92
Archaeoglobus fulgidus
AE000782 2178400
48.58
Pyrobaculum aerophilum
AE009441 2222430
51.36
Общая длина последовательностей выборки — 12 242 336 нуклеотидов.
В описании результатов идентификатор этой выборки — normal_ar.
5. Выборка GC-бедных геномов бактерий
В данную выборку входят несколько полных геномов бактерий, последовательно-
сти которых имеют относительно низкий GC-состав. Полный список этих геномов
и их характеристики приведены в таблице (2.5).
Таблица 2.5
Геномы, входящие в выборку GC-бедных геномов бактерий
Организм
EMBL ID
Длина
GC-состав
Mycoplasma bovis
CP002188 1003404
29.31
Buchnera aphidicola
AE013218
641454
25.33
Mycoplasma leachii
CP002108 1008951
23.75
Acholeplasma laidlawii CP000896 1496992
31.93
Sulcia muellery
CP002163
276511
21.13
Blochmannia floridanus BX248583 705557
27.38
Общая длина последовательностей выборки — 5 132 869 нуклеотидов.
В описании результатов идентификатор этой выборки — low_bact.
6. Выборка GC-бедных геномов архей
В данную выборку входят несколько полных геномов архей, последовательности
которых имеют относительно низкий GC-состав. Полный список этих геномов и
их характеристики приведены в таблице (2.6).
23
Таблица 2.6
Геномы, входящие в выборку GC-бедных геномов архей
Организм
EMBL ID
Длина
GC-состав
Nanoarchaeum equitans
AE017199
490885
31.56
Methanocaldococcus vulcanius
CP001787 1746329
31.49
Methanobrevibacter ruminantium CP001719 2937203
32.64
Methanothermococcus okinawensis CP002792 1662525
29.30
Sulfolobus tokodaii
BA000023 2694756
32.79
Methanosphaera stadtmanae
CP000102 1767403
27.63
Methanococcus voltae
CP002057 1936387
28.59
Общая длина последовательностей выборки — 13 235 488 нуклеотидов.
В описании результатов идентификатор этой выборки — low_ar.
Отметим, что для всех перечисленных полных геномов аннотированы последо-
вательности как транспортных, так и рибосомальных РНК.
24
Глава
3
РЕЗУЛЬТАТЫ
И
ОБСУЖДЕНИЕ
3.1
Сравнение программ поиска гомологов
Сравнение различных способов корректировки матрицы замен нуклеотидов под-
разумевает сравнение результатов работы одной или нескольких программ поиска го-
мологов нуклеотидных последовательностей, запущенных с различными вариантами
корректировки или вообще без таковой. Таким образом, необходимо понять, каким об-
разом оценивать качество работы программы поиска гомологов.
На сегодняшний день не существует общепризнанного критерия, по которому
можно судить о качестве результатов работы программы поиска гомологов. Сравнение
количества или суммарной длины выравниваний из выдачи программы, имеющих вели-
чину 𝐸𝑣𝑎𝑙𝑢𝑒 меньше заданного порога, в общем случае не подходит на роль критерия,
так как сама величина 𝐸𝑣𝑎𝑙𝑢𝑒 в разных программах вычисляется по-разному. Кроме
того, сравнение любых численных свойств выданных выравниваний, таких как вес вы-
равнивания или его 𝐸𝑣𝑎𝑙𝑢𝑒, не дает нам представления о биологической значимости
найденных гомологов.
3.1.1
Сценарии для сравнения качества программ поиска
гомологов
Для оценки качества работы программ поиска гомологов нуклеотидных после-
довательностей нами был разработан набор сценариев на языке программирования
Python.
Эти сценарии сопоставляют координаты участков банка данных, найденных про-
граммой поска, с аннотацией того же банка данных. Для этого на вход сценариям, по-
мимо последовательностей запроса и аннотированных последовательностей банка дан-
ных, даётся информация о том, какие именно участки банка данных аннотированны
25
как гомологи последовательностей запроса. Такие участки считаются верными наход-
ками, а все остальные — ложными. На основании анализа аннотаций находок выдаются
графические и численные характеристики качества поиска.
Перед началом работы с указанным набором сценариев необходимо подготовить
и поместить в соответствующие директории следующие данные:
1. Файл, содержащий одну или несколько последовательностей запроса в формате
FASTA.
2. Файл, содержащий одну или несколько последовательностей базы данных в фор-
мате FASTA.
3. Файлы, содержащие информацию об участках, которые содержатся в последова-
тельностях банка данных и гомологичны любой из последовательности запроса.
Эта информация включает в себя координаты участка и направление цепи, на
которой он лежит. Каждой последовательности банка данных соответствует один
такой файл.
В ходе наших исследований последовательностями запроса являлись последова-
тельности транспортных РНК или фрагменты последовательностей рибосомаль-
ных РНК различных организмов. Последовательности банка данных обычно пред-
ставляли из себя полные геномы прокариот.
Информация о гомологичных участках, содержащихся в последовательностях бан-
ка данных, автоматически извлекается из аннотаций геномов, конкретно — из
файлов в формате EMBL или GenBank.
4. Файл конфигурации, содержащий информацию о том, какие программы и с каки-
ми параметрами будут запущены для проверки качества выданных ими результа-
тов. Здесь же указывается, на каком наборе входных данных (последовательности
запроса и банка данных) должен быть произведен поиск гомологов.
Затем все указанные пользователем программы с заданными им наборами пара-
метров начинают проводить поиск гомологов входных последовательтей в банке дан-
ных. Координаты найденных гомологов сравниваются с подготовленными координата-
ми гомологичных участков (в нашем случае гомологи тРНК запроса сравнивались с ан-
нотированными тРНК банка данных, а гомологи рРНК — с аннотированными рРНК).
Для каждого найденного выравнивания подсчитывается длина “правильного” фрагмен-
та (то есть в нашем случае число нуклеотидов найденного программой гомолога, лежа-
щих в пределах координат любой аннотированной тРНК/рРНК) и длина “неправиль-
ного” фрагмента (то есть в нашем случае число нуклеотидов найденного программой
гомолога, лежащих вне предела координат любой аннотированной тРНК/рРНК).
26
Затем для каждого выравнивания рассчитывается две величины. Первая вели-
чина — 𝑅𝑣𝑎𝑙𝑢𝑒 — равна суммарной длине “правильных” фрагментов для всех выравни-
ваний, имеющих вес больший, чем данное выравнивание. Вторая величина — 𝑊𝑣𝑎𝑙𝑢𝑒 —
равна суммарной длине “неправильных” фрагментов для всех выравниваний, имеющих
вес больший, чем данное выравнивание.
Все выравнивания, полученные каждой программой (набором параметров) в от-
дельности, сортируются по величине веса, выданной программой (набором параметров)
для каждого выравнивания. На основании этих данных строятся двумерные графики,
на которых каждому найденному выравниванию соответствует точка определенного
цвета. Каждая совокупность точек (выравниваний), имеющих одинаковый цвет, соот-
ветствует одной программе (набору параметров).
Значение ординаты каждой точки (выравнивания) численно равно величине 𝑅𝑣𝑎𝑙𝑢𝑒
для данного выравнивания. Значение абсциссы численно равно величине 𝑊𝑣𝑎𝑙𝑢𝑒 для
данного выравнивания. Точки, обозначающие самые лучшие (имеющие наибольший
вес) выравнивания, находятся в начале координатной сетки. Чем выше на координат-
ной плоскости расположена получившаяся совокупность точек (выравниваний), тем
лучше работает поиск.
Описанный выше метод позволяет наглядно оценить биологическую значимость
выравниваний, найденных программой поиска гомологов. Кроме того, он также позво-
ляет сравнить между собой качество работы нескольких программ (наборов парамет-
ров) на одинаковых наборах входных данных.
3.1.2
Примеры результатов сравнения
Для тестирования описанного выше метода сравнения программ поиска гомоло-
гов мы использовали различные программы поиска гомологов нуклеотидных последо-
вательностей, такие как BLASTN, discontiguous MEGABLAST, FASTA, YASS и Nhunt.
В некоторых случаях результаты, выдаваемые одной программой, сильно различались
в зависимости от набора параметров программы. Поэтому мы также сравнивали эф-
фективность поиска гомологов для различных наборов параметров одной программы.
Сравнение программ семейства BLASTN
Сравнение между собой программ семейства BLASTN представляет большой ин-
терес. Во-первых, на сегодняшний день это наиболее популярные программы для поис-
ка гомологов нуклеотидных последовательностей. Во-вторых, именно для программы
BLASTN наблюдается наибольшая зависимость качества поиска гомологов от набора
входных параметров.
27
Во всех случаях применения программ семейства BLASTN использовался пакет
BLAST+ версии 2.2.25.
Ниже перечислены программы и наборы параметров, которые были взяты для
сравнения программ семейства BLASTN. Для каждой программы также указан ко-
роткий идентификатор на латинице, который используется в графике 3.1. При этом
опущена та область графика, в которой характер поведения кривых (наборов вырав-
ниваний) для каждой программы (или набора параметров) не меняется относительно
видимой части графика (то есть за пределами видимой части графика новые “правиль-
ные” выравнивания практически на встречаются).
1. Программа BLASTN с наиболее чувствительными параметрами
Программа BLASTN со следующим набором параметров: –dust no –word_size 4
–penalty –4 –reward 5 –gapopen 10 –gapextend 6. По нашим наблюдением, именно
такой набор параметров позволяет добиться наибольшей чувствительности для
программы BLASTN. Данный набор параметров доступен только в локальной
версии программы BLASTN. Короткий идентификатор — blast_best.
2. Программа BLASTN с наиболее чувствительными параметрами, до-
ступными на ее веб-интерфейсе
Программа BLASTN со следующим набором параметров: –dust no –word_size 7
–penalty –1 –reward 1 –gapopen 5 –gapextend 2. При запуске программы BLASTN с
сайта веб-интерфейса, доступного по адресу [15], пользователю на выбор предо-
ставляется лишь несколько заранее предустановленных наборов параметров. По
нашим наблюдениям, из этих доступных наборов параметров именно указанный
набор параметров позволяет добиться наибольшей чувствительности для програм-
мы BLASTN. Короткий идентификатор — blast_web.
3. Программа BLASTN с параметрами по умолчанию
Программа BLASTN со следующим набором параметров: –dust no –word_size 7 –
penalty –1 –reward 1 –gapopen 5 –gapextend 2. Данный набор параметров установлен
по умолчанию как в локальной версии программы BLASTN, так и на сайте ее веб-
интерфейса. Короткий идентификатор — blast_def.
4. Программа discontiguous MEGABLAST
Программа discontiguous MEGABLAST со следующим набором параметров:
–template_length 16 –dust no –template_type optimal –penalty –4 –reward 5 –gapopen
10 –gapextend 6. Данный набор выставлен по аналогии с наиболее чувствительны-
ми параметрами для программы BLASTN. Тем не менее, по нашим наблюдениям,
выбор различных наборов параметров программы практически не влияет на ка-
чество поиска гомологов. Короткий идентификатор — megablast.
28
Все вышеперечисленный программы с соответствующими наборами параметров
были использованы для поиска гомологов транспортных РНК из генома плазмодия
(выборка plasm_tRNA) в геномах различных бактерий с GC-составом, близким к нор-
мальному (выборка normal_bact). В результате обработки полученных выравниваний
и координат гомологов в соответствии с процедурой, описанной в разделе 3.1.1, был по-
лучен график 3.1, отражающий качество поиска гомологов различными программами
(наборами параметров).
Рис. 3.1
Сравнение программ семейства BLASTN. blast_best — результаты работы программы
BLASTN с наиболее чувствительными параметрами, время работы — 84 сек. blast_def
— результаты работы программы BLASTN с параметрами по умолчанию, время работы
— 6.9 сек. blast_web — результаты работы программы BLASTN с наиболее чувстви-
тельными параметрами, доступными на ее веб-интерфейсе, время работы — 9.5 сек.
megablast — результаты работы программы discontiguous MEGABLAST, время рабо-
ты — 13.2 сек. По оси ординат для каждой точки (выравнивания) отложена соответ-
ствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс — величина 𝑊𝑣𝑎𝑙𝑢𝑒.
Данный график является типичным для всех возможных пар выборок входных
последовательностей и последовательностей банка данных. Можно заключить, что па-
29
раметры, выставленные в программе BLASTN по умолчению, являются оптимальными
с точки зрения времени работы программы, но не с точки зрения качества поика гомо-
логов. Аналогично, с точки зрения качества поиска гомологов не оптимальны наборы
параметров, доступные на сайте веб-интерфейса программы BLASTN. Таким образом,
лучшее качество поиска гомологов нуклеотидных последовательностей доступно поль-
зователю только при установке локального варианта программы BLASTN.
Необъясненным остается тот факт, что программа discontiguous MEGABLAST,
призванная улучшить чувствительность поиска гомологов по отношению к программе
BLASTN, на деле выдает чрезвычайно малое количество находок — как правильных,
так и неправильных. Возможно, это следствие других эвристик программы, или же в
данном случае сказывается какая-либо неисправленная ошибка.
Сравнение различных программ поиска гомологов
Кроме того, между собой были сравнены различные существующие программы
поиска гомологов. Параметры в них подбирались таким образом, чтобы обеспечить
максимальную чувствительность программы (в программе Nhunt — чтобы обеспечить
оптимальное соотношение между чувствительностью и скоростью работы). При ука-
зании параметров упомянуты лишь те из них, которые отличаются от параметров по
умолчанию.
Ниже перечислены программы и наборы параметров, которые были взяты для
сравнения различных программ поиска гомологов. Для каждой программы также ука-
зан короткий идентификатор на латинице, который используется в графике 3.2.
1. Программа BLASTN с наиболее чувствительными параметрами
Программа BLASTN со следующим набором параметров: –dust no –word_size 4
–penalty –4 –reward 5 –gapopen 10 –gapextend 6. По нашим наблюдением, именно
такой набор параметров позволяет добиться наибольшей чувствительности для
программы BLASTN. Короткий идентификатор — blast_best.
2. Программа FASTA
Программа FASTA с параметром ktup = 3 (см. раздел 1.1.1). По нашим наблюде-
нием, именно такой набор параметров позволяет добиться наибольшей чувстви-
тельности для программы FASTA. Короткий идентификатор — fasta.
3. Программа Nhunt
Программа Nhunt взята с параметрами по умолчанию. Короткий идентификатор
— nhunt.
30
4. Программа YASS
Программа YASS взята с параметрами по умолчанию. По нашим наблюдениям,
выбор других наборов параметров программы практически не влияет на качество
поиска гомологов. Короткий идентификатор — yass.
Все вышеперечисленные программы с соответствующими наборами параметров
были использованы для поиска гомологов транспортных РНК из генома плазмодия
(выборка plasm_tRNA) в геномах различных бактерий с GC-составом, близким к нор-
мальному (выборка normal_bact). В результате обработки полученных выравниваний
и координат гомологов в соответствии с процедурой, описанной в разделе 3.1.1, был по-
лучен график 3.2, отражающий качество поиска гомологов различными программами
(наборами параметров). При этом опущена та область графика, в которой характер по-
ведения кривых (наборов выравниваний) для каждой программы (или набора парамет-
ров) не меняется относительно видимой части графика (то есть за пределами видимой
части графика новые “правильные” выравнивания практически на встречаются).
31
Рис. 3.2
Сравнение различных программ поиска гомологов. blast_best — результаты работы
программы BLASTN с наиболее чувствительными параметрами, время работы — 84
сек. fasta — результаты работы программы FASTA, время работы — 73 сек. yass —
результаты работы программы YASS, время работы — 310 сек. nhunt — результаты
работы программы Nhunt, время работы — 340 сек. По оси ординат для каждой точ-
ки (выравнивания) отложена соответствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс —
величина 𝑊𝑣𝑎𝑙𝑢𝑒.
Данный график является типичным для всех возможных пар выборок вход-
ных последовательностей и последовательностей банка данных. Можно заключить, что
программа Nhunt, хотя и не является оптимальной с точки зрения времени работы,
позволяет заметно улучшить качество поиска гомологов по сравнению с программами
BLASTN и YASS. Тот факт, что программа FASTA выдает чрезвычайно малое коли-
чество правильных находок, требует отдельного анализа, выходящего за рамки данной
работы.
32
3.2
Адаптация метода Ю – Альтшуля для матрицы
замен нуклеотидов
Для аминокислотных последовательностей можно использовать так называемый
метод Ю – Альтшуля [9]. Этот метод основан на минимизации величины относительной
энтропии старой и новой матриц частот аминокислотных остатков (см. раздел 1.2.2).
Однако все математические выкладки, применяемые для расчета скорректированной
матрицы замен аминокислотных остатков, формально могут быть применены и для
расчета скорректированной матрицы замен нуклеотидов.
В случае, если вероятности появления любого нуклеотида в каждой позиции
последовательности постоянны и равны, произвольная исходная матрица замен нук-
леотидов может быть записана в виде:
𝑠
𝑖𝑗
=
1
𝜆
ln(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
),
(3.1)
где 𝑝
𝑖
= 𝑝
𝑗
=
1
4
— частоты встречаемости нуклеотидов 𝑖 и 𝑗, 𝜆 — параметр нормиров-
ки, а 𝑞
𝑖𝑗
(элемент матрицы частот 𝑞) — предполагаемая частота позиций правильного
выравнивания, в которых нуклеотид 𝑖 выровнен с нуклеотидом 𝑗.
В отличие от аминокислотных последовательностей, для нуклеотидных после-
довательностей не существует общепринятых эталонных выравниваний, на основании
которых можно было бы напрямую расчитать значения 𝑞
𝑖𝑗
. Поэтому традиционно эти
величины рассчитываются на основе 𝑠
𝑖𝑗
— элементов матрицы замен нуклеотидов 𝑠. В
соответствии с формулой 1.6, элемент матрицы частот 𝑞
𝑖𝑗
может быть вычислен как
𝑞
𝑖𝑗
= 𝑝
𝑖
𝑝
𝑗
exp(𝜆𝑠
𝑖𝑗
),
(3.2)
а 𝜆 — восстановлено как единственное положительное решение уравнения ∑︀𝑞
𝑖𝑗
= 1.
Пусть 𝑃
𝑖
— частоты нуклеотидов в первой последовательности выравнивания,
а 𝑃
′
𝑗
— во второй. Пусть также 𝑄
𝑖𝑗
— элемент новой, скорректированной матрицы ча-
стот 𝑄, которую мы будем использовать вместо 𝑞. На 𝑄 накладываются естественные
ограничения, следующие из биологического смысла матрицы частот:
∑︁
𝑗
𝑄
𝑖𝑗
= 𝑃
𝑖
(3.3)
и
∑︁
𝑖
𝑄
𝑖𝑗
= 𝑃
′
𝑗
(3.4)
Это 2𝑁 − 1 независимое ограничение. В частности, из них следует, что ∑︀𝑄
𝑖𝑗
= 1.
33
Для упрощения и ускорения вычисления скорректированной матрицы замен нук-
леотидов нами был несколько изменен способ получения скорректированной матрицы
частот 𝑄, описанный в [9]. Вместо того, чтобы применять к системе уравнений 1.13
многомерный метод Ньютона, искомые значения 𝐴
𝑖
и 𝐵
𝑗
могут быть получены мето-
дом последовательных приближений.
Для реализации метода последовательных приближений всем переменным 𝐵
𝑗
(кроме 𝐵
𝑁
) сперва присваивается начальное значение 1. Затем вычисляются начальные
значения для переменных 𝐴
𝑖
:
𝐴
𝑖
=
𝑃
𝑖
∑︀
𝑗
𝐵
𝑗
𝑞
𝑖𝑗
(3.5)
После этого вычисляются новые значения для переменных 𝐵
𝑗
(кроме 𝐵
𝑁
):
𝐵
𝑗
=
𝑃
′
𝑗
∑︀
𝑖
𝐴
𝑖
𝑞
𝑖𝑗
(3.6)
Затем значения 𝐴
𝑖
и 𝐵
𝑗
поочередно вычисляются по формулам 3.5 и 3.6 соответственно.
Итерации прекращаются, когда все новые значения переменных отличаются от значе-
ний предыдущего шага не более, чем на наперед заданную величину, много меньшую
единицы. На практике сходимость достигается достаточно быстро, в пределах несколь-
ких десятков итераций.
После решения системы 1.13 значения скорректированной матрицы замен нахо-
дятся аналогично тому, как было описано выше: элементы 𝑄
𝑖𝑗
матрицы частот нуклео-
тидов выводятся в соответствии с 1.12, а новая матрица замен вычисляется с помощью
1.6.
Для сравнения относительной скорости исходного и модифицированного метода
корректировки матрицы замен оба метода были запущены для расчета матрицы частот
нуклеотидов при одинаковых исходных данных. Каждый метод корректировки был
запущен 100 000 раз при различных частотах нуклеотидов последовательностей запроса
и банка данных.
Результатом каждого запуска было время, за которое получено решение — зна-
чения матрицы частот нуклеотидов. Оба метода являются итерационными, то есть чем
большее количество итераций проводится, тем ближе получаемые значения частот к
точному решению. Итерации каждого из методов в тестах прекращались, как только
все текущие значения матрицы частот отличались от значений, полученных при преды-
дущей итерации, не более чем на 10
−5
.
Для каждой пары методов и наборов исходных данных запуск и подсчет времени
работы проводились 10 раз, после чего высчитывалось среднее время работы. Это время
представлено в таблице (3.1). Исходные частоты последовательностей запроса и банка
данных для краткости выражены в виде GC-состава.
34
Таблица 3.1
Сравнение скорости работы исходного и модифицированного методов Ю – Альтшуля
(100 000 запусков)
Запрос Банк Модифицированный метод Исходный метод
50%
50%
1.53 сек
3.2 сек
70%
50%
1.5 сек
4.2 сек
70%
70%
2.0 сек
3.7 сек
70%
30%
2.0 сек
4.4 сек
Запрос — GC-состав последовательности запроса, Банк — GC-состав
последовательности банка данных, Модифицированный метод — время работы
модифицированного метода корректировки, Исходный метод — время работы
исходного метода корректировки.
Таким образом, модифицированный метод корректировки во всех протестиро-
ванных случаях ускоряет процесс нахождения решения примерно в два раза.
3.3
Новые методы корректировки матрицы замен
Нами также разработано два варианта оригинального метода корректировки
матрицы замен. Оба они основаны на рассмотрении не всего пространства матриц ча-
стот, а лишь матриц специального вида, в которых элемент матрицы 𝑞
𝑖𝑗
явным образом
зависит от соответствующих частот нуклеотидов 𝑝
𝑖
и 𝑝
𝑗
. Отличие двух способов друг
от друга состоит в том, что в первом случае скорректированная матрица замен сохра-
няет ожидаемый счет исходной матрицы замен, а во втором случае скорректированная
матрица частот сохраняет относительную энтропию (информационную наполненность)
исходной матрицы частот.
Положим 𝑟
𝑖𝑗
= 𝑞
𝑖𝑗
− 𝑝
𝑖
𝑝
𝑗
, где 𝑞
𝑖𝑗
, 𝑝
𝑖
и 𝑝
𝑗
— те же, что в формуле (1.6). Из этого,
в частности, следует, что ∑︀
𝑖
𝑟
𝑖𝑗
= ∑︀
𝑗
𝑟
𝑖𝑗
= 0.
Пусть 𝑃
𝑖
— частота встречаемости нуклеотидов в первой выравниваемой после-
довательности, а 𝑃
′
𝑗
— во второй. Положим также
𝑄
𝑖𝑗
= 𝑃
𝑖
𝑃
′
𝑗
+ 𝑐𝑟
𝑖𝑗
где 𝑐 > 0 — параметр, который зависит от 𝑃
𝑖
и 𝑃
′
𝑗
. Такой вид матрицы удобен тем,
что при любом 𝑐 выполняются условия ∑︀
𝑗
𝑄
𝑖𝑗
= 𝑃
𝑖
и ∑︀
𝑖
𝑄
𝑖𝑗
= 𝑃
′
𝑗
, то есть естественные
условия, налагаемых на матрицу 𝑄. Замена 𝑝
𝑖
на 𝑃
𝑖
, 𝑝
𝑗
на 𝑃
′
𝑗
и 𝑞
𝑖𝑗
на 𝑄
𝑖𝑗
в формуле (1.6)
позволяет сразу получить элементы скорректированной матрицы замен нуклеотидов
𝑆
𝑖𝑗
.
35
3.3.1
Сохранение ожидаемого счета
Ниже описан способ корректировки матрицы замен нуклеотидов, основанный на
сохранении ожидаемого счета новой матрицы.
Одним из наиболее важных параметров, определяющих матрицу замен нуклео-
тидов, является ее ожидаемый счет 𝐸 (среднее значение веса сопоставления двух слу-
чайно взятых букв). Ожидаемый счет вычисляется по формуле
𝐸 =
∑︁
𝑖,𝑗
𝑝
𝑖
𝑝
𝑗
𝑠
𝑖𝑗
=
∑︁
𝑖,𝑗
𝑝
𝑖
𝑝
𝑗
1
𝜆
ln(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
).
Формула 1.1 выведена при предположении отрицательных значений величины
𝐸 [1]. Кроме того, очевидно, что при положительных значених 𝐸 выравнивания по
случайным причинам будут раширяться на негомологичные участки.
В практике выравниваний нуклеотидных последовательностей матрицы замен
с различными значениями отрицательного ожидаемого счета обычно трактуются сле-
дующим образом: чем меньше ожидаемый счет матрицы, тем менее алгоритм Смита
– Ватермана склонен расширять выравнивания в обе стороны, на слабогомологичные
участки последовательностей.
Вычисление параметра 𝑐 основано на сохранении ожидаемого счета в новой мат-
рице по сравнению с исходной:
∑︁
𝑖,𝑗
𝑝
𝑖
𝑝
𝑗
1
𝜆
ln(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
) =
∑︁
𝑖,𝑗
𝑃
𝑖
𝑃
′
𝑗
1
𝜆
ln(
𝑄
𝑖𝑗
𝑃
𝑖
𝑃
′
𝑗
).
Это уравнение имеет единственное решение относительно 𝑐, которое может быть
найдено численными методами. В случае, если все 𝑃
𝑖
и 𝑃
′
𝑗
равны 1/4 (то есть частоты
появления всех нуклеотидов в обоих последовательностях одинаковы), получаем 𝑐 = 1,
𝑄
𝑖𝑗
= 𝑞
𝑖𝑗
и 𝑆
𝑖𝑗
= 𝑠
𝑖𝑗
.
3.3.2
Сохранение относительной энтропии
Другим важным параметром, который характеризует матрицу частот (и, соот-
ветственно, матрицу замен) нуклеотидов, является ее относительная энтропия (инфор-
мационная наполненность) [9]. Ниже описан способ корректировки матрицы замен
нуклеотидов, основанный на сохранении относительной энтропии новой матрицы.
Относительная энтропия вычисляется по формуле
𝐻 =
∑︁
𝑖,𝑗
𝑞
𝑖𝑗
𝑙𝑛(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
)
36
Теперь вычисление параметра 𝑐 основано на сохранении относительной энтропии
новой матрицы по сравнению с исходной:
∑︁
𝑖,𝑗
𝑞
𝑖𝑗
𝑙𝑛(
𝑞
𝑖𝑗
𝑝
𝑖
𝑝
𝑗
) =
∑︁
𝑖,𝑗
𝑄
𝑖𝑗
𝑙𝑛(
𝑄
𝑖𝑗
𝑃
𝑖
𝑃
′
𝑗
).
Это уравнение имеет единственное решение относительно 𝑐, которое может быть
найдено численными методами. В случае, если все 𝑃
𝑖
и 𝑃
′
𝑗
равны 1/4 (то есть частоты
появления всех нуклеотидов в обоих последовательностях одинаковы), получаем 𝑐 = 1,
𝑄
𝑖𝑗
= 𝑞
𝑖𝑗
и 𝑆
𝑖𝑗
= 𝑠
𝑖𝑗
.
3.4
Учет частот нуклеотидов в программе Nhunt
Важным этапом работы большинства программ поска гомологов нуклеотидных
последовательностей является этап отбора лучших диагоналей. Для этого этапа нами
был разработан способ, который позволяет более точно высчитывать 𝑃
𝑐𝑜𝑖𝑛𝑐
— вероят-
ность случайного совпадения слова из последовательности запроса и слова из банка
данных (см. формулу 1.3 для исходной версии программы Nhunt).
Пусть 𝑃
𝑖
— частоты остатков в последовательности запроса, а 𝑃
′
𝑗
— в последо-
вательности банка данных. Вероятность 𝐻
1
совпадения двух случайных букв в одной
позиции может быть рассчитана как 𝐻
1
= ∑︀
𝑖
𝑃
𝑖
𝑃
′
𝑗
. Тогда вероятность 𝐻
𝑤+1
того, что
все 𝑤 +1 букв слова из последовательности запроса совпадут со всеми 𝑤 +1 букв слова
из последовательности банка данных, рассчитывается как 𝐻
𝑤+1
= (𝐻
1
)
𝑤+1
. Аналогич-
но вероятность 𝐻
𝑤
того, что два слова не совпадают лишь в одной букве в конкретной
позиции, рассчитывается как 𝐻
𝑤
= (𝐻
1
)
𝑤
· (1 − 𝐻
1
). В итоге 𝑃
𝑐𝑜𝑖𝑛𝑐
— вероятность того,
что два слова совпадут по крайней мере в 𝑤 из 𝑤 + 1 позиций, — может быть рассчи-
тана как 𝑃
𝑐𝑜𝑖𝑛𝑐
= 𝐻
𝑤+1
+ 𝑤 · 𝐻
𝑤
. Именно это значение применяется в ходе расчета 𝐸
𝑞
—
ожидаемого счета диагонали (см. формулу 1.4).
Все три вышеперечисленных метода корректировки матрицы замен также были
внедрены в программу поиска гомологов нуклеотидных последовательностей Nhunt,
созданную в ходе предварительной работы. Пользователь может выбрать, на каких
этапах поиска гомологов использовать корректировку матрицы, а также какой именно
способ корректировки матрицы следует применять.
За выбор этапа поиска, на котором применяется корректировка матрицы, отвеча-
ет опция −𝑎 программы Nhunt. Данная опция может принимать следующие значения:
1. −𝑎 = 0. В этом случае корректировка матрицы замен вообще не применяется. По-
иск гомологов проходит в точности так, как описано в разделе 1.1.4. Это значение
является значением по умолчанию.
37
2. −𝑎 = 1. Корректировка матрицы замен происходит на этапе локального выравни-
вания для лучших диагоналей и на этапе слияния участков локального сходства
(см. одноименные части в разделе 1.1.4). Вместо параметров 𝑚𝑎𝑡𝑐ℎ и 𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ
для каждой пары нуклеотидов используется соответствующий элемент матрицы
замен. Сама же матрица рассчитывается на основе частот встречаемости нуклео-
тидов в последовательности запроса и в последовательности банка данных.
3. −𝑎 = 2. Корректировка матрицы замен применяется только на этапе оконча-
тельного подсчета веса выравнивания. Для каждого выданного выравнивания вы-
числяется скорректированная матрица замен нуклеотидов. Эта матрица рассчи-
тывается на основе частот встречаемости нуклеотидов в тех фрагментах после-
довательностей запроса и банка данных, которые вошли в выравнивание. Затем
для каждой пары последовательностей из одного выравнивания вновь проводит-
ся локальное выравнивание с помощью алгоритма Смита – Ватермана, при этом
в качестве матрицы замен нуклеотидов применяется скорректированная матри-
ца замен. Если выравнивание изменилось, то на основании новых частот встре-
чаемости нуклеотидов рассчитывается новая скорректированная матрица замен,
проводится новое выравнивание с помощью алгоритма Смита – Ватермана и так
далее. Таким образом, выравнивание как бы “уточняется” до тех пор, пока оно не
перестанет изменяться после очередной итерации. Затем для полученных вырав-
ниваний на основании новых весов вычисляются и новые значения 𝐸𝑣𝑎𝑙𝑢𝑒.
4. −𝑎 = 3. Корректировка матрицы замен применяется на всех возможных этапах,
перечисленных в пунктах 2 и 3.
За выбор одного из трех способов корректировки матрицы замен нуклеотидов
отвечает опция −𝐴 программы Nhunt. Данная опция может принимать следующие зна-
чения:
1. −𝐴 = 1. Для корректировки матрицы замен нуклеотидов применяется метод Ю
– Альтшуля (см. раздел 3.2). Это значение является значением по умолчанию.
2. −𝐴 = 2. Для корректировки матрицы замен нуклеотидов применяется метод,
описанный в разделе 3.3.1 и основанный на сохранении ожидаемого счета матрицы
замен.
3. −𝐴 = 3. Для корректировки матрицы замен нуклеотидов применяется метод,
описанный в разделе 3.3.2 и основанный на сохранении относительной энтропии
матрицы частот.
38
Отметим, что выбор одного из значений данной опции имеет эффект лишь в том
случае, если значение опции −𝑎 не равно 0.
Новую версию программы Nhunt, в которую включены различные варианты
корректировки матрицы замен нуклеотидов, можно скачать в Интернете по адресу
http://mouse.genebee.msu.ru/∼bennigsen/nhunt.
3.5
Сравнение разных методов корректировки
матрицы замен
Для того чтобы выяснить, насколько эффективна корректировка матрицы замен
при поиске гомологов нуклеотидных последовательностей, было проведено сравнение
качества поиска гомологов между программой BLASTN и программой Nhunt с раз-
личными наборами параметров. Ниже перечислены программы и наборы параметров,
которые были взяты для выяснения эффективности корректировки матрицы замен на
различных примерах. При указании параметров упомянуты лишь те из них, которые
отличаются от параметров по умолчанию. Для каждой программы также указаны ко-
роткие идентификаторы на латинице, которые используются в графиках 3.4, 3.5, 3.6 и
3.3 и в таблице 3.2.
1. Программа BLASTN с наиболее чувствительными параметрами
Программа BLASTN со следующим набором параметров: –dust no –word_size 4
–penalty –4 –reward 5 –gapopen 10 –gapextend 6. По нашим наблюдениям, именно
такой набор параметров позволяет добиться наибольшей чувствительности для
программы BLASTN. Программа BLASTN была взята для сравнения в силу то-
го, что она обладает наибольшей чувствительностью при поиске гомологов после
программы Nhunt. Отметим, что в программе BLASTN не используется метод кор-
ректировки матрицы замен нуклеотидов. Короткий идентификатор в графиках —
blast_best, в таблице 3.2 — A.
2. Программа Nhunt без корректировки матрицы
Программа Nhunt со следующим набором параметров: –z 7 –p 7 (см. раздел 1.1.4).
Отметим, что при данном наборе параметров в программе Nhunt не используется
метод корректировки матрицы замен нуклеотидов. Тем не менее, частоты нук-
леотидов учитываются при расчете ожидаемого счета диагонали (см. раздел 3.4).
Короткий идентификатор в графиках — nhunt_0, в таблице 3.2 — B.
3. Программа Nhunt с корректировкой матрицы с помощью модифициро-
ванного метода Ю – Альтшуля
39
Программа Nhunt со следующим набором параметров: –z 7 –p 7 -a 3 -A 1 (см.
разделы 1.1.4 и 3.4). При данном наборе параметров в программе Nhunt исполь-
зуется метод корректировки матрицы замен нуклеотидов, основанный на моди-
фицированном методе Ю – Альтшуля (см. раздел 3.2). Короткий идентификатор
в графиках — nhunt_1, в таблице 3.2 — C.
4. Программа Nhunt с корректировкой матрицы с сохранением ожидаемо-
го счета
Программа Nhunt со следующим набором параметров: –z 7 –p 7 -a 3 -A 2 (см. раз-
делы 1.1.4 и 3.4). При данном наборе параметров в программе Nhunt используется
метод корректировки матрицы замен нуклеотидов, основанный на рассмотрении
матриц специального вида с сохранением ожидаемого счета (см. раздел 3.3.1).
Программа Nhunt с таким набором параметров не отражена в графиках, так как
в подавляющем большинстве случаев кривая (совокупность точек), соответствую-
щая этому набору параметров, почти не отличается от кривой, соответствующей
методу с рассмотрением матриц специального вида с сохранением относительной
энтропии (см. следующий пункт). Короткий идентификатор в таблице 3.2 — D.
5. Программа Nhunt с корректировкой матрицы с сохранением относи-
тельной энтропии
Программа Nhunt со следующим набором параметров: –z 7 –p 7 -a 3 -A 3 (см. раз-
делы 1.1.4 и 3.4). При данном наборе параметров в программе Nhunt использует-
ся метод корректировки матрицы замен нуклеотидов, основанный на рассмотре-
нии матриц специального вида с сохранением относительной энтропии (см. раздел
3.3.2). Короткий идентификатор в графиках — nhunt_3, в таблице 3.2 — E.
Все вышеперечисленные программы с соответствующими наборами параметров
были использованы для поиска гомологов последовательностей из каждой выборки за-
просов в геномах каждой выборки банков данных (см. раздел 3.1.2). Для каждой пары
“запрос – банк данных” был получен набор выравниваний (результатов поиска гомо-
логов), которые были выданы каждой из пяти программ (наборов параметров). Всего
было получено 24 таких набора выравниваний (что соответствует четырем различным
выборкам последовательностей запроса и шести различным выборкам последователь-
ностей банка данных).
Для того, чтобы в численном виде выразить качество работы каждой програм-
мы поиска гомологов, нами был применен способ, описанный ниже. Он применяется
отдельно для каждой программы (набора параметров) поиска гомологов в каждой из
24 пар “запрос – банк данных”.
40
1. Координаты найденных гомологов сравнивались с подготовленными координата-
ми гомологичных участков — аннотированных тРНК/рРНК. Для каждого най-
денного выравнивания была подсчитана длина “правильного” фрагмента (то есть
число нуклеотидов найденного программой гомолога, лежащих в пределах коор-
динат любой аннотированной тРНК/рРНК) и длина “неправильного” фрагмента
(то есть число нуклеотидов найденного программой гомолога, лежащих вне пре-
дела координат любой аннотированной тРНК/рРНК).
2. Выравнивания сортируются по величине веса, выданной программой (набором
параметров) для каждого выравнивания.
3. Для каждого выравнивания рассчитывается две величины. Первая величина —
𝑅𝑣𝑎𝑙𝑢𝑒 — равна суммарной длине “правильных” фрагментов для всех выравнива-
ний, имеющих вес больший, чем данное выравнивание. Вторая величина — 𝑊𝑣𝑎𝑙𝑢𝑒
— равна суммарной длине “неправильных” фрагментов для всех выравниваний,
имеющих вес больший, чем данное выравнивание.
4. Для каждого выравнивания рассчитывается величина чувствительности 𝑆𝑒𝑛𝑠, ко-
торая была достигнута программой поиска гомологов на всех выравниваниях, име-
ющих вес больший, чем данное выравнивание. Эта величина рассчитывается по
следующей формуле:
𝑆𝑒𝑛𝑠 =
𝑅𝑣𝑎𝑙𝑢𝑒
𝑇𝑎𝑟𝑔𝑒𝑡𝐿𝑒𝑛
,
где 𝑇𝑎𝑟𝑔𝑒𝑡𝐿𝑒𝑛 — сумарная длина всех аннотированных гомологов последователь-
ностей запроса в геномах банка данных.
Для случая поиска гомологов транспортных РНК величина 𝑇𝑎𝑟𝑔𝑒𝑡𝐿𝑒𝑛 дополни-
тельно домножалась на общее число последовательностей запроса, так как после-
довательность любой транспортной РНК из геномов банка данных потенциаль-
но может быть обнаружена как гомолог последовательности любой транспорт-
ной РНК из выборки запроса. Для случая поиска гомологов рибосомальных РНК
(выборка последовательностей запроса plasm_rRNA) такая поправка не применя-
лась. Как выяснилось, между некоторыми фрагментами рРНК, использованными
в качестве запросов, имеется существенное сходство. В результате рассчитанное
значение чувствительности в этом случае иногда превышает 1.
5. Для каждого выравнивания рассчитывается величина избирательности 𝑆𝑒𝑙, кото-
рая была достигнута программой поиска гомологов на всех выравниваниях, име-
ющих вес больший, чем данное выравнивание. Эта величина рассчитывается по
следующей формуле:
𝑆𝑒𝑙 =
𝑅𝑣𝑎𝑙𝑢𝑒
𝑅𝑣𝑎𝑙𝑢𝑒 + 𝑊𝑣𝑎𝑙𝑢𝑒
.
41
6. Для каждого выравнивания рассчитывается величина качества поиска 𝑄𝑢𝑎𝑙, кото-
рая была достигнута программой поиска гомологов на всех выравниваниях, име-
ющих вес больший, чем данное выравнивание. Эта величина рассчитывается по
следующей формуле:
𝑄𝑢𝑎𝑙 = 𝑆𝑒𝑛𝑠 · 𝑆𝑒𝑙.
7. Качество работы каждой программы (набора параметров) оценивается как макси-
мальное значение качества поиска 𝑄𝑢𝑎𝑙 по всем выравниваниям, выданным дан-
ной программой (набором параметров).
Данные по численному выражению качества работы каждой протестированной
программы (набора параметров) для каждой из 24 пар “запрос – банк данных” собраны
в таблице (4.1) Приложения. Соответствующее время работы для каждой программы
(набора параметров) приведены в таблице (4.2) Приложения.
В таблице (3.2) приведены данные о том, какая программа (набор параметров)
обладает наилучшим (самым высоким) качеством поиска для каждой из 24 пар “запрос
– банк данных”.
Таблица 3.2
Программа, показавшая самое высокое качество поиска
Банк данных / запрос ecoli_tRNA mit_tRNA plasm_tRNA plasm_rRNA
rich_bact
E
C
B
C
rich_ar
E
E
E
C
normal_bact
B
C
B
C
normal_ar
B
C
B
C
low_bact
B
E
B
E
low_ar
B
E
B
C
A — программа BLASTN с наиболее чувствительными параметрами. B — программа
Nhunt без корректировки матрицы. C — программа Nhunt с корректировкой матрицы
с помощью модифицированного метода Ю – Альтшуля. D — программа Nhunt с
корректировкой матрицы с сохранением ожидаемого счета. E — программа Nhunt с
корректировкой матрицы с сохранением относительной энтропии.
Из таблицы видно, что варианты A (BLAST с наиболее чувствительными пара-
метрами) и D (метод корректировки матрицы замен с сохранением ожидаемого счета)
ни разу не стали лучшими. Сравнение вариантов С и E показывает полезность их обоих.
Достаточно часто лучшим оказывается вариант с полным отсутсвием корректировки;
причины такого поведения пока окончательно не выяснены.
Кроме вычисления такой достаточно грубой характеристики, как описанное вы-
ше “качество”, все наборы выравниваний, полученные разными программами для каж-
42
дой пары “запрос – банк данных”, были обработаны в соответствии с процедурой, опи-
санной в разделе 3.1.1. Затем для каждой пары “запрос – банк данных” был построен
график, отражающий качество поиска для каждой программы (набора параметров).
Четыре наиболее типичных графика приведены ниже. При этом приведена только та
область графика, в которой для каждой программы достигается максимальная вели-
чина 𝑄𝑢𝑎𝑙 среди всех выравниваний, выданных данной программой.
Поиск гомологов тРНК E. coli в геномах архей
В данном примере программы поиска использовались для поиска гомологов транс-
портных РНК E. coli (выборка ecoli_tRNA) в геномах различных архей с GC-составом,
близким к нормальному (выборка normal_ar). В результате обработки полученных вы-
равниваний и координат гомологов в соответствии с процедурой, описанной в разделе
3.1.1, был получен график 3.3, отражающий качество поиска гомологов различными
программами (наборами параметров).
43
Рис. 3.3
Сравнение различных программ поиска гомологов. blast_best — результаты работы
программы BLASTN с наиболее чувствительными параметрами. nhunt_0 — резуль-
таты работы программы Nhunt без корректировки матрицы. nhunt_1 — результаты
работы программы Nhunt с корректировкой матрицы с помощью модифицированного
метода Ю – Альтшуля. nhunt_3 — результаты работы программы Nhunt с корректи-
ровкой матрицы с сохранением относительной энтропии. По оси ординат для каждой
точки (выравнивания) отложена соответствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс
— величина 𝑊𝑣𝑎𝑙𝑢𝑒.
В данном случае по неясным пока причинам результат применения набора па-
раметров, при котором не используется корректировка матрицы замен нуклеотидов,
превосходит по качеству результаты применения набора параметров, при котором ис-
пользуются различные методы корректировки матрицы замен.
Поиск гомологов тРНК митохондрии человека в GC-бедных геномах бакте-
рий
В данном примере программы поиска использовались для поиска гомологов транс-
портных РНК митохондрии человека (выборка mit_tRNA) в геномах различных бак-
44
терий с низким GC-составом (выборка low_bact). В результате обработки полученных
выравниваний и координат гомологов в соответствии с процедурой, описанной в разде-
ле 3.1.1, был получен график 3.4, отражающий качество поиска гомологов различными
программами (наборами параметров).
Рис. 3.4
Сравнение различных программ поиска гомологов. blast_best — результаты работы
программы BLASTN с наиболее чувствительными параметрами. nhunt_0 — резуль-
таты работы программы Nhunt без корректировки матрицы. nhunt_1 — результаты
работы программы Nhunt с корректировкой матрицы с помощью модифицированного
метода Ю – Альтшуля. nhunt_3 — результаты работы программы Nhunt с корректи-
ровкой матрицы с сохранением относительной энтропии. По оси ординат для каждой
точки (выравнивания) отложена соответствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс
— величина 𝑊𝑣𝑎𝑙𝑢𝑒.
Из-за малого сходства последовательностей тРНК митохондрии человека и ге-
номов архейных микроорганизмов все программы плохо разделяют ложные и правиль-
ные находки. Тем не менее, даже в этом случае применение корректировки матрицы
замен позволяет найти значительное число выравниваний гомологичных последова-
тельностей, имеющих наивысший счет среди остальных выравниваний.
45
Поиск гомологов рРНК плазмодия в GC-бедных геномах бактерий
В данном примере программы поиска использовались для поиска гомологов ри-
босомальных РНК плазмодия (выборка plasm_rRNA) в геномах различных бактерий с
низким GC-составом (выборка low_bact). Напомним, что рибосомальная РНК плазмо-
дия также имеет экстремально низкий GC-состав. В результате обработки полученных
выравниваний и координат гомологов в соответствии с процедурой, описанной в разде-
ле 3.1.1, был получен график 3.5, отражающий качество поиска гомологов различными
программами (наборами параметров).
Рис. 3.5
Сравнение различных программ поиска гомологов. blast_best — результаты работы
программы BLASTN с наиболее чувствительными параметрами. nhunt_0 — резуль-
таты работы программы Nhunt без корректировки матрицы. nhunt_1 — результаты
работы программы Nhunt с корректировкой матрицы с помощью модифицированного
метода Ю – Альтшуля. nhunt_3 — результаты работы программы Nhunt с корректи-
ровкой матрицы с сохранением относительной энтропии. По оси ординат для каждой
точки (выравнивания) отложена соответствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс
— величина 𝑊𝑣𝑎𝑙𝑢𝑒.
Здесь применение корректировки матрицы замен нуклеотидов заметно улучшает
46
качество поиска гомологов; фактически, применение любого из методов корректиров-
ки матрицы замен позволяет найти подавляющее большинство гомологов, при этом
негомологичные выравнивания практически не встречаются.
Поиск гомологов рРНК плазмодия в GC-богатых геномах архей
В данном примере программы поиска использовались для поиска гомологов ри-
босомальных РНК плазмодия (выборка plasm_rRNA) в геномах различных архей с
высоким GC-составом (выборка rich_ar). В результате обработки полученных вырав-
ниваний и координат гомологов в соответствии с процедурой, описанной в разделе 3.1.1,
был получен график 3.6, отражающий качество поиска гомологов различными програм-
мами (наборами параметров).
47
Рис. 3.6
Сравнение различных программ поиска гомологов. blast_best — результаты работы
программы BLASTN с наиболее чувствительными параметрами. nhunt_0 — резуль-
таты работы программы Nhunt без корректировки матрицы. nhunt_1 — результаты
работы программы Nhunt с корректировкой матрицы с помощью модифицированного
метода Ю – Альтшуля. nhunt_3 — результаты работы программы Nhunt с корректи-
ровкой матрицы с сохранением относительной энтропии. По оси ординат для каждой
точки (выравнивания) отложена соответствующая ей величина 𝑅𝑣𝑎𝑙𝑢𝑒, по оси абсцисс
— величина 𝑊𝑣𝑎𝑙𝑢𝑒.
Здесь применение корректировки матрицы замен нуклеотидов также заметно
улучшает качество поиска гомологов.
Резюмируя, можно предположить, что в случаях, когда последовательности за-
проса сами имеют сильно сдвинутый нуклеотидный состав, поиск их гомологов замет-
но более эффективен с применением корректировки матрицы замен нуклеотидов. При
этом неважно, насколько сдвинутый нуклеотидный состав имеют последовательности
банка данных. В том случае, когда нуклеотидный состав последовательностей запро-
са не сильно отличается от равномерного, для поиска гомологов более эффективными
могут оказаться как те алгоритмы, в которых применяется корректировка матрицы
замен, так и те, в которых она не применяется.
48
В тех примерах, когда наиболее эффективными являются методы, включающие
в себя корректировку матрицы замен нуклеотидов, примерно в половине случаев са-
мым эффективным методом корректировки является модифицированный метод Ю –
Альтшуля. В остальных случаях самым эффективным методом является метод коррек-
тировки, оснванный на сохранении относительной энтропии. При этом метод коррек-
тировки, основанный на сохранении ожидаемого счета, почти всегда очень близок по
качеству к методу, оснванному на сохранении относительной энтропии.
49
Глава
4
ЗАКЛЮЧЕНИЕ
В данной работе продемонстрирована эффективность применения методов кор-
ректировки матрицы замен нуклеотидов для решения проблемы малой сложности при
поиске гомологов нуклеотидных последовательностей.
Создан пакет компьютерных сценариев, позволяющий оценивать качество рабо-
ты различных программ поиска гомологов нуклеотидных последовательностей, исходя
из биологического смысла найденных выравниваний. Данная система успешно проте-
стирована при сравнении качества работы нескольких существующих программ поиска
гомологов нуклеотидных последовательностей. В ходе теста было выяснено, что для
наиболее популярной на сегодняшний день программы поиска гомологов BLASTN па-
раметры, выставленные по умолчанию, а также доступные на сайте веб-интерфейса
к данной программе, в большинстве случаев не являются оптимальными для поиска
гомологов нуклеотидных последовательностей.
Было проведено исследование существующих подходов к решению проблемы ма-
лой сложности как для аминокислотных, так и для нуклеотидных последовательностей.
Оказалось, что для нуклеотидных последовательностей применяется лишь метод мас-
кировки участков малой сложности, тогда как для аминокислотных последовательно-
стей успешно применяется метод корректировки матрицы замен, который лишен зна-
чительных недостатков, присущих методу маскировки. Тем не менее, существующие
способы корректировки матрицы замен аминокислотных остатков формально полно-
стью применимы к матрице замен нуклеотидов.
Исходя из этого, была предложена модификация наиболее простого и эффек-
тивного метода корректировки матрицы замен, применяемого в программе BLASTP.
Также было разработано два оригинальных метода корректировки матрицы замен нук-
леотидов. Все три метода были встроены в программу поиска гомологов нуклеотидных
последовательностей Nhunt, созданную в ходе предварительной работы.
Для тестирования методов корректировки матрицы замен нуклеотидов была ис-
50
пользована созданная нами система проверки качества поиска гомологов нуклеотидных
последовательностей. На ряде примеров показано, что использование корректировки
матрицы замен чаще всего заметно улучшает качество работы программы поиска го-
мологов нуклеотидных последовательностей без заметного увеличения времени работы
программы. Частные случаи, в которых применение корректировки матрицы замен
ухудшало качество поиска, требуют более детального рассмотрения.
Полученные результаты говорят о том, что методы корректировки матрицы за-
мен, которые сейчас широко используются лишь для аминокислотных последователь-
ностях, могут успешно применяться и для поиска гомологов нуклеотидных последо-
вательностей. Это позволит значительно улучшить качество поиска гомологов в тех
случаях, когда последовательности запроса и / или банка данных содержат протяжен-
ные участки малой сложности.
51
ВЫВОДЫ
1. Разработана методика оценки и сравнения качества работы различных программ
поиска гомологов нуклеотидных последовательностей. Большая часть компонен-
тов данной методики реализована в виде программного пакета.
2. Метод корректировки матрицы замен аминокислотных остатков, лежащий в ос-
нове борьбы с малой сложностью в программе BLASTP, адаптирован для случая
нуклеотидных последовательностей, при этом лежащий в его основе алгоритм мо-
дицифирован и ускорен.
3. Предложены два оригинальных метода корректировки матрицы замен для случая
нуклеотидных последовательностей, основанные на использовании матриц замен
нуклеотидов специального вида.
4. Все три вышеописанных метода корректировки матрицы замен были встроены в
Nhunt — программу поиска гомологов нуклеотидных последовательностей.
5. Применение методов корректировки матрицы замен для решения проблемы малой
сложности в нуклеотидных последовательностях показало свою эффективность и
позволило заметно улучшить качество поиска гомологов на большей части иссле-
дованных примеров. Вместе с тем, в некоторых случаях применение корректиров-
ки матрицы замен по неясным пока причинам незначительно ухудшает качество
поиска гомологов.
52
Литература
1. Karlin, S. and Altschul, S. F. 1990. Methods for assessing the statistical significance
of molecular sequence features by using general scoring schemes. Proceedings of the
National Academy of Sciences of the USA 87:2264-2268.
2. Noe, L. and Kucherov, G. 2004. Improved hit criteria for DNA local alignment. BMC
Bioinformatics, 5:149.
3. Noe, L. and Kucherov, G. 2005. YASS: enhancing the sensitivity of DNA similarity
search. Nucleic Acids Research, 33 (Web Server issue):W540–W543
4. Pearson, W. R. and Lipman, D. J. 1988. Improved tools for biological sequence
comparison. Proceedings of the National Academy of Sciences of the USA 4:2444-2448.
5. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. and Lipman, D. J. 1990. Basic local
alignment search tool. Journal of Molecular Biology 215:403-410.
6. Autschul, S. F., Madden, T. L., Schaffer, A. A., Zhang, J., Zhang, Z., Miller, W. and
Lipman, D. J. 1997. Gapped BLAST and PSI-BLAST: a new generation of protein
database search programs. Nucleic Acids Research 25:3389-3402.
7. Smith, T. F. and Waterman, M. S. 1981. Identification of common molecular
subsequences. Journal of Molecular Biology 147:195-197.
8. Yu Y.-K., Wootton J. C. and Altschul S. F. 2003. The compositional adjustment of
amino acid substitution matrices. Proceedings of the National Academy of Sciences of
the United States of America 87:15688-15693.
9. Yu Y.-K. and Altschul S. F. 2005. The construction of amino acid substitution matrices
for the comparison of proteins with non-standard compositions. Bioinformatics 21:902-
911
10. Altschul S. F., Wootton J.C., Gertz E.M., Agarwala R., Morgulis A., Schaffer A.A.,
Yu Y.-K. 2005. Protein database searches using compositionally adjusted substitution
matrices. FEBS J. 272:5101-5109
53
11. Yu Y.-K., Gertz E.M., Agarwala R., Schaffer A.A., Altschul S.F. 2006. Retrieval
accuracy, statistical significance and compositional similarity in protein sequence
database searches. Nucleic Acids Res. 34:5966-5973
12. Morgulis A., Gertz E.M., Schaffer A.A., Agarwala R. 2006. A fast and symmetric DUST
implementation to mask low-complexity DNA sequences. J Comput Biol. 2006 13:1028-
1040.
13. http://faculty.virginia.edu/wrpearson/fasta/
14. ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/user_manual.pdf
15. http://blast.ncbi.nlm.nih.gov/Blast.cgi
16. http://www.ncbi.nlm.nih.gov/blast/discontiguous.shtml
17. http://mouse.belozersky.msu.ru/∼bennigsen/nhunt
54
ПРИЛОЖЕНИЯ
55
Приложение 1
Таблица 4.1
Качество работы различных программ поиска гомологов
Запрос
Банк данных blastn_b
nhunt_0
nhunt_1
nhunt_2
nhunt_3
ecoli_tRNA
rich_bact
0.418
0.471
0.475
0.467
0.483
ecoli_tRNA
rich_ar
0.147
0.161
0.172
0.170
0.177
ecoli_tRNA
normal_bact
0.451
0.513
0.487
0.485
0.494
ecoli_tRNA
normal_ar
0.195
0.220
0.188
0.176
0.182
ecoli_tRNA
low_bact
0.427
0.507
0.401
0.318
0.341
ecoli_tRNA
low_ar
0.193
0.259
0.166
0.116
0.131
mit_tRNA
rich_bact
1.99 · 10
−3
1.96 · 10
−3
2.46 · 10
−3
1.48 · 10
−3
1.78 · 10
−3
mit_tRNA
rich_ar
6.19 · 10
−4
9.24 · 10
−4
1.07 · 10
−3
1.66 · 10
−3
1.82 · 10
−3
mit_tRNA
normal_bact 5.87 · 10
−4
5.47 · 10
−4
2.40 · 10
−3
1.38 · 10
−3
1.62 · 10
−3
mit_tRNA
normal_ar
2.71 · 10
−4
3.19 · 10
−4
6.55 · 10
−4
6.39 · 10
−4
6.39 · 10
−4
mit_tRNA
low_bact
4.14 · 10
−4
4.14 · 10
−4
2.06 · 10
−3
3.28 · 10
−3
3.56 · 10
−3
mit_tRNA
low_ar
3.45 · 10
−4
1.59 · 10
−4
2.83 · 10
−3
3.03 · 10
−3
3.35 · 10
−3
plasm_tRNA
rich_bact
0.224
0.288
0.272
0.244
0.256
plasm_tRNA
rich_ar
0.124
0.151
0.147
0.151
0.158
plasm_tRNA normal_bact
0.196
0.294
0.267
0.249
0.257
plasm_tRNA
normal_ar
0.167
0.237
0.206
0.173
0.180
plasm_tRNA
low_bact
0.243
0.323
0.239
0.189
0.203
plasm_tRNA
low_ar
0.174
0.257
0.165
0.100
0.116
plasm_rRNA
rich_bact
1.009
0.954
1.423
1.292
1.292
plasm_rRNA
rich_ar
0.449
0.406
0.823
0.782
0.786
plasm_rRNA normal_bact
0.700
0.655
1.200
1.123
1.125
plasm_rRNA
normal_ar
0.279
0.262
0.742
0.657
0.679
plasm_rRNA
low_bact
0.539
0.466
1.253
1.254
1.259
plasm_rRNA
low_ar
0.194
0.164
0.714
0.699
0.700
blast_best — качество работы программы BLASTN с наиболее чувствительными
параметрами. nhunt_0 — качество работы программы Nhunt без корректировки
матрицы. nhunt_1 — качество работы программы Nhunt с корректировкой матрицы
с помощью модифицированного метода Ю – Альтшуля. nhunt_2 — качество работы
программы Nhunt с корректировкой матрицы с сохранением ожидаемого счета.
nhunt_3 — качество работы программы Nhunt с корректировкой матрицы с
сохранением относительной энтропии.
56
Приложение 2
Таблица 4.2
Время работы различных программ поиска гомологов
Запрос
Банк данных blastn_b nhunt_0 nhunt_1 nhunt_2 nhunt_3
ecoli_tRNA
rich_bact
196
889
835
1232
1170
ecoli_tRNA
rich_ar
92
507
455
620
582
ecoli_tRNA
normal_bact
111
892
987
1501
1362
ecoli_tRNA
normal_ar
67
626
711
1075
1031
ecoli_tRNA
low_bact
28
217
349
458
430
ecoli_tRNA
low_ar
55
423
603
888
782
mit_tRNA
rich_bact
40
202
294
457
406
mit_tRNA
rich_ar
18
91
134
196
176
mit_tRNA
normal_bact
42
268
292
458
413
mit_tRNA
normal_ar
30
206
212
323
294
mit_tRNA
low_bact
24
183
118
115
103
mit_tRNA
low_ar
45
370
214
259
231
plasm_tRNA
rich_bact
142
650
676
1019
905
plasm_tRNA
rich_ar
64
359
365
507
460
plasm_tRNA normal_bact
91
513
607
977
894
plasm_tRNA
normal_ar
61
369
439
701
648
plasm_tRNA
low_bact
29
140
213
287
259
plasm_tRNA
low_ar
53
268
381
576
521
plasm_rRNA
rich_bact
65
346
439
541
512
plasm_rRNA
rich_ar
28
226
199
267
251
plasm_rRNA normal_bact
111
508
580
638
646
plasm_rRNA
normal_ar
79
470
423
473
502
plasm_rRNA
low_bact
472
207
147
112
108
plasm_rRNA
low_ar
837
2311
646
437
510
blast_best — время работы программы BLASTN с наиболее чувствительными
параметрами. nhunt_0 — время работы программы Nhunt без корректировки
матрицы. nhunt_1 — время работы программы Nhunt с корректировкой матрицы с
помощью модифицированного метода Ю – Альтшуля. nhunt_2 — время работы
программы Nhunt с корректировкой матрицы с сохранением ожидаемого счета.
nhunt_3 — время работы программы Nhunt с корректировкой матрицы с
сохранением относительной энтропии.
57
Информация о работе Оптимизация поиска гомологов нуклеотидных последовательностей