import prody
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
from IPython.display import Image
struc = prody.parsePDB('3q3y')
mean_betas = [(r, np.mean(r.getBetas())) for r in struc.iterResidues() if 'CA' in r.getNames()]
len(mean_betas)
mean_betas.sort(key=lambda x: x[1])
print(mean_betas[0], mean_betas[-1], sep='\n')
r_min, r_max = mean_betas[0][0], mean_betas[-1][0]
for r in r_min, r_max:
betas = r.getBetas()
print(r.getResname(),
f'std: {round(np.var(betas), 6)}',
f'var: {round(np.std(betas), 6)}\n{"-"*10}', sep='\n')
Таким образом, остаток с наименьшим средним B-фактором -- аланин цепи А, с наибольшим - знакомый нам по второму практикуму лизин 0.
Разброс значений B-фактора в случае лизина намного больше, чем аналогичный для аланина.
В выводе выше представлены стандартное отклонение и дисперсия для аланина и лизина соответственно.
Image('1_1.png')
Рисунок 1: Структура 3Q3Y. Остатки LYS-0 и ALA-172 с наибольшим и наименьшим средними значениями B-фактора. Окрашены в зеленый и розовый соответственно.
Как видно из рисунка 1, лизину предположительно может соответствовать высокая подпижность, так как он "выходит на поверхность" своим боковым радикалом, а вот аланин находится внутри глобулы, таким образом и подвижность мы должны наблюдать минимальную.
Локализация остатков относительно поверхности лучше видна на рисунке 2.
Image('1_2.png')
Рисунок 2: Структура 3Q3Y. Остатки LYS-0 и ALA-172 с наибольшим и наименьшим средними значениями B-фактора. Окрашены в зеленый и розовый соответственно. Добавлена surface для наглядности расположения относительно поверхности.
Боковой радикал Lys-0 стабилизирован лишь потенциальными водородными связями с плохой геометрией (рисунок 3), а Ala-172 находится в бета складчатом слое и связан с соседними листами водородными связями (рисунок 4). Образование водородных связей для аланина через боковой радикал не характерно.
Image('1_3.png')
Рисунок 3: Остаток LYS-0 и его окружение.
Image('1_4.png')
Рисунок 4 Остаток ALA-172 и его окружение.
Prody и B-факторы часть 2¶
protein_mass_center = prody.calcCenter(struc, weights=struc.getMasses())
mean_betas, distances = [], []
for r in struc.iterResidues():
if 'CA' in r.getNames():
mean_betas.append(np.mean(r.getBetas()))
mass_center = prody.calcCenter(r, weights=r.getMasses())
distances.append(prody.calcDistance(protein_mass_center, mass_center))
fig, ax = plt.subplots(figsize=(12, 8))
sns.scatterplot(distances, mean_betas, color='darkblue')
ax.set_title('Зависимость B-фатора от расстояния между центрами масс', )
ax.set_xlabel('Расстояния между центрацми масс остатка и белка, Å')
ax.set_ylabel('B-фактор')
fig.show()
import statsmodels.formula.api as smf
model1 = smf.ols(formula = 'B_factor ~ Distances',
data=pd.DataFrame({'Distances':distances, 'B_factor':mean_betas}))
print(model1.fit().summary())
Предположим, что регрессор в виде расстояния между центрами масс линейно описиывает зависимую переменную, то бишь B-фактор.
Как видно из модели выше, линейно зависимость не описывается, R-squared никакой и т.д.
Единственное, что интересное, - коэффициент Дурбина-Уотсона. Если он стремится к 0, то мы отвергаем нулевую гипотезу о независимости отклонений. Здесь значение где-то посередине, но автокорреляция отклонений B-факторов в рамках белка звучит крайне логично, так как есть у нас идёт, например, череда бета слоёв, то подвижность в рамках них будет минимальна. Получается, автокорреляция.
from scipy.stats import spearmanr, pearsonr
print(spearmanr(mean_betas, distances), pearsonr(mean_betas, distances), sep='\n')
Ни монотонная, ни линейная корреляция также не обнаружены.
Если о монотонной можно как-то аккуратно поговорить, то вот линейной точно не видно, в чём мы уже убедились выше, потерпев фиаско в построении модели.
Задание 3. Как работает восстановление функции электронной плотности по экспериментальным данным¶
run compile-func.py -g 30,3,5+24,3,4+40,2.5,12+30,4,10.5+25,3,11.2
Рисунок 4: Модель электронной плотности одномерной системы атомов из двух молекул, состоящих из 2 и 3 атомов соответственно.
Отныне и впредь все несовместимости исходных скриптов с python3 будут фикситься, потому что добавлять новый интерпретатор в kernel сложнее, чем добавить скобки к print.
С помощью скрипта получаем файл с амплитудами и фазами¶
run func2fourier.py -i func.txt -o out.txt
Image('screen1.png')
Мною были выбраны следующие критерии для классификации восстановления по качеству:
- Плохое: возможно только угадывать со свечкой в руках
- Среднее: определение положения атомов возможно частично
- Хорошее: возможно предположить положение всех максимумов, зная число атомов, но максимумы при этом не отличимы от шума
- Отличное: можно определить местоположение всех атомов
Сначала рассмотрим несколько полных наборов гармоник без добавления шума.
Для интервалов взял степени двойки: 0-2, 0-4, 0-8, 0-16, 0-32, 0-64, 0-128.
Результаты дальнейшего моделирования занесены в таблицу.
run fourier-filter.py -i out.txt -r 0-2 -o 11.txt
run fourier2func.py -f func.txt -i 11.txt -o 12.txt
Рисунок 5: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-2).
Качество электронной плотности в данном случае однозначно плохое, так как понятно, что ничего непонятно. Можно лишь предположить, что где-то в этой области были атомы.
run fourier-filter.py -i out.txt -r 0-4 -o 21.txt
run fourier2func.py -f func.txt -i 21.txt -o 22.txt
Рисунок 6: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-4).
Пока ещё всё плохо, но уже можно вполне уверенно предполагать, что мы имеем дело с двумя молекулами. Количество атомов остаётся неизвестным.
run fourier-filter.py -i out.txt -r 0-8 -o 31.txt
run fourier2func.py -f func.txt -i 31.txt -o 32.txt
Рисунок 7: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-8).
Стало лучше, но стоит попробовать ещё раз увеличить набор.
run fourier-filter.py -i out.txt -r 0-16 -o 41.txt
run fourier2func.py -f func.txt -i 41.txt -o 42.txt
Рисунок 8: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-16).
Стоит ещё раз увеличить набор.
run fourier-filter.py -i out.txt -r 0-32 -o 51.txt
run fourier2func.py -f func.txt -i 51.txt -o 52.txt
Рисунок 9: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-32).
Теперь мы уже можем предположить, где находятся максимумы, а для первой молекулы всё видно однозначно. Хорошее восстановление.
run fourier-filter.py -i out.txt -r 0-64 -o 61.txt
run fourier2func.py -f func.txt -i 61.txt -o 62.txt
Рисунок 10: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64).
Первая восстановление, которое я классифицирую как отличное.
run fourier-filter.py -i out.txt -r 0-128 -o 71.txt
run fourier2func.py -f func.txt -i 71.txt -o 72.txt
Рисунок 11: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-128).
Видно, что теперь увеличение числа гармоник не сказывается на качестве восстановления ЭП и мы всё так же точно определяем максимумы гауссовых слагаемых.
На основании проведенных экспериментов можно сделать вывод, что качество восстановления ЭП увеличивается с длиной интервала гармоник, что логично, так как мы с большей точностью способы апроксимировать нашу функцию тригонометрией.
Теперь попробуем восстановление ЭП с дополнительным шумом амплитуды и/или фазы в размере 25, 50, 75%.
Буду производить восстановление по полному набору гармоник на интервале 0-64.
run func2fourier.py -i func.txt -o F11.txt -F 25
run fourier-filter.py -i F11.txt -r 0-64 -o F12.txt
run fourier2func.py -f func.txt -i F12.txt -o F13.txt
Рисунок 12: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум амплитуды равен 25%.
Даже при 25% уровне зашумления качество восстановления всё ещё можно назвать отличным.
run func2fourier.py -i func.txt -o F21.txt -F 50
run fourier-filter.py -i F21.txt -r 0-64 -o F22.txt
run fourier2func.py -f func.txt -i F22.txt -o F23.txt
Рисунок 13: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум амплитуды равен 50%.
При 50% уровне зашумления всё еще можно точно определить положения всех атомов, а значит модель всё еще отличная, хоть и, действительно, шумненькая.
run func2fourier.py -i func.txt -o F31.txt -F 75
run fourier-filter.py -i F31.txt -r 0-64 -o F32.txt
run fourier2func.py -f func.txt -i F32.txt -o F32.txt
Рисунок 14: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум амплитуды равен 75%.
Положения исходных атомов всё еще можно точно определить, но уровень шума настолько велик, что можно предположить и существования других атомов.
Если знать состав, то можно с большой долей вероятности угадать правильную локализацию всех атомов. А значит и модель хорошая.
Думаю, восстановление получается неплохим из-за большого интервала гармоник. Уменьшим его до 48.
run func2fourier.py -i func.txt -o F41.txt -F 75
run fourier-filter.py -i F41.txt -r 0-48 -o F42.txt
run fourier2func.py -f func.txt -i F42.txt -o F43.txt
Рисунок 15: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-48), шум амплитуды равен 75%.
Получаем хорошее качество, но хуже, чем в предыдущем эксперименте
run func2fourier.py -i func.txt -o P11.txt -P 25
run fourier-filter.py -i P11.txt -r 0-64 -o P12.txt
run fourier2func.py -f func.txt -i P12.txt -o P13.txt
Рисунок 16: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум фазы равен 25%.
Отличная модель.
run func2fourier.py -i func.txt -o P21.txt -P 50
run fourier-filter.py -i P21.txt -r 0-64 -o P22.txt
run fourier2func.py -f func.txt -i P22.txt -o P23.txt
Рисунок 17: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум фазы равен 50%.
Средняя модель.
run func2fourier.py -i func.txt -o P31.txt -P 75
run fourier-filter.py -i P31.txt -r 0-64 -o P32.txt
run fourier2func.py -f func.txt -i P32.txt -o P33.txt
Рисунок 18: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум фазы равен 75%.
Средняя, но намного более зашумленная модель. По ней, мне кажется, нельзя уверенно что-либо сказать, так что качество восстановления очень близко к плохому.
Пока можно предположить, что зашумленность фазы критичнее для восстановления электронной плотности, чем шум амплитуды.
Попробуем теперь скомбинировать шумы.
run func2fourier.py -i func.txt -o T11.txt -P 10 -F 40
run fourier-filter.py -i T11.txt -r 0-64 -o T12.txt
run fourier2func.py -f func.txt -i T12.txt -o T13.txt
Рисунок 19: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум фазы равен 10%, амплитуды -- 40%.
Восстановление можно считать отличным. Из рисунка можно делать однозначные выводы.
run func2fourier.py -i func.txt -o T21.txt -P 40 -F 10
run fourier-filter.py -i T21.txt -r 0-64 -o T22.txt
run fourier2func.py -f func.txt -i T22.txt -o T23.txt
Рисунок 20: Восстановление функции электронной плотности по полному набору гармоник (интервал гармоник: 0-64), шум фазы равен 40%, амплитуды -- 10%.
Восстановление можно считать хорошим. При этом кажется, что теперь шум фазы уже менее критичен, потому что восстановление выглядит более похожим на истину.
Теперь попробуем восстановить ЭП без шумов, но по неполным наборам гармоник: 8-56, 6-64, 0-50.
run fourier-filter.py -i out.txt -r 8-56 -o N11.txt
run fourier2func.py -f func.txt -i N11.txt -o N12.txt
Рисунок 21: Восстановление функции электронной плотности по неполному набору гармоник (интервал гармоник: 8-56).
Модель средняя. Потеря начала и конца критична.
run fourier-filter.py -i out.txt -r 6-64 -o N21.txt
run fourier2func.py -f func.txt -i N21.txt -o N22.txt
Рисунок 22: Восстановление функции электронной плотности по неполному набору гармоник (интервал гармоник: 6-64).
Большая беда с baseline, но картинку всё ещё можно успешно интерпретировать. Оценим модель как хорошую.
run fourier-filter.py -i out.txt -r 0-50 -o N31.txt
run fourier2func.py -f func.txt -i N31.txt -o N32.txt
Рисунок 23: Восстановление функции электронной плотности по неполному набору гармоник (интервал гармоник: 0-50).
Получаем отличную картинку.
ps. Я часто классифицировал модели как отилчные в виду отсутствия некоторого атома с маленьким пиком. В нашем случае, например, водорода, а то, что я моделировал, больше похоже на неорганику.