Статистическая оценка значимости значений оценочной функции при выравнивании последовательностей

Когда в 1980;х гг. появились алгоритмы выравнивания последовательностей, возникла необходимость оценивать статистическую достоверность получаемых результатов. Первые методы оценки достоверности результатов исходили из предпосылки, что значения оценочной функции при выравнивании случайных последовательностей распределены согласно нормальному закону. Обычно процедура такой оценки включала:

• 1) получение случайных последовательностей путем перестановки местами остатков внутри выравниваемых последовательностей и повторное выравнивание получившихся последовательностей;
• 2) накопление статистики результатов выравнивания посредством множественного повторения пункта 1. Обычно повторность составляла несколько сотен итераций;
• 3) оценка среднего выборки и стандартного отклонения на основе полученных результатов;
• 4) расчет Z-индекса (Z-score) по формуле

где р и, а — среднее выборки и стандартное отклонение для распределения случайных значений оценочной функции; s (a, b) — значение оценочной функции для выравнивания реальных последовательностей А и В.

Фактически Z-индекс показывает, как далеко, если считать в стандартных отклонениях, отстоит значение оценочной функции для реальных последовательностей от среднего значения для выравниваний случайных п ос лед о вате л ь н осте й.

Хотя последующие исследования показали, что для локальных выравниваний предположение о нормальном характере распределения некорректно и приводит к переоценке степени достоверности результатов (рис. 3.9), подходы, оперирующие Z-индексами, развиваются и в наши дни. Так, теоретически было доказано, что при значениях Z-индекса больше 8 выравнивание можно считать достоверным вне зависимости от закона распределения. Однако использование заведомо высоких ограничений для статистической оценки результатов делает метод менее чувствительным относительно аналогов, апеллирующих к реально наблюдаемому закону распределения. Таким образом, применение Z-иидекса оправдано только в тех случаях, когда применение других методов не представляется возможным: глобальные выравнивания, некоторые алгоритмы выравнивания пространственных структур и тредиига, оценка значимости заведомо не случайных выравниваний, например, при выравнивании гомологичных генов или белков.

Рис. 3.9. Сравнение значений функций плотности распределения экстремальных значений (а) и нормального распределения (б):

видно, что вероятность встретить выравнивание со значением оценочной функции меньше X (площадь графика слева от вертикальной черты) будет больше при использовании нормального распределения. Поскольку реально статистика значений оценочной функции описывается распределением экстремальных значений, это приводит к завышению оценки значимости выравнивания Теоретические основы расчета статистической значимости для локальных выравниваний последовательностей были разработаны в конце 1980;х — начале 1990;х гг. на примере теории анализа фрагментов одиночных последовательностей.

