analisis_traj.rmd

---
title: "Análisis de Trayectorias"
author: "Joel Ricci López"
output: html_document
---

# Manejo de trayectorias con *pytraj*

```{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) 
```

Como primer paso vamos a cargar los datos de la trayectoria de la DM, considerando cada una de las fases realizadas.  
Para ello, comenzamos importando el módulo `pytraj`.

```{python}
import pytraj as pyt
```

A continuación, cargamos la estructura inicial de la porteína, que corresponde a la conformación lineal a partir de la cual realizamos la minimización y la dinámica molecular.  
Observa que la función `pyt.load()` recibe los parámetros `filename` y `top` correspondientes al archivo de coordenadas y de topología, respectivamente.

```{python}
# Cargamos la restructura inicial, es decir, la estructura lineal del péptido
lineal_strc = pyt.load(
    filename = 'dm_sources_1L2Y/1-topologia/tc5b.pdb',
    top = 'dm_sources_1L2Y/1-topologia/tc5b.psf')
```

Ahora, si llamamos a la variable *`cristal_strc`*, veremos que ésta corresponde a un objeto *`pytraj.Trajectory`* con los siguientes atributos:

```{python}
lineal_strc
```

## Visualización de la esturctura
Recordemos que podemos utilizar la librería *nglview* para visualizar la estructura de la proteína:

```{python eval=FALSE, collapse=TRUE}
# Importamos nglview
import nglview as nv
import warnings
lineal_strc_view = nv.show_pytraj(lineal_strc)
lineal_strc_view.background = '#303030'
lineal_strc_view
```

<div style="height: 400px; width: 100%; position: relative;" class='viewer_3Dmoljs justify-content-center border' data-href='https://raw.githubusercontent.com/jRicciL/MD_namd_python/master/dm_sources_1L2Y/1-topologia/tc5b.pdb' 
data-backgroundcolor='0x303030' 
data-style1='cartoon:color=spectrum' 
data-surface1='opacity:.5;color:white'>
</div> 

## Carga de la trayectoria
Ahora, vamos a proceder a cargar en memoria los archivos de trayectoria `dcd` de todas las fases de la dinámica. Para ello, es necesario primero definir el archivo de topología que utilizaremos.

```{python}
f_topology = 'dm_sources_1L2Y/4-run_dm_2/tc5b_wb.psf'
top_psf = pyt.load_topology(f_topology)
top_psf
```

Este archivo de topoligía incluye a todas las molécuals del sistema, sin embargo, por ahora sólo nos interesa analizar la trayectoria de la proteína, por lo que redefiniremos el archivo de topología para que sólo incluya los átomos de la misma (los primeros 20 residuos del sistema).

```{python}
top_prot_psf = top_psf[':1-20']
top_prot_psf 
```

Ahora procederemos a cargar los archivos dcd. Sin embargo, para este análisis cargaremos archivos dcd a los cuales se les ha removido el solvente, con el objetivo de reducir el tamaño de los archivos y poder distribuirlos en línea.  
A continuación mostramos el ejemplo de cómo se generaron estos archivos "reducidos" a partir de los dcd originales.

```{python}
# Para optimizar el análisis, previamente se han extraido de las trayectorias
# únicamente los átomos correspondientes a la proteína y guardado en la carpeta
# 5_traj_analysis

# Un ejemplo de cómo realizar dicha extracción es el siguiente:
_min = pyt.load(filename = 'dm_sources_1L2Y/4-run_dm_2/1_min/tc5b_wb_min.dcd', 
                top = 'dm_sources_1L2Y/4-run_dm_2/tc5b_wb.psf')
print(_min)
# Ahora guardamos un nuevo archivo de trayectoria:
pyt.write_traj(filename = 'dm_sources_1L2Y/5-traj_analysis/tc5b_PROT_MIN.dcd', 
               traj = _min[':1-20'], 
               overwrite=True)
print(F'Archivo guardado:\n{_min[":1-20"]}')
```

### Archivos *.dcd*

Ahora ubicamos el directorio y nombre de cada uno de los archivos dcd. Esto incluye también al archivo .pdb con la estructura inicial.

```{python}
# PDB inicial
f_inicial = 'dm_sources_1L2Y/2-solvatar_wt/tc5b_wb.pdb'
# Archivos de dinámica
dir_traj = 'dm_sources_1L2Y/5-traj_analysis' # Directorio de los archivos
f_min = F'{dir_traj}/tc5b_PROT_MIN.dcd'
f_heat = F'{dir_traj}/tc5b_PROT_SA.dcd'
f_eq = F'{dir_traj}/tc5b_PROT_EQ.dcd'
f_prod = F'{dir_traj}/tc5b_PROT_PROD.dcd'
```

Ahora cargamos cada uno de los archivos utilizando el objeto de topología `top_psf`. Como ejemplo puedes ver que también es posible usar el parámetro mask para determinar qué átomos deberán ser considerados al cargarse en memoria. Sin embargo, podemos omitirlo en las fases de minimización en adelante, ya que esta selección de los átomso de laproteína ya la hicimos al cargar el archivo de topología.

