Был выбран домен Acetone_carb_G, основные параметры которого приведены в таблице:

ID Accession Type Seed Full Uniprot Average length Average %id Average coverage Change status Description
Acetone_carb_G PF08882 Domain 27 307 1790 106.8 36 48.03 Changed Acetone carboxylase gamma subunit

Для работы была выбрана двухдоменная архитектура A0A1H6FYJ5_THEAL Hydantoinase_B - Acetone_carb_G (в порядке следования доменов). Такая архитектура состоит из 38 последовательностей.

По данным pfam, длина профиля для данного домена равна 114.

In [1]:
import pandas as pd
import plotly.express as px
import plotly.io as pio

from Bio import SeqIO
from Bio.Seq import Seq
from collections import defaultdict
from sklearn.metrics import auc, f1_score

pio.renderers.default = 'iframe'

Для начала был получен список AC, содержащих двухдоменную архитектуру.

In [2]:
ac_2d = []

with open('2dac.txt') as records:
    for line in records:
        ac_2d.append(line.split()[0])

Далее были получены последовательности всех белков содержащих Acetone_carb_G, а после, уже из него, были получены последовательности двухдоменных белков. Как и ожидалось, вытащились 38 последовательностей.

In [3]:
full = SeqIO.parse('full.fasta', 'fasta')
records = [SeqIO.SeqRecord(seq = i.seq, id = i.id, description = '') for i in full if i.id in ac_2d]
SeqIO.write(records, '2dac.fasta', 'fasta')
Out[3]:
38

Далее посчитаем длину каждой последовательности в выборке full, а также получим список их AC.

In [4]:
len_full = [len(i) for i in SeqIO.parse('full.fasta', 'fasta')]
ac_full = [i.id for i in SeqIO.parse('full.fasta', 'fasta')]
In [5]:
with open('ac_full.txt', 'w') as out:
    for i in ac_full:
        print(i, file = out)

Результат выравнивания (muscle) и pезультат ревизии (jalview). Во время ревизии удалялись слишком длинные/слишком короткие последовательности, колонки, состоящие из гепов, а также небольшие участки до первого домена и участки после второго домена. После ревизии осталось 26 последовательностей. Список AC был также получен.

In [6]:
alig = SeqIO.parse('alig.fasta', 'fasta')
ac_alig = [i.id.split('/')[0] for i in alig]
In [7]:
with open('ac_alig.txt', 'w') as out:
    for i in ac_alig:
        print(i, file = out)

Теперь всё готово к постройке HMM. Команды ниже - это команды для создания профиля (его длина 702 а.о.), его калибровки, а также результата работы профиля по выборке full.

hmm2build profile alig.fasta
hmm2calibrate profile
hmmsearch --cpu=1 profile full.fasta &> hmmsearch_log.txt

Вычада hmmsearh была вручную переделана в нужную табличку.

In [8]:
hmm = pd.read_csv('hmm_table.csv')
hmm.head()
Out[8]:
e-value score bias seq
0 0.0 1393.9 20.7 A0A2H5VYN0_9BACT
1 0.0 1295.8 0.0 A0A562QB61_9BACI
2 0.0 1289.4 0.0 A0A2N5MCK9_9BACI
3 0.0 1288.8 33.1 A0A5B8UD28_9ACTN
4 0.0 1277.4 46.5 A0A5B8U2X4_9ACTN

Далее был получен список находок по профилю.

In [9]:
ac_hmm = list(hmm.seq)
In [10]:
with open('ac_hmm.txt', 'w') as out:
    for i in ac_hmm:
        print(i, file = out)

Теперь попробуем сделать табличку, в которой в качестве индексов были бы названия последовательностей, а в качестве колонок: содержит ли нужную архитектуру? - architecture; использовалась ли в выравнивании? - alig; нашлась ли в hmm? - hmm. Часть таблицы представлена ниже, но она получилась неслишком информативной.

In [11]:
def checker(ac_sm):
    result = {}
    
    for i in ac_full:
        if i in ac_sm:
            result[i] = 1
        else:
            result[i] = 0
    
    return result

dd = defaultdict(list)

for d in (checker(ac_2d), checker(ac_alig), checker(ac_hmm)):
    for key, value in d.items():
        dd[key].append(value)
        
df = pd.DataFrame(dd, index = ['architecture', 'alig', 'hmm']).transpose()
In [12]:
df.head()
Out[12]:
architecture alig hmm
A0A1I3KG61_9BACL 0 0 1
A0A133XH84_9RHOO 0 0 1
A0A0Q8KZJ6_9RHIZ 0 0 1
Q0RWI2_RHOJR 0 0 1
A0A4U0GLW4_9GAMM 0 0 1

