-
Notifications
You must be signed in to change notification settings - Fork 4
/
rmsd_y_rmsf.rmd
198 lines (158 loc) · 7.47 KB
/
rmsd_y_rmsf.rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
---
title: "Análisis de Trayectorias"
author: "Joel Ricci López"
output: html_document
---
```{r echo=FALSE, include = FALSE}
library(reticulate)
use_virtualenv('mds')
pyt <- import('pytraj')
plt <- import('matplotlib')
np <- import('numpy')
mda <- import('MDAnalysis')
nv <- import('nglview')
knitr::opts_chunk$set(fig.width=15, fig.height=5)
```
# RMSD y RMSF
## Carga de la tryectoria total
Bien, ahora podemos comenzar con el análisis de la trayectoria. Si bien, en general debemos enfocarnos en la fase de producción, para este ejemplo comenzaremos observando la evolución conformacional de la proteína a través de todas las fases de la trayectoria.
Comenzamos cargando las librerías y archivos necesarios:
```{python}
import pytraj as pyt
# cargamos el archivo de topología seleccionando sólo
# los primeros 20 residuos
f_topology = 'dm_sources_1L2Y/4-run_dm_2/tc5b_wb.psf'
top_prot_psf = pyt.load_topology(f_topology)[':1-20']
# Cargamos la trayectoria total de la dinámica, que guardamos
# en la sección anterior
name_file_full_traj = 'dm_sources_1L2Y/5-traj_analysis/tc5b_PROT_FULL_TRAJ.dcd'
full_traj = pyt.load(filename = name_file_full_traj,
top = top_prot_psf)
full_traj
```
## Cálculo del RMSD
Ahora procederemos a calcular el rmsd de todos los frames de la trayectoria con respecto a una estructura de referencia, que en este caso será la estructura lineal de la proteína; correspondiente al archivo `tc5b_wb.pdb`.
```{python}
# Cargamos la restructura lineal/inicial
init_struc = pyt.load(
filename = 'dm_sources_1L2Y/2-solvatar_wt/tc5b_wb.pdb',
top = 'dm_sources_1L2Y/2-solvatar_wt/tc5b_wb.psf', mask = ':1-20')
```
Ahora calculamos el RMSD de toda la trayectoria usando `lineal_strc` como estructura de referencia:
```{python}
# Cálculo del RMSD de la trayectoria total
# tomando como referencia la estructura inicial (lineal)
rmsd_full_REF_lineal = pyt.rmsd(traj = full_traj,
mask = ':4-17,@CA',
ref = init_struc)
```
### Grafica de RMSD
Para graficar el RMSD importamos las librerías matplotlib y seaborn y establecemos algunos parámetros por default para las gráficas.
```{python }
# Importamos las librerías y parámetros para generar las gráficas
from matplotlib import pyplot as plt
import seaborn as sns
import pylab
sns.set( context = 'talk', style = 'ticks', palette = "Spectral")
pylab.rcParams['figure.figsize'] = (15, 5)
```
Ahora generamos la gráfica de la trayectoria total:
```{python collapse=TRUE}
# Creamos la figura del RMSD de la trayectoria total
plt.scatter( x = range(full_traj.n_frames),
y = rmsd_full_REF_lineal, s = 20, alpha = 0.3 )
plt.title("RMSD trayectoria TOTAL. Ref = conformación lineal")
plt.xlabel("# Frame")
plt.show()
```
Como es de esperarse, conforme el tiempo de simulación avanza hay un aumento en el RMSD de cada frame con respecto a la conformación inicial.
Ahora bien, para un análisis más interesante usaremos como referencia la estructura por NMR de la proteína, con el objetivo de saber si, a lo largo de la simulación, las estructura de la proteína se aproxima a la estructura experimental.
Para ello cargamos la estructura NMR de la proteína.
```{python}
# Cargamos la estructura NMR
nmr_strc = pyt.load(filename = 'dm_sources_1L2Y/5-traj_analysis/1l2y_Hs.pdb',
top = 'dm_sources_1L2Y/5-traj_analysis/1l2y_Hs.psf')
```
Calculamos el RMSD de la trayectoria usando como referencia a nmr_strc. Seguido de ello, generamos la gráfica.
```{python}
rmsd_full_REF_nmr = pyt.rmsd( traj = full_traj,
mask = ':4-17,@CA',
ref = nmr_strc)
```
```{python}
plt.scatter( x = range(full_traj.n_frames), y = rmsd_full_REF_nmr,
s = 20, alpha = 0.3, c='orange' )
plt.title("RMSD trayectoria TOTAL. Ref = conformación NMR 1l2y")
plt.xlabel("# Frame")
plt.show()
```
Como ves, al principio de la simulación la conformación inicial y la estructura por NMR son muy distintas (entre 6 a 8 A), sin embargo, a lo largo de la dinámica y sobre todo en la producción, la conformación de la porteína tiende a ser más parecida a la estructura NMR.
Finalmente, visualicemos únicamente la fase de producción con ambas referencias:
```{python}
import pandas as pd
# Creamos una tabla con pandas
# Dado que la producción corresponde a los últimos 5000 frames,
# usamos el slicing [-5000:] para extraer esos últimos 5000 datos
rmsd_prod = pd.DataFrame({"Ref_lineal": rmsd_full_REF_lineal[-5000:],
"Ref_NMR": rmsd_full_REF_nmr[-5000:]})
# Sólo por variar, vamos a usar seaborn para generar el gráfico:
ax = sns.lineplot(data = rmsd_prod)
ax.set_title("RMSD producción")
```
## Cálculo del RMSF
Vamos a finalizar este análisis evaluando el RMSF de la proteína a lo largo de la **fase de producción**.
Carguemos primero dicha trayectoria (Aunque bien podriamos usar los últimos frames del objeto `full_traj`).
```{python}
f_prod = 'dm_sources_1L2Y/5-traj_analysis/tc5b_PROT_PROD.dcd'
prod_traj = pyt.load(filename = f_prod, top = top_prot_psf)
```
#### Un pequeño paréntesis:
**_¿Es `full_traj[:-5000]` igual a `prod_traj`?_**
Es decir, ¿los últimos 5000 frames del objeto `full_traj` son iguales a todos frames de `prod_traj`?
```{python}
print(prod_traj)
print(full_traj[-5000:])
prod_traj == full_traj[-5000:]
```
Si bien ambos objetos, `` ``, tienen la misma información, *python* nos dice que `prod_traj == full_traj[-5000:]` son distintos (`False`). Esto se debe a que las variables **por _referencia_** `prod_traj` y `full_traj` apuntan a objetos distintos en memoria, pero esto no significa que ambos objetos puedan tener información similar. Esto es importante y lo debemos de tener en cuenta al momento de querer copiar o manipular secciones de trayectorias. Finálicemos este paréntesis con el siguiente código:
```{python}
prod_traj_NO_copia = prod_traj # Apuntan al mismo objeto
prod_traj_copia = prod_traj.copy() # Apuntan a diferentes objetos
```
```{python}
prod_traj == prod_traj_NO_copia
prod_traj == prod_traj_copia
```
### Obtención del frame promedio y superposición de la trayectoria
Ahora volvamos al ejercicio para proceder a superponer toda la trayectoria a un frame promedio.
```{python}
# Obtenemos el frame promedio
avg_prod_frame = pyt.mean_structure(traj = prod_traj)
# Y superponemos la trayectoria a dicho frame
prod_traj.superpose(ref = avg_prod_frame, mask = ':3-18@CA')
```
Ahora calculamos el RMSF de toda la producción:
```{python}
rmsf_prod_traj = pyt.rmsf( prod_traj, mask = '@CA' )
# Graficamos
plt.plot( rmsf_prod_traj.T[1], c = 'red')
plt.title("RMSF producción total")
plt.xlabel("# Residuo")
plt.show()
```
PAra terminar esta sección, un análisis interesante sería ver cómo el RMSF de la proteína va cambiando tomando ciertos intervalos de la producción.
```{python}
# Ahora calculamos el RMSF a diferentes intervalos de la producción
num_interv = 5
interv_size = int(prod_traj.n_frames / num_interv)
df_rmsf = pd.DataFrame(
{F'Interv. {i}': pyt.rmsf( prod_traj[i*interv_size : (i+1)*interv_size],
mask = '@CA').T[1] for i in range(num_interv)} )
```
Graficamos usando `seaborn`.
```{python}
sns.set_palette("husl")
ax = sns.lineplot(data = df_rmsf, dashes = False)
ax.set_title("RMSF producción")
```
En general, podremos observar que hacia la última parte de la trayectoria, el backbone de la proteína tiende a ser más estable. Continuaremos con más análisis de la trayectoria en la siguiente sección.