prac 4

Сначала подключим необходимые библиотеки для расчетов.

import matplotlib.pyplot as plt
import prody as pd

А теперь загрузим информацию о пространственном расположении атомов нашего белка - гипотетического активатора секреции в возбудителе бруциллеза. С помощью пакета Prody определим какой остаток будет иметь максимальное, а какой минимальное, среди всех остатков, среднее по всем его атомам значение B-фактора (mean_beta). Нам интересено узнать, в каком участке белка каждый из этих остатков расположен и какие остаки составляют их окружение. Также посмотрим, есть ли разброс в значениях B-фактора для атомов интересующих нас остатков.

myProtein=pd.parsePDB('7dnp') # инфа о координатах записалась в объекте myProtein

@> Connecting wwPDB FTP server RCSB PDB (USA).
@> 7dnp downloaded (C:\Users\kiril\...\7dnp.pdb.gz)
@> PDB download via FTP completed (1 downloaded, 0 failed).
@> 1431 atoms and 1 coordinate set(s) were parsed in 0.06s.

mean_betas = [] # заведем список, чтобы хранить в нем сведения о mean_beta каждого из остатков нашего белка 
for residue in myProtein.iterResidues():
    mean_beta = np.mean(residue.getBetas())
    mean_betas.append([residue, mean_beta])

# У нас всего 245 остатков в белке
a=0
for residue in myProtein.iterResidues():
    a+=1
print(a)

245

# отсортируем список mean_betas поубыванию значений, затем по их возрастанию и узнаем какие соответсвуют им остатки
print(sorted(mean_betas, key=lambda x:x[1])[0])

sorted(mean_betas, key=lambda x:x[1], reverse=True)[0]

[<Residue: PHE 89 from Chain A from 7dnp (11 atoms)>, 15.123636363636365]

[<Residue: SER 173 from Chain A from 7dnp (6 atoms)>, 46.975]

#Определим разброс значений B-фактора для атомов остатка PHE 89 (у него минимальная mean-beta)
PHE89=myProtein.select('resnum 89 and chain A')
SER173=myProtein.select('resnum 173 and chain A')

PHE_index,SER_index=[],[]

for i in PHE89:
    PHE_index.append(i.getIndex())

for i in SER173:
    SER_index.append(i.getIndex())
    
PHE_betas=PHE89.getBetas()
PHE_sigma=np.sum((PHE89.getBetas()-np.mean(PHE89.getBetas()))**2)/len(PHE89)

SER_betas=SER173.getBetas()
SER_sigma=np.sum((SER173.getBetas()-np.mean(SER173.getBetas()))**2)/len(SER173)

fig, axs = plt.subplots(1, 2, sharey=False)
axs[0].plot(PHE_index,
         PHE_betas,
         color='#000000',
         marker='o',
         mfc='#2ec0f9',
         mec='black',
         ls='--',
         label = 'PHE_89')
axs[1].plot(SER_index,
            SER_betas,
            color="black",
            marker='o',
            mfc='pink',
            mec='black',
            ls='--',
            label = 'SER_173')
axs[0].set_xlabel('resi')
axs[0].set_ylabel('B-factor')
axs[1].set_xlabel('resi')
axs[1].set_ylabel('B-factor')
axs[0].set_title("B-factor of PHE_89'atoms")
axs[1].set_title("B-factor of SER_173'atoms")
fig.tight_layout()
axs[0].set_xticks(np.arange(min(PHE_index),max(PHE_index)+1,1))
axs[0].set_yticks(np.arange(12,19+1,1))
axs[1].set_xticks(np.arange(min(SER_index),max(SER_index)+1, 1))
axs[1].set_yticks(np.arange(44,49+1,1))
axs[0].annotate(r'$\sigma=$'+str(PHE_sigma)[:5], xytext=(656, 18),xy=(656, 18),)
axs[1].annotate(r'$\sigma=$'+str(SER_sigma)[:5], xytext=(1294.5, 48),xy=(1294.5, 48),)

fig = matplotlib.pyplot.gcf()
fig.set_size_inches(12.5, 5.5)

from IPython.display import Image
Image(filename='7dnp.png')

Image(filename='phe89.png')

Image(filename='ser173.png')

Image(filename='b7dnp.png')

2 TASK. А сейчас представим на графике зависимость среднего B-factorа по атомам остатка нашего белка от расстояния между центром масс каждого остатка и центром масс всего белка.

protein_mass_center=pd.calcCenter(myProtein, weights=myProtein.getMasses())

distances = []
betas = []
for residue in myProtein.iterResidues():
    if 'CA' in residue.getNames():
        beta = np.mean(residue.getBetas())
        betas.append(beta)
        
        mass_center = pd.calcCenter(residue, weights=residue.getMasses())
        d = pd.calcDistance(protein_mass_center, mass_center)
        distances.append(d)
        

import numpy.polynomial.polynomial as poly
def scatter_hist(x, y, ax, ax_histx, ax_histy):
    # no labels
    ax_histx.tick_params(axis="x", labelbottom=False)
    ax_histy.tick_params(axis="y", labelleft=False)
    # the scatter plot:
    ax.scatter(x, y,fc='lightpink',ec='black')
   

    # now determine nice limits by hand:
    binwidth = 0.5
    xymax = max(np.max(np.abs(x)), np.max(np.abs(y)))
    lim = (int(xymax/binwidth) + 1) * binwidth
    
    #plot polyfit
    x_new = np.linspace(2, 26, num=len(distances)*10)
    coefs = poly.polyfit(distances, betas, 2)
    ffit = poly.polyval(x_new, coefs)
    ax.plot(x_new, ffit,color='black')
     
    bins = np.arange(-lim, lim + binwidth, binwidth)
    ax_histx.hist(x, bins=bins,facecolor='#FBB102',edgecolor='black')
    ax_histy.hist(y, bins=bins, orientation='horizontal',facecolor='lightgreen',edgecolor='black')
    
    
# definitions for the axes
left, width = 0.21, 0.65
bottom, height = 0.1, 0.65
spacing = 0.005


rect_scatter = [left, bottom, width, height]
rect_histx = [left, bottom + height + spacing, width, 0.2]
rect_histy = [left + width + spacing, bottom, 0.2, height]

# start with a square Figure
fig = plt.figure(figsize=(8, 8))

ax = fig.add_axes(rect_scatter)
ax_histx = fig.add_axes(rect_histx, sharex=ax)
ax_histy = fig.add_axes(rect_histy, sharey=ax)
ax.set_xlim(2,26)
ax.set_ylim(14,46)
ax.set_xlabel("Disatance between mass centers, A")
ax.set_ylabel('Beta factor')
# use the previously defined function
scatter_hist(distances, betas, ax, ax_histx, ax_histy)

plt.show()

3.Восстановление функции ЭП

Image(filename='full.png')

Image(filename='full2.png')

Image(filename='partial.png')