```{python}
#Estructura inicial
init_struc = pyt.load(filename = f_inicial, 
                      top = top_psf, mask = ':1-20')
print( F'Estructura inicial: {init_struc.n_frames} frames')

# Minimización
min_traj = pyt.load(filename = f_min, top = top_prot_psf)
print( F'Minimización: {min_traj.n_frames} frames')
# Calentamiento
heat_traj = pyt.load(filename = f_heat, top = top_prot_psf)
print( F'Calentamiento: {heat_traj.n_frames} frames')
# Equilibrado
eq_traj = pyt.load(filename = f_eq, top = top_prot_psf)
print( F'Equilibrado: {eq_traj.n_frames} frames')
# Producción
prod_traj = pyt.load(filename = f_prod, top = top_prot_psf)
print(F'Producción: {prod_traj.n_frames} frames')
```

Además de cargar los archivos, hemos usado el atributo n_frames que nos dice cuántos frames posee cada tryectoria. En python es posible utilizar las funciones `type()` y `dir()` para conocer qué métodos y atributos posee un objeto según la clase a que pertenece.

```{python}
print(type( prod_traj ))
print(dir( prod_traj ))
```

## Concatenación de las trayectorias
Ahora, con el objetivo de ver la evolución conformacional de la proteína a lo largo de todas las fases, vamos a unir todas las trayectorias en una sola. Esto también nos permitirá guardar una única trayectoria para posteriores análisis.  
Primero generamos una lista con los objetos de trayectoria:

```{python collapse=TRUE}
# Combinamos todas las etapas para visualizar una única trayectoria
trajs_list = [min_traj, heat_traj, eq_traj, prod_traj]
```

El siguiente paso es usar la lista traj_list para iterar sobre sus elementos.  
Para ello vamos  utilizar una expresión conosida como list comprehension. Por ahora las usaremos para uardar en una nueva lista el número de frames de cada fase.

```{python}
# Número total de frames
n_frames_list = [ traj.n_frames for traj in trajs_list ] # Observa la sintaxis
n_frames_list
```

```{python}
n_full_frames = sum(n_frames_list)
print( F'Número total de frames: {n_full_frames}')
```

Ahora crearemos un nuevo objeto `Trajectory` inicializándolo como una copia del objeto `init_struct`, al que iremos añadiendo las coordenadas de los objetos de la lista `trajs_list` usando el método `append_xyz`:

```{python}
# Inicializamos la trayectoria con una copia de la estructura inicial, 
# que es un objeto de la clase Trajectory de pytraj
full_traj = init_struc.copy()

# Iteramos sobre la lista de trayecorias, añadiendo las coordenadas
# al objeto full_traj
for traj in trajs_list:
    full_traj.append_xyz(traj.xyz)

full_traj
```

### Superposición de las estructuras
  
Primero vamos a superponer todas las conformaciones de la proteína. Para ello, el método superpose realiza moviemientos translacionales y rotacionales de forma iterativa minimizando el RMSD entre todas las confromaciones.
Observa que el alineamiento lo realizaremos considerando los carbonos alfa únicamente de los residuos 4 a 17, es decir, hemos omitido los primeros y los últimos tres reiduos de la proteína. Si lo deseas puedes reducir o extender esta selección.

```{python}
# Superposición de los frames
full_traj.superpose(mask = ':4-17@CA', ref = 0) 
# Por default el primer frame es tomado como referencia
```

### Guardando el archivo de tryectoria
Como vimos en el ejemplo de la tayectoria de la minimización, podemos usar la función `write_traj()`. Observa qué parámetros se requieren para guardar el nuevo archivo.

```{python eval=FALSE}
# Podemos guardar el archivo para trabajar con él en un posterior análisis
name_file_full_traj = F'{dir_traj}/tc5b_PROT_FULL_TRAJ.dcd'
pyt.write_traj(filename = name_file_full_traj, 
               traj = full_traj, overwrite = True)
```


### Visualización de la trayectoria
Podemos usar *nglview* para visualizar la tryectoria con los frames superpuestos.
```{python  eval=FALSE}
full_traj_sup_view = nv.show_pytraj(full_traj)
full_traj_sup_view.add_representation("licorice")
full_traj_sup_view
```

<div class='container text-center'>
<script src="https://unpkg.com/ngl@0.10.4/dist/ngl.js"></script>
  <script>
    document.addEventListener("DOMContentLoaded", function () {
       var stage = new NGL.Stage( "viewport", {backgroundColor:'white'}  );
       stage.loadFile("https://raw.githubusercontent.com/jRicciL/MD_namd_python/master/dm_sources_1L2Y/1-topologia/tc5b.psf").then(function (o) {
          NGL.autoLoad("https://raw.githubusercontent.com/jRicciL/MD_namd_python/master/dm_sources_1L2Y/5-traj_analysis/tc5b_PROT_SA.dcd").then(function (frames) {
            o.addTrajectory(frames, {
              initialFrame: 0,
              superpose: false
            })
            var sele = "not backbone or .CA or (PRO and .N)"
            o.addRepresentation("licorice")
            o.addRepresentation("cartoon", { color: "residueindex", aspectRatio: 4, scale: 0.5 })
            var traj = o.trajList[0].trajectory;
            var player = new NGL.TrajectoryPlayer( traj, { } );
            traj.setPlayer( player );
            traj.player.play();
            stage.autoView();
        })
      })
    });
  </script>
<div id="viewport" style="height: 300px; width: 100%; position: relative;" class='justify-content-center'></div>
</div>