Для получения хорошей информативности - объединим обе таблички.

In [13]:
merged = pd.merge(df, hmm, how = 'right', left_on = df.index, right_on = 'seq')
In [14]:
merged.to_csv('merged.csv')

Теперь таблица выглядит куда более информативной. С ней и будем работать далее.

In [15]:
merged
Out[15]:
architecture alig hmm e-value score bias seq
0 1 1 1 0.000000 1393.9 20.7 A0A2H5VYN0_9BACT
1 1 1 1 0.000000 1295.8 0.0 A0A562QB61_9BACI
2 1 1 1 0.000000 1289.4 0.0 A0A2N5MCK9_9BACI
3 1 1 1 0.000000 1288.8 33.1 A0A5B8UD28_9ACTN
4 1 1 1 0.000000 1277.4 46.5 A0A5B8U2X4_9ACTN
... ... ... ... ... ... ... ...
302 0 0 1 0.000012 16.0 0.1 A0A1M7FN25_9RHOB
303 0 0 1 0.000030 14.7 0.0 D4YQC9_9MICO
304 0 0 1 0.000340 11.2 0.2 A0A1R1SBK6_9ACTN
305 0 0 1 NaN 0.0 NaN A0A143B505_9DELT
306 0 0 1 NaN 0.0 NaN A0A511VFX0_9BACL

307 rows × 7 columns

Расчитает False Positive Rate = 1 - spec.

In [16]:
fpr = []
for i in range(len(merged.architecture)):
    
    try:
        zeros = merged.architecture[:i].value_counts()[0]
    except KeyError:
        zeros = 0
    
    fpr.append(zeros/269) # 269 - число последовательностей из full минус двухдоменные последовательности

Теперь посчитает True Positive Rate = sens.

In [17]:
tpr = []
for i in range(len(merged.architecture)):
    
    try:
        ones = merged.architecture[:i].value_counts()[1]
    except KeyError:
        ones = 0
    
    tpr.append(ones/38) # 38 - двухдоменные последовательности

Далее визуализируем всё получившееся - получим ROC кривую. В нашем случае получился очень сильно почти идеальный классификатор (AUC = 0.9963). Возможно, это связано с резким падением score. Для идеальной классификации используется вес 95.9.

(график интерактивный, при наведении мышкой на линию сверху высвечивается score, а снизу значения параметров)

In [18]:
fig = px.area(x = fpr, y = tpr,
    title = f'ROC кривая (AUC={auc(fpr, tpr):.4f})',
    labels = dict(x = '1-spec', y = 'sens'),
    width=700, height=500,
    hover_name = merged.score
)

fig.add_shape(
    type = 'line', line = dict(dash='dash'),
    x0 = 0, x1 = 1, y0 = 0, y1 = 1
)

fig.update_yaxes(scaleanchor="x", scaleratio=1)
fig.update_xaxes(constrain='domain')
fig.update_layout(margin=dict(l=0, r=0, t=50, b=20))
fig.show()

Также был произведен расчет F1-score:

In [19]:
f1_score(merged.hmm, merged.architecture)
Out[19]:
0.22028985507246376

Как и предполагалось выше, наблюдается действительно резкий скачок в в распределении весов.

In [20]:
fig = px.scatter(y = merged.score, hover_name = merged.seq)

fig.update_layout(
    margin=dict(l=10, r=10, t=10, b=10),
    showlegend = False,
    yaxis_title = 'score',
    xaxis_title = '',
    hovermode = 'closest'
)
fig.show()

Также была построена гистограмма распредления длин белков, входящих в выборку full, а также маленький боксплотик сверху.

Видно, что белки образует два кластера. Собственно, те, которые подлиннее, состоят из двудоменной архитектуры (> 300 а.к.), а те, которые короче (их ~ 90% от full), состоят только из одного домена.

In [21]:
fig = px.histogram(x = len_full, marginal = 'box', histnorm = 'percent')

fig.update_layout(
    margin=dict(l=10, r=10, t=10, b=10),
    showlegend = False,
    yaxis_title = '%',
    xaxis_title = 'длина',
    hovermode = 'closest',
    height=500
)

fig.show()

Итого, получился идеальный профиль, для предсказания исследумой двухдоменной структуры.

In [24]:
!jupyter nbconvert --to html hmm.ipynb

[NbConvertApp] Converting notebook hmm.ipynb to html
[NbConvertApp] Writing 626249 bytes to hmm.html