Из базы данных Pfam выбрал GAGA (PF09237). Домен имеет 2 белка в seed и 186 белков в full, длина профиля HMM - 46 аминокислот.
Всего имеется 26 архитектур. Выбрал архитектуру из трех доменов (BTB, GAGAx2), встречается в 15 последовательностях. Схема ниже:
Сначала скачал все последовательности в FASTA формате [1].После чего отобрал уникальные AC всех последовательностей.[2] Далее с сайта PFAM скопировал AC последовательностей соотвествующих архитектуре.[3] После чего было необходимо соотнести данный список со всеми последовательностями [4] и получить новый сокращенный список. [5]
Результат я загрузил в Jalview для выравнивания. Последовательностей 6, их колличество решил не сокращать. Удалил фрагменты выравнивания после 3 домена. Общая длина профиля составила 441 а/к. Визуализация ниже
К полученному отредактированному выравниваю применил серию команд:
hmm2build profile fin.fasta
hmm2calibrate profile
hmmsearch --cpu=1 profile full.fasta &> log.txt
hmm2search --cpu=1 profile full.fasta &> log2.txt
hmm2build строит HMM-профиль.[6] hmm2calibrate калибрует его. hmmsearch и hmm2search производит поиск изучаемой двудоменной архитектуры в переданных последовательностях с порогом E-value 0.01. [7][8] Из логов достал таблицу результатов.[9] удалил строки с p-value выше 0.01, также в 3 случаях профиль не был найден, эти строки также исключил.
Далее по данным выше полученных файлов(таблица лог, список всех последоваательностей, список их АС, список АС трехдоменной структуры) была получена общая таблица (колонки: АС, есть ли наша структура, участвовал ли в создании профиля, hmm2search, score, e-value, длина белка)[10][11]
Из второго лога получил таблицу с 129 находками.[12] Поместив результат в excel [13], посчитал specificity, 1 - sensitivity и F1. В последствии получил два графика: сначала гарфик распредееления весов. Падение наблюдается не слишком резкое, визуально можно оценить порог(300-400). Далее визуализировал ROC кривую. Для определния порога оптимальными величинами специфичности и чувствиетльности посчитал F1 (для ROC и F1 добавил стоблец architect таблицы говорящий о принадлежности белка к архитектуре). Максимальное возможное значение F1 наблюдалось при пороге 723.6. При оценке ROC кривой значение будет приблизительно тоже (точка 0.03 по оХ и 0.83 по оУ), сопотавляя значения с таблицей получаем совпадение. Что сильно отличается от определения на глаз. На графике распределния весов в этой точке начинается самый резкий скачок из имеющихся.