После работы на сервере были получены следующие файлы:
dna_pbc_1.pdb
rms_1.xvg
rms_2.xvg
sas_dna.xvg
hbond_dna.xvg
hbond_dna_sol.xvg
Затем они были скачаны на домашний компьютер и дальнейший анализ проводился на нем.
Прежде всего файл с фрагментом ДНК был модифицирован для возможности дальнейшей работы. Для последующего моделирования (задача №1647686) была проведена оптимизация геометрии системы для удаления "плохих" контактов в молекуле. Начальное и конечное значение максимальной силы равно соответственно: 3045 для атома 50 и 593 для атома 32.
Также необходимо было добавить в систему воду и произвести утряску. После утряски вода располагается более равномерно и хаотично (справа), мы не наблюдаем никаких структур.
display(Image('/Users/anna/Desktop/with_water.png'))
Визуальный анализ движения молекулы ДНК показывает, что примерно на 5 кадре происходит переход в В-форму (справа): можно заметить, как пары азотистых оснований образуют параллельные "стопки", характерные для В-формы, меняется изгиб цепи.
Примерное время - 800 пс
(Улучшение визуализации не сильно помогло, все равно удобнее просматривать по кадрам)
display(Image('/Users/anna/Desktop/transformation.png'))
def video(fname, mimetype):
from IPython.display import HTML
video_encoded = open(fname, "rb").read().encode("base64")
video_tag = '<video controls alt="test" src="data:video/{0};base64,{1}">'.format(mimetype, video_encoded)
return HTML(data=video_tag)
video("/Users/anna/Downloads/pymol.mov", "mp4")
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image, FileLink, FileLinks
Изменение среднеквадратичного отклонения от начального положения в зависимости от времени.
x, y = [], []
with open('/Users/anna/Downloads/struct/rms_1.txt') as f:
for line in f:
cols = line.split()
x.append(float(cols[0]))
y.append(float(cols[1]))
plt.plot(x, y, 'b.')
plt.xlabel('time')
plt.ylabel('RMSD')
plt.show()
Изменение среднеквадратичного отклонения относительно каждой предыдущей структуры на растоянии 400 кадров в зависимости от времени. Как мы видим, графики очень похожи для обоих случаев. Отклонение к концу не уменьшается, видимо, конформационный переход еще не завершен полностью.
x, y = [], []
with open('/Users/anna/Downloads/struct/rms_2.txt') as f:
for line in f:
cols = line.split()
x.append(float(cols[0]))
y.append(float(cols[1]))
plt.plot(x, y, 'b.')
plt.xlabel('time')
plt.ylabel('RMSD')
plt.show()
Зависимость изменения доступных растворителю гидрофобной и гидрофильной поверхностей от времени.
x, y1, y2 = [], [], []
with open('/Users/anna/Downloads/struct/sas_dna.txt') as f:
for line in f:
cols = line.split()
x.append(float(cols[0]))
y1.append(float(cols[1]))
y2.append(float(cols[2]))
plt.plot(x, y1,label='hydrophobic')
plt.plot(x,y2,label='hydrophilic')
plt.xlabel('time')
plt.legend()
fitfunc = lambda p, x: p[0]*pow(p[1]-x,1) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1,1,0] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x[100:], y1[100:]))
p2, success = optimize.leastsq(errfunc, p0[:], args=(x[100:], y2[100:]))
plt.plot( x,fitfunc(p1,x),"r-",c='blue',alpha=0.5)
plt.plot( x,fitfunc(p2,x),"r-",c='red',alpha=0.5)
plt.show()
В целом гидрофобная и гидрофильная поверхности, доступные растворителю, изменяются в некоторых пределах на всем времени симуляции, но в среднем остаются примерно постоянными. Возможно, это связано с малым размером участка ДНК, т.е. на таком масштабе изменение не видно. Рассмотрим подробнее начало: здесь наблюдается локальное увеличение гидрофильной поверхности и уменьшение гидрофобной. Вспомним, что конформационный переход приходится примерно на 800. Это хорошо соотносится с изменением гидрофильной/гидрофобной поверхности. Т.о. возможно, переход из А-формы в B происходит из-за того, что уменьшение свободной энергии соответствует уменьшению гидрофобной поверхности, доступной растворителю.
Водородные связи внутри ДНК. Их количество меняется на протяжении моделирования, преобладают значения 13-14-15. Т.е. в среднем 14, что хорошо, поскольку совпадает со стандартным количеством водородных связей для утсон-криковских пар в нашей молекуле.
x, y1 = [], []
with open('/Users/anna/Downloads/struct/hbond_dna.txt') as f:
for line in f:
cols = line.split()
x.append(float(cols[0]))
y1.append(float(cols[1]))
plt.plot(x, y1,'ro')
plt.xlabel('time')
plt.ylabel('# hhbond DNA-DNA')
plt.show()
Зависимость количества водородных связей между водой и ДНК не была построена, так как файл выдавал нули либо скопление связей в одной временной точке, изменение параметра cut-off не помогло. Видимо, где-то закралась ошибка!
Известно, что А-формы ДНК не бывает в воде (но бывает двухцепочечная РНК). Ожидаемо было бы увидеть небольшое увеличение количества связей, т.к. гидратация ДНК играет важную роль в превращении А-формы в В-форму, вода усиливает стэкинг-взаимодействие.