Для написания отчета по данному практикуму я выбрал семейство Ribosome recycling factor (AC: PF01765, ID: RRF).
Белки с данным доменом катализируют освобождение рибосомы с мРНК после терминации трансляции.
Последовательностей в seed и full 631 и 12825 соответственно.
Подсемейство было выбрано по доменной архитектуре. Оно содержало два домена: HSP70 и RRF.
Всего в подсемейство содержится 26 белков. Эти последовательности использовались в выравнивании для создания HMM профиля.
Рисунок 1. Выравнивание подсемейства с доменом RRF из 26 белков. Выделение домена показано красной полосой сверху.
При помощи MAFFT и Jalview было создано выравнивание, из которого был вырезан домен:
hmmbuild domen.hmm domenRRF.fa
Полученный профиль был использован для поиска по базе всех белков:
Для анализа использовался порог score = 200, выбранный по резкому разрыву между соседними значениями.
Для обработки результатов был написан Python-скрипт:
#!/usr/bin/env python3
THRESHOLD = 200.0
def fasta_ids(path):
ids = set()
with open(path) as f:
for line in f:
if line.startswith(">"):
ids.add(line[1:].split()[0])
return ids
all_ids = fasta_ids("all_proteins.fasta")
subfamily_ids = fasta_ids("protein-sequences-3.fasta")
scores = {}
with open("out.tbl") as f:
for line in f:
if line.startswith("#") or not line.strip():
continue
parts = line.split(maxsplit=6)
target = parts[0]
score = float(parts[5])
scores[target] = score
TP = FP = FN = TN = 0
for prot in all_ids:
score = scores.get(prot, float("-inf"))
in_subfamily = prot in subfamily_ids
above = score >= THRESHOLD
if in_subfamily and above:
TP += 1
elif (not in_subfamily) and above:
FP += 1
elif in_subfamily and (not above):
FN += 1
else:
TN += 1
print(TP, FP, FN, TN)
Результаты:
TP
FP
FN
TN
26
156
0
30303
Все белки из подсемейства прошли порог, также его прошли некоторые другие белки, не относящиеся к подсемейству.
Всего белков, которые прошли порог - 182.
