Практикум 9. Домены и профили
В ходе этого практикума для некоторого семейства белков (в данном случае - для белков определённой доменной архитектуры) был построен HMM-профиль, который затем тестировался в качестве правила отнесения белка к этому семейству.
Построение входного выравнивания
Для изучения были выбраны белки двухдоменной архитектуры NAD_binding_1-Fer2. Домен Fer2 (2Fe-2S iron-sulfur cluster binding domain, Pfam AC: PF00111) ранее рассматривался в ходе одного из практикумов; почитать о нём подробнее можно в отчёте.
Построение входного seed-выравнивания для профиля требует наличия представительной выборки белков семейства. Чтобы получить её, был проведён поиск белков, содержащих домены NAD_binding_1 и Fer2, в UniProt (без каких-либо ограничений по таксону, запрос: database:(type:pfam pf00175) database:(type:pfam pf00111)); выдача была сохранена в файл .xlsx, который можно скачать по ссылке. Затем из выдачи были удалены белки, содержащие другие домены, помимо изучаемых. В таблице 1 приведённого файла представлен полученный список белков с нужной доменной архитектурой. Далее была построена гистограмма распределения длин этих белков (рисунок 1). Можно видеть, что большинство из них (более 70% от общего числа) имеют длину 310 - 330 а. о., значит, такая длина является типичной для семейства.
С помощью скрипта из таблицы 1 были отобраны 120 белков указанной длины; для того, чтобы в выборке были представлены максимально различные последовательности, скрипт обязательно включал в неё белки из организмов, принадлежащих к различным отделам, а также обеспечивал наибольшее возможное разнообразие на уровне семейств. Список AC выбранных белков приведён ниже.
Просмотреть список
A0A3M1S5K4 A0A3M2D4B9 A0A2G7T4Z1 A0A559NRJ7 A0A6L9ZHG5 A0A6P0R177 W2SVQ9 A0A7I7ZY41 A0A5D0TPQ3 A0A3N5AWQ5 A0A516WW27 A0A653RUR9 A0A1B2HTM6 A0A1I3NKV1 A0A544Y3D4 A0A7J5X4X0 A0A285EDE4 A0A0B2AM47 A0A4V2EYB0 A0A544Y3D4 A0A1X0ANR9 A0A255GXK4 A0A6N9USK9 K6V8Q4 A0A5B8U9E1 A0A3Q9C5Y9 A0A135GFP7 A0A6L3VM62 A0QZU4 A0A3B0FLH2 A0A149ZJ26 A0A4Z0HL62 A0A128ATK7 A0A3S9MIA4 X0Q7U1 A0A1R0KXP7 A0A178X3H5 A0A1H4IIY6 A0A2A3FUL6 I4EU71 A0A4R3IY64 A0A2G7FBA4 A0A417XRX7 K8XM39 A0A0Q6HFN9 A0A433RXR6 A0A0K9YQA0 A0A098ERC7 A0A561DMU0 V6T6H6 A0A1B3XK57 A0A561DMU0 A0A7H8QCZ6 A0A3M8D937 A0A1B3XK57 A0A270A1V8 A0A2A4BPQ9 A0A0Q9VW89 A0A6I1FBI0 A0A426GRL0 A0A2A8MEZ9 A0A3M8D937 A0A3S4PW54 A0A544TFC3 A0A2A8TPP7 W4RQ31 A0A4U2LYF3 A0A433RXR6 A0A544TFC3 A0A3N9U978 A0A0Q9VW89 A0A2A8FJX9 A0A6I1FBI0 A0A4Y8UIF7 A0A6I1FBI0 A0A109MS45 A0A7H8QCZ6 A0A2A8MEZ9 A0A561DMU0 A0A223EMP0 A0A0U2U4E3 A0A4Q7FGE3 W6V0S3 A0A017HL72 A0A0B1YD87 A0A5N7TYJ3 A0A537R606 A0A3E1SEQ2 A0A210VXX3 A0A1Y0LFF2 A0A3S2BB31 Q1R0Q3 H8MV48 A0A198YI16 A0A4U9HX30 A0A261RQX6 A0A370U9J0 A0A323UPM0 A0A2S6Z517 A0A375D9A3 A0A0C2E584 A0A5R9QEY8 A0A2G6X6W3 A0A2K4TXS2 A0A1N6YYZ2 A0A3S0AKV3 A0A1G9Z268 A0A285ZEZ8 A0A6G6EQE2 A0A2W5QFB1 A0A4Q4GZW7 A0A2U1WWL1 A0A244ELA6 A0A6I6TTR5 A0A0W0JT40 A0A0W0HN51
Затем последовательности выбранных белков были скачаны из UniProt и выровнены программой muscle. С помощью Jalview выравнивание было отредактировано: удалены неконсервативные колонки на N и C-концах, убраны последовательности, значительно отличающиеся от других на участках высокой консервативности (вероятно, их наличие объясняется тем, что в выборку попали белки, в которых домены NAD_binding_1 и Fer2 присутствуют, но расположены в обратном порядке); также были удалены избыточные последовательности (порог идентичности - 85%). Полученное выравнивание, в котором осталось 70 последовательностей, доступно по ссылке.
Построение HMM-профиля семейства белков
Задача построения профиля была решена с помощью пакета HMMER. На вход программе hmm2build было подано описанное выше seed-выравнивание, на выходе получен HMM-профиль:
hmm2build NAD_binding_1_Fer2_profile.hmm NAD_binding_1_Fer2.fasta
Далее профиль был откалиброван программой hmm2calibrate, которая на основе применения профиля к случайным последовательностям рассчитывает коэффициенты нормализации весов:
hmm2calibrate NAD_binding_1_Fer2_profile.hmm
Файл с профилем доступен по ссылке.
Тестирование профиля
Для того, чтобы проверить, позволяет ли полученный профиль распознавать белки данной доменной архитектуры, из UniProt были скачаны последовательности всех белков, содержащих домен Fer_2 (без ограничений по таксону), после чего среди них был проведён поиск по профилю:
hmm2search --cpu 1 -E 0.01 NAD_binding_1_Fer2_profile.hmm Fer_2.fasta > search_res.txt
Результаты поиска приведены в таблице 2 файла .xlsx по ссылке. Далее было визуализировано распределение весов находок; полученный график изображён на рисунке 2, на нём хорошо заметна ступенька веса.
Чтобы обосновать выбор порога веса для отнесения белка к семейству, был применён метод автоматического перебора порогов: для каждого возможного значения порога (т. е. для каждой находки) были рассчитаны чувствительность, специфичность и значение F1, затем был найден максимум последнего; он составил примерно 0.5026 при значении веса 434.1, оно и было принято за пороговое. Все вычисления также приведены в таблице 2. Кроме того, по полученным данным была построена ROC-кривая, изображённая на рисунке 3.
Можно отметить некоторую странность этой кривой, а именно то, что ни при каком значении порога чувствительность теста не достигает 100%. Возможно, это связано с тем, что среди белков, для которых должна быть известна их принадлежность к семейству, значительную долю составляют белки с доменами, расположенными в обратном порядке (они попали в их число из-за ограниченных возможностей поиска UniProt, который не учитывает порядок доменов). Поиск по профилю не обнаруживает такие белки, но они ошибочно учитываются в числе принадлежащих к семейству, что сказывается при расчёте чувствительности.
Для оценки работы профиля как правила отнесения белков к семейству была составлена таблица «предсказание против истины» (см. таблицу 3 ниже):
| Истинно + | Истинно - | |
|---|---|---|
| Предсказано + | 4360 | 4822 |
| Предсказано - | 1892 | 42378 |
Исходя из этой таблицы, можно сделать вывод, что правило работает плохо, выдавая неприлично много ложноположительных результатов. Однако, если проверить, какие белки вызывают эти ложноположительные результаты, можно заметить, что значительная их часть действительно имеет нужную доменную архитектуру (по аннотации UniProt), но вот в Pfam по какой-то причине обозначен только один домен. Это наводит ны мысль о том, что не вполне корректно использовать ссылки на Pfam в записи белка в качестве критерия его принадлежности к семейству. В связи с этим была предпринята попытка искать по доменам в поле FT (запрос: annotation:(type:"positional domain" "fad-binding fr-type") annotation:(type:"positional domain" "2fe-2s ferredoxin-type")), выдача была сохранена, белки с другими доменами - убраны из неё, и теперь этот список белков был принят за семейство. Все описанные выше расчёты по оценке работы правила были проведены ещё раз для нового критерия «истины» (ссылка на файл), новый порог веса составил 258.9 при максимуме F1 0.7503. Ниже приедена новая таблица «предсказание против истины»:
| Истинно + | Истинно - | |
|---|---|---|
| Предсказано + | 10373 | 5468 |
| Предсказано - | 1445 | 36166 |
Можно видеть, что, если судить по новой проверке, правило работает несколько лучше. К тому же, теперь пороговое значение веса соответствует ступеньке в распределении весов, что также свидетельствует в пользу большей корректности новой оценки. На рисунке 4 приведена новая ROC-кривая:
Как видно из графика, проблемы с чувствительностью не наблюдается.