Моделирование вращения нуклеотидов в РНК

Постановка задачи

Предположим, что мы моделируем вращение нуклеотидов при сворачивании шпильки GCAA tetraloop (GGGCGCAAGCCU). Нужно реализовать вращение случайно выбранных нуклеотидов.

В программном обеспечении Rosetta

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

a 0 -157.884 174.248 87.1186 -115.37 -62.0781 53.6519 34.9201 90.5193 0 0 0 0
u -75.1395 173.083 51.5993 78.3736 -130.967 -49.3271 -106.588 39.5321 92.0595 0 0 1 0
g -53.9271 172.606 60.4137 76.7298 -166.974 -80.0958 69.6379 35.731 104.345 0 1 0 0
c -71.8925 -170.195 55.8123 84.2024 -146.216 -75.9283 -97.0879 35.2839 95.6198 0 1 0 0
c -58.663 174.597 44.3207 82.0836 -156.653 -70.2968 -98.5083 35.8355 96.9331 0 1 0 0
a -74.1112 179.128 55.1376 83.5659 -154.388 -72.7722 81.6324 35.6556 95.5311 0 1 0 0
g -67.829 164.242 69.9991 80.48 -160.055 -66.653 72.7704 35.6534 99.6295 0 1 0 0
c -74.5377 -176.783 54.773 84.2962 -142.208 -81.7383 -100.013 36.082 95.0934 0 1 0 0
u -62.2138 173.078 43.5885 81.2923 -154.694 -65.1851 -98.1924 38.1824 90.2845 0 1 0 0

Формат таблицы: 1 значение вид нуклеотида (vall_current_sequence), далее 9 углов соответствующих атомов, далее признак chainbreak, далее 3 признака типа углов.

Допустим была случайно выбран нуклеотид на 8-ой позиции для осуществления вращения. Вырежем последовательность из 3-х нуклеотидов, начиная с этой позиции - agc. Заменим g и a на r, а c и u на y, согласно следующему алгоритму

for (int i = 0; i < 3; i++) 
{
    if (RNA_string[ i ] == 'g' || RNA_string[ i ] == 'a' )
    {
        RNA_string_local[ i ] = 'r';
    }
    else
    {
        RNA_string_local[ i ] = 'y';
    }
}

В нашем случае получим строку RNA_string = rry. Из всего множества нужно отобрать только те, которые подходят для нашей тройки нуклеотидов. Это можно сделать в соответствии с алгоритмом:

for (int i = 1; i <= 2705 - 3 + 1; i++ )
{
    bool match = true;

    for (int offset = 0; offset < 3; offset++ )
    {
        if ( vall_is_chainbreak_( i + offset ) ||
                !compare_RNA_char( vall_current_sequence[offset], RNA_string[ offset ] ) )
        {
            match = false;
            break;
        }
    }

    if (match) 
    {
        fragment_library_p->add_torsion( vall_torsions_, i, size, true /*flip_chi_purines*/,  vall_current_sequence, vall_name_ );
    }
}

bool compare_RNA_char(char char1, char char2) 
{
    if (char1 == char2) return true;
    if (char1 == 'r' && (char2 == 'a' || char2 == 'g')) return true;
    if (char1 == 'y' && (char2 == 'c' || char2 == 'u')) return true;
    if (char2 == 'r' && (char1 == 'a' || char1 == 'g')) return true;
    if (char2 == 'y' && (char1 == 'c' || char1 == 'u')) return true;
    return false;
}

В результате будет найдено 334 потенциальных фрагмента для возможного вращения.