RNAFoldingAI

Внимание: Открытый исследовательский проект, к которому может присоединится каждый ...

RNAFoldingAI - программное обеспечение молекулярного моделирования РНК, с элементами искусственного интеллекта. Разрабатывается участниками Викиверситета для проведения исследований в трех направлениях: геномика, искусственный интеллект, гибкое программное обеспечение высокого качества.

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

Ревизии

Программное обеспечение находится в процессе разработки поэтому его описание будет появляться практически в интерактивном режиме.

Исходники данного ПО находятся в репозитории SVN, т.е. используется контроль версий. Для ознакомления список ревизий можно посмотреть здесь Далее следует список ревизий относящихся к данному проекту:

1. №22 - Первая черновая версия проекта RNAFoldingAI. Создает .pdb файл стартового состояния по первичной структуре РНК.
2. №23 - Черновая версия модуля RNARotate „Вращение 3-х мерной цепи РНК“
3. №24 - Черновая версия модуля RNAEnergyСalc „Расчет энергетических характеристик РНК“
4. №25 - Черновая версия модуля FoldingSearchRnd „Поиск структуры с минимальной оценкой энергии цепи случайным перебором“.

• Итого на данный момент выполнен реинжиниринг РНК расчетов из ПО Rosseta. Функциональность воспроизведена, но требует тестирования. Можно начинать думать об архитектуре и написании модуля ИИ. Но уже сейчас можно отметить улучшение ряда архитектурных решений: ПО разделено на предметные модули, убраны лишние ненужные связи между классами.

Причины создания

Предметные причины создания описываются в статье «Геномика бросает вызов искусственному интеллекту». Здесь же мы рассмотрим причины только с точки зрения создания нового программного обеспечения RNAFoldingAI 0.1.

При написании данного ПО, был исследован исходных код ПО Rosetta 2.3 (см. Молекулярное моделирование), который написан преимущественно на Си++ с элементами Фортрана и при этом ориентирован на ОС Linux, и дальнейшее его развитие пошло исключительно для ОС Linux. При ближайшем рассмотрении оказалось:

