В качестве исследуемого семейства 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.
Предсказано +
Предсказано -
Фактически +
TP = 9
FN = 4
Фактически -
FP = 19
TN = 19687