Итак, рассмотрим простой случай одиночной последовательности. Допустим, у нас имеется последовательность, составленная из п типов символов, которые встречаются с вероятностями р_{, р₂, р_п (для аминокислотных последовательностей п = 20). Также допустим, что каждому типу символов соответствует определенное значение параметров оценочной функции s_v s₂, •••> V Как ^{и в} случае требований к оценочной функции для локальных выравниваний, потребуем, чтобы ожидаемое значение суммы параметров оценочной функции для случайной последовательности было меньше нуля:

но при этом хотя бы одно из значений s, было больше нуля¹. В случае выполнения перечисленных выше требований можно утверждать, что существует единственный положительный корень^[1][2] X. В уравнении.

X в данном случае можно рассматривать как натуральный шкалирующий фактор для используемой оценочной функции.

Рассмотрим необычные фрагменты последовательностей, для которых значение оценочной функции больше нуля^[3]. Очевидно, что вероятность получить значение оценочной функции S для фрагментов больше некоторой величины x (P (S > х)) зависит как от параметров функции Sp s₂, s_JV так и от длины исследуемой последовательности. Теория показывает, что для сравнения между собой фрагментов, полученных с использованием различных параметров оценочной функции и из последовательностей различной длины и состава, значения оценочной функции для этих фрагментов необходимо нормализовать по формуле.

где У — нормализованное значение оценочной функции для фрагмента; X — введенный в формуле (3.7) шкалирующий фактор; К — фактор, характеризующий параметры оценочной функции и частоты встречаемости различных типов символов; N — длина исходной последовательности.

Поскольку ожидаемое значение суммы параметров оценочной функции для случайной последовательности меньше нуля, то вероятность обнаружить в некоторой последовательности определенное количество фрагментов со значением оценочной функции больше некоторого положительного значения распределена по закону Пуассона (закон распределения редких событий):

где к — количество фрагментов, для которых значение оценочной функции больше некоторой положительной величины х р — математическое ожидание для распределения Пуассона, в нашем случае.

Далее мы можем определить вероятность того, что ни один из фрагментов не превысит заданный порог отсечения:

Следовательно, вероятность найти хотя бы один сегмент в случайной последовательности, для которого оценочная функция достигнет или превысит пороговое значение х, хорошо аппроксимируется формулой.

Для упомянутых выше нормализованных значений У формула вероятности упрощается до.

В случае локального выравнивания без пропусков пар последовательностей с длинами пит описанные выше формулы изменяются с точностью замены длины последовательности Nwa произведение длин последовательностей т • п. Такая замена связана с тем, что поиск фрагментов происходит на двумерной матрице размерностью т на п.

Методы оценки параметров распределения значений оценочной функции. Теория статистического анализа распределения значений оценочной функции для фрагментов последовательностей оперирует двумя параметрами X и К, определяющими форму распределения. Для задачи анализа фрагментов одиночной последовательности и выравнивания последовательностей без пропусков изложенная выше теория дает формулы для расчета этих параметров, например, так это делается при статистической оценке результатов алгоритма BLAST. Однако для описания локального выравнивания последовательностей с пропусками и более сложных задач на данный момент никакой теории нет, хотя эмпирически было замечено, что при наличии определенной зависимости между распределением значений оценочной функции и длинами выравниваемых последовательностей это распределение также следует закону EVD.

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

Поскольку для локальных выравниваний с пропусками отсутствует возможность рассчитать значения основных параметров распределения значений оценочной функции, эти параметры оцениваются на основе результатов симуляции выравниваний для случайных последовательностей. Случайные последовательности для построения выравниваний получают либо методом перестановки символов, так же как и в случае расчета Z-индексов, либо случайной генерацией на основании частот символов, встречаемых в реальных последовательностях.

В случае обычных выравниваний такой расчет производится каждый раз при изменении системы [оценочная функция с конкретным набором параметров] — |последовательности с определенным соотношением частот встречаемости символов], поскольку любое изменение параметров оценочной функции или использование для выравнивания последовательностей с существенно отличающимися от «модельных» частотами символов требует переоценки характеристик распределения. Поскольку обычные выравнивания последовательностей не требуют переоценки параметров распределения для каждого выравнивания, этот расчет можно произвести простым методом максимального правдоподобия (maximum likelihood estimation) на основании кумулятивной функции распределения (3.11) или функции плотности распределения:

Метод максимального правдоподобия, предложенный Р. Фишером в 1912 г., основан на исследовании вероятности получения выборки наблюдений случайной величины X (x_v х₂, …, х_п) с известным законом распределения /(х, Т). В данном случае требуется получить оценку 9 параметра распределения Т, в общем случае векторного. Согласно методу максимального правдоподобия вероятность получить выборку значений (xj, х₂,…, х_п) равна f (x_v Т) • /(х₂, Т) • …/(х,₇, Т) dx₁ dx₂… dx_n. Совместная плотность вероятности