1. Архитектурное качество ПО Rosetta 2.3 находится на очень низком уровне
2. Смешан объектно-ориентированный и процедурный стиль, часто используются глобальные переменные
3. Использование Си++ утяжеляет программные конструкции
4. Использование Фортрана сильно затрудняет процесс отладки (debug`ing)
5. Объем кода большой и он не может быть поделен на предметные части, все смысловые части сильно переплетены между собой

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

Планы

Поэтому было решено полностью переписать данное ПО, а точнее создать заново на современном языке C#. Но при этом как в плане изучения, так и тестирования ориентироваться на функциональность ПО Rosetta 2.3. При этом, первейшей задачей является создание пяти независимых модулей:

1. Формирование 3-х мерной цепи РНК по ее первичной структуре Шаблон:Сделано
2. Вращение 3-х мерной цепи РНК Шаблон:Сделано
3. Расчет энергетических характеристик РНК Шаблон:Сделано
4. Поиск структуры с минимальной оценкой энергии цепи случайным перебором Шаблон:Сделано
5. Поиск структуры с минимальной оценкой энергии цепи применяя методы ИИ

В дальнейшем они могут быть расширены и для экспериментов с ДНК и белками.

Описание функциональности

Модуль RNACreate „Формирование 3-х мерной цепи РНК по ее первичной структуре“

Входные, выходные данные

Считывает файл с расширением .seq, содержащий первичную структуру РНК.

Пример входного файла .seq  

>test.pdb
gggcgcaagccu

На выходе формирует два файла: промежуточный выходной файл с расширением .out и файл 3-х мерного расположения молекулы с расширением .pdb (стандартный файл, который можно просмотреть программами визуализации молекулярной структуры, например, VMD).

Пример промежуточного файла .out  

SEQUENCE: gggcgcaagccu
SCORE: description
SCORE: X_1
   1 0.000 -178.226 54.267 84.464 -152.000 -76.410 89.836 36.841 95.283 179.911 X 1.509 -1.478 -1.022 X X_1
   2 -28.330 -178.226 54.267 84.464 -152.000 -76.410 89.836 36.841 95.283 179.911 X -4.068 -1.073 -3.676 X X_1
   3 -28.330 -178.226 54.267 84.464 -152.000 -76.410 89.836 36.841 95.283 179.911 X -8.530 -1.404 -7.954 X X_1
   4 -28.330 -178.226 54.267 84.464 -152.000 -76.410 -87.682 36.841 93.507 179.911 X -12.167 -3.661 -12.425 X X_1
   5 -28.330 -178.226 54.267 84.464 -152.000 -76.410 89.836 36.841 95.283 179.911 X -15.933 -7.489 -15.504 X X_1
   6 -28.330 -178.226 54.267 84.464 -152.000 -76.410 -87.682 36.841 93.507 179.911 X -20.677 -11.274 -16.721 X X_1
   7 -28.330 -178.226 54.267 84.464 -152.000 -76.411 90.623 36.841 94.812 179.911 X -26.464 -13.437 -17.100 X X_1
   8 -28.330 -178.226 54.267 84.464 -152.000 -76.410 90.622 36.841 94.812 179.911 X -32.524 -13.699 -18.336 X X_1
   9 -28.330 -178.226 54.267 84.464 -152.000 -76.410 89.836 36.841 95.283 179.911 X -37.866 -13.304 -21.437 X X_1
  10 -28.330 -178.226 54.267 84.464 -152.000 -76.410 -87.682 36.841 93.507 179.911 X -42.077 -14.025 -25.916 X X_1
  11 -28.330 -178.226 54.267 84.464 -152.000 -76.410 -87.682 36.841 93.507 179.911 X -45.649 -16.738 -30.182 X X_1
  12 -28.330 -178.226 54.267 84.464 0.000 0.000 -89.325 36.841 93.744 179.911 X -49.587 -20.724 -32.813 X X_1

В заголовке файла .out приводится последовательность цепи. Далее следуют две строки SCORE - описание и значения энергетических характеристик РНК - они рассчитываются модулем RNAEnergyСalc и описаны ниже.

Далее для каждого нуклеотида приводятся 10 углов базовых атомов нуклеотида и координаты x, y, z нуклеотида. Дополнительные предметные детали см. Энергия связей аденин-урацил в РНК.

Основное содержимое .pdb файла это координаты x, y, z всех атомов цепи.

Пример файла .pdb для цепи из двух нуклеотидов au  

ATOM 1 P rA 1 4.445 -0.215 -1.094 1.00 51.50
ATOM 2 O1P rA 1 5.131 1.101 -1.160 1.00 49.22
ATOM 3 O2P rA 1 3.920 -0.832 -2.338 1.00 49.22
ATOM 4 O5* rA 1 3.248 -0.096 -0.051 1.00 49.43
ATOM 5 C5* rA 1 2.954 -1.168 0.864 1.00 49.92
ATOM 6 C4* rA 1 1.495 -1.541 0.748 1.00 50.16
ATOM 7 O4* rA 1 0.680 -0.365 1.001 1.00 49.50
ATOM 8 C3* rA 1 1.052 -2.022 -0.625 1.00 50.01
ATOM 9 O3* rA 1 0.078 -3.066 -0.534 1.00 49.43
ATOM 10 C2* rA 1 0.480 -0.752 -1.242 1.00 50.14
ATOM 11 O2* rA 1 -0.614 -1.001 -2.103 1.00 49.34
ATOM 12 C1* rA 1 0.000 0.000 0.000 1.00 50.16
ATOM 13 N1 rA 1 -2.115 4.525 1.186 1.00 49.26
ATOM 14 C2 rA 1 -2.540 3.265 1.332 1.00 50.50
ATOM 15 N3 rA 1 -1.942 2.127 0.987 1.00 49.25
ATOM 16 C4 rA 1 -0.749 2.374 0.420 1.00 50.43
ATOM 17 C5 rA 1 -0.176 3.614 0.201 1.00 50.28
ATOM 18 C6 rA 1 -0.911 4.739 0.612 1.00 50.46
ATOM 19 N6 rA 1 -0.486 5.996 0.465 1.00 49.23
ATOM 20 N7 rA 1 1.065 3.481 -0.407 1.00 49.29
ATOM 21 C8 rA 1 1.213 2.186 -0.541 1.00 50.34
ATOM 22 N9 rA 1 0.153 1.455 -0.061 1.00 49.95
ATOM 23 1H5* rA 1 3.184 -0.860 1.894 1.00 50.09
ATOM 24 2H5* rA 1 3.586 -2.040 0.640 1.00 50.09
ATOM 25 H4* rA 1 1.368 -2.361 1.470 1.00 50.09
ATOM 26 H3* rA 1 1.863 -2.467 -1.220 1.00 50.09
ATOM 27 1H2* rA 1 1.211 -0.219 -1.868 1.00 50.09
ATOM 28 2HO* rA 1 -0.949 -0.136 -2.477 1.00 50.43
ATOM 29 H1* rA 1 -1.082 -0.182 0.086 1.00 50.09
ATOM 30 H2 rA 1 -3.526 3.154 1.804 1.00 50.13
ATOM 31 1H6 rA 1 -1.056 6.757 0.778 1.00 50.38
ATOM 32 2H6 rA 1 0.401 6.178 0.043 1.00 50.38
ATOM 33 H8 rA 1 2.104 1.726 -0.995 1.00 50.12
ATOM 34 P rU 2 0.333 -4.460 -1.292 1.00 51.50
ATOM 35 O1P rU 2 -0.140 -4.188 -2.673 1.00 49.22
ATOM 36 O2P rU 2 1.745 -4.863 -1.066 1.00 49.22
ATOM 37 O5* rU 2 -0.612 -5.575 -0.657 1.00 49.43
ATOM 38 C5* rU 2 -0.111 -6.494 0.331 1.00 49.92
ATOM 39 C4* rU 2 -1.254 -6.977 1.192 1.00 50.16
ATOM 40 O4* rU 2 -1.628 -5.927 2.122 1.00 49.50
ATOM 41 C3* rU 2 -2.540 -7.301 0.447 1.00 50.01
ATOM 42 O3* rU 2 -2.519 -8.623 -0.100 1.00 49.43
ATOM 43 C2* rU 2 -3.592 -7.143 1.537 1.00 50.14
ATOM 44 O2* rU 2 -3.708 -8.286 2.361 1.00 49.34
ATOM 45 C1* rU 2 -3.020 -5.978 2.345 1.00 50.16
ATOM 46 N1 rU 2 -3.520 -4.655 1.942 1.00 49.66
ATOM 47 C2 rU 2 -4.788 -4.289 2.351 1.00 50.55
ATOM 48 O2 rU 2 -5.510 -5.025 3.000 1.00 49.55
ATOM 49 N3 rU 2 -5.176 -3.028 1.971 1.00 49.54
ATOM 50 C4 rU 2 -4.441 -2.118 1.239 1.00 50.53
ATOM 51 O4 rU 2 -4.930 -1.017 0.976 1.00 49.52
ATOM 52 C5 rU 2 -3.141 -2.571 0.854 1.00 49.85
ATOM 53 C6 rU 2 -2.734 -3.795 1.209 1.00 50.20
ATOM 54 1H5* rU 2 0.650 -6.003 0.954 1.00 50.09
ATOM 55 2H5* rU 2 0.378 -7.348 -0.161 1.00 50.09
ATOM 56 H4* rU 2 -0.871 -7.898 1.656 1.00 50.09
ATOM 57 H3* rU 2 -2.719 -6.656 -0.427 1.00 50.09
ATOM 58 1H2* rU 2 -4.602 -6.989 1.130 1.00 50.09
ATOM 59 2HO* rU 2 -4.409 -8.128 3.056 1.00 50.43
ATOM 60 H1* rU 2 -3.332 -6.125 3.390 1.00 50.09
ATOM 61 H3 rU 2 -6.090 -2.739 2.254 1.00 50.36
ATOM 62 H5 rU 2 -2.478 -1.918 0.266 1.00 50.10
ATOM 63 H6 rU 2 -1.731 -4.127 0.904 1.00 50.14
TER

complete: res ss phi psi omega frag rama seq
   1 0.000 -127.005 60.379 0.00 a
   2 -98.421 154.297 180.000 0.00 u

Формирование цепи

1. Рассчитать координаты x, y, z каждого атома системы
2. Построить дерево связанных атомов
3. Рассчитать углы phi, theta, d для каждого атома системы

Модуль RNARotate „Вращение 3-х мерной цепи РНК“

1. Сформировать множество потенциально возможных поворотов для n (=3, но может быть упрощенно =1) нуклеотидов
2. Повернуть 9 атомов n нуклеотидов на угол phi
3. На основе измененных углов пересчитать координаты x, y, z атомов системы

Модуль RNAEnergyСalc „Расчет энергетических характеристик РНК“

В оценке энергии цепи РНК участвуют следующие составляющие:

SCORE = (VDW + RG) + (RNA_BS + RNA_BP_W + RNA_BP_H + RNA_BP_S) + (RNA_NONB + RNA_O2ST + RNA_PHOS) + (RNA_AXIS*0.2 + RNA_STAG)

Модуль FoldingSearchRnd „Поиск структуры с минимальной оценкой энергии цепи случайным перебором“

1. Установить начальное значение BestScore = LowScore = CurrentScore равным первому значению оценки энергии
2. На основе оценки энергии и случайных поворотов отобрать лучшею позицию
   1. Если CurrentScore > BestScore
      1. Расчитать некоторой уровень возможности Probability такого поворота
      2. Если случайное число R > Probability, то отклонить такой поворот (rejected) и закончить
      3. Иначе, условно принять (accepted by chance)
   2. Если CurrentScore < LowScore
      1. Найден правильный поворот в траектории сворачивания, записать BestScore = LowScore = CurrentScore
      2. Иначе (CurrentScore > LowScore and CurrentScore < BestScore) поворот принять, но не считать удовлетворительным (accepted), записать BestScore = CurrentScore
   3. Случайно выбрать следующею позицию и вид поворота, рассчитать оценку энергии, перейти к пункту 2.1

Пример протокола поиска траектории сворачивания цепи (выполнено 100 итераций поиска)  

Position: 2 FragID: 1280 Change: 8.539 LowScore: 8.539
Position: 4 FragID: 2093 Change: -0.076 LowScore: 8.463
Position: 2 FragID: 1187 Change: -1.328 LowScore: 7.135
Position: 4 FragID: 2550 Change: -0.290 LowScore: 6.845
Position: 3 FragID: 1836 Change: -1.786 LowScore: 5.059
Position: 3 FragID: 1890 Change: -0.376 LowScore: 4.683
Position: 1 FragID: 246 Change: -0.342 LowScore: 4.340
Position: 1 FragID: 169 Change: -2.107 LowScore: 2.233
Position: 3 FragID: 1858 Change: -1.897 LowScore: 0.337
Position: 2 FragID: 1074 Change: -0.313 LowScore: 0.024
Position: 1 FragID: 644 Change: -0.183 LowScore: -0.159
Position: 1 FragID: 453 Change: -0.130 LowScore: -0.289
Position: 1 FragID: 418 Change: -0.270 LowScore: -0.559
Position: 4 FragID: 2630 Change: -0.123 LowScore: -0.682
Position: 1 FragID: 332 Change: -0.156 LowScore: -0.838
Position: 2 FragID: 993 Change: -0.064 LowScore: -0.902
Position: 4 FragID: 2271 Change: -0.484 LowScore: -1.386

Соответствующие конечное состояние цепи (.out файл)  

SEQUENCE: auau
SCORE: SCORE VDW RG RNA_BS RNA_BP_W RNA_BP_H RNA_BP_S RNA_AXIS RNA_STAGGER RNA_BULGE RNA_NONBASEBASE RNA_O2STAR RNA_PHOSPHATE description
SCORE: -1.39 4.42 4.60 0.00 -3.93 -2.85 0.00 0.00 -3.61 0.00 0.00 0.00 0.00 S_1
0 L 0.000 147.040 87.894 83.546 -164.798 -58.486 62.258 35.770 98.030 179.911 X 0.803 -1.500 -1.600 X 332 S_1
1 L -75.525 -176.571 57.560 80.425 -147.629 -76.798 -101.354 36.941 95.812 179.911 X -5.396 -0.599 -1.081 X 993 S_1
2 L 85.934 -163.647 -170.190 82.697 -153.296 -79.199 96.800 36.959 93.513 179.911 X -9.295 3.277 -3.705 X 1858 S_1
3 L -59.665 171.277 54.897 82.011 0.000 0.000 -112.583 35.516 97.364 179.912 X -10.521 9.269 -5.309 X 2271 S_1

Модуль FoldingSearchAI „Поиск структуры с минимальной оценкой энергии цепи применяя методы ИИ“

• В основу будет положено Обучение с подкреплением. Начнем с Q-learning (хотя потенциально возможно использовать и генетические алгоритмы и искусственные нейронные сети ).

Постановка задачи в терминах Q-learning

В качестве агента будет выступать собственно цепь РНК. Средой для агента будет энергетическая поверхность окружающая цепь РНК и формирующая ее 3-х мерную поверхность.

Действием агента будем считать все возможные повороты элементов цепи, а вознаграждением от среды значение оценки энергии этой цепи.

Действие определяется двумя числами: (1) выбранной позицией - нуклеотид который поворачивается и (2) идентификатором поворота - определяет 9 углов поворота с последующим пересчетом позиций атомов нуклеотида. Тогда множество доступных агенту действий - это открытое множество не запрещенных для выбранной позиции поворотов. Учитывая, что действие не задается простым перечислением, оно должно быть реализовано, в отличии от классического версии, динамическим массивом. Элементами массива, будет объект действие (а не просто число), характеризующийся рядом параметров - в данном случае двумя позиция и идентификатор поворота.

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

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

0 - 332 // читается как положение начального (нулевого) нуклеотида в цепи в положении 332
1 - 993 
2 - 1858
3 - 2271

и такое положение однозначно определяется значением энергии равной -1.39 .

Множество возможных положений не превышает 3000. Но множество состояний, в таком случае, зависит от длины цепи. Число состояний можно оценить как $n^{k}$, где n - длина цепи, а k - число возможных положений. То есть в нашем примере $4^{3000}$ - огромное число состояний для нашей маленькой цепи. Нам же нужно исследовать как минимум цепи с длиной около 100 нуклеотидов, т.е. с числом состояний $100^{3000}$. Конечно многие из них вырожденные, но какие именно это еще нужно выяснить.

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

Поэтому от такого способа мы должны отказаться и рассмотреть другой способ задания множества состояний.

Состояния на базе динамического определения.