Для работы возьмём доменную архитектуру IPR000477 - IPR002711, включающую 195 белков, согласно базе данных InterPro. Первый домен в этой архитектуре - обратная транскриптаза, а второй - HNH-эндонуклеаза.
На сайте InterPro я сгенерировал fasta-файл с последовательностями этих 195 белков. В программе Jalview я провёл Mafft-выравнивание этих последовательностей (см. рис. 1).
Вот ссылка на файл с выравниванием.
Далее возьмём только те координаты из этого выравнивания, которые содержат интересующие нас домены. Для этого посмотрим на координаты этих доменов в репрезентативном для выбранной архитектуры белке с АС A0A0A0HXF5. Эти координаты: 87 - 576. Вырежем именно эти координаты в выравнивании и далее работать будем с этим уменьшённым файлом. К тому же удалим некоторые сиквенсы, которые выровнялись совсем плохо (на глаз). Вот получившийся файл, состоящий из 172 последовательностей.
С помощью несложного скрипта на Python3 я записал в файл accession numbers для всех последовательностей в полученном выравнивании.
На сервере Kodomo я запустил программы из пакета HMMER 2.3.2.
Чтобы создать HMM-профиль для выравнивания:
hmm2build hmmout align_ed.faЧтобы откалибровать полученный файл:
hmm2calibrate hmmoutРезультат: файл hmmout.
Теперь выполним поиск по HMM-профилю в файле со всеми белками, содержащими данную архитектуру (kiselev-195-full.fasta). Это файл с исходными сиквенсами, скачанными из InterPro:
hmm2search --cpu=1 hmmout kiselev-195-full.fasta > hmm_results.txtВ этом файле представлено несколько таблиц, первая из которых особенно интересна. В ней указаны веса и значения E-value для всех сиквенсов в kiselev-195-full.fasta. Построим диаграмму распределения весов последовательностей:
За пороговое значение веса можно взять 600.
В Python3 построили пару графиков, показывающих достоверность полученных данных:
