import numpy as np, pandas as pd, prody as prd

Prody - это библиотека, позволяющая быстро анализировать структуры биомолекул. Нам было предложено использовать Prody для краткого анализа B-фактора структуры 7avz, соответствующей тирозинкиназному домену Mer-киназы.

Задание 1

s1 = prd.parsePDB("7AVZ")
mean_betas = []
for residue in s1.iterResidues():
    if "CA" in residue.getNames(): 
        mean_beta = np.mean(residue.getBetas())
        mean_betas.append([residue, mean_beta])
betas = sorted(mean_betas, key = lambda x: x[1])

@> PDB file is found in working directory (7avz.pdb.gz).
@> 2214 atoms and 1 coordinate set(s) were parsed in 0.03s.

B-фактор - параметр атома, характеризующий неопределённость его электронной плотности. В файле PDB содержится информация о B-факторах атомов. Нам требуется найти аминокислотные остатки с самым низким и самым высоким средним B-фактором.

betas[0]

[<Residue: ALA 708 from Chain A from 7AVZ (5 atoms)>, 40.33]

Вот остаток с наименьшим B-фактором. Ala-708 располагается внутри глобулы в окружении гирофобных остатков, сложенных в альфа-спирали. Такая структура очень устойчива.

Благодаря плотной упаковке окружения и его гидрофобности остаток лишен не только свободы передвижения, но и доступа воды.

betas[-1]

[<Residue: SER 661 from Chain A from 7AVZ (6 atoms)>, 103.745]

Остаток Ser-661 обладает самым высоким средним B-фактором. Он располагается снаружи глобулы и смотрит в раствор, то есть полностью сольватирован. ЭП радикала этого остатка не отображается уже на уровне подрезки 1. Видно, что у соседей остатка ЭП радикалов тоже весьма слабая.

Интересно посмотреть на разброс значений B-факторов атомов этих остатков. Как показано ниже, у остатка с высоким B-фактором (Ser-661) очень высокое стандартное отклонение B-фактора. У атомов этого остатка B-фактор сильно отличается: самые низкие значения - у азота аминогруппы и кислорода CО-группы остова, а у остальных - достаточно высокие.

res_max = s1.select("resnum 661")
np.std(res_max.getBetas())

16.095086983300217

Напротив, у остатка с самым низким значением среднего B-фактора (Ala-708) разброс его значений небольшой.

res_min = s1.select("resnum 708")
np.std(res_min.getBetas())

2.5889457313740665

Задание 2

Нам было предложено построить график зависимости среднего B-фактора остатков от его расстояния до центра массы белковой структуры.

prot_mass_center = prd.calcCenter(s1, weights=s1.getMasses())

betas_new = []
distances = []
for residue in s1.iterResidues():
    if "CA" in residue.getNames():
        beta = np.mean(residue.getBetas())
        betas_new.append(beta)
        
        mass_center = prd.calcCenter(residue, weights=residue.getMasses())
        distance = prd.calcDistance(prot_mass_center, mass_center)               
        distances.append(distance)

import seaborn as sns, matplotlib.pyplot as plt
sns.set(rc={'figure.figsize':(9,6)}) 
plot = sns.scatterplot(x=distances, y=betas_new)
plot.set_xlabel('Расстояние до центра масс, А')
plot.set_ylabel('B-фактор')
plot.set_title("B-фактор = f(d(остаток, центр масс))")

Text(0.5, 1.0, 'B-фактор = f(d(остаток, центр масс))')

На графике выстраивается отчётливый тренд к росту B-фактора остатков с удалением от центра масс. Вероятно, есть смысл говорить о такой зависимости: чем дальше от центра масс, тем, если представить белковую глобулу в очень грубом приближении шаром, ближе к поверхности глобулы, контактирующей с раствором. Кроме того, что эти остатки могут принимать множественные положения из-за действия растворителя, они еще и не стабилизированы окружением. Кроме того, как мы выяснили в практикуме 2, у остаткаов с высоким B-фактором есть склонность входить в состав неструктурированных фрагментов. Однако не все так просто: некоторые остатки, образующие вторичную структуру, обладают аномально высоким B-фактором; в свою очередь и не все остатки, расположенные близко к центру масс, имеют низкие значения параметра. В связи с этим предполагать сферическую модель сферы было бы крайне неосмотрительно.

Ниже представлены коэффициенты корреляции Спирмана (мера ранговой зависимости) и Пирсона (мера линейной зависимости) для зависимости B-фактора от рассояния до центра масс. Коэффициент Пирсона в 0.6 даёт некоторый повод подумать о линейной зависимости, но по графику для меня лично это неочевидно.

