HMM профиль, ROC-кривая

Составление списка белков целевого семейства

В качестве целевого семейства было выбрано семейство белков RF-1 (release factor 1, фактор терминации трансляции 1), содержащее домены RF-1 (PF00472). Согласно указаниям по выполнению практикума, на профиль были наложены дополнительные условия. В данном практикуме мною будет построен профиль для класса Clostridia. Для поиска аннотированных последовательностей представителей класса Clostridia, содержащих RF-1 домен в Uniprot был введен следующий поисковый запрос:

database:(type:pfam id:PF00472) taxonomy:clostridia AND reviewed:yes

Колонки выходной таблицы были обработаны таким образом, чтобы в ней были ID записи (Entry name), Fragments, Protein length, Protein name, таксономия нужного уровня (class), домены Pfam, после чего файл, содержащий последовательности, был сохранен в формате tab-separated и fasta. Таким образом была получена информация о 42 находках, отвечающих заданным требованиям, и сами последовательности. Содержание файла формата tab-separated было добавлено в итоговую excel таблицу (лист «Лист1») а fasta-файл экспортирован в JalView для дальнейшего редактирования. Множественное выравнивание было выровнено при помощи сервиса Muscle и покрашено Clustal, by conservation 50. После визуальной оценки качества выравнивания мною были удалены 17 из 42 находок, содержащих слишком много гэпов (неинформативные) или плохо выровненных относительно большинства последовательностей. Файл с оставшимися последовательностями можно скачать по ссылке.

Рисунок 1. Выравнивание факторов терминации трансляции 1 (RF-1) класса Clostridia

Построение и калибровка профиля

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

hmm2build -g profile.out cleaned.fasta

hmm2calibrate profile.out

Результат работы программы: profile.out

После чего проведен поиск с помощью полученного профиля по базе данных uniport, команда:

hmm2search --domE 1000 --domT -50 profile.out /srv/databases/emboss/data/uniprot/uniprot_sprot.fasta > result.txt

Результат работы программы: result.txt
domE - порог E-value для выравнивания белковых доменов, domT - пороговое значение score для выравнивания белковых доменов. Результат работы программы, записанный в файл, включает список параметров запуска, таблицу найденных последовательностей, таблицу найденных доменов, выравнивания находок относительно профиля (показан консенсус профиля) и гистограмму весов находок. Информация из полученного файла также была добавлена в итоговый excel файл (лист «Лист2»). Из таблицы доменов, содержащей ID последовательности, координаты выравненного участка в последовательности и в профиле, E-value и вес T., на другой лист excel документа («Лист3») были импортированы данные об ID последовательностей, E-value и score, после чего был проведен поиск совпадений со списком ID находок из результатов поиска по базе данных Uniprot с дальнейшим выравниванием в JalView (т.е. 25 находок, оставшихся после чистки выравнивания, информация о которых содержится на первом листе файла excel).

На рисунках 2 и 3 соответственно изображены гистограмма и диаграмма весов (score) находок.

Рисунок 2. Гистограмма весов находок.	Рисунок 3. Диаграмма весов находок.

Затем для полученных данных были построена ROC-кривая (Receiver Operator Characteristic) — кривая, которая наиболее часто используется для анализа качества моделей (показывает зависимость количества верно классифицированных положительных примеров от количества неверно классифицированных отрицательных примеров). Для построения кривой были определены:

TP (True Positives) — верно классифицированные положительные примеры (так называемые истинно положительные случаи). Это количество последовательностей, расположенных выше порога и достоверно содержащих искомый домен;

TN (True Negatives) — верно классифицированные отрицательные примеры (истинно отрицательные случаи). Это количество последовательностей, расположенных ниже порога и не содержащих искомый домен;

FP (False Positives) — отрицательные примеры, классифицированные как положительные (ошибка II рода, ложно положительные случаи). Это количество последовательностей, расположенных выше порога, но не содержащих домен;

FN (False Negatives) — положительные примеры, классифицированные как отрицательные (ошибка I рода, ложно отрицательные примеры). Это количество последовательностей, расположенных ниже порога и содержащих домен;

(SP) — специфичность, доля истинно отрицательных случаев, которые были правильно идентифицированы моделью (т.е. доля предсказанных белков, не содержащих домен, от общего количества последовательностей, известно не содержащих этот домен): SP = TN/(TN+FP);

(SE) — чувствительность, доля истинно положительных случаев (т.е. доля предсказанных белков, содержащих домен, от общего количества последовательностей, известно содержащих домен). Вышеупомянутые показатели считаются по следующим формулам: SE = TP/(TP+FN).

Таблица с данными, полученными для построения ROC-кривой нашего профиля приведена ниже.

Таблица 1.

Рисунок 4. ROC-кривая

Заключение

В случае, если мы хотим найти всех представителей выбранного таксона (класс Clostridia) с целевым доменом RF-1 (PF00472), включенных в итоговое выравнивание (25 последовательностей), и допускаем наличие ложно положительных находок, то нам придется принять значение cutoff = 758,9. И если мы не хотим обнаружить ни одной ложной находки, в таком случае принимаем значение Cutoff = 851,9.

Таблица 2.

На мой взгляд, данный профиль является неэффективным для поиска: если по нему находятся все целевые последовательности, список находок включает в себя 37,5% (15/(15+25)) побочных находок; если список находок на 100% состоит из достоверных находок, то он находит лишь 52% (13/25) целевых последовательностей.

К семестрам