import numpy as np
import prody as pd
import pandas
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter

Загрузим данные¶

Предварительно структуры были подрезаны с концов, чтобы наблюдалось соответствие длин (работали с биологическими сборками)

opened = pd.parsePDB('/home/julybel/study/3d/prak-07/task3/2fff_clean_v2.pdb')
opened = pd.parseDSSP('/home/julybel/study/3d/prak-07/task3/2fff.dssp', opened)
opened_sol_acc = opened.getData('dssp_acc')

@> 3492 atoms and 1 coordinate set(s) were parsed in 0.07s.

closed = pd.parsePDB('/home/julybel/study/3d/prak-07/task3/2bg1_clean_v2.pdb')
closed = pd.parseDSSP('/home/julybel/study/3d/prak-07/task3/2bg1.dssp', closed)
closed_sol_acc = closed.getData('dssp_acc')

@> 3488 atoms and 1 coordinate set(s) were parsed in 0.06s.

Cопоставление "номер остатка: экспонированность"¶

sp_open = []
for atom, acc in zip(opened, opened_sol_acc):
    if atom.getName() == 'N':
        sp_open.append((atom.getResname() + '_' + str(atom.getResnum()), acc))
    
sp_close = []
for atom, acc in zip(closed, closed_sol_acc):
    if atom.getName() == 'N':
        sp_close.append((atom.getResname() + '_' + str(atom.getResnum()), acc))

Проверим, есть ли такие остатки, типы которых для определенной позиции не совпадают:

print(*[(i, j) for i, j in zip(sp_open, sp_close) if i[0] != j[0]], sep='\n')

(('ARG_686', 159.0), ('GLN_686', 126.0))
(('ARG_687', 154.0), ('GLN_687', 90.0))

Оказалось, что есть - это аминокислоты в позициях 686 и 687. Их изменения в экспонированности не будут учитываться при соответствующем анализе.

x, y = [], []

for i, j in zip(sp_open, sp_close):
    x.append(int(i[0].split('_')[1]))
    y.append(i[1] - j[1])
    
fig, ax = plt.subplots(figsize=(15,5))
ax.set_title("Differences in residues exposure\n", fontdict={'fontsize' : 18})
ax.set_xlabel('Residue number', fontdict={'fontsize' : 16})
ax.set_ylabel('Exposure(open) - Exposure(close)', fontdict={'fontsize' : 16})
plt.plot(x, y, color='forestgreen')

[<matplotlib.lines.Line2D at 0x7f177c27fb10>]

x, y = [], []

for i, j in zip(sp_open, sp_close):
    x.append(int(i[0].split('_')[1]))
    y.append(abs(i[1] - j[1]))
    
fig, ax = plt.subplots(figsize=(15,5))
ax.set_title("Abs differences in residues exposure\n", fontdict={'fontsize' : 18})
ax.set_xlabel('Residue number', fontdict={'fontsize' : 16})
ax.set_ylabel('Exposure(open) - Exposure(close)', fontdict={'fontsize' : 16})
plt.plot(x, y, color='forestgreen')

[<matplotlib.lines.Line2D at 0x7f177ba099d0>]

print('Residue with max change in exposure:')
print(*[(i, j) for i, j in zip(x, y) if j == max(y)])

Residue with max change in exposure:
(361, 95.0)

Мы видим, что наибольшее изменение в экспонированности (по модулю) наблюдается для остатка 361 - Val.

NB1¶

Приведем эмпирические значения SASA из литературы:
https://sci-hub.do/10.1371/journal.pone.0080635 (Table 1, p. 3)

dat = pandas.read_csv('/home/julybel/study/3d/prak-07/task3/table.csv', sep='\t', header=None)[[0,3]].set_index(0).T
dat

reopen = []
for atom, acc in zip(opened, opened_sol_acc):
    if atom.getName() == 'N':
        #recalc = round(acc*10/float(dat[atom.getResname()]), 2)
        recalc = round(acc/float(dat[atom.getResname()]), 2)
        reopen.append((atom.getResname() + '_' + str(atom.getResnum()), recalc))
    
reclose = []
for atom, acc in zip(closed, closed_sol_acc): 
    if atom.getName() == 'N':
        #recalc = round(acc*10/float(dat[atom.getResname()]), 2)
        recalc = round(acc/float(dat[atom.getResname()]), 2)
        reclose.append((atom.getResname() + '_' + str(atom.getResnum()), recalc))

x, y = [], []

xo, yo = [], []
xc, yc = [], []

for i, j in zip(reopen, reclose):
    x.append(int(i[0].split('_')[1]))
    y.append(abs(i[1] - j[1]))
    
    xo.append(int(i[0].split('_')[1]))
    yo.append(i[1])
    
    xc.append(int(j[0].split('_')[1]))
    yc.append(j[1])
    
fig, ax = plt.subplots(figsize=(15,5))
ax.set_title("Abs differences in recalculated residues exposure\n", fontdict={'fontsize' : 18})
ax.set_xlabel('Residue number', fontdict={'fontsize' : 16})
ax.set_ylabel('Exposure(open) - Exposure(close)', fontdict={'fontsize' : 16})
plt.plot(x, y, color='forestgreen')

[<matplotlib.lines.Line2D at 0x7f177b9a6b50>]

x, y = [], []

for ox, oy, cx, cy in zip(xo, yo, xc, yc):
    x.append(ox)
    #print(oy, '\t', cy, '\t', round(oy/cy, 2))
    y.append(round(oy/cy, 2))

fig, ax = plt.subplots(figsize=(15,5))
ax.set_title("Ratio of exposures\n", fontdict={'fontsize' : 18})
ax.set_xlabel('Residue number', fontdict={'fontsize' : 16})
ax.set_ylabel('Ratio (open/close)', fontdict={'fontsize' : 16})
plt.plot(x, y, color='forestgreen')

/home/julybel/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in double_scalars
  
/home/julybel/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in double_scalars

[<matplotlib.lines.Line2D at 0x7f177967e910>]

Как мы видим, при выборе отношения экспонированностей открытой формы к закрытой возникают значения nan (прерывания в графике), что может объясняться делением на 0 (на самом деле, делением 0 на 0, если посмотреть численные значения). Соответственно, в таком случае, при использовании отношения, а не разности, мы не можем дать оценку изменению экспонированности остатка.

При делении числа на 0, например, 0.01 на 0.0 возникает значение inf, что также не дает корректную оценку интересующей нас величины

Вероятно, корректнее было бы использовать разность величин (модуль разности) для оценки изменения экспонированности остатка)

print('Residue with max change in exposure:')
print(*[(i, j) for i, j in zip(x, y) if j == max(y)])

Residue with max change in exposure:
(361, 0.5900000000000001)

Максимальную по модулю экспонированность все еще имеет 361 остаток (валин).
В целом, можно сказать, что значительных изменений в разнице экспонированности остатков открытой и закрытой форм при использовании для расчета двух разных формул мы не увидели