1. Составление списка белков целевого семейства из `SwissProt`
В качестве целевого семейства были выбраны белки Млекопитающих, содержащих домен TIR (PF01582), который ранее анализировался в
практикуме 7 . Для этого в Uniprot был введен следующий поисковый запрос:
database:(type:pfam id:PF01582) taxonomy:"Mammalia [40674]" AND reviewed:yes
Полученные данные были собраны в таблице excel. Скачать файл по
ссылке, страница "all data". Для анализа было выбрано подсемейство белков, содержащее домены PF13895;PF01582 и PF13855;PF01582; (Первый содержит помимо TIR еще иммуноглобулиновый домен, а второй - лейцин-богатый домен). Список отобранных белков можно найти на второй странице ("selected_subfamily") того же файла.
2. Построение профиля
Затем по одной и той же схеме было построено два профиля. Один - рандомная выборка белков, содержащих домен PF01582 (в Jalview PF01582 seed), другая - на основе выравнивания 22 белков, принадлежащих к выбранному подсемейству (выравнивание получено с помощью
скрипта как в практикуме 7). C помощью пакета HMM был получен профиль, откалиброван, а также проведен поиск с помощью этого профиля по базе данных uniprot.
Сводная таблица, содержащая команды и результаты.
3. Анализ полученных результатов
На 3 и 4 листе файла
excel находится информация о находках доменов для двух построенных профилей. Профиль, составленный из рандомной выборки нашел 15 из 30 последовательностей указанного подсемейства и 47 из 99 последовательностей, принадлежащих млекопитающим. Профиль, составленный из белков выбранного подсемейства нашел 29 из 30 последовательностей указанного подсемейства и 83 из 99 последовательностей белков млекопитающих. Уже при таком грубом анализе можно сказать, что профиль, построенный на основе подсемейства ищет белки подсемейства и млекопитающих лучше, чем профиль, построенный на рандомной выборке. Гистограммы весов находок, для данных профилей приведены на рисунке 1.
Рисунок 1. Гистограммы весов для находок для профиля, составленного на основе рандомной выборки (выше) и на основе выбранного подсемейства (снизу). Красной стрелочкой отмечена ступенька (вероятностный порог), который я изначально выбрала бы "на глаз", зеленой стрелочкой - реально установленный порог, полученный после построения ROС-кривой.
Затем для полученных данных были получены ROC-кривые (Receiver Operator Characteristic) (показывают зависимость количества верно классифицированных положительных примеров от количества неверно классифицированных отрицательных примеров). Для построения кривой были определены:
TP (True Positives) — верно классифицированные положительные примеры (так называемые истинно положительные случаи). Это количество последовательностей, расположенных выше порога и содержащие искомый домен;
TN (True Negatives) — верно классифицированные отрицательные примеры (истинно отрицательные случаи). Это количество последовательностей, расположенных ниже порога и не содержащих искомый домен;
FP (False Positives) — отрицательные примеры, классифицированные как положительные (ошибка II рода, ложно положительные случаи). Это количество последовательностей, расположенных выше порога, но не содержащих домен;
FN (False Negatives) — положительные примеры, классифицированные как отрицательные (ошибка I рода, ложно отрицательные примеры). Это количество последовательностей, расположенных ниже порога и содержащих домен.
Кроме того необходимо определить еще несколько параметров: специфичность (SP) — доля истинно отрицательных случаев, которые были правильно идентифицированы моделью (т.е. доля предсказанных белков, не содержащих домен, от общего количества последовательностей, известно не содержащих этот домен) и чувствительность (SE) — доля истинно положительных случаев (т.е. доля предсказанных белков, содержащих домен, от общего количества последовательностей, известно содержащих домен). Эти показатели считаются по следующим формулам: SP = TN/(TN+FP), SE = TP/(TP+FN).
ROC-кривая представляет собой зависимость SE от 1-SP. Приведенные вычисления проведены на 5 и 6 страницах указанного выше excel-файла. Полученные ROC-кривые приведены на рисунке 2. Для выбора порога необходимо максимизировать SE и SP, то есть выбрать такую пару, для которой максимально значение (SE+SP-1). Для неспецифичного профиля это оказалось 0,229, соответствующее значение порога для веса: 143,6; для специфичного профиля - 0,585, вес - 213,8.
Рисунок 2. ROC-кривые для неспецифичного профиля (выше) и специфичного (ниже)
Таблица 2. Результаты поиска для выбранного порога
|
TP |
TN |
FP |
FN |
profile1 |
12 |
27 |
36 |
3 |
profile2 |
26 |
44 |
20 |
3 |
На основе полученных результатов можно сказать, что профиль, построенный по выборке белков из подсмейста выделяет это подсемемейство лучше, чем профиль, выбранный из рандомных последовательностей, содержащих этот домен. Кроме того полученный профиль обладает специфичностью 68,8% и чувствительностью 89,7%, что позволяет его использовать для поиска белков млекопитающих, содежащих доменные архитектуры PF13895;PF01582 и PF13855;PF01582.
4. Поиск blast по SwissProt
Для оценки качества полученного профиля был проведен референсный поиск гомологов белка, входящего в выбранное семейство (Q9HB29, ILRL2_HUMAN). Поиск проводился в swissprot с помощью алгоритма blastp, wordsize = 6, E-value = 100. Список находок приведен на листе 7 (blast). Бласт нашел 102 находки, однако их параметры значительно хуже обоих профилей. Для выбранного подсемейства FP = 11, TN = 91, FP =19. Для млекопитающих FP = 26, TN = 76, FP = 73.
5. Источники
[1]
http://www.uniprot.org