рассматриваемая как функция параметра Т, и называется функцией правдоподобия. В качестве оценки 0 параметра Т следует взять то значение, которое обращает функцию правдоподобия в максимум. Для нахождения оценки необходимо заменить в функции правдоподобия Г на 0 и решить уравнение д L/d 0 = 0. В целях упрощения вычислений переходят от функции правдоподобия к ее логарифму In L. Такое преобразование допустимо, так как функция правдоподобия — положительная функция и она достигает максимума в той же точке, что и ее логарифм. Если параметр распределения — векторная величина 0 = (0_t, 0₂, ОД то оценки максимального правдоподобия находят из системы уравнений.

В нашем случае (при использовании закона распределения (12)) необходимо решить уравнения.

где Хи К — параметры распределения, подлежащие оценке; п — размер выборки для оценки параметров; N — размер области поиска (произведение длин выравниваемых последовательностей), при дальнейшей оптимизации системы уравнений этот параметр сокращается.

После необходимых преобразований системы уравнений (3.17) параметр К можно выразить через X, а само значение X легко находится численными методами¹. Исследования показали, что 10 000 итераций симуляции (это фактически размер выборки для оценки параметров) достаточно, чтобы возможная ошибка для параметра X составляла около 1%^[4][5].

Рис. 3.10. Иллюстрация к «методу островков» для оценки параметров распределения значений оценочной функции:

а — изображен граф, на котором выделены все возможные выравнивания. Окружности обозначают вершины, являющиеся якорями для выравниваний;

6 — трехмерное изображение поверхности, формируемой значениями оценочной функции. «Метод островков» назван так, поскольку локальные выравнивания формируют «островки» положительных значений в «море» отрицательных значений оценочной функции Более сложным случаем является оценка статистической достоверности при выравниваниях типа профиль — последовательность и профиль — профиль, а также при применении трединга (см. подпараграф 3.4.2). Отличие состоит в том, что в этих случаях параметры оценочной функции изменяются для каждого конкретного выравнивания, так как сам профиль или композиция структурного остова в сочетании с энергетической функцией и являются параметрами оценочной функции. Поскольку проведение для каждого реального выравнивания тысячи дополнительных выравниваний «случайных» последовательностей делает процесс прямой оценки вычислительно крайне невыгодным^[6], был найден способ более эффективной оценки параметров распределения.

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

Также можно отметить, что эти субоптимальные выравнивания не пересекаются с основным и друг с другом, таким образом, представляя независимые, а точнее — слабо зависимые события. Поскольку островки не являются строго независимыми, оценка параметров распределения методом максимального правдоподобия производится на основании формулы для распределения Пуассона, а точнее — формулы математического ожидания распределения (3.10), показывающего ожидаемое количество «островков», для которых значение оценочной функции больше некоторого критического значения. В этом случае параметры распределения рассчитываются по формулам:

где с — минимальное значение оценочной функции для отбора островков (см. пояснения ниже); R_c — количество островков, имеющих значение оценочной функции больше с; Л — область поиска, т. е. если было проанализировано В выравниваний последовательностей с длинами п и т, то А = Впт:

где 5(f) — значение оценочной функции для i-го островка, а суммирование идет по всем элементам множества 1_СУ включающего все островки, для которых значение оценочной функции 5 больше некоторого значения с.

Отслеживание всех островков — довольно простая задача, требующая минимальной модификации алгоритма выравнивания последовательностей. Для этого на каждой вершине, в которой мы начинаем новое выравнивание, мы должны делать особую отметку. Подобные вершины называют якорями выравниваний. В дальнейшем для каждого шага алгоритма нам нужно передавать отметку якоря данного выравнивания «следующим» вершинам. Алгоритм выравнивания, реализующий эту логику, имеет такую же вычислительную сложность (0(/г²)), как и обычный алгоритм, но при этом каждое выравнивание дает нам более одного значения для последующей оценки параметров распределения.

Исследования показали, что для более точной оценки параметров этим способом необходимо производить некоторую фильтрацию субоптимальных выравниваний, отбирая только те из них, для которых значение оценочной функции больше определенного значения с (см. формулы (3.18) и (3.19)). Это происходит потому, что островки с низкими значениями оценочной функции практически не содержат пропусков, в результате чего оцениваемые значения параметров распределения смещаются к характерным для выравнивания без пропусков. Еще один эффект, требующий внимания^[7], называется «граничный эффект». Он связан с тем, что выравниваемые последовательности не бесконечны, в результате чего оценочная функция для пограничных островков не достигает максимальных значений, что вносит ошибку в оценку параметров. Для борьбы с этим эффектом используются лишь островки, якорь которых находится в «центральной части» выравнивания (рис. 3.11).

Рис. 3.11. Коррекция эффекта границы введением дополнительных областей

вокруг выравнивания

• [1] Для параметров оценочной функции, основанных на логарифмах отношений частот, таких как матрицы BLOSSUM, оба требования выполняются автоматически. Все остальныетипы параметров можно изменить, чтобы они удовлетворяли этому условию.
• [2] Если рассматривать все возможные корни, то указанное уравнение всегда имеет ещеодно решение X = 0.
• [3] В качестве примера задач, соответствующих данной модели, можно привести поискгидрофобных или заряженных кластеров в аминокислотной последовательности. В этом случае параметрами оценочной функции .v, s2, s2o будут индексы гидрофобности аминокислотили их заряд соответственно. Преимущество данного подхода относительно рассмотренныхв предыдущей главе методов, усредняющих значения в пределах некоторого «окна», является отсутствие необходимости вводить само искусственное понятие «окна» и выбирать егоразмер.
• [4] Отличие численных методов решения уравнений от аналитических легко пояснитьна примере решения определенного интеграла, известного из школьной программы. В случаеаналитического метода для подынтегрального выражения находится первообразная, послечего в формулу подставляются значения для границ интервала, на котором производитсяинтегрирование. Численные методы используются при отсутствии возможности нахожденияпервообразной, в этом случае график подынтегрального выражения разбивается на фрагменты методом прямоугольников, и рассчитывается их суммарная площадь. Контролируястепень разбиения (величину шага аргумента), достигается требуемая точность решения.
• [5] Поскольку параметр X входит в формулу распределения EVD как показатель степениэкспоненты, то даже отклонение в 4% дает серьезную ошибку в 2,7 раза при оценке статистической достоверности выравниваний.
• [6] Следует напомнить, что сам по себе поиск в базе данных может потребовать проведениядесятков тысяч реальных выравниваний.
• [7] Этот эффект проявляется и при других способах оценки параметров.