Практикум 10. HMM-профили и эволюционные домены


В качестве исследуемого семейства PFAM был выбран C-концевой домен альфа-субъединиц F-АТФаз (PF00306). Трансмембранные АТФазы — мембранные ферменты, которые либо используют гидролиз АТФ для переноса протонов, либо синтезируют АТФ за счёт протонного градиента; F-, V- и A-типы состоят из каталитического (F1/V1/A1) и мембранного (F0/V0/A0) комплексов и работают как роторные моторы. У F-АТФаз каталитический центр F1 образован субъединицами α и β (β — каталитические, α — регуляторные), а вращение γ-субъединицы, движимое протонами через F0, индуцирует синтез АТФ. У V- и A-АТФаз, напротив, каталитической является субъединица α, а β — регуляторной; данная запись представляет C-концевой домен α-субъединицы.

В выбранном семействе было выбрано подсемейство по доменной архитектуре PF00137 - PF02874 - PF00006 - PF00306, представленное 13 белками.

Последовательности белков выбранного подсемейства были скачаны и выравнены с помощью MAFFT с параметрами по умолчанию. Далее из этого выравнивания были выделены последовательности выбранных доменов и по ним был построен HMM-профиль, используя следующую команду:

hmmbuild --amino hmm.output domains_pr10.fa

Далее нужно провести поиск по всем последовательностям выбранного семейства с помощью программы hmmsearch. Суммарно в данном семействе описано 49.000 белков, что слишком много. Поэтому было решено отобрать белки только из группы Viridiplantae (зеленые растения) и Metazoa (многоклеточные животные), поскольку в выбранном мной подсемействе находятся белки из этих организмов. Суммарно было скачано 19.725 последовательностей (16.958 с Viridiplantae и 2.767 с Metazoa).

С помощью следующей команды был осуществлен поиск по выделенным последовательностям с помощью программы hmmsearch:

hmmsearch -o hmmsearch.output --tblout tableoutput.txt hmm.output allseq.fasta

Опция --tblout позволяет создать удобный файл для парсинга

Далее определим порог, по которому мы наиболее эффективно будем различать белки нашего семейства. Сначала соберем в один файл информацию о белках и их скорах при поиске hmmsearch. Была использована следующая команда:

grep -v '#' tableoutput.txt | tr -s ' ' | cut -d' ' -f1,6 > allscores.txt

Был составлен скрипт на Python, который выдает метрики для каждого score по белкам, позволяя определить наиболее оптимальное значение скора. Выдача программы:

score = 318.1, TP = 13, FP = 19047, FN = 0, TN = 659

score = 685.6, TP = 12, FP = 13354, FN = 1, TN = 6352

score = 689.1, TP = 11, FP = 13342, FN = 2, TN = 6364

score = 769.4, TP = 10, FP = 8343, FN = 3, TN = 11363

score = 788.4, TP = 9, FP = 19, FN = 4, TN = 19687

score = 791.7, TP = 8, FP = 18, FN = 5, TN = 19688

score = 797.9, TP = 7, FP = 18, FN = 6, TN = 19688

score = 799.5, TP = 6, FP = 18, FN = 7, TN = 19688

score = 806.2, TP = 5, FP = 15, FN = 8, TN = 19691

score = 808.3, TP = 4, FP = 14, FN = 9, TN = 19692

score = 810.7, TP = 3, FP = 12, FN = 10, TN = 19694

score = 814.4, TP = 2, FP = 6, FN = 11, TN = 19700

score = 821.4, TP = 1, FP = 5, FN = 12, TN = 19701

Как мне кажется, самым лучшим порогом будет являться тот, при котором наблюдается максимально возможное количество True Positives, минимальное количество False Positives и максимальное количество True Neagtives. Как мне кажется, наиболее оптимальным порогом является 788.4. Метрики для этого порога изображены в таблице 1.

Таблица 1. Метрики для наиболее оптимального найденного скора.

Предсказано + Предсказано -
Фактически + TP = 9 FN = 4
Фактически - FP = 19 TN = 19687