HMM-профиль подсемейства дельта эндотоксинов.

Дельта-токсины.

Для анализа было выбрано pfam-семейство дельта эндотоксинов (PF03944). Дельта эндотоксины это белковые токсины, вырабатывающиеся бактериями (например Bacillus thuringiensis) и убивающие некоторых насекомых (благодаря чему применяются в генной модификации растений).

Подсемейство для анализа было выбрано по доменной архитектуре: PF03945 - PF00555 - PF03944 - PF17997 - PF21463. Оно содержит в себе 266 белков и является вторым подсемейством по количеству (первое содержит 268 белков). В самом семействе содержится 1116 последовательностей белков.

HMM-профиль: построение и поиск.

Для построения профиля использовался домен PF03944 из последовательностей белков выбранного подсемейства. Построение профиля проводилось программой hmmbuild:
hmmbuild domain.hmm domain.fasta

После этого, используя полученный профиль и программу hmmsearch, был совершен поиск по нему во всех последовательностях семейства:
hmmsearch --tblout found_domain.txt domain.hmm protein-matching-PF03944.fasta

Анализ находок.

Для анализа выходного файла был написан python-скрипт:

#!/usr/bin/env python3

import argparse

expected = set()
with open("domain.fasta") as f:
    for line in f:
        if line.startswith(">"):
            prot_id = line[1:].split()[0]
            expected.add(prot_id.rsplit("/", 1)[0] if "/" in prot_id else prot_id)

parser = argparse.ArgumentParser()
parser.add_argument("threshold", type=float)
parser.add_argument("input_file", type=argparse.FileType("r"))
args = parser.parse_args()

tp = fp = fn = tn = 0
for line in args.input_file:
    if line.startswith("#"):
        continue
    parts = line.split()
    if len(parts) < 6:
        continue
    score = float(parts[5])
    in_expected = parts[0] in expected
    if score >= args.threshold:
        if in_expected:
            tp += 1
        else:
            fp += 1
    else:
        if in_expected:
            fn += 1
        else:
            tn += 1

print(tp, fp, fn, tn)
	

python script.py 250 found_domain.txt

Вывод программы (tp, fp, fn, tn): 213 271 53 461

Также результат приведен в таблице 1. К TN добавлены последовательности белков, не попавшие в выдачу hmmsearch (118).