from scipy.stats import spearmanr
spearmanr(distances, betas_new)

SpearmanrResult(correlation=0.6411401910293151, pvalue=1.9916967142226675e-32)

from scipy.stats import pearsonr
pearsonr(distances, betas_new)

(0.5985658106681595, 1.9319328149618057e-27)

Нам было интересно придумать другую модель для нашей молекулы (которая и не молекула никакая на самом деле, а просто вырезанный домен тирозин-киназы). Можно упрощенно представить нашу белковую структуру (одна цепь) как две условно одинаковые суб-глобулы, соединённые перешейком примерно в центре масс. B-фактор, характеризующий подвижность, конечно, зависит от расстояния до центра масс, однако есть что-то еще, что влияет на него. Нам показалось, что имеет смысл такая же зависимость B-фактора от расстояния до центра масс глобулы, который принадлежит остаток. Это разумно, так как в центре суб-глобулы на остаток действют такие же стабилизирующие факторы, как и в центре глобулы. Как мы видели, линейная зависимость B-фактора от расстояния до центра масс немного сомнительна. Еще более сомнительна линенйная зависимость от расстояния до центра глобулы. Мы разбили структуру на две УСЛОВНЫЕ суб-глобулы "на глаз".

dom1 = s1.select("resnum 562 to 677 731 to 742")
dom1_mass_center = prd.calcCenter(dom1, weights=dom1.getMasses())
r1 = list(range(562, 678))
r2 = list(range(731, 742))

betas_new_dom1 = []
distances_dom1 = []
for residue in s1.iterResidues():
    if "CA" in residue.getNames() and ((residue.getResnum() in r1) or(residue.getResnum() in r2)):
        beta = np.mean(residue.getBetas())
        betas_new_dom1.append(beta)
        
        mass_center = prd.calcCenter(residue, weights=residue.getMasses())
        distance = prd.calcDistance(dom1_mass_center, mass_center)               
        distances_dom1.append(distance)
        
        
dom2 = s1.select("resnum 678 to 730 763 to 862")
dom2_mass_center = prd.calcCenter(dom2, weights=dom2.getMasses())
r1 = list(range(678, 731))
r2 = list(range(763, 863))

betas_new_dom2 = []
distances_dom2 = []
for residue in s1.iterResidues():
    if "CA" in residue.getNames() and ((residue.getResnum() in r1) or(residue.getResnum() in r2)):
        beta = np.mean(residue.getBetas())
        betas_new_dom2.append(beta)
        
        mass_center = prd.calcCenter(residue, weights=residue.getMasses())
        distance = prd.calcDistance(dom2_mass_center, mass_center)               
        distances_dom2.append(distance)
        
sns.set(rc={'figure.figsize':(14,6)})       
fig, ax =plt.subplots(1,2)
plot = sns.scatterplot(x=distances_dom1, y=betas_new_dom1, ax=ax[0])
plot.set_xlabel('Расстояние до центра с/д 1')
plot.set_ylabel('B-фактор')
plot.set_title("B-фактор в остатках субдомена 1")
plot = sns.scatterplot(x=distances_dom2, y=betas_new_dom2, ax=ax[1])
plot.set_xlabel('Расстояние до центра с/д 2')
plot.set_ylabel('B-фактор')
plot.set_title("B-фактор в остатках субдомена 2")

Text(0.5, 1.0, 'B-фактор в остатках субдомена 2')

Имеет смысл в этом случае рассмотреть зависимость B-фактора от расстояния в какой-то степени. Кроме того, есть какой-то базовый уровень B-фактора из-за собственной подвижности белка как единицы. На небольших расстояниях мы можем считать расстояния до центра масс и центра суб-глобулы независимыми.

r1_1 = list(range(562, 678))
r2_1 = list(range(731, 743))
r1_2 = list(range(678, 731))
r2_2 = list(range(763, 863))

betas_all = []
ds_to_mc = []
ds_to_sdc = []
for residue in s1.iterResidues():
    if "CA" in residue.getNames():
        beta = np.mean(residue.getBetas())
        betas_all.append(beta)

        mass_center = prd.calcCenter(residue, weights=residue.getMasses())
        distance = prd.calcDistance(prot_mass_center, mass_center)               
        ds_to_mc.append(distance)

        if (residue.getResnum() in r1_1) or(residue.getResnum() in r2_1):
            distance1 = prd.calcDistance(dom1_mass_center, mass_center) 
        elif (residue.getResnum() in r1_2) or(residue.getResnum() in r2_2):
            distance1 = prd.calcDistance(dom2_mass_center, mass_center)
        else:
            distance1 = "NA"
            print(residue.getResnum())
        ds_to_sdc.append(distance1)

dom2_mass_center

array([ 6.66366584, 22.96348023, -3.94233646])

dom1_mass_center

array([ 27.16646747,  36.55479335, -12.95993441])

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

from scipy.optimize import curve_fit

def func(X, a, b, c, n1, n2):
    d1,d2 = X
    return a + b*np.power(d1, n1) + c*np.power(d2, n2)

d1 = np.array(ds_to_mc)
d2 = np.array(ds_to_sdc)

p0 = 40., .05, .05, 1, 1
popt, pcov = curve_fit(func, (d1,d2), betas_all, p0, maxfev=1200)
popt

array([34.31995044,  0.07013715,  1.09585979,  1.70806096,  0.99494132])

perr = np.sqrt(np.diag(pcov))
perr

array([6.81991929, 0.13535577, 1.80637028, 0.5472898 , 0.46368291])

По нашей модели B-фактор зависит от расстояния до центра масс все-таки нелинейно (в степени 1.7), а от расстояния до центра суб-глобулы - линейно (степень 1). Однако доверять такой модели, наверное, стоит с осторожностью, так как стандартные ошибки параметров очень большие.

sns.set(rc={'figure.figsize':(14,6)})       
fig, ax =plt.subplots(1,2)
plot = sns.scatterplot(x=distances, y=betas_new, ax=ax[0])
plot.set_xlabel('Расстояние до центра масс')
plot.set_ylabel('B-фактор')
plot.set_title("Данные")
plot = sns.scatterplot(x=d1, y=func((d1,d2), *popt), ax=ax[1])
plot.set_xlabel('Расстояние до центра')
plot.set_ylabel('B-фактор')
plot.set_ylim(38, 105)
plot.set_title("Модель")

Text(0.5, 1.0, 'Модель')

Модель пыталась объяснить дисперсию при средних значениях расстояний до центра масс, но у нее все равно не получилось. Какие могут быть проблемы?

1. Некорректное разделение на глобулы - более того, размеры этих глобул не равны и это не было учтено.

2. Можно было бы учесть локальную плотность остатков при фиксированных расстояниях d1 и d2 (кольцо, плоскость которого параллельна плоскости, разделяющей глобулы).

Задание 3

Предлагалось придумать одномерную молекулярную систему и смоделировать эксперимент по расшифровке этой структуры. Необходимо было получить реалистичную функцию реальной электронной плотности системы. Смоделировать ЭП, получаемую в ходе эксперимента можно, разложив функцию на большое количество гармоник, и выбрав из них какое-то ограниченное их количество. Из фаз и амплитуд этих гармоник получается "экспериментально полученная" функция ЭП.

Я решила сгенерировать функцию для системы "O2 + N3-". В азид-анионе заряд делокализован: электроны предпочитают находиться в области крайних атомов иона.

Функция ЭП была сгенерирована командой "python compile-func.py -g 8,2.9,6+8,2.9,7.21+7.9,3,17+6.2,3,18.15+7.9,3,19.3".

Разложение функции ЭП на гармоники было сделано командой "python func2fourier.py -i func_ox_and_az.txt" с ситуативным добалением параметров -F и -P - шум амплитуд и фаз.

Нужные гармоники выбирались вручную копипастой.

Получилось так:

Далее были смоделированы эксперименты. Интервалы - это гармоники. Если гармоники не указаны, то интервалом считать 0-25.

summ = pd.read_csv("bioinf_pr4/summary.csv", delimiter=";")
summ

	Harmonics set	Resolution (A)	Data integrity (%)	F noise (% of phi)	Phase noise (% of phi)	Quality
0	0_10	3.0 A	100%	0	0	bad
1	0_20	1.5 A	100%	0	0	medium
2	0_25	1.2 A	100%	0	0	perfect
3	0_30	1.0 A	100%	0	0	perfect
4	0_25	1.2 A	100%	10	10	perfect
5	0_25	1.2 A	100%	0	20	good
6	0_25	1.2 A	100%	20	0	perfect
7	1_10, 15_25	1.2 A	80%	10	10	perfect
8	4_25	1.2 A	84%	10	10	medium
9	3_7, 10_25	1.2 A	80%	10	10	good
10	2_7, 10_17, 20_25	1.2 A	76%	10	10	perfect

Какие выводы можно сформулировать из выше сделанного?

1. При более высоком разрешении качество структуры лучше. (опыты 0 - 3)

2. Шум фазы гораздо сильнее затрудняет расшифровку ЭП, чем шум амплитуды. (опыты 4 - 6)

3. Необязательно бОльшая полнота данных означает лучшее качество: гораздо важнее присутствие первых гармоник (самого низкого разрешения). (опыты 7 - 10)