prac 4 Python для структур и не только

ЗАДАНИЕ1.

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

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]

PHE 89 from Chain A имеет наименьшее среднее по атомам значение b-фактора. Оно равно 15.12. SER 173 from Chain A имеет наибольшее значение b-фактора. Оно равно 46.975.

from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')

Определим разброс значений B-фактора для атомов остатка PHE 89 (у него минимальная mean-beta-factor) и SER 173. Разброс для каждого остатка не превышает 5 единиц B-factor.

#Определим разброс значений 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)

В контексте всей структуры белка, PHE89 находится между альфа-спиралями (см. рисунок ниже), в ядре белка. SER173 находится на поверхности белка и вероятно за счет полярного гидроксила образует водородные связи с окружением, увеличивая таким образом свою подвижность - а это причина иметь высокий B-factor.

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

В окружении PHE89 вся вода связана, поэтому она не вызывает колебаний этого остатка. Напротив, SER173 свободен для воды и может вступать с ней в связи и колебаться.

Image(filename='phe89.png')

Image(filename='ser173.png')

Приводим еще раз изображение белка, но теперь окрашенного по значению B-factor его остатков. Тут видно что серин лежит в петле - в самом подвижном из всех возможных элементов структуры.

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)
        

Зависимости строгой мы не наблюдаем: остатки, которые расположены относительно далеко от центра масс белка, имеют всевозможные значения b-factor. Это значит, что есть другие факторы, нежели удаленность от центра, которые более прямо определяют B-фактор. Одним из них может быть тип вторичной структуры, на котором сидит рассматриваемый остаток. Если остаток будет в петле, то его B-фактор будет больше, чем у того, кто сидит на альфа-спирали.

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.Восстановление функции ЭП

Посмотрим на умозрительной одномерной модели, как потеря данных, полученных из эксперимента РСА, 

снижает качество восстановленной по ним функции ЭП и усложняет процесс ее интерпретации. Наша модель состоит из 2х молекул, в одной 3 атома, в другой 2. В идеальном случае мы бы восстанавливали функцию ЭП по 3000 рефлексам (точнее по их модулям и фазам структурных факторов). Но в ходе эксперимента мы не можем избавиться от ошибок в значениях модуля(F) и фазы (P) или от их потери, и имеем дело с меньшим числом рефлексов. На картинке ниже приведена серия восстановления функции ЭП по модулям и фазам, чьи значения отличаются от истинного не более чем на 5, 15, 30 процентов. Непрерывная желтая кривая здесь и далее соответсвует истинной функции. Темный пунктир есть экспериментально полученная функция. Вы можете видеть, что фазовый шум более критичен для качества функции, чем модульный шум. Удивительно, но шум от фазы гасится шумом модуля. По-моему это обстоятельство не было оставлено без внимания при реализации комбинированного синтеза Фурье.

Image(filename='full.png')

Мы рассматриваем простую модель, что нам хватает 40+1(если нулевую тоже учитывать) гармоник, чтобы успешно восстановить функцию ЭП. Разрешение равно (30)/(40+1)=0.732 A. Надо сказать, что уже 10+1 гармоник достаточно, чтобы назвать восстановлением среднего качества. Хорошим оно будет при 20+1 использованных гармониках.

Image(filename='full2.png')

Узнаем, как восстановится функция по неполному набор гармоник. 
Левая картинка получилась по набору, где не хватало

первых 2 гармоник. Сигналы видны хорошо, а широкий колокол на 20 A нам не мешает интерпретировать функцию. Средняя картинка получится, если из полного набора из сорока гармоник пропадут 10% гармоник из середины. Теперь на 20A мы видим мелкие колоколы, которые можно спутать с сигналом. Если знать число молекул и число атомов в них, то восстановим мы функцию с хорошей оценкой, но без этого знания можно допустить ошибку. Если к этому неполному набору добавить 50ую гармонику, то получим правую картинку. Она похожа на среднюю, значит не имеют большого значения в восстановлении гармоники с номерами, большими номера, при котором восстановление по полному набору имеет отличное качество.

Image(filename='partial.png')

Приведем для рассмотренных примеров таблицу, где суммируем наши результаты.

Набор гармоник	Разрешение (А)	Полнота данных (%)	Шум амплитуды (% от величины F)	Шум фазы (% от величины P)	Качество восстановления
Полный набор гармоник
0-40	0.732	100	5	0	отл
0-40	0.732	100	0	5	отл
0-40	0.732	100	5	5	отл
0-40	0.732	100	15	0	отл
0-40	0.732	100	0	15	отл
0-40	0.732	100	15	15	отл
0-40	0.732	100	30	0	отл
0-40	0.732	100	0	30	удовл
0-40	0.732	100	30	30	отл
0	30	100	0	0	неуд
0-1	15	100	0	0	неуд
0-5	5	100	0	0	неуд
0-10	2.7	100	0	0	удовл
0-20	1.43	100	0	0	хор
0-40	0.73	100	0	0	отл
Неполный набор гармоник
2-40	0.73	95.12	0	0	отл
0-18&23-40	0.73	90.24	0	0	хор
0-18&23-40&50	0.73	90.24	0	0	хор

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