HMM-профили и эволюционные домены

Выбор семейства и подсемейства Pfam

Мное было выбранно семейство Ribosomal protein S8e. Информация об этом семействе:

Белки, содержащие домен PF01201, являются структурными компонентами малой (40S) субъединицы рибосомы. Они играют ключевую роль в биогенезе рибосом и процессе трансляции, участвуя в связывании с рибосомной РНК (рРНК) и формировании функционального центра рибосомы. Кроме того, некоторые белки этого семейства (например, NSA2) участвуют в биогенезе рибосом на более ранних этапах.

Далее было выбранно подсемейство: PF01201 - PF00240 - PF01599, опираясь на доменную архитектуру (чтобы было хотя бы 3 домена), а также на количество белков в этом подсемействе (чтобы их не было слишком мало). В этом подсемействе 18 белков, и последовательности всех использовались в выравнивании для построения HMM-профиля.

Построение HMM-профиля

Я скачал последовательности 18 белков с Pfam, оставил в них только домены и выровнял программой Muscle. Далее это выравнивание использовалось для построения HMM-профиля.

# Построение HMM-профиля
hmmbuild my_profile.hmm for_hmm.fa

Настройка HMM-профиля

Для настройки порога на оптимальный вес находки, построенный мною профиль был применён ко всем белкам из выбранного семейства. В этом семействе всего 13 423 белка.

# Применение HMM-профиля к белкам из выбранного семейства
hmmsearch --tblout search_results.txt my_profile.hmm protein-matching-PF01201.fasta.gz

Опция --tblout нужна для того, чтобы полученный файл был в виде удобной таблицы.

Далее были посмотрены результаты для всех белков и только из подсемейства (для подсемейства представленно ниже).

target name                                 score
A0A423TU39|unreviewed|Ubiquitin-ribosomal   472.4
A0AAW2HZS4|unreviewed|Ribosome              568.2
A0A9J6GXU3|unreviewed|Ribosome              572.5
A0A1D2N322|unreviewed|Ribosome              591.4
A0A9N9S041|unreviewed|Ribosome              613.7
A0AAE1NWN3|unreviewed|Ribosome              615.1
A0ACB7RS43|unreviewed|Uncharacterized       617.7
A0A232ES75|unreviewed|Ribosome              622.0
A0A4Y7M7C8|unreviewed|Ribosome              629.6
A0A4Y7NKT8|unreviewed|Ribosome              631.5
A0A9J6E8G4|unreviewed|Ribosome              635.9
A0A212FHW0|unreviewed|Ribosome              639.8
A0A1A9ZVP2|unreviewed|Ribosome              642.4
A0A8S4REJ0|unreviewed|Ribosome              645.1
A0A182TUI2|unreviewed|Ribosome              647.2
A0A182SGG4|unreviewed|Ribosome              648.4
A0AAV2SG52|unreviewed|Ribosome              648.5
W5J6B5|unreviewed|Ribosome                  648.9

Мной был выбран порог score >= 568.2, так как при пороге равном минимальному весу подходящего белка, количество ложноположительных находок равно 30, а при повышении порога до указанного мы теряем всего один подходящий белок (не находим его), но при этом снижаем количество ложноположительных результатов в 3 раза (с 30 до 10) и получаем что количество истиноположительных белков оказывается найдённо больше, чем ложноположительных. Ниже приведена таблица (табл. 1) с параметрами TP (True Positive), FP (False Positive), FN (False Negative), TN (True Negative).

Таблица 1. Численные характеристики выделения подсемейства профилем
Positive Negative
True 17 13 395
False 10 1

К сожалению у меня не достаточно знаний, что бы оценить какой порог был бы наилучшим. Для приведённого порога чувствительность равна 0,944 и специфичность равна 0,999.