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

Выбор семейства

Для выполнения работы, я выбрала семейство:

Pfam AC: PF04068
Pfam ID: Fer4_RLI
Name: Possible Fer4-like domain in RNase L inhibitor, RLI
В выравнивании seed: 65 последовательностей
В выравнивании full: 5926 последовательностей

Это семейство включает возможные металл-связывающие домены белков в ингибиторе эндорибонуклеазы РНКазы L. Этот домен был обнаружен на N-конце RLI-белков, рядом с 4Fe-4S-связывающим доменом. РНКаза L является компонентом интерферон-зависимого ответа: она участвует в противовирусной защите, а также, возможно, может влиять на стабильность клеточных РНК и процессы.

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

Подсемейство было выбрано до домменной архитектуре PF04068 - PF00037, в нем 38 белков(рис 1).
Репрезентативным белком данной архитектуры является Q01512 — белок с 4Fe-4S ферредоксин-подобным доменом, предположительно участвующим в связывании железо-серного кластера.

Рис 1. Доменная архитектура выбранного подсемейства.

Построение профиля hmm для подсемейства. Поиск по полученному профилю по белкам семейства.

После скачивания последовательностей белков подсемейства, они были выровнены программой Muscle, и был выделен наш домен PF04068 (рис 2).

Рис 2. Выравнивание белков выбранного подсемейства с выделенным доменом PF04068.

Домен был вырезан в файл PF04068.fa и запущена программа построения hmm профиля:
hmmbuild PF04068.hmm PF04068.fa

Скачали последовательности всех белков семейства (protein-matching-PF04068.fasta), и запустили поиск полученным профилем по этим белкам:
hmmsearch --tblout output_hmmsearch.txt PF04068.hmm protein-matching-PF04068.fasta

Поиск оптимального порога на вес (score)

Для начала, с помощью таблички из выдачи hmmsearch создадим файл, где будут только имена белков и их веса:
grep -v '^#' output_hmmsearch.txt | tr -s ' ' | cut -d' ' -f1,6 > scores.txt
Также из файла со всеми белками подсемейства достанем их названия и тоже запишем в отдельный файл
grep '^>' PF04068.fa | sed 's/^>//' | awk '{print $1}' >domain.txt

Далее полученные файлы анализировала в скрипте на python.
Было выяснено, что всего hmmsearch обнаружил 6050 находки, и все 38 белков подсемейства было найдено.
Минимальный score белков подсемейства — 34.8, а максимальный — 69.0.

Для поиска порога будем использовать следующие метрики:

TP — True positives, число находок, принадлежащих подсемейству с весом выше порога
FP — False positives, число находок не из подсемейства, но в весом выше порога
FN — False negatives, число находок из подсемейства, но в весом ниже порога
TN — True negatives, число находок не из подсемейства и с весом ниже порога

Для удобства сравнения порогов дополнительно рассчитывался F1-cкор. Он учитывает одновременно два требования: сохранить как можно больше представителей подсемейства (со скором больше порога), и уменьшить число ложноположительных находок. Поэтому, F1 зависит от других двух метрик: precision и recall. Precision показывает, какая доля белков выше порога действительно относится к подсемейству, а recall показывает, какая доля белков подсемейства была найдена при данном пороге.
Precision = TP / (TP + FP)
Recall = TP / (TP + FN)
F1 = 2 * Precision * Recall / (Precision + Recall)

Оптимальный порог выбирался не только по максимальному F1, но и с учётом числа TP, чтобы TP был не меньше 20. Это важно, потому что слишком строгий порог может давать мало ложноположительных находок, но при этом терять большую часть настоящего подсемейства, а это грустно. При переборе разных порогов, был найден оптимальный: 57.7

По итогом поиска порога хочется сказать, что профиль оказался далеко не специфичным для выбранной архитектуры: при мягких порогах число FP было большим, а при строгих порогах резко возрастало число FN. Возможно это связано с тем, что подсемейство выделялось по полной доменной архитектуре, тогда как HMM-профиль строился только по одному домену PF04068.

Численные характеристики выделения подсемейства профилем

Для порога 57.7 получились следующие параметры (таблица 1).

Таблица 1. TP, FP, TN, FN для порога 57.7

	True	False
Positives	30	4037
Negatives	1975	8