Предсказание вторичной структуры молекул РНК

Вопрос предсказания вторичной структуры молекул РНК оформился в 1960—1970 гг. после появления работы Фреско и др^[1]., в которой было показано, что одноцепочечная молекула РНК способна «замыкаться» сама на себя с формированием водородных связей между азотистыми основаниями.

В 1980;х гг. получили развитие первые алгоритмы предсказания вторичной структуры, основанные на методиках динамического программирования. Данные алгоритмы производили поиск оптимальных, а чуть позже — и субоптимальных структур на основании определения наиболее энергетически выгодного паттерна спаривания нуклеотидов.

Приведем пример такого алгоритма — алгоритм Зукера, доступный в программе mfold (http://mfold.bioinfo.rpi.edu/cgi-bin/rna-forml.cgi). Данный алгоритм использует стандартное определение вторичной структуры РНК (см. подпараграф 1.2.6), т. е. вторичная структура РНК определена как список нуклеотидных пар такой, что для пары с индексами г, j, где i < j, и другой пары k, I (г < к) выполняются следующие условия:

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

• 1) заполнения квадратной матрицы динамического программирования;
• 2) обратного прохода по матрице в направлении оптимальных значений.

В ходе выполнения первого этапа алгоритма производится заполнение квадратной матрицы, строки и столбцы которой соответствуют нуклеотидам анализируемой последовательности, т. е. каждая ячейка такой матрицы (г, j) соответствует наличию взаимодействия между г-м и j-м нуклеотидами. Алгоритм последовательно, по диагоналям, начиная от центральной, строит все более протяженные элементы вторичной структуры, обладающие минимальной энергией. Информация об энергии обнаруженных фрагментов структуры сохраняется в нескольких элементах массива, описывающего вершины графа. Используется несколько типов значений: V (i, j), W (i, j), FL (i, j), FH (i, j).

V (i, j) — содержит значения минимальной энергии всех структур, в которых нуклеотиды 5, и 5_; формируют пару (для нуклеотидов, не способных формировать пару, V (i, j) ставится равным бесконечности). Данная энергия рассчитывается на основании значений, рассчитанных для внутренних по отношению к паре (i, j) структур (3.29). Возможные варианты «внутренних» структур приведены на рис. 3.27, а. Вычисления данного значения, выполняемые на каждом шагу первого этапа алгоритма, описываются следующей формулой:

Рис. 3.26. Элементы вторичной структуры РНК, моделируемые алгоритмом Зукера при расчете промежуточных значений.

W (i, j) — содержит значения минимальной энергии всех возможных типов структур, вне зависимости от непосредственного спаривания нуклеотидов 5, и S: и вычисляется для каждой пары согласно следующей формуле:

FH (iJ) — характеризует энергию петли шпилечной структуры, замыкаемой остатками 5, — и Sj, и рассчитывается на основании табличных значений, известных из экспериментов по «плавлению» РНК.

FL (i, j, k, /) — характеризует энергию внутреннего участка структуры, располагающегося между парами S_h Sj и S_k, 5/. Это может быть:

• • приращение стеблевой структуры (стабилизированное стекинговыми взаимодействиями);
• • одностороннее выпячивание;
• • внутренняя петля.

Так же как и FH (i, j), данная величина рассчитывается на основании табличных значений.

Выбор, сделанный при вычислении формул (3.30) и (3.31), запоминается (как в случае алгоритма выравнивания последовательностей). Конформация с наименьшей энергией находится путем обратной трассировки (второй этап алгоритма) из ячейки матрицы W (l, п) по значениям W и V, заполненным в результате первого прохода. Вычислительная сложность данного алгоритма оценивается как 0(п^А), однако если ограничить длину анализируемых FL участков, то возможно снижение вычислительной сложности до 0(п³).

По мере исследования реальных структур молекул РНК было выяснено, что структура с минимальной энергией, получаемая в результате работы описанного выше алгоритма, далеко не всегда соответствует биологически активной форме. Другими словами, реальные биологические структуры молекул РНК являются субоптимальными с точки зрения минимальной энергии вторичной структуры. Это потребовало разработки алгоритмов, позволяющих исследовать ансамбли структур с энергией, несколько превышающей минимальную.

Для понимания логики работы алгоритма анализа субоптимальных структур следует обратиться к вопросу предсказания структуры кольцевой (замкнутой) РНК. Напомним, что алгоритм поиска оптимальной структуры линейной молекулы постепенно рассматривает все более и более протяженные участки молекулы, заканчивая заполнение матрицы расчетом энергии W ( 1, л), включающей всю молекулу. В случае кольцевой молекулы РНК при использовании подобного алгоритма выбор любых двух нуклеотидов делит последовательность на «внутренний» и «внешний» фрагменты, а сумма значений V (i, j) и V (j, г) для всех комбинаций i и j дает нам энергии всех вариантов пространственных структур, формируемых молекулой. Применение данной логики к линейным молекулам РНК возможно, если обрабатывать «внешний» фрагмент (включающий переход от S_n к $) как особый вариант «разорванной» петли. В итоге после заполнения матрицы алгоритмом динамического программирования и обнаружения структуры с минимальной энергией (Е_тin) становится возможен анализ субоптимальных структур. Пусть Р будет числом от 0 до 100, тогда Р-оптимальной парой нуклеотидов S_j9 Sj будем называть те нуклеотиды, для которых выполняется следующее неравенство:

Р-оптимальные пары можно либо отложить на точечном графике (рис. 3.27), что позволит оценить, насколько детерминированными являются различные фрагменты пространственной структуры, либо проведя обратную трассировку от интересующей пары — выделить содержащую ее субоптимальную структуру. Данный алгоритм не гарантирует нахождения всех субоптимальных структур, но позволяет выделить все возможные фрагменты таких структур.

Рис. 3.27. Точечные графики структуры с минимальной энергией (а),.

5% (б) и 10% (в) Р-оптимальных структур последовательности самосплайсируемого нитрона одной из рибосомальных РНК Tetrahymena thermophila

Дальнейшее развитие методов моделирования вторичной структуры молекул РНК заключалось в анализе вероятности встретить пару нуклеотидов Sj, Sj в ансамбле всех возможных структур. Расчет производится посредством определения нормализующего фактора Z, основанного на формуле Больцмана:

где Q — множество всех возможных вторичных структур молекулы РНК; S — вариант вторичной структуры молекулы; E (S) — энергия структуры S; R — константа Больцмана; Т — температура.

Вероятность конкретного варианта 5₀ структуры в таком случае будет равна.

Отметим, что нормализующую константу Z можно эффективно рассчитать с помощью модифицированной версии приведенного выше алгоритма поиска структуры с минимальной энергией:

• • заменив в формулах (3.30) и (3.31) расчет энергий компонентов структуры (/77, FL, E_pnin E_unpair) на их больцмановские вероятности;
• • используя вместо знаков сложения энергии — умножение вероятностей (комбинирование независимых фрагментов структуры), а вместо выбора минимального значения энергии (min) — суммирование вероятностей;
• • добавив ряд ограничений в рекурсию, чтобы гарантировать однократное вхождение каждого варианта фрагмента структуры в вычисления.

После расчета нормализующего фактора появляется возможность определить вероятность р_{} встретить пару S_v Sj в ансамбле всех возможных структур как суммарную вероятность встретить эту пару в каждой конфигурации молекулы, отнесенную к нормализующему фактору.

Знание индивидуальных вероятностей р_{; для каждой пары позволило создать вероятностные схемы формирования субоптимальных структур (обратный обход матрицы с Больцмановской вероятностью выбора состояний в противовес детерминистическому «обратному» проходу по матрице алгоритма, когда с каждым шагом строго выбирается оптимальный вариант), формирующих термодинамический ансамбль вторичных структур.

Дальнейшее развитие предсказания пространственных структур РНК требует учета энергетического вклада от формирования элементов третичной структуры алгоритмами динамического программирования с более высокой сложностью (0(п⁶)) либо построением выравниваний молекул РНК с известной структурой (сравнительное моделирование).

• [1] Fresco J. R., Alberts В. М., Doty Р. Some molecular details of the secondary structure ofribonucleic acid //Nature. 1960. Oct 8; 188:98—101.