In [ ]:
%matplotlib inline
import numpy as np
import pylab as pl
import matplotlib.patches as mpatches
import matplotlib.ticker as ticker
import os
import shutil
from IPython.display import Image
from matplotlib.ticker import FormatStrFormatter
In [ ]:
ruta=os.getcwd()
c=input('Nombre de la trayectoria para realizar el análisis... Ejemplo: run001....')
if os.path.isdir(c):
indir = '/'+c
print (indir)
ruta_old_traj=ruta+indir
print (ruta)
print (ruta_old_traj)
else:
print ('La carpetac'+c+' no existe...')
In [ ]:
#
ruta_scripts=ruta+'/scripts_fimda'
print (ruta_scripts)
if os.path.exists(ruta_scripts):
print ('Ruta identificada para búsqueda de scripst adicionales ===>',ruta_scripts)
else:
print ('La carpeta de scripst adicionales no existe, copiar en '+ruta_scripts+' ..!!!')
In [ ]:
#Verificando que exista la nueva carpeta para la conversión de trayectorias
#nuevaruta = ruta+'/'+indir+'_XTC'
nuevaruta = ruta+indir+'_Dinamica'
print ( nuevaruta )
if not os.path.exists(nuevaruta):
os.makedirs(nuevaruta)
print ('Se ha creado la ruta ===>',nuevaruta)
else:
print ("La ruta "+nuevaruta+" existe..!!!")
In [ ]:
print ('Obtenemos los archivos a convertir')
#Buscamos el archivo DCD, PDB y PSF para realizar las operaciones
for filename in os.listdir(ruta_old_traj):
if filename.endswith('.dcd'):
dcd_file=filename
if filename.endswith('.psf'):
psf_file=filename
if filename.endswith('.pdb'):
pdb_file=filename
print ('pdb file =>', pdb_file)
print ('psf file =>', psf_file)
print ('dcd file =>', dcd_file)
print ( 'Nos vemos a ....', ruta_old_traj )
os.chdir( ruta_old_traj )
print ('\nEjecutando CATDCD para convertir la trayectoria....')
output_catdcd=!catdcd -otype trr -o output.trr $dcd_file
print (output_catdcd.n)
print ('\nEjecutando TRJCONV para convertir la trayectoria....')
output_trjconv=!trjconv -f output.trr -o output.xtc -timestep 20
#print (output_trjconv.n)
print ('\nBorrando archivos temporales de conversión...')
output_rm=!rm output.trr
print ('\nMoviendo los archivos de salida al directorio '+nuevaruta)
source_file=ruta_old_traj+'/output.xtc'
dest_file=nuevaruta+'/output.xtc'
shutil.move(source_file,dest_file)
print ('\Copiando el archivo ionized.pdb a '+nuevaruta)
source_file=ruta_old_traj+'/ionized.pdb'
dest_file=nuevaruta+'/ionized.pdb'
shutil.copy(source_file,dest_file)
print ('\nCopiando el archivo ionized.psf a '+nuevaruta)
source_file=ruta_old_traj+'/ionized.psf'
dest_file=nuevaruta+'/ionized.psf'
shutil.copy(source_file,dest_file)
print('\nTrayectoria convertida, regresando a '+ruta)
os.chdir( ruta )
In [ ]:
print ('Visualizando la nueva trayectoria')
file_psf=nuevaruta+'/'+psf_file
traj = nuevaruta+'/output.xtc'
!vmd $file_psf $traj
In [ ]:
### Creando el directorio para el análisis del RMSD
#Verificando que exista la nueva carpeta para la conversión de trayectorias
#nuevaruta = ruta+'/'+indir+'_XTC'
ruta_rmsd = nuevaruta+'/rmsd'
print ( ruta_rmsd )
if not os.path.exists(ruta_rmsd):
os.makedirs(ruta_rmsd)
print ('Se ha creado la ruta ===>',ruta_rmsd)
else:
print ("La ruta "+ruta_rmsd+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rmsd )
os.chdir( ruta_rmsd )
In [ ]:
print ('Ejecutando el análisis de rmsd...')
!echo 3 3 | g_rms -f ../output.xtc -s ../ionized.pdb -a avgrp.xvg
In [ ]:
#Inicializando vector
rmsd=[]
try:
archivo = open( 'rmsd.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in archivo.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
num=float(sl[0])
#num2=float(sl[1])
num=num/1000
rmsd.append(repr(num)+'\t'+sl[1]+'\n')
i=i+1
#Escribiendo el archivo RMSD
f = open('rmsd.dat', 'w')
#f.write('@ title "RMSD" \n')
f.write('@ xaxis label " Time (ns)" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label " RMSD (nm)" \n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 1.5\n')
f.write('@TYPE xy \n')
#f.write('@ subtitle "C-alpha after lsq fit to C-alpha" \n')
f.write("".join(rmsd))
f.close()
In [ ]:
#Cargando el archivo para visualizar en xmgrace
!xmgrace rmsd.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='rmsd.png')
In [ ]:
#Inicializando vector
rmsd_residue=[]
try:
archivo_rmsd = open( 'aver.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=1
for linea in archivo_rmsd.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
num=int(sl[0])
print ('Residuo =>',num+1)
rmsd_residue.append(repr(num+1)+'\t'+sl[1]+'\n')
i=i+1
#Escribiendo el archivo RMSD_RESIDUE
f = open('rmsd_residue.dat', 'w')
#f.write('@ title "C-alpha" \n')
f.write('@ xaxis label "Residue" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label " RMSD (nm)" \n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 2.5\n')
f.write('@ s0 symbol 1\n')
f.write('@ s0 symbol size 1.000000\n')
f.write('@ s0 symbol color 1\n')
f.write('@ s0 symbol pattern 1\n')
f.write('@ s0 symbol fill color 2\n')
f.write('@ s0 symbol fill pattern 1\n')
f.write('@ s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(rmsd_residue))
f.close()
In [ ]:
!xmgrace rmsd_residue.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='rmsd_residue.png')
In [ ]:
data_rmsd=np.loadtxt('rmsd.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsd[:,0]/1000, data_rmsd[:,1], linewidth = 2, markeredgewidth=3, color='black')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.suptitle('RMSD', fontsize=50)
#pl.title('C-alpha after lsq fit to C-alpha', fontsize=30)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
In [ ]:
data_rmsd_res=np.loadtxt('aver.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsd_res[:,0]+1, data_rmsd_res[:,1], '-o', color='black', markersize=25,
markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.title('C-alpha', fontsize=40)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
pl.xlim(0, len(data_rmsd_res[:,1]))
In [ ]:
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_rmsf = nuevaruta+'/rmsf'
print ( ruta_rmsf )
if not os.path.exists(ruta_rmsf):
os.makedirs(ruta_rmsf)
print ('Se ha creado la ruta ===>',ruta_rmsf)
else:
print ("La ruta "+ruta_rmsf+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rmsf )
os.chdir( ruta_rmsf )
In [ ]:
print ('Ejecutando el análisis de rmsf...')
!echo 3 | g_rmsf -f ../output.xtc -s ../ionized.pdb -oq bfac.pdb -o rmsf.xvg -res
In [ ]:
#Inicializando vector
rmsf=[]
rmsf_x=[]
rmsf_y=[]
try:
file_rmsf = open( 'rmsf.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_rmsf.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
print ('Residue =>',cadena)
rmsf.append(sl[0]+'\t'+sl[1]+'\n')
rmsf_x.append(int(sl[0]))
rmsf_y.append(float(sl[1]))
i=i+1
file_rmsf.close()
#Escribiendo el archivo RMSD
f = open('rmsf.dat', 'w')
#f.write('@ title "RMSF fluctuation" \n')
f.write('@ xaxis label " Residue" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "RMSF (nm)" \n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 2.5\n')
f.write('@ s0 symbol 1\n')
f.write('@ s0 symbol size 1.000000\n')
f.write('@ s0 symbol color 1\n')
f.write('@ s0 symbol pattern 1\n')
f.write('@ s0 symbol fill color 2\n')
f.write('@ s0 symbol fill pattern 1\n')
f.write('@ s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(rmsf))
f.close()
In [ ]:
!xmgrace rmsf.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='rmsf.png')
In [ ]:
data_rmsf=np.loadtxt('rmsf.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rmsf[:,0], data_rmsf[:,1], '-o', color='black', markersize=25,
markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 40)
pl.ylabel('RMSF (nm)', fontsize = 40)
#pl.title('RMSF Fluctuation', fontsize=40)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
pl.xlim(0, len(data_rmsf[:,1]))
In [ ]:
#Inicializando vector
bfactors=[]
try:
file_bfactor = open( 'bfac.pdb' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_bfactor.readlines():
fila = linea.strip()
sl = fila.split()
if (sl[0]=='ATOM'):
#print (sl[0])
idresidue=fila[23:26]
bfactor=fila[60:66]
print (idresidue + '\t'+bfactor)
bfactors.append(idresidue+'\t'+bfactor+'\n')
#i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('bfactor.dat', 'w')
#f.write('@ title "B-factors" \n')
foo = 'baz "\\"'
f.write('@ xaxis label " Residue" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "B-factors (' +"\\"+'cE'+"\\"+'C)"\n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 2.5\n')
f.write('@ s0 symbol 1\n')
f.write('@ s0 symbol size 1.000000\n')
f.write('@ s0 symbol color 1\n')
f.write('@ s0 symbol pattern 1\n')
f.write('@ s0 symbol fill color 2\n')
f.write('@ s0 symbol fill pattern 1\n')
f.write('@ s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(bfactors))
f.close()
In [ ]:
!xmgrace bfactor.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='bfactor.png')
In [ ]:
#Inicializando vector
bfactors=[]
try:
file_bfactor = open( 'bfac.pdb' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
print ('Residuo' + '\t'+'bfactor')
for linea in file_bfactor.readlines():
fila = linea.strip()
sl = fila.split()
if (sl[0]=='ATOM'):
#print (sl[0])
idresidue=fila[23:26]
bfactor=fila[60:66]
print (idresidue + '\t'+bfactor)
bfactors.append(idresidue+'\t'+bfactor+'\n')
#i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('bfactor.dat', 'w')
f.write("".join(bfactors))
f.close()
data_bfactor=np.loadtxt('bfactor.dat',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
#ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_bfactor[:,0], data_bfactor[:,1], '-o', color='black', markersize=25,
markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel('Residue', fontsize = 40)
pl.ylabel('B-factors ('+ r'$\AA$'+')' , fontsize = 40)
#pl.title('B-Factors', fontsize=40)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
pl.xlim(0, len(data_bfactor[:,1]))
In [ ]:
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_ss = nuevaruta+'/estructura'
print ( ruta_ss )
if not os.path.exists(ruta_ss):
os.makedirs(ruta_ss)
print ('Se ha creado la ruta ===>',ruta_ss)
else:
print ("La ruta "+ruta_ss+" existe..!!!")
print ( 'Nos vamos a ....', ruta_ss )
os.chdir( ruta_ss )
In [ ]:
print ('Ejecutando el análisis de esctructura secundaria...')
!echo 5 | do_dssp -f ../output.xtc -s ../ionized.pdb -o sec_est.xpm -tu ns
In [ ]:
print ('\n Convirtiendo el archivo a ps...')
!xpm2ps -f sec_est.xpm -by 6 -bx .1 -o est_sec.eps
In [ ]:
print('\nConvirtiendo a png...')
!convert -density 600 est_sec.eps -resize 1024x1024 est_sec.png
In [ ]:
print ('Cargando el archivo...')
Image(filename='est_sec.png', width=1024)
In [ ]:
### Creando el directorio para el análisis del r-gyro
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_rgyro = nuevaruta+'/rgyro'
print ( ruta_rgyro )
if not os.path.exists(ruta_rgyro):
os.makedirs(ruta_rgyro)
print ('Se ha creado la ruta ===>',ruta_rgyro)
else:
print ("La ruta "+ruta_rgyro+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rgyro)
os.chdir( ruta_rgyro )
In [ ]:
print ('Ejecutando el análisis de rgyro...')
!echo 3 | g_gyrate -f ../output.xtc -s ../ionized.pdb -o gyrate.xvg
In [ ]:
#Inicializando vector
rgyro=[]
try:
file_rmsf = open( 'gyrate.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_rmsf.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
num=float(sl[0])
#num2=float(sl[1])
num=num/1000
rgyro.append(repr(num)+'\t'+sl[1]+'\n')
i=i+1
#Escribiendo el archivo RGYRO.DAT
f = open('rgyro.dat', 'w')
#f.write('@ title "Radius of gyration" \n')
f.write('@ xaxis label " Time (ns)" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "Rg (nm)" \n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 2.5\n')
f.write('@TYPE xy \n')
f.write("".join(rgyro))
f.close()
In [ ]:
!xmgrace rgyro.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='rgyro.png')
In [ ]:
data_rgyro=np.loadtxt('gyrate.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
pl.plot(data_rgyro[:,0]/1000, data_rgyro[:,1], linewidth = 2, color='black')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('Rg (nm)', fontsize = 40)
#pl.suptitle('Radius of gyration', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
Para realizar este análisis se debe cargar el pdb original de la proteina que se encuentra en la carpeta 01_BUILD.
Cargarlo con VMD y dirigirse al Menú EXTENSIONS -> ANALYSIS -> SEQUENCE VIEWER, en la cual se tomará el rango de átomos del campo Struct (H), el cual se proporcionará de la forma "resid X1 to X2" donde X1 es primer átomo de la helix y X2 el último átomo de la helix.
In [ ]:
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_helix = nuevaruta+'/rmsd_helix'
print ( ruta_helix )
if not os.path.exists(ruta_helix):
os.makedirs(ruta_helix)
print ('Se ha creado la ruta ===>',ruta_helix)
else:
print ("La ruta "+ruta_helix+" existe..!!!")
print ( 'Nos vamos a ....', ruta_helix)
os.chdir( ruta_helix )
In [ ]:
num=input('Número de hélices con las que cuenta la proteína:')
print (num)
In [ ]:
if (int(num)==1):
indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
print (indices_ha1)
r_helix_1=1
r_helix_2=0
r_helix_3=0
r_helix_4=0
if (int(num)==2):
indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
print (indices_ha1)
indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
print (indices_ha2)
r_helix_1=1
r_helix_2=1
r_helix_3=0
r_helix_4=0
if (int(num)==3):
indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
print (indices_ha1)
indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
print (indices_ha2)
indices_ha3=input('Proporciona el rango de índices de la Hélice 3:')
print (indices_ha3)
r_helix_1=1
r_helix_2=1
r_helix_3=1
r_helix_4=0
if (int(num)==4):
indices_ha1=input('Proporciona el rango de índices de la Hélice 1:')
print (indices_ha1)
indices_ha2=input('Proporciona el rango de índices de la Hélice 2:')
print (indices_ha2)
indices_ha3=input('Proporciona el rango de índices de la Hélice 3:')
print (indices_ha3)
indices_ha4=input('Proporciona el rango de índices de la Hélice 4:')
print (indices_ha4)
r_helix_1=1
r_helix_2=1
r_helix_3=1
r_helix_4=1
In [ ]:
#Script para vmd de la Hélice Alfa 2
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
if (r_helix_1==1):
f = open('ha1.tcl', 'w')
print(f)
f.write('set psfFile '+ psf+' \n')
f.write('set dcdFile '+ dcd+' \n')
f.write('\nmol load psf $psfFile dcd $dcdFile\n')
f.write('set outfile ' +'[open ' +'rmsd_ha1.dat'+' w]\n')
f.write('set nf [molinfo top get numframes]\n')
f.write('\n#RMSD calculation loop\n')
f.write('set f1 [atomselect top "'+indices_ha1+' " frame 0]\n')
f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f.write(' set sel [atomselect top "'+indices_ha1+' " frame $i]\n')
f.write(' $sel move [measure fit $sel $f1]\n')
f.write(' set time [expr {$i +1}]\n')
f.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f.write(' puts $outfile "$time $time"\n')
f.write('}\n')
f.write('close $outfile')
f.close()
if (r_helix_2==1):
f = open('ha2.tcl', 'w')
print(f)
f.write('set psfFile '+ psf+' \n')
f.write('set dcdFile '+ dcd+' \n')
f.write('\nmol load psf $psfFile dcd $dcdFile\n')
f.write('set outfile ' +'[open ' +'rmsd_ha2.dat'+' w]\n')
f.write('set nf [molinfo top get numframes]\n')
f.write('\n#RMSD calculation loop\n')
f.write('set f1 [atomselect top "'+indices_ha2+' " frame 0]\n')
f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f.write(' set sel [atomselect top "'+indices_ha2+' " frame $i]\n')
f.write(' $sel move [measure fit $sel $f1]\n')
f.write(' set time [expr {$i +1}]\n')
f.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f.write(' puts $outfile "$time $time"\n')
f.write('}\n')
f.write('close $outfile')
f.close()
if (r_helix_3==1):
f = open('ha3.tcl', 'w')
print(f)
f.write('set psfFile '+ psf+' \n')
f.write('set dcdFile '+ dcd+' \n')
f.write('\nmol load psf $psfFile dcd $dcdFile\n')
f.write('set outfile ' +'[open ' +'rmsd_ha3.dat'+' w]\n')
f.write('set nf [molinfo top get numframes]\n')
f.write('\n#RMSD calculation loop\n')
f.write('set f1 [atomselect top "'+indices_ha3+' " frame 0]\n')
f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f.write(' set sel [atomselect top "'+indices_ha3+' " frame $i]\n')
f.write(' $sel move [measure fit $sel $f1]\n')
f.write(' set time [expr {$i +1}]\n')
f.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f.write(' puts $outfile "$time $time"\n')
f.write('}\n')
f.write('close $outfile')
f.close()
if (r_helix_4==1):
f = open('ha4.tcl', 'w')
print(f)
f.write('set psfFile '+ psf+' \n')
f.write('set dcdFile '+ dcd+' \n')
f.write('\nmol load psf $psfFile dcd $dcdFile\n')
f.write('set outfile ' +'[open ' +'rmsd_ha4.dat'+' w]\n')
f.write('set nf [molinfo top get numframes]\n')
f.write('\n#RMSD calculation loop\n')
f.write('set f1 [atomselect top "'+indices_ha4+' " frame 0]\n')
f.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f.write(' set sel [atomselect top "'+indices_ha4+' " frame $i]\n')
f.write(' $sel move [measure fit $sel $f1]\n')
f.write(' set time [expr {$i +1}]\n')
f.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f.write(' puts $outfile "$time $time"\n')
f.write('}\n')
f.write('close $outfile')
f.close()
In [ ]:
if (r_helix_1==1):
#Calculando con VMD hélice 1
!vmd -dispdev text < ha1.tcl
if (r_helix_2==1):
#Calculando con VMD hélice 2
!vmd -dispdev text < ha2.tcl
if (r_helix_3==1):
#Calculando con VMD hélice 3
!vmd -dispdev text < ha3.tcl
if (r_helix_4==1):
#Calculando con VMD hélice 4
!vmd -dispdev text < ha4.tcl
In [ ]:
if (int(num)==1):
#Graficando
data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
#pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#pl.title('RMSD Helix Alfa', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
if (int(num)==2):
#Graficando
data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
#pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#pl.title('RMSD Helix Alfa', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
if (int(num)==3):
#Graficando
data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
data_ha3=np.loadtxt('rmsd_ha3.dat',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
#pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
pl.plot(data_ha3[:,1]*0.02, data_ha3[:,0]/10, linewidth = 3, color='green')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (nm)', fontsize = 40)
#pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#pl.title('RMSD Helix Alfa', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
if (int(num)==4):
#Graficando
data_ha1=np.loadtxt('rmsd_ha1.dat',comments=['#', '@'])
data_ha2=np.loadtxt('rmsd_ha2.dat',comments=['#', '@'])
data_ha3=np.loadtxt('rmsd_ha3.dat',comments=['#', '@'])
data_ha4=np.loadtxt('rmsd_ha4.dat',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f'))
#pl.plot(data_ha1[:,0], data_ha1[:,1], linewidth = 3)
pl.plot(data_ha1[:,1]*0.02, data_ha1[:,0]/10, linewidth = 3, color='black')
pl.plot(data_ha2[:,1]*0.02, data_ha2[:,0]/10, linewidth = 3, color='red')
pl.plot(data_ha3[:,1]*0.02, data_ha3[:,0]/10, linewidth = 3, color='green')
pl.plot(data_ha4[:,1]*0.02, data_ha4[:,0]/10, linewidth = 3, color='blue')
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('RMSD (A)', fontsize = 40)
#pl.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#pl.title('RMSD Helix Alfa', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
In [ ]:
### Creando el directorio para el análisis del SASA
### NOTA: se calcula con gromacs4 ya que arroja bien los resultados comparado con gromacs5
ruta_sasa = nuevaruta+'/sasa'
print ( ruta_sasa )
if not os.path.exists(ruta_sasa):
os.makedirs(ruta_sasa)
print ('Se ha creado la ruta ===>',ruta_sasa)
else:
print ("La ruta "+ruta_sasa+" existe..!!!")
print ( 'Nos vamos a ....', ruta_sasa )
os.chdir( ruta_sasa )
In [ ]:
print ('Ejecutando el análisis de sasa con Gromacs 4 utilizando la opción 1 (protein)...')
!echo 1 1 | /opt/gromacs4/bin/g_sas -f ../output.xtc -s ../ionized.pdb -o solven-accessible-surface.xvg -oa atomic-sas.xvg -or residue-sas.xvg
In [ ]:
#Inicializando vector
sasa_residuo=[]
try:
residue_sas = open( 'residue-sas.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in residue_sas.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
print ('Residue =>',cadena)
sasa_residuo.append(sl[0]+'\t'+sl[1]+'\n')
i=i+1
#Escribiendo el archivo RMSD
f = open('sasa-residuo.dat', 'w')
#f.write('@ title "Area per residue over the trajectory" \n')
f.write('@ xaxis label " Residue " \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "Area (nm' +"\\"+'S2'+"\\N"+')"\n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 2.5\n')
f.write('@ s0 symbol 1\n')
f.write('@ s0 symbol size 1.000000\n')
f.write('@ s0 symbol color 1\n')
f.write('@ s0 symbol pattern 1\n')
f.write('@ s0 symbol fill color 2\n')
f.write('@ s0 symbol fill pattern 1\n')
f.write('@ s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(sasa_residuo))
f.close()
In [ ]:
!xmgrace sasa-residuo.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='sasa-residuo.png')
In [ ]:
data_sasa_residue=np.loadtxt('residue-sas.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
pl.plot(data_sasa_residue[:,0], data_sasa_residue[:,1], '-o', color='black', markersize=25,
markerfacecolor='red',markeredgecolor='black',markeredgewidth=3, linewidth = 4, )
pl.xlabel("Residue", fontsize = 30)
#pl.ylabel('Area (nm2)', fontsize = 30)
pl.ylabel('Area ( nm'+ r'$\ ^2$'+')' , fontsize = 40)
#pl.title('Area per residue over the trajectory', fontsize=40)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
pl.xlim(0, len(data_sasa_residue[:,1]))
In [ ]:
#Inicializando vector
sasa=[]
try:
sasafile = open( 'solven-accessible-surface.xvg' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in sasafile.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (not '#' in cadena) and (not '@' in cadena):
#print (cadena)
num=float(sl[0])
num=num/1000
sasa.append(repr(num)+'\t'+sl[1]+'\t'+sl[2]+'\t'+sl[3]+'\n')
i=i+1
cel2=float(sl[2])
print(cel2)
#Escribiendo el archivo RMSD
f = open('sasa.dat', 'w')
#f.write('@ title "Solven Accessible Surface" \n')
f.write('@ xaxis label " Time (ns) " \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 3.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "Area (nm' +"\\"+'S2'+"\\N"+')"\n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 3.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
#f.write('@ s0 legend "Hydrophobic"\n')
#if (cel2>0):
#f.write('@ s1 legend "Hydrophilic"\n')
f.write('@TYPE xy \n')
f.write("".join(sasa))
f.close()
In [ ]:
!xmgrace sasa.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='sasa.png')
In [ ]:
data_sasa=np.loadtxt('solven-accessible-surface.xvg',comments=['#', '@'])
#Engrosar marco
fig=pl.figure(figsize=(20, 12), dpi=100, linewidth=3.0)
ax = fig.add_subplot(111)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
#ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
pl.xlabel("Time (ns)", fontsize = 40)
pl.ylabel('Area ( nm'+ r'$\ ^2$'+')' , fontsize = 40)
#pl.title('Solvent Accessible Surface', fontsize=50)
pl.xticks(fontsize=30)
pl.yticks(fontsize=30)
dato=data_sasa[:,2]
dato2=dato[0]
if (dato2>0):
pl.plot(data_sasa[:,0]/1000, data_sasa[:,1], linewidth = 2, color='black')
pl.plot(data_sasa[:,0]/1000, data_sasa[:,2], linewidth = 2, color='red')
else:
pl.plot(data_sasa[:,0]/1000, data_sasa[:,1], linewidth = 2, color='black')
In [ ]:
### Creando el directorio para el análisis del SASA
### NOTA: se calcula con gromacs4 ya que arroja bien los resultados comparado con gromacs5
ruta_m_rmsd = nuevaruta+'/matriz'
print ( ruta_m_rmsd )
if not os.path.exists(ruta_m_rmsd):
os.makedirs(ruta_m_rmsd)
print ('Se ha creado la ruta ===>',ruta_m_rmsd)
else:
print ("La ruta "+ruta_m_rmsd+" existe..!!!")
print ( 'Nos vamos a ....', ruta_m_rmsd )
os.chdir( ruta_m_rmsd )
In [ ]:
print ('\nCopiando el archivo rmsd_matrix.tcl a '+ruta_m_rmsd)
source_file=ruta_scripts+'/rmsd_matriz/rmsd_matrix.tcl'
dest_file=ruta_m_rmsd+'/rmsd_matrix.tcl'
shutil.copy(source_file,dest_file)
#print ( 'Nos vemos a ....', ruta_old_traj )
#os.chdir( ruta_old_traj )
file_dcd=ruta_old_traj+'/'+dcd_file
file_psf=ruta_old_traj+'/'+psf_file
print (file_dcd)
print ('\nEjecutando CATDCD para obtener 100 frames de la trayectoria original....')
output_catdcd=!catdcd -o 100.dcd -stride 50 $file_dcd
print (output_catdcd.n)
In [ ]:
#Arrancando VMD para cargar el script rmsd_matrix.tcl
!vmd 100.dcd $file_psf
In [ ]:
ruta_matriz=os.getcwd()
if os.path.isfile('salida.dat'):
print ('El archivo salida.dat existe')
else:
print ('El archivo salida.dat no existe.. ejecutar desde MATRIZ DE RMSD...')
In [ ]:
#Creando el gráfico
data_matriz=np.loadtxt('salida.dat',comments=['#', '@'])
print(data_matriz.shape)
pl.figure(figsize=(20, 12), dpi=100)
imgplot = pl.imshow(data_matriz, origin='lower', cmap=pl.cm.Greens, interpolation='nearest')
#imgplot = pl.imshow(data_matriz, origin='lower', cmap=pl.cm.coolwarm, interpolation='nearest')
pl.xlabel("Time (ns)", fontsize = 30)
pl.ylabel('Time (ns)', fontsize = 30)
#pl.suptitle('RMSD', fontsize=50)
#pl.title('C-Alpha RMSD matrix', fontsize=40)
pl.xticks(fontsize=20)
pl.yticks(fontsize=20)
pl.xlim(0, 100)
pl.ylim(0, 100)
pl.colorbar()
In [ ]:
### Creando el directorio para el análisis del RMSF
#Verificando que exista la nueva carpeta para la conversión de trayectorias
ruta_matriz_dm = nuevaruta+'/matriz_dm'
print ( ruta_matriz_dm )
if not os.path.exists(ruta_matriz_dm):
os.makedirs(ruta_matriz_dm)
print ('Se ha creado la ruta ===>',ruta_matriz_dm)
else:
print ("La ruta "+ruta_matriz_dm+" existe..!!!")
print ( 'Nos vamos a ....', ruta_matriz_dm )
os.chdir( ruta_matriz_dm )
In [ ]:
!echo 4 | g_mdmat -f ../output.xtc -s ../ionized.pdb -mean average -frames frames -dt 10000
In [ ]:
!xpm2ps -f frames.xpm -o frames.eps
!xpm2ps -f average.xpm -o average.eps
print('\nConvirtiendo a png...')
!convert -density 600 frames.eps -resize 1024x1024 frames.png
!convert -density 600 average.eps -resize 1024x1024 average.png
In [ ]:
print ('Cargando el archivo average...')
Image(filename='average.png', width=800)
Para el cálculo de la energía libre se requiere el valor mínimo y máximo del RMSD y del radio de gyro, así como el valor de la temperatura a la cual se realizó la simulación. Estos datos son de entrada para el script del cálculo del mismo.
In [ ]:
### Creando el directorio para el análisis de la libre energía
ruta_f_energy = nuevaruta+'/free_energy'
print ( ruta_f_energy )
if not os.path.exists(ruta_f_energy):
os.makedirs(ruta_f_energy)
print ('Se ha creado la ruta ===>',ruta_f_energy)
else:
print ("La ruta "+ruta_f_energy+" existe..!!!")
print ( 'Nos vamos a ....', ruta_f_energy )
os.chdir( ruta_f_energy )
In [ ]:
#Solicita la temperatura
t=input('Temperatura a la cual se realizó la simulación:')
temperatura=int(t)
print ('Temperatura=>',temperatura)
In [ ]:
print ('Ejecutando el análisis de rmsd...')
!echo 3 3 | g_rms -f ../output.xtc -s ../ionized.pdb -a avgrp.xvg
print ('Ejecutando el análisis de rgyro...')
!echo 3 | g_gyrate -f ../output.xtc -s ../ionized.pdb -o gyrate.xvg
In [ ]:
print ('\nCopiando el archivo generateFES.py a '+ruta_f_energy)
source_file=ruta_scripts+'/free_energy/generateFES.py'
dest_file=ruta_f_energy+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
In [ ]:
#Cargando valores del RMSD
data_rmsd=np.loadtxt('rmsd.xvg',comments=['#', '@'])
#Cargnaod valores del R-GYRO
data_rgyro=np.loadtxt('gyrate.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del rmsd
min_rmsd=np.amin(data_rmsd[:,1])
max_rmsd=np.amax(data_rmsd[:,1])
print ('Minimo RMSD=>',min_rmsd)
print ('Máximo RMSD=>',max_rmsd)
#Obteniendo los valores máximo y mínimo del r-gyro
min_rgyro=np.amin(data_rgyro[:,1])
max_rgyro=np.amax(data_rgyro[:,1])
print ('Minimo RGYRO=>',min_rgyro)
print ('Máximo RGYRO=>',max_rgyro)
#Creando los archivos de entrada para el script
np.savetxt('rmsd.dat',data_rmsd[:,1], fmt='%1.7f')
np.savetxt('rgyro.dat',data_rgyro[:,1], fmt='%1.7f')
!paste rgyro.dat rmsd.dat > fes.dat
#Ejecutando el script de FES
!python generateFES.py fes.dat $min_rgyro $max_rgyro $min_rmsd $max_rmsd 200 200 $temperatura FEES.dat
#Cargando el archivo generado para plotear con matplotlib
data_fes=np.loadtxt('FEES.dat',comments=['#', '@'])
In [ ]:
# This loads the magics for gnuplot
%load_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "free_energy.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "Rg (nm)
set ylabel "RMSD (nm)"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "FEES.dat" with pm3d
In [ ]:
### Creando el directorio para el análisis del PCA
ruta_pca = nuevaruta+'/pca'
print ( ruta_pca )
if not os.path.exists(ruta_pca):
os.makedirs(ruta_pca)
print ('Se ha creado la ruta ===>',ruta_pca)
else:
print ("La ruta "+ruta_pca+" existe..!!!")
print ( 'Nos vamos a ....', ruta_pca )
os.chdir( ruta_pca )
In [ ]:
#Calculando matriz de covarianza
!echo 1 1 | g_covar -s ../ionized.pdb -f ../output.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm
In [ ]:
!echo 1 1 | g_anaeig -s ../ionized.pdb -f ../output.xtc -v eigenvectors.trr -eig eigenvalues.xvg -first 1 -last 2 -2d 2dproj_1_2.xvg
In [ ]:
#pcaX, pcaY=np.loadtxt('2dproj_1_2.xvg',comments=['#', '@'], unpack=True)
data_pca=np.loadtxt('2dproj_1_2.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del pca
min_pcaX=np.amin(data_pca[:,0])
max_pcaX=np.amax(data_pca[:,0])
print ('Minimo PCA_X=>',min_pcaX)
print ('Máximo PCA_X=>',max_pcaX)
min_pcaY=np.amin(data_pca[:,1])
max_pcaY=np.amax(data_pca[:,1])
print ('Minimo PCA_Y=>',min_pcaY)
print ('Máximo PCA_Y=>',max_pcaY)
#Creando los archivos de entrada para el script
np.savetxt('PCA.dat',data_pca, fmt='%1.5f')
#Copiando el script generateFES de la carpeta Free_energy
print ('\nCopiando el archivo generateFES.py a '+ruta_pca+ ' desde '+ ruta_f_energy)
source_file=ruta_f_energy+'/generateFES.py'
dest_file=ruta_pca+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Ejecutando el script de FES
!python generateFES.py PCA.dat $min_pcaX $max_pcaX $min_pcaY $max_pcaY 200 200 $temperatura FEES_PCA.dat
In [ ]:
#Volver a cargar el kernel de gnuplot para limpiar su buffer
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "pca.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "projection on eigenvector 1 (nm)"
set ylabel "projection on eigenvector 2 (nm)"
set title " "
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "FEES_PCA.dat" with pm3d
In [ ]:
from htmd import *
In [ ]:
### Creando el directorio para el análisis de los RMSD de los puentes
ruta_rmsd_diedros = nuevaruta+'/rmsd_diedros'
print ( ruta_rmsd_diedros )
if not os.path.exists(ruta_rmsd_diedros):
os.makedirs(ruta_rmsd_diedros)
print ('Se ha creado la ruta ===>',ruta_rmsd_diedros)
else:
print ("La ruta "+ruta_rmsd_diedros+" existe..!!!")
print ( 'Nos vamos a ....', ruta_rmsd_diedros)
os.chdir( ruta_rmsd_diedros )
Para este análisis se deberá revisar el archivo psf_charmm.tcl de la carpeta 01_BUILD, en el cual se tiene la definición de los puentes como la siguiente:
patch DISU A:4 A:22
patch DISU A:8 A:18
El número del puente se determinará de acuerdo al orden en que se encuentran definidos en este archivo, por ejemplo, la nota anterior:
DB1 4-22
DB2 8-18
La entrada de datos será por los índices del lado izquierdo y derecho respectivamente, con los cuales se creará la estructura completa de cada uno de ellos tomando los valores de los indices para su respectivo análisis.
In [ ]:
# Cargando la molécula
mol = Molecule('../ionized.pdb')
# Solicitando los datos de entrada
px1l=input('Índice del DB1 izquierdo:')
px1r=input('Índice del DB1 derecho:')
px2l=input('Índice del DB2 izquierdo:')
px2r=input('Índice del DB2 derecho:')
revisa1=1
revisa2=1
In [ ]:
if (revisa1>0):
#Obteniendo lado izquierdo del DB1
x1l_name=mol.get('name','resname CYS and noh and resid '+px1l)
x1l_index=mol.get('index','resname CYS and noh and resid '+px1l)
x1l_resid=mol.get('resid','resname CYS and noh and resid '+px1l)
#Obteniendo lado derecho del DB1
x1r_name=mol.get('name','resname CYS and noh and resid '+px1r)
x1r_index=mol.get('index','resname CYS and noh and resid '+px1r)
x1r_resid=mol.get('resid','resname CYS and noh and resid '+px1r)
if (revisa2>0):
#Obteniendo el lado izquierdo del DB2
x2l_name=mol.get('name','resname CYS and noh and resid '+px2l)
x2l_index=mol.get('index','resname CYS and noh and resid '+px2l)
x2l_resid=mol.get('resid','resname CYS and noh and resid '+px2l)
#Obteniendo el lado derecho del DB2
x2r_name=mol.get('name','resname CYS and noh and resid '+px2r)
x2r_index=mol.get('index','resname CYS and noh and resid '+px2r)
x2r_resid=mol.get('resid','resname CYS and noh and resid '+px2r)
#Obteniendo la lista de índices de los puentes
print ('Generando la lista de los índices para enviarlos')
db1x1l=[]
db1x2l=[]
db1x3m=[]
db1x2r=[]
db1x1r=[]
db1l_name_l=[]
db1l_index_l=[]
db1r_name_l=[]
db1r_index_l=[]
db2l_name_l=[]
db2l_index_l=[]
db2r_name_l=[]
db2r_index_l=[]
db3l_name_l=[]
db3l_index_l=[]
db3r_name_l=[]
db3r_index_l=[]
if (revisa1>0):
#Obteniendo los índices del DB1
for i in range(len(x1l_name)):
if (x1l_name[i]=='N' or x1l_name[i]=='CA' or x1l_name[i]=='CB' or x1l_name[i]=='SG'):
db1l_name_l.append(str(x1l_name[i]))
db1l_index_l.append(str(x1l_index[i]))
for i in range(len(x1r_name)):
if (x1r_name[i]=='N' or x1r_name[i]=='CA' or x1r_name[i]=='CB' or x1r_name[i]=='SG'):
db1r_name_l.append(str(x1r_name[i]))
db1r_index_l.append(str(x1r_index[i]))
print ('DB1 X1L =>',db1l_name_l)
print (db1l_index_l)
print ('DB1 X1R =>',db1r_name_l)
print (db1r_index_l)
if (revisa2>0):
#Obteniendo los índices del DB2
for i in range(len(x2l_name)):
if (x2l_name[i]=='N' or x2l_name[i]=='CA' or x2l_name[i]=='CB' or x2l_name[i]=='SG'):
db2l_name_l.append(str(x2l_name[i]))
db2l_index_l.append(str(x2l_index[i]))
for i in range(len(x2r_name)):
if (x2r_name[i]=='N' or x2r_name[i]=='CA' or x2r_name[i]=='CB' or x2r_name[i]=='SG'):
db2r_name_l.append(str(x2r_name[i]))
db2r_index_l.append(str(x2r_index[i]))
print ('DB2 X1L =>',db2l_name_l)
print (db2l_index_l)
print ('DB2 X1R =>',db2r_name_l)
print (db2r_index_l)
In [ ]:
#Generando el DB1 completo ordenado
filas=8
col=2
DB1_i=[]
DB1_N=[]
DB2_i=[]
DB2_N=[]
DB3_i=[]
DB3_N=[]
for i in range(0,filas):
DB1_N.append([' '])
DB1_i.append(['0'])
DB2_N.append([' '])
DB2_i.append(['0'])
DB3_N.append([' '])
DB3_i.append(['0'])
if (revisa1>0):
#Cargando índices para el puente 1
for i in range(len(db1l_name_l)):
if db1l_name_l[i]=='N':
DB1_N[0] = db1l_name_l[i]
DB1_i[0]='index '+db1l_index_l[i]
if db1l_name_l[i]=='CA':
DB1_N[1] = db1l_name_l[i]
DB1_i[1]='index '+db1l_index_l[i]
if db1l_name_l[i]=='CB':
DB1_N[2] = db1l_name_l[i]
DB1_i[2]='index '+db1l_index_l[i]
if db1l_name_l[i]=='SG':
DB1_N[3] = db1l_name_l[i]
DB1_i[3]='index '+db1l_index_l[i]
for i in range(len(db1r_name_l)):
if db1r_name_l[i]=='SG':
DB1_N[4] = db1r_name_l[i]
DB1_i[4]='index '+db1r_index_l[i]
if db1r_name_l[i]=='CB':
DB1_N[5] = db1r_name_l[i]
DB1_i[5]='index '+db1r_index_l[i]
if db1r_name_l[i]=='CA':
DB1_N[6] = db1r_name_l[i]
DB1_i[6]='index '+db1r_index_l[i]
if db1r_name_l[i]=='N':
DB1_N[7] = db1r_name_l[i]
DB1_i[7]='index '+db1r_index_l[i]
print ('Puente DB1 = resid '+px1l+':'+px1r)
print ('Names DB1=>',DB1_i)
print ('Index DB1=>',DB1_N)
print ('\n')
if (revisa2>0):
#Cargando índices para el puente 2
for i in range(len(db2l_name_l)):
if db2l_name_l[i]=='N':
DB2_N[0] = db2l_name_l[i]
DB2_i[0]='index '+db2l_index_l[i]
if db2l_name_l[i]=='CA':
DB2_N[1] = db2l_name_l[i]
DB2_i[1]='index '+db2l_index_l[i]
if db2l_name_l[i]=='CB':
DB2_N[2] = db2l_name_l[i]
DB2_i[2]='index '+db2l_index_l[i]
if db2l_name_l[i]=='SG':
DB2_N[3] = db2l_name_l[i]
DB2_i[3]='index '+db2l_index_l[i]
for i in range(len(db2r_name_l)):
if db2r_name_l[i]=='SG':
DB2_N[4] = db2r_name_l[i]
DB2_i[4]='index '+db2r_index_l[i]
if db2r_name_l[i]=='CB':
DB2_N[5] = db2r_name_l[i]
DB2_i[5]='index '+db2r_index_l[i]
if db2r_name_l[i]=='CA':
DB2_N[6] = db2r_name_l[i]
DB2_i[6]='index '+db2r_index_l[i]
if db2r_name_l[i]=='N':
DB2_N[7] = db2r_name_l[i]
DB2_i[7]='index '+db2r_index_l[i]
print ('Puente DB2 = resid '+px2l+':'+px2r)
print ('Names DB2=>',DB2_i)
print ('Index DB2=>',DB2_N)
print ('\n')
In [ ]:
if (revisa1>0):
#Creando script para DB1_x1l
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f1 = open('DB1_x1l.tcl', 'w')
print(f1)
f1.write('set psfFile '+ psf+' \n')
f1.write('set dcdFile '+ dcd+' \n')
f1.write('\nmol load psf $psfFile dcd $dcdFile\n')
f1.write('set outfile ' +'[open ' +'db1_x1l.dat'+' w]\n')
f1.write('set nf [molinfo top get numframes]\n')
f1.write('\n#RMSD calculation loop\n')
f1.write('set f1 [atomselect top "'+DB1_i[0]+' or '+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' " frame 0]\n')
f1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f1.write(' set sel [atomselect top "'+DB1_i[0]+' or '+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' " frame $i]\n')
f1.write(' $sel move [measure fit $sel $f1]\n')
f1.write(' set time [expr {$i +1}]\n')
f1.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f1.write(' puts $outfile " $time"\n')
f1.write('}\n')
f1.write('close $outfile')
f1.close()
#Creando script para DB1_x2l
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f2 = open('DB1_x2l.tcl', 'w')
print(f2)
f2.write('set psfFile '+ psf+' \n')
f2.write('set dcdFile '+ dcd+' \n')
f2.write('\nmol load psf $psfFile dcd $dcdFile\n')
f2.write('set outfile ' +'[open ' +'db1_x2l.dat'+' w]\n')
f2.write('set nf [molinfo top get numframes]\n')
f2.write('\n#RMSD calculation loop\n')
f2.write('set f1 [atomselect top "'+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' " frame 0]\n')
f2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f2.write(' set sel [atomselect top "'+DB1_i[1]+' or '+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' " frame $i]\n')
f2.write(' $sel move [measure fit $sel $f1]\n')
f2.write(' set time [expr {$i +1}]\n')
f2.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f2.write(' puts $outfile " $time"\n')
f2.write('}\n')
f2.write('close $outfile')
f2.close()
#Creando script para DB1_x3m
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f3 = open('DB1_x3m.tcl', 'w')
print(f3)
f3.write('set psfFile '+ psf+' \n')
f3.write('set dcdFile '+ dcd+' \n')
f3.write('\nmol load psf $psfFile dcd $dcdFile\n')
f3.write('set outfile ' +'[open ' +'db1_x3m.dat'+' w]\n')
f3.write('set nf [molinfo top get numframes]\n')
f3.write('\n#RMSD calculation loop\n')
f3.write('set f1 [atomselect top "'+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' " frame 0]\n')
f3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f3.write(' set sel [atomselect top "'+DB1_i[2]+' or '+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' " frame $i]\n')
f3.write(' $sel move [measure fit $sel $f1]\n')
f3.write(' set time [expr {$i +1}]\n')
f3.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f3.write(' puts $outfile " $time"\n')
f3.write('}\n')
f3.write('close $outfile')
f3.close()
#Creando script para DB1_x2r
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f4 = open('DB1_x2r.tcl', 'w')
print(f4)
f4.write('set psfFile '+ psf+' \n')
f4.write('set dcdFile '+ dcd+' \n')
f4.write('\nmol load psf $psfFile dcd $dcdFile\n')
f4.write('set outfile ' +'[open ' +'db1_x2r.dat'+' w]\n')
f4.write('set nf [molinfo top get numframes]\n')
f4.write('\n#RMSD calculation loop\n')
f4.write('set f1 [atomselect top "'+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' " frame 0]\n')
f4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f4.write(' set sel [atomselect top "'+DB1_i[3]+' or '+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' " frame $i]\n')
f4.write(' $sel move [measure fit $sel $f1]\n')
f4.write(' set time [expr {$i +1}]\n')
f4.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f4.write(' puts $outfile " $time"\n')
f4.write('}\n')
f4.write('close $outfile')
f4.close()
#Creando script para DB1_x1r
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f5 = open('DB1_x1r.tcl', 'w')
print(f5)
f5.write('set psfFile '+ psf+' \n')
f5.write('set dcdFile '+ dcd+' \n')
f5.write('\nmol load psf $psfFile dcd $dcdFile\n')
f5.write('set outfile ' +'[open ' +'db1_x1r.dat'+' w]\n')
f5.write('set nf [molinfo top get numframes]\n')
f5.write('\n#RMSD calculation loop\n')
f5.write('set f1 [atomselect top "'+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' or '+DB1_i[7]+' " frame 0]\n')
f5.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f5.write(' set sel [atomselect top "'+DB1_i[4]+' or '+DB1_i[5]+' or '+DB1_i[6]+' or '+DB1_i[7]+' " frame $i]\n')
f5.write(' $sel move [measure fit $sel $f1]\n')
f5.write(' set time [expr {$i +1}]\n')
f5.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f5.write(' puts $outfile " $time"\n')
f5.write('}\n')
f5.write('close $outfile')
f5.close()
if (revisa2>0):
##########################################################################################
## Creando los archivos para DB2
#######################################################################################
#Creando script para DB2_x1l
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f6 = open('DB2_x1l.tcl', 'w')
print(f6)
f6.write('set psfFile '+ psf+' \n')
f6.write('set dcdFile '+ dcd+' \n')
f6.write('\nmol load psf $psfFile dcd $dcdFile\n')
f6.write('set outfile ' +'[open ' +'db2_x1l.dat'+' w]\n')
f6.write('set nf [molinfo top get numframes]\n')
f6.write('\n#RMSD calculation loop\n')
f6.write('set f1 [atomselect top "'+DB2_i[0]+' or '+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' " frame 0]\n')
f6.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f6.write(' set sel [atomselect top "'+DB2_i[0]+' or '+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' " frame $i]\n')
f6.write(' $sel move [measure fit $sel $f1]\n')
f6.write(' set time [expr {$i +1}]\n')
f6.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f6.write(' puts $outfile " $time"\n')
f6.write('}\n')
f6.write('close $outfile')
f6.close()
#Creando script para DB1_x2l
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f7 = open('DB2_x2l.tcl', 'w')
print(f7)
f7.write('set psfFile '+ psf+' \n')
f7.write('set dcdFile '+ dcd+' \n')
f7.write('\nmol load psf $psfFile dcd $dcdFile\n')
f7.write('set outfile ' +'[open ' +'db2_x2l.dat'+' w]\n')
f7.write('set nf [molinfo top get numframes]\n')
f7.write('\n#RMSD calculation loop\n')
f7.write('set f1 [atomselect top "'+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' " frame 0]\n')
f7.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f7.write(' set sel [atomselect top "'+DB2_i[1]+' or '+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' " frame $i]\n')
f7.write(' $sel move [measure fit $sel $f1]\n')
f7.write(' set time [expr {$i +1}]\n')
f7.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f7.write(' puts $outfile " $time"\n')
f7.write('}\n')
f7.write('close $outfile')
f7.close()
#Creando script para DB1_x3m
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f8 = open('DB2_x3m.tcl', 'w')
print(f8)
f8.write('set psfFile '+ psf+' \n')
f8.write('set dcdFile '+ dcd+' \n')
f8.write('\nmol load psf $psfFile dcd $dcdFile\n')
f8.write('set outfile ' +'[open ' +'db2_x3m.dat'+' w]\n')
f8.write('set nf [molinfo top get numframes]\n')
f8.write('\n#RMSD calculation loop\n')
f8.write('set f1 [atomselect top "'+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' " frame 0]\n')
f8.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f8.write(' set sel [atomselect top "'+DB2_i[2]+' or '+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' " frame $i]\n')
f8.write(' $sel move [measure fit $sel $f1]\n')
f8.write(' set time [expr {$i +1}]\n')
f8.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f8.write(' puts $outfile " $time"\n')
f8.write('}\n')
f8.write('close $outfile')
f8.close()
#Creando script para DB1_x2r
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f9 = open('DB2_x2r.tcl', 'w')
print(f9)
f9.write('set psfFile '+ psf+' \n')
f9.write('set dcdFile '+ dcd+' \n')
f9.write('\nmol load psf $psfFile dcd $dcdFile\n')
f9.write('set outfile ' +'[open ' +'db2_x2r.dat'+' w]\n')
f9.write('set nf [molinfo top get numframes]\n')
f9.write('\n#RMSD calculation loop\n')
f9.write('set f1 [atomselect top "'+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' " frame 0]\n')
f9.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f9.write(' set sel [atomselect top "'+DB2_i[3]+' or '+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' " frame $i]\n')
f9.write(' $sel move [measure fit $sel $f1]\n')
f9.write(' set time [expr {$i +1}]\n')
f9.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f9.write(' puts $outfile " $time"\n')
f9.write('}\n')
f9.write('close $outfile')
f9.close()
#Creando script para DB1_x1r
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print(psf)
f10 = open('DB2_x1r.tcl', 'w')
print(f10)
f10.write('set psfFile '+ psf+' \n')
f10.write('set dcdFile '+ dcd+' \n')
f10.write('\nmol load psf $psfFile dcd $dcdFile\n')
f10.write('set outfile ' +'[open ' +'db2_x1r.dat'+' w]\n')
f10.write('set nf [molinfo top get numframes]\n')
f10.write('\n#RMSD calculation loop\n')
f10.write('set f1 [atomselect top "'+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' or '+DB2_i[7]+' " frame 0]\n')
f10.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
f10.write(' set sel [atomselect top "'+DB2_i[4]+' or '+DB2_i[5]+' or '+DB2_i[6]+' or '+DB2_i[7]+' " frame $i]\n')
f10.write(' $sel move [measure fit $sel $f1]\n')
f10.write(' set time [expr {$i +1}]\n')
f10.write(' puts -nonewline $outfile "[measure rmsd $sel $f1]"\n')
f10.write(' puts $outfile " $time"\n')
f10.write('}\n')
f10.write('close $outfile')
f10.close()
In [ ]:
if (revisa1>0):
#Calculando con VMD rmsd DB1 X1L
!vmd -dispdev text < DB1_x1l.tcl
#Calculando con VMD DB1 X2L
!vmd -dispdev text < DB1_x2l.tcl
#Calculando con VMD DB1 X3M
!vmd -dispdev text < DB1_x3m.tcl
#Calculando con VMD DB1 X2R
!vmd -dispdev text < DB1_x2r.tcl
#Calculando con VMD DB1 X1R
!vmd -dispdev text < DB1_x1r.tcl
if (revisa2>0):
#Calculando con VMD rmsd DB2 X1L
!vmd -dispdev text < DB2_x1l.tcl
#Calculando con VMD DB2 X2L
!vmd -dispdev text < DB2_x2l.tcl
#Calculando con VMD DB2 X3M
!vmd -dispdev text < DB2_x3m.tcl
#Calculando con VMD DB2 X2R
!vmd -dispdev text < DB2_x2r.tcl
#Calculando con VMD DB2 X1R
!vmd -dispdev text < DB2_x1r.tcl
In [ ]:
escale_y=[]
fig = pl.figure(figsize=(25,8))
fig.subplots_adjust(hspace=.4, wspace=0.3)
#Formateando los valores de los ejes
#Engrosando marcos
ax = fig.add_subplot(2,5,1)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,2)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,3)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,4)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,5)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,6)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,7)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,8)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,9)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax = fig.add_subplot(2,5,10)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
if (revisa1>0):
#Datos de DB1
data_db1_x1l=np.loadtxt('db1_x1l.dat',comments=['#', '@'])
data_db1_x2l=np.loadtxt('db1_x2l.dat',comments=['#', '@'])
data_db1_x3m=np.loadtxt('db1_x3m.dat',comments=['#', '@'])
data_db1_x2r=np.loadtxt('db1_x2r.dat',comments=['#', '@'])
data_db1_x1r=np.loadtxt('db1_x1r.dat',comments=['#', '@'])
sub1 = fig.add_subplot(251) # instead of plt.subplot(2, 2, 1)
#sub1.set_title('DB1_X1L')
sub1.set_xlabel('Time (ns)')
sub1.set_ylabel('RMSD (nm)')
sub1.plot(data_db1_x1l[:,1]*0.02, data_db1_x1l[:,0]/10, color='black', linewidth = 1, label='DB1_X1L')
x1,x2,y1,y2=sub1.axis()
escale_y.append(y2)
sub2 = fig.add_subplot(252)
#sub2.set_title('DB1_X2L')
sub2.set_xlabel('Time (ns)')
sub2.set_ylabel('RMSD (nm)')
sub2.plot(data_db1_x2l[:,1]*0.02, data_db1_x2l[:,0]/10, color='black', linewidth = 1, label='DB1_X2L')
x1,x2,y1,y2=sub2.axis()
escale_y.append(y2)
sub3 = fig.add_subplot(253)
#sub3.set_title('DB1_X3M')
sub3.set_xlabel('Time (ns)')
sub3.set_ylabel('RMSD (nm)')
sub3.plot(data_db1_x3m[:,1]*0.02, data_db1_x3m[:,0]/10, color='black', linewidth = 1, label='DB1_X3M')
x1,x2,y1,y2=sub3.axis()
escale_y.append(y2)
sub4 = fig.add_subplot(254)
#sub4.set_title('DB1_X2R')
sub4.set_xlabel('Time (ns)')
sub4.set_ylabel('RMSD (nm)')
sub4.plot(data_db1_x2r[:,1]*0.02, data_db1_x2r[:,0]/10, color='black', linewidth = 1, label='DB1_X2R')
x1,x2,y1,y2=sub4.axis()
escale_y.append(y2)
sub5 = fig.add_subplot(255)
#sub5.set_title('DB1_X1R')
sub5.set_xlabel('Time (ns)')
sub5.set_ylabel('RMSD (nm)')
sub5.plot(data_db1_x1r[:,1]*0.02, data_db1_x1r[:,0]/10, color='black', linewidth = 1, label='DB1_X1R')
x1,x2,y1,y2=sub5.axis()
escale_y.append(y2)
if (revisa2>0):
#DAtos de DB2
data_db2_x1l=np.loadtxt('db2_x1l.dat',comments=['#', '@'])
data_db2_x2l=np.loadtxt('db2_x2l.dat',comments=['#', '@'])
data_db2_x3m=np.loadtxt('db2_x3m.dat',comments=['#', '@'])
data_db2_x2r=np.loadtxt('db2_x2r.dat',comments=['#', '@'])
data_db2_x1r=np.loadtxt('db2_x1r.dat',comments=['#', '@'])
#Ploteando DB2
sub6 = fig.add_subplot(256)
#sub6.set_title('DB2_X1L')
sub6.set_xlabel('Time (ns)')
sub6.set_ylabel('RMSD (nm)')
sub6.plot(data_db2_x1l[:,1]*0.02, data_db2_x1l[:,0]/10, color='red', linewidth = 1, label='DB2_X1L')
x1,x2,y1,y2=sub6.axis()
escale_y.append(y2)
sub7 = fig.add_subplot(257)
#sub7.set_title('DB2_X2L')
sub7.set_xlabel('Time (ns)')
sub7.set_ylabel('RMSD (nm)')
sub7.plot(data_db2_x2l[:,1]*0.02, data_db2_x2l[:,0]/10, color='red', linewidth = 1, label='DB2_X2L')
x1,x2,y1,y2=sub7.axis()
escale_y.append(y2)
sub8 = fig.add_subplot(258)
#sub8.set_title('DB2_X3M')
sub8.set_xlabel('Time (ns)')
sub8.set_ylabel('RMSD (nm)')
sub8.plot(data_db2_x3m[:,1]*0.02, data_db2_x3m[:,0]/10, color='red', linewidth = 1, label='DB2_X3M')
x1,x2,y1,y2=sub8.axis()
escale_y.append(y2)
sub9 = fig.add_subplot(259)
#sub9.set_title('DB2_X2R')
sub9.set_xlabel('Time (ns)')
sub9.set_ylabel('RMSD (nm)')
sub9.plot(data_db2_x2r[:,1]*0.02, data_db2_x2r[:,0]/10, color='red', linewidth = 1, label='DB2_X2R')
x1,x2,y1,y2=sub9.axis()
escale_y.append(y2)
sub10 = fig.add_subplot(2,5,10)
#sub10.set_title('DB2_X1R')
sub10.set_xlabel('Time (ns)')
sub10.set_ylabel('RMSD (nm)')
sub10.plot(data_db2_x1r[:,1]*0.02, data_db2_x1r[:,0]/10, color='red', linewidth = 1, label='DB2_X1R')
x1,x2,y1,y2=sub10.axis()
escale_y.append(y2)
#escale_y
escale_y.sort(reverse=True)
escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,escale_y[0]))
sub2.axis((x1,x2,y1,escale_y[0]))
sub3.axis((x1,x2,y1,escale_y[0]))
sub4.axis((x1,x2,y1,escale_y[0]))
sub5.axis((x1,x2,y1,escale_y[0]))
sub6.axis((x1,x2,y1,escale_y[0]))
sub7.axis((x1,x2,y1,escale_y[0]))
sub8.axis((x1,x2,y1,escale_y[0]))
sub9.axis((x1,x2,y1,escale_y[0]))
sub10.axis((x1,x2,y1,escale_y[0]))
In [ ]:
### Creando el directorio para el análisis de las distancias de enlace de los puentes
ruta_diedros = nuevaruta+'/diedros_intra'
print ( ruta_diedros )
if not os.path.exists(ruta_diedros):
os.makedirs(ruta_diedros)
print ('Se ha creado la ruta ===>',ruta_diedros)
else:
print ("La ruta "+ruta_diedros+" existe..!!!")
print ( 'Nos vamos a ....', ruta_diedros)
os.chdir( ruta_diedros )
In [ ]:
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
if (revisa1>0):
#Creando script para DB1_x1l
d1 = open('dihed_DB1_x1l.tcl', 'w')
print(d1)
d1.write('set psfFile '+ psf+' \n')
d1.write('set dcdFile '+ dcd+' \n')
d1.write('\nmol load psf $psfFile dcd $dcdFile\n')
d1.write('set outfile ' +'[open ' +'dihed_db1_x1l.dat'+' w]\n')
d1.write('set nf [molinfo top get numframes]\n')
d1.write(' \n')
d1.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[0]+'"] get index]\n')
d1.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
d1.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
d1.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
d1.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d1.write(' set x [measure dihed $dihed frame $i]\n')
d1.write(' set time [expr {$i +1}]\n')
d1.write(' puts $outfile "$time $x"\n')
d1.write('}\n')
d1.close()
#Creando script para DB1_x2l
d2 = open('dihed_DB1_x2l.tcl', 'w')
print(d2)
d2.write('set psfFile '+ psf+' \n')
d2.write('set dcdFile '+ dcd+' \n')
d2.write('\nmol load psf $psfFile dcd $dcdFile\n')
d2.write('set outfile ' +'[open ' +'dihed_db1_x2l.dat'+' w]\n')
d2.write('set nf [molinfo top get numframes]\n')
d2.write(' \n')
d2.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
d2.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
d2.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
d2.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
d2.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d2.write(' set x [measure dihed $dihed frame $i]\n')
d2.write(' set time [expr {$i +1}]\n')
d2.write(' puts $outfile "$time $x"\n')
d2.write('}\n')
d2.close()
#Creando script para DB1_x3m
d3 = open('dihed_DB1_x3m.tcl', 'w')
print(d3)
d3.write('set psfFile '+ psf+' \n')
d3.write('set dcdFile '+ dcd+' \n')
d3.write('\nmol load psf $psfFile dcd $dcdFile\n')
d3.write('set outfile ' +'[open ' +'dihed_db1_x3m.dat'+' w]\n')
d3.write('set nf [molinfo top get numframes]\n')
d3.write(' \n')
d3.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
d3.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
d3.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
d3.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
d3.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d3.write(' set x [measure dihed $dihed frame $i]\n')
d3.write(' set time [expr {$i +1}]\n')
d3.write(' puts $outfile "$time $x"\n')
d3.write('}\n')
d3.close()
#Creando script para DB1_x2r
d4 = open('dihed_DB1_x2r.tcl', 'w')
print(d4)
d4.write('set psfFile '+ psf+' \n')
d4.write('set dcdFile '+ dcd+' \n')
d4.write('\nmol load psf $psfFile dcd $dcdFile\n')
d4.write('set outfile ' +'[open ' +'dihed_db1_x2r.dat'+' w]\n')
d4.write('set nf [molinfo top get numframes]\n')
d4.write(' \n')
d4.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
d4.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
d4.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
d4.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
d4.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d4.write(' set x [measure dihed $dihed frame $i]\n')
d4.write(' set time [expr {$i +1}]\n')
d4.write(' puts $outfile "$time $x"\n')
d4.write('}\n')
d4.close()
#Creando script para DB1_x1r
d5 = open('dihed_DB1_x1r.tcl', 'w')
print(d5)
d5.write('set psfFile '+ psf+' \n')
d5.write('set dcdFile '+ dcd+' \n')
d5.write('\nmol load psf $psfFile dcd $dcdFile\n')
d5.write('set outfile ' +'[open ' +'dihed_db1_x1r.dat'+' w]\n')
d5.write('set nf [molinfo top get numframes]\n')
d5.write(' \n')
d5.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
d5.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
d5.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
d5.write('set selatoms4 [[atomselect top "protein and chain A and '+DB1_i[7]+'"] get index]\n')
d5.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d5.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d5.write(' set x [measure dihed $dihed frame $i]\n')
d5.write(' set time [expr {$i +1}]\n')
d5.write(' puts $outfile "$time $x"\n')
d5.write('}\n')
d5.close()
if (revisa2>0):
#####################################################################
########## Puente 2
##########################################3
#Creando script para DB2_x1l
d6 = open('dihed_DB2_x1l.tcl', 'w')
print(d6)
d6.write('set psfFile '+ psf+' \n')
d6.write('set dcdFile '+ dcd+' \n')
d6.write('\nmol load psf $psfFile dcd $dcdFile\n')
d6.write('set outfile ' +'[open ' +'dihed_db2_x1l.dat'+' w]\n')
d6.write('set nf [molinfo top get numframes]\n')
d6.write(' \n')
d6.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[0]+'"] get index]\n')
d6.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
d6.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
d6.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
d6.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d6.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d6.write(' set x [measure dihed $dihed frame $i]\n')
d6.write(' set time [expr {$i +1}]\n')
d6.write(' puts $outfile "$time $x"\n')
d6.write('}\n')
d6.close()
#Creando script para DB2_x2l
d7 = open('dihed_DB2_x2l.tcl', 'w')
print(d7)
d7.write('set psfFile '+ psf+' \n')
d7.write('set dcdFile '+ dcd+' \n')
d7.write('\nmol load psf $psfFile dcd $dcdFile\n')
d7.write('set outfile ' +'[open ' +'dihed_db2_x2l.dat'+' w]\n')
d7.write('set nf [molinfo top get numframes]\n')
d7.write(' \n')
d7.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
d7.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
d7.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
d7.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
d7.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d7.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d7.write(' set x [measure dihed $dihed frame $i]\n')
d7.write(' set time [expr {$i +1}]\n')
d7.write(' puts $outfile "$time $x"\n')
d7.write('}\n')
d7.close()
#Creando script para DB2_x3m
d8 = open('dihed_DB2_x3m.tcl', 'w')
print(d8)
d8.write('set psfFile '+ psf+' \n')
d8.write('set dcdFile '+ dcd+' \n')
d8.write('\nmol load psf $psfFile dcd $dcdFile\n')
d8.write('set outfile ' +'[open ' +'dihed_db2_x3m.dat'+' w]\n')
d8.write('set nf [molinfo top get numframes]\n')
d8.write(' \n')
d8.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
d8.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
d8.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
d8.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
d8.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d8.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d8.write(' set x [measure dihed $dihed frame $i]\n')
d8.write(' set time [expr {$i +1}]\n')
d8.write(' puts $outfile "$time $x"\n')
d8.write('}\n')
d8.close()
#Creando script para DB2_x2r
d9 = open('dihed_DB2_x2r.tcl', 'w')
print(d9)
d9.write('set psfFile '+ psf+' \n')
d9.write('set dcdFile '+ dcd+' \n')
d9.write('\nmol load psf $psfFile dcd $dcdFile\n')
d9.write('set outfile ' +'[open ' +'dihed_db2_x2r.dat'+' w]\n')
d9.write('set nf [molinfo top get numframes]\n')
d9.write(' \n')
d9.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
d9.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
d9.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
d9.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
d9.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d9.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d9.write(' set x [measure dihed $dihed frame $i]\n')
d9.write(' set time [expr {$i +1}]\n')
d9.write(' puts $outfile "$time $x"\n')
d9.write('}\n')
d9.close()
#Creando script para DB2_x1r
d10 = open('dihed_DB2_x1r.tcl', 'w')
print(d10)
d10.write('set psfFile '+ psf+' \n')
d10.write('set dcdFile '+ dcd+' \n')
d10.write('\nmol load psf $psfFile dcd $dcdFile\n')
d10.write('set outfile ' +'[open ' +'dihed_db2_x1r.dat'+' w]\n')
d10.write('set nf [molinfo top get numframes]\n')
d10.write(' \n')
d10.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
d10.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
d10.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
d10.write('set selatoms4 [[atomselect top "protein and chain A and '+DB2_i[7]+'"] get index]\n')
d10.write('set dihed [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] [lindex $selatoms4] ]\n')
d10.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
d10.write(' set x [measure dihed $dihed frame $i]\n')
d10.write(' set time [expr {$i +1}]\n')
d10.write(' puts $outfile "$time $x"\n')
d10.write('}\n')
d10.close()
In [ ]:
if (revisa1>0):
#Calculando con VMD rmsd DB1 X1L
!vmd -dispdev text < dihed_DB1_x1l.tcl
#Calculando con VMD DB1 X2L
!vmd -dispdev text < dihed_DB1_x2l.tcl
#Calculando con VMD DB1 X3M
!vmd -dispdev text < dihed_DB1_x3m.tcl
#Calculando con VMD DB1 X2R
!vmd -dispdev text < dihed_DB1_x2r.tcl
#Calculando con VMD DB1 X1R
!vmd -dispdev text < dihed_DB1_x1r.tcl
if (revisa2>0):
#Calculando con VMD rmsd DB2 X1L
!vmd -dispdev text < dihed_DB2_x1l.tcl
#Calculando con VMD DB2 X2L
!vmd -dispdev text < dihed_DB2_x2l.tcl
#Calculando con VMD DB2 X3M
!vmd -dispdev text < dihed_DB2_x3m.tcl
#Calculando con VMD DB2 X2R
!vmd -dispdev text < dihed_DB2_x2r.tcl
#Calculando con VMD DB2 X1R
!vmd -dispdev text < dihed_DB2_x1r.tcl
In [ ]:
print ('\nCopiando el archivo generateFES.py a '+ruta_diedros)
source_file=ruta_f_energy+'/generateFES.py'
dest_file=ruta_diedros+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
In [ ]:
if (revisa1>0):
#Cargando valores del DB1_X1L
data_db1_x1l=np.loadtxt('dihed_db1_x1l.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_db1_x1r=np.loadtxt('dihed_db1_x1r.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB1_X1L
min_x1l=np.amin(data_db1_x1l[:,1])
max_x1l=np.amax(data_db1_x1l[:,1])
print ('Minimo DB1_X1L=>',min_x1l)
print ('Máximo DB1_X1L=>',max_x1l)
#Obteniendo los valores máximo y mínimo del DB1_X1R
min_x1r=np.amin(data_db1_x1r[:,1])
max_x1r=np.amax(data_db1_x1r[:,1])
print ('Minimo DB1_X1R=>',min_x1r)
print ('Máximo DB1_X1R=>',max_x1r)
#Creando los archivos de entrada para el script
np.savetxt('db1_x1l.dat',data_db1_x1l[:,1], fmt='%1.14f')
np.savetxt('db1_x1r.dat',data_db1_x1r[:,1], fmt='%1.14f')
!paste db1_x1l.dat db1_x1r.dat > DB1_x1_lr.dat
#Ejecutando el script de FES
!python generateFES.py DB1_x1_lr.dat $min_x1l $max_x1l $min_x1r $max_x1r 200 200 $temperatura XL1_XR1.dat
###################################################################
#Cargando valores del DB1_X2l
data_db1_x2l=np.loadtxt('dihed_db1_x2l.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_db1_x2r=np.loadtxt('dihed_db1_x2r.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB1_X1L
min_x2l=np.amin(data_db1_x2l[:,1])
max_x2l=np.amax(data_db1_x2l[:,1])
print ('Minimo DB1_X2L=>',min_x2l)
print ('Máximo DB1_X2L=>',max_x2l)
#Obteniendo los valores máximo y mínimo del DB1_X1R
min_x2r=np.amin(data_db1_x2r[:,1])
max_x2r=np.amax(data_db1_x2r[:,1])
print ('Minimo DB1_X2R=>',min_x2r)
print ('Máximo DB1_X2R=>',max_x2r)
#Creando los archivos de entrada para el script
np.savetxt('db1_x2l.dat',data_db1_x2l[:,1], fmt='%1.14f')
np.savetxt('db1_x2r.dat',data_db1_x2r[:,1], fmt='%1.14f')
!paste db1_x2l.dat db1_x2r.dat > DB1_x2_lr.dat
#Ejecutando el script de FES
!python generateFES.py DB1_x2_lr.dat $min_x2l $max_x2l $min_x2r $max_x2r 200 200 $temperatura XL2_XR2.dat
######################################################################################
#Generando los archivos para X3M
data_db1_x3m=np.loadtxt('dihed_db1_x3m.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB1_X1L
min_x3m=np.amin(data_db1_x3m[:,1])
max_x3m=np.amax(data_db1_x3m[:,1])
print ('Minimo DB1_X3M=>',min_x3m)
print ('Máximo DB1_X3M=>',max_x3m)
print ('Minimo DB1_X1L=>',min_x1l)
print ('Máximo DB1_X1L=>',max_x1l)
print ('Minimo DB1_X2L=>',min_x2l)
print ('Máximo DB1_X2L=>',max_x2l)
print ('Minimo DB1_X1R=>',min_x1r)
print ('Máximo DB1_X1R=>',max_x1r)
print ('Minimo DB1_X2R=>',min_x2r)
print ('Máximo DB1_X2R=>',max_x2r)
#Creando los archivos de entrada para el script
np.savetxt('db1_x3m.dat',data_db1_x3m[:,1], fmt='%1.14f')
!paste db1_x3m.dat db1_x1l.dat > DB1_x3m_x1l.dat
!paste db1_x3m.dat db1_x2l.dat > DB1_x3m_x2l.dat
!paste db1_x3m.dat db1_x1r.dat > DB1_x3m_x1r.dat
!paste db1_x3m.dat db1_x2r.dat > DB1_x3m_x2r.dat
#Ejecutando el script de FES
!python generateFES.py DB1_x3m_x1l.dat $min_x3m $max_x3m $min_x1l $max_x1l 200 200 $temperatura XM3_XL1.dat
!python generateFES.py DB1_x3m_x2l.dat $min_x3m $max_x3m $min_x2l $max_x2l 200 200 $temperatura XM3_XL2.dat
!python generateFES.py DB1_x3m_x1r.dat $min_x3m $max_x3m $min_x1r $max_x1r 200 200 $temperatura XM3_XR1.dat
!python generateFES.py DB1_x3m_x2r.dat $min_x3m $max_x3m $min_x2r $max_x2r 200 200 $temperatura XM3_XR2.dat
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xl1_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_1}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XL1_XR1.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xl2_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_2}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XL2_XR2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xm3_vs_xl1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XL1.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xm3_vs_xl2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XL2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xm3_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XR2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_xm3_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "XM3_XR1.dat" with pm3d
In [ ]:
if (revisa2>0):
#Cargando valores del DB2_X1L
data_db2_x1l=np.loadtxt('dihed_db2_x1l.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_db2_x1r=np.loadtxt('dihed_db2_x1r.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB2_X1L
min_db2_x1l=np.amin(data_db2_x1l[:,1])
max_db2_x1l=np.amax(data_db2_x1l[:,1])
print ('Minimo DB2_X1L=>',min_db2_x1l)
print ('Máximo DB2_X1L=>',max_db2_x1l)
#Obteniendo los valores máximo y mínimo del DB2_X1R
min_db2_x1r=np.amin(data_db2_x1r[:,1])
max_db2_x1r=np.amax(data_db2_x1r[:,1])
print ('Minimo DB2_X1R=>',min_db2_x1r)
print ('Máximo DB2_X1R=>',max_db2_x1r)
#Creando los archivos de entrada para el script
np.savetxt('db2_x1l.dat',data_db2_x1l[:,1], fmt='%1.14f')
np.savetxt('db2_x1r.dat',data_db2_x1r[:,1], fmt='%1.14f')
!paste db2_x1l.dat db2_x1r.dat > DB2_x1_lr.dat
#Ejecutando el script de FES
!python generateFES.py DB2_x1_lr.dat $min_db2_x1l $max_db2_x1l $min_db2_x1r $max_db2_x1r 200 200 $temperatura DB2_XL1_XR1.dat
###################################################################
#Cargando valores del DB2_X2l
data_db2_x2l=np.loadtxt('dihed_db2_x2l.dat',comments=['#', '@'])
#Cargando valores del DB2_X1R
data_db2_x2r=np.loadtxt('dihed_db2_x2r.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB2_X1L
min_db2_x2l=np.amin(data_db2_x2l[:,1])
max_db2_x2l=np.amax(data_db2_x2l[:,1])
print ('Minimo DB2_X2L=>',min_db2_x2l)
print ('Máximo DB2_X2L=>',max_db2_x2l)
#Obteniendo los valores máximo y mínimo del DB2_X1R
min_db2_x2r=np.amin(data_db2_x2r[:,1])
max_db2_x2r=np.amax(data_db2_x2r[:,1])
print ('Minimo DB2_X2R=>',min_db2_x2r)
print ('Máximo DB2_X2R=>',max_db2_x2r)
#Creando los archivos de entrada para el script
np.savetxt('db2_x2l.dat',data_db2_x2l[:,1], fmt='%1.14f')
np.savetxt('db2_x2r.dat',data_db2_x2r[:,1], fmt='%1.14f')
!paste db2_x2l.dat db2_x2r.dat > DB2_x2_lr.dat
#Ejecutando el script de FES
!python generateFES.py DB2_x2_lr.dat $min_db2_x2l $max_db2_x2l $min_db2_x2r $max_db2_x2r 200 200 $temperatura DB2_XL2_XR2.dat
######################################################################################
#Cargando valores del DB2_X3M
data_db2_x3m=np.loadtxt('dihed_db2_x3m.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB2_X3M
min_db2_x3m=np.amin(data_db2_x3m[:,1])
max_db2_x3m=np.amax(data_db2_x3m[:,1])
print ('Minimo DB2_X3M=>',min_db2_x3m)
print ('Máximo DB2_X3M=>',max_db2_x3m)
print ('Minimo DB2_X1R=>',min_db2_x1r)
print ('Máximo DB2_X1R=>',max_db2_x1r)
print ('Minimo DB2_X2R=>',min_db2_x2r)
print ('Máximo DB2_X2R=>',max_db2_x2r)
print ('Minimo DB2_X1L=>',min_db2_x1l)
print ('Máximo DB2_X1L=>',max_db2_x1l)
print ('Minimo DB2_X2L=>',min_db2_x2l)
print ('Máximo DB2_X2L=>',max_db2_x2l)
#Creando los archivos de entrada para el script
np.savetxt('db2_x3m.dat',data_db2_x3m[:,1], fmt='%1.14f')
!paste db2_x3m.dat db2_x1r.dat > DB2_x3m_x1r.dat
!paste db2_x3m.dat db2_x2r.dat > DB2_x3m_x2r.dat
!paste db2_x3m.dat db2_x1l.dat > DB2_x3m_x1l.dat
!paste db2_x3m.dat db2_x2l.dat > DB2_x3m_x2l.dat
#Ejecutando el script de FES
!python generateFES.py DB2_x3m_x1r.dat $min_db2_x3m $max_db2_x3m $min_db2_x1r $max_db2_x1r 200 200 $temperatura DB2_XM3_XR1.dat
!python generateFES.py DB2_x3m_x2r.dat $min_db2_x3m $max_db2_x3m $min_db2_x2r $max_db2_x2r 200 200 $temperatura DB2_XM3_XR2.dat
!python generateFES.py DB2_x3m_x1l.dat $min_db2_x3m $max_db2_x3m $min_db2_x1l $max_db2_x1l 200 200 $temperatura DB2_XM3_XL1.dat
!python generateFES.py DB2_x3m_x2l.dat $min_db2_x3m $max_db2_x3m $min_db2_x2l $max_db2_x2l 200 200 $temperatura DB2_XM3_XL2.dat
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xl1_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_1}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XL1_XR1.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xl2_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^L_2}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XL2_XR2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xm3_vs_xl1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XL1.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xm3_vs_xl2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^L_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XL2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xm3_vs_xr2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_2}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XR2.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_xm3_vs_xr1.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 X@^M_3}"
set ylabel "{/=30 X@^R_1}"
set title "Free Energy Surface Intramolecular DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB2_XM3_XR1.dat" with pm3d
In [ ]:
############################################
#### Intermolecular DB1- DB2 - X1L
############################################
#Creando el DB1-DB2-X1L
!paste db1_x1l.dat db2_x1l.dat > DB1_DB2_x1l.dat
print('Minimo DB1-X1L=>',min_x1l)
print('Máximo DB1-X1L=>',max_x1l)
print('Minimo DB2-X1L=>',min_db2_x1l)
print('Máximo DB2-X1L=>',max_db2_x1l)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x1l.dat $min_x1l $max_x1l $min_db2_x1l $max_db2_x1l 200 200 $temperatura DB1_DB2_X1L.dat
#########################################
#### Intermolecular DB1- DB2 - X2L
############################################
#Creando el DB1-DB2-X2L
!paste db1_x2l.dat db2_x2l.dat > DB1_DB2_x2l.dat
print('Minimo DB1-X2L=>',min_x2l)
print('Máximo DB1-X2L=>',max_x2l)
print('Minimo DB2-X2L=>',min_db2_x2l)
print('Máximo DB2-X2L=>',max_db2_x2l)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x2l.dat $min_x2l $max_x2l $min_db2_x2l $max_db2_x2l 200 200 $temperatura DB1_DB2_X2L.dat
############################################
#### Intermolecular DB1- DB2 - X3M
############################################
#Creando el DB1-DB2-X3M
!paste db1_x3m.dat db2_x3m.dat > DB1_DB2_x3m.dat
print('Minimo DB1-X3M=>',min_x3m)
print('Máximo DB1-X3M=>',max_x3m)
print('Minimo DB2-X3M=>',min_db2_x3m)
print('Máximo DB2-X3M=>',max_db2_x3m)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x3m.dat $min_x3m $max_x3m $min_db2_x3m $max_db2_x3m 200 200 $temperatura DB1_DB2_X3M.dat
############################################
#### Intermolecular DB1- DB2 - X2R
############################################
#Creando el DB1-DB2-X2R
!paste db1_x2r.dat db2_x2r.dat > DB1_DB2_x2r.dat
print('Minimo DB1-X2R=>',min_x2r)
print('Máximo DB1-X2R=>',max_x2r)
print('Minimo DB2-X2R=>',min_db2_x2r)
print('Máximo DB2-X2R=>',max_db2_x2r)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x2r.dat $min_x2r $max_x2r $min_db2_x2r $max_db2_x2r 200 200 $temperatura DB1_DB2_X2R.dat
############################################
#### Intermolecular DB1- DB2 - X1R
############################################
#Creando el DB1-DB2-X1R
!paste db1_x1r.dat db2_x1r.dat > DB1_DB2_x1r.dat
print('Minimo DB1-X1R=>',min_x1r)
print('Máximo DB1-X1R=>',max_x1r)
print('Minimo DB2-X1R=>',min_db2_x1r)
print('Máximo DB2-X1R=>',max_db2_x1r)
#Ejecutando el script de FES
!python generateFES.py DB1_DB2_x1r.dat $min_x1r $max_x1r $min_db2_x1r $max_db2_x1r 200 200 $temperatura DB1_DB2_X1R.dat
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "DB1_DB2_X1L.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^L_1}"
set ylabel "{/=30 DB2 X@^L_1}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X1L.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "DB1_DB2_X2L.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^L_2}"
set ylabel "{/=30 DB2 X@^L_2}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X2L.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "DB1_DB2_X3M.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^M_3}"
set ylabel "{/=30 DB2 X@^M_3}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X3M.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "DB1_DB2_X2R.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set xyplane 0
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^R_2}"
set ylabel "{/=30 DB2 X@^R_2}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "DB1_DB2_X2R.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "DB1_DB2_X1R.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 DB1 X@^R_1}"
set ylabel "{/=30 DB2 X@^R_1}"
set title "Free Energy Surface Intermolecular DB1-DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
set cbrange[8:10]
splot "DB1_DB2_X1R.dat" with pm3d
In [ ]:
hist_escale_y=[]
fig = pl.figure(figsize=(25,8))
fig.subplots_adjust(hspace=.4, wspace=.3)
#subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
#left = 0.125 # the left side of the subplots of the figure
#right = 0.9 # the right side of the subplots of the figure
#bottom = 0.1 # the bottom of the subplots of the figure
#top = 0.9 # the top of the subplots of the figure
#wspace = 0.2 # the amount of width reserved for blank space between subplots
#hspace = 0.2 # the amount of height reserved for white space between subplots
#Formateando los valores de los ejes
#Engrosando marcos
ax = fig.add_subplot(2,5,1)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,2)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,3)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,4)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
ax = fig.add_subplot(2,5,5)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(3)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
#Cargando valores del DB1
data_h_db1_x1l=np.loadtxt('db1_x1l.dat',comments=['#', '@'])
data_h_db1_x2l=np.loadtxt('db1_x2l.dat',comments=['#', '@'])
data_h_db1_x3m=np.loadtxt('db1_x3m.dat',comments=['#', '@'])
data_h_db1_x2r=np.loadtxt('db1_x2r.dat',comments=['#', '@'])
data_h_db1_x1r=np.loadtxt('db1_x1r.dat',comments=['#', '@'])
#Cargando valores del DB2
data_h_db2_x1l=np.loadtxt('db2_x1l.dat',comments=['#', '@'])
data_h_db2_x2l=np.loadtxt('db2_x2l.dat',comments=['#', '@'])
data_h_db2_x3m=np.loadtxt('db2_x3m.dat',comments=['#', '@'])
data_h_db2_x2r=np.loadtxt('db2_x2r.dat',comments=['#', '@'])
data_h_db2_x1r=np.loadtxt('db2_x1r.dat',comments=['#', '@'])
sub1 = fig.add_subplot(251) # instead of plt.subplot(2, 2, 1)
sub1.set_xlabel('Angle (Degree) ', fontsize=10)
sub1.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub1.hist(data_h_db1_x1l,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub1.hist(data_h_db2_x1l,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub1.axis()
hist_escale_y.append(y2)
sub2 = fig.add_subplot(252) # instead of plt.subplot(2, 2, 1)
sub2.set_xlabel('Angle (Degree) ', fontsize=10)
sub2.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub2.hist(data_h_db1_x2l,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub2.hist(data_h_db2_x2l,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub2.axis()
hist_escale_y.append(y2)
sub3 = fig.add_subplot(253) # instead of plt.subplot(2, 2, 1)
sub3.set_xlabel('Angle (Degree) ', fontsize=10)
sub3.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub3.hist(data_h_db1_x3m,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub3.hist(data_h_db2_x3m,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub3.axis()
hist_escale_y.append(y2)
sub4 = fig.add_subplot(254) # instead of plt.subplot(2, 2, 1)
sub4.set_xlabel('Angle (Degree) ', fontsize=10)
sub4.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub4.hist(data_h_db1_x2r,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub4.hist(data_h_db2_x2r,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub4.axis()
hist_escale_y.append(y2)
sub5 = fig.add_subplot(255) # instead of plt.subplot(2, 2, 1)
sub5.set_xlabel('Angle (Degree) ', fontsize=10)
sub5.set_ylabel('P(Angle)')
n1, bins1, rectangles1 = sub5.hist(data_h_db1_x1r,100, normed=True, color='black',histtype='step', linewidth=3)
n2, bins2, rectangles2 = sub5.hist(data_h_db2_x1r,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=sub5.axis()
hist_escale_y.append(y2)
#escale_y
hist_escale_y.sort(reverse=True)
hist_escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,hist_escale_y[0]))
sub2.axis((x1,x2,y1,hist_escale_y[0]))
sub3.axis((x1,x2,y1,hist_escale_y[0]))
sub4.axis((x1,x2,y1,hist_escale_y[0]))
sub5.axis((x1,x2,y1,hist_escale_y[0]))
In [ ]:
### Creando el directorio para el análisis de las distancias de enlace de los puentes INTERMOLECULAR
ruta_bonds_puentes = nuevaruta+'/bonds_puentes'
print ( ruta_bonds_puentes )
if not os.path.exists(ruta_bonds_puentes):
os.makedirs(ruta_bonds_puentes)
print ('Se ha creado la ruta ===>',ruta_bonds_puentes)
else:
print ("La ruta "+ruta_bonds_puentes+" existe..!!!")
print ( 'Nos vamos a ....', ruta_bonds_puentes)
os.chdir( ruta_bonds_puentes )
In [ ]:
print ('\nCopiando el archivo generateFES.py a '+ruta_bonds_puentes)
source_file=ruta_scripts+'/free_energy/generateFES.py'
dest_file=ruta_bonds_puentes+'/generateFES.py'
shutil.copy(source_file,dest_file)
#Cambiando permisos de ejecución
!chmod +x generateFES.py
In [ ]:
psf=ruta_old_traj+'/'+psf_file
dcd=ruta_old_traj+'/'+dcd_file
print ('Puente DB1=>',DB1_N)
print ('Puente DB1=>',DB1_i)
print ('Puente DB2=>',DB2_N)
print ('Puente DB2=>',DB2_i)
puente=2
if (int(puente)==2):
#Creando script para Bond X1 Left
b1 = open('bond_DB1_left.tcl', 'w')
print(b1)
b1.write('set psfFile '+ psf+' \n')
b1.write('set dcdFile '+ dcd+' \n')
b1.write('\nmol load psf $psfFile dcd $dcdFile\n')
b1.write('set outfile ' +'[open ' +'bond_db1_left.dat'+' w]\n')
b1.write('set nf [molinfo top get numframes]\n')
b1.write(' \n')
b1.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[1]+'"] get index]\n')
b1.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[2]+'"] get index]\n')
b1.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[3]+'"] get index]\n')
b1.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
b1.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
b1.write(' set x [measure angle $angle frame $i]\n')
b1.write(' set time [expr {$i +1}]\n')
b1.write(' puts $outfile "$time $x"\n')
b1.write('}\n')
b1.close()
#Creando script para Bond X1 Right
b2 = open('bond_DB1_right.tcl', 'w')
print(b2)
b2.write('set psfFile '+ psf+' \n')
b2.write('set dcdFile '+ dcd+' \n')
b2.write('\nmol load psf $psfFile dcd $dcdFile\n')
b2.write('set outfile ' +'[open ' +'bond_db1_right.dat'+' w]\n')
b2.write('set nf [molinfo top get numframes]\n')
b2.write(' \n')
b2.write('set selatoms1 [[atomselect top "protein and chain A and '+DB1_i[4]+'"] get index]\n')
b2.write('set selatoms2 [[atomselect top "protein and chain A and '+DB1_i[5]+'"] get index]\n')
b2.write('set selatoms3 [[atomselect top "protein and chain A and '+DB1_i[6]+'"] get index]\n')
b2.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
b2.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
b2.write(' set x [measure angle $angle frame $i]\n')
b2.write(' set time [expr {$i +1}]\n')
b2.write(' puts $outfile "$time $x"\n')
b2.write('}\n')
b2.close()
#Creando script para Bond DB2 X1 Left
b3 = open('bond_DB2_left.tcl', 'w')
print(b3)
b3.write('set psfFile '+ psf+' \n')
b3.write('set dcdFile '+ dcd+' \n')
b3.write('\nmol load psf $psfFile dcd $dcdFile\n')
b3.write('set outfile ' +'[open ' +'bond_db2_left.dat'+' w]\n')
b3.write('set nf [molinfo top get numframes]\n')
b3.write(' \n')
b3.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[1]+'"] get index]\n')
b3.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[2]+'"] get index]\n')
b3.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[3]+'"] get index]\n')
b3.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
b3.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
b3.write(' set x [measure angle $angle frame $i]\n')
b3.write(' set time [expr {$i +1}]\n')
b3.write(' puts $outfile "$time $x"\n')
b3.write('}\n')
b3.close()
#Creando script para Bond DB2 X1 Right
b4 = open('bond_DB2_right.tcl', 'w')
print(b4)
b4.write('set psfFile '+ psf+' \n')
b4.write('set dcdFile '+ dcd+' \n')
b4.write('\nmol load psf $psfFile dcd $dcdFile\n')
b4.write('set outfile ' +'[open ' +'bond_db2_right.dat'+' w]\n')
b4.write('set nf [molinfo top get numframes]\n')
b4.write(' \n')
b4.write('set selatoms1 [[atomselect top "protein and chain A and '+DB2_i[4]+'"] get index]\n')
b4.write('set selatoms2 [[atomselect top "protein and chain A and '+DB2_i[5]+'"] get index]\n')
b4.write('set selatoms3 [[atomselect top "protein and chain A and '+DB2_i[6]+'"] get index]\n')
b4.write('set angle [list [lindex $selatoms1] [lindex $selatoms2] [lindex $selatoms3] ]\n')
b4.write('for {set i 0} {$i < $nf} {incr i 1} {\n')
b4.write(' set x [measure angle $angle frame $i]\n')
b4.write(' set time [expr {$i +1}]\n')
b4.write(' puts $outfile "$time $x"\n')
b4.write('}\n')
b4.close()
In [ ]:
#Calculando con VMD bond DB1 Left
!vmd -dispdev text < bond_DB1_left.tcl
#Calculando con VMD bond DB1 Right
!vmd -dispdev text < bond_DB1_right.tcl
#Calculando con VMD bond DB2 Left
!vmd -dispdev text < bond_DB2_left.tcl
#Calculando con VMD bond DB2 Right
!vmd -dispdev text < bond_DB2_right.tcl
In [ ]:
#Cargando valores del DB1
data_bond_db1_left=np.loadtxt('bond_db1_left.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_bond_db1_right=np.loadtxt('bond_db1_right.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB1 Left
min_bond1_left=np.amin(data_bond_db1_left[:,1])
max_bond1_left=np.amax(data_bond_db1_left[:,1])
print ('Minimo DB1_Left=>',min_bond1_left)
print ('Máximo DB1_Left=>',max_bond1_left)
#Obteniendo los valores máximo y mínimo del DB1 Right
min_bond1_right=np.amin(data_bond_db1_right[:,1])
max_bond1_right=np.amax(data_bond_db1_right[:,1])
print ('Minimo DB1_Right=>',min_bond1_right)
print ('Máximo DB1_Right=>',max_bond1_right)
#Creando los archivos de entrada para el script
np.savetxt('bond_DB1_left.dat',data_bond_db1_left[:,1], fmt='%1.14f')
np.savetxt('bond_DB1_right.dat',data_bond_db1_right[:,1], fmt='%1.14f')
!paste bond_DB1_left.dat bond_DB1_right.dat > angles_DB1.dat
#Ejecutando el script de FES
!python generateFES.py angles_DB1.dat $min_bond1_left $max_bond1_left $min_bond1_right $max_bond1_right 200 200 $temperatura Angles_DB1.dat
###################################################################3
#Cargando valores del DB2
data_bond_db2_left=np.loadtxt('bond_db2_left.dat',comments=['#', '@'])
#Cargando valores del DB1_X1R
data_bond_db2_right=np.loadtxt('bond_db2_right.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del DB2 Left
min_bond2_left=np.amin(data_bond_db2_left[:,1])
max_bond2_left=np.amax(data_bond_db2_left[:,1])
print ('Minimo DB2_Left=>',min_bond2_left)
print ('Máximo DB2_Left=>',max_bond2_left)
#Obteniendo los valores máximo y mínimo del DB2 Right
min_bond2_right=np.amin(data_bond_db2_right[:,1])
max_bond2_right=np.amax(data_bond_db2_right[:,1])
print ('Minimo DB2_Right=>',min_bond2_right)
print ('Máximo DB2_Right=>',max_bond2_right)
#Creando los archivos de entrada para el script
np.savetxt('bond_DB2_left.dat',data_bond_db2_left[:,1], fmt='%1.14f')
np.savetxt('bond_DB2_right.dat',data_bond_db2_right[:,1], fmt='%1.14f')
!paste bond_DB2_left.dat bond_DB2_right.dat > angles_DB2.dat
#Ejecutando el script de FES
!python generateFES.py angles_DB2.dat $min_bond2_left $max_bond2_left $min_bond2_right $max_bond2_right 200 200 $temperatura Angles_DB2.dat
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db1_a1_a2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 C@^1_{/Symbol a}}-{/=30 C@^1_{/Symbol b}}-{/=30 S@^1_{/Symbol g}}"
set ylabel "{/=30 C@^2_{/Symbol a}}-{/=30 C@^2_{/Symbol b}}-{/=30 S@^2_{/Symbol g}}"
set title "Free Energy Surface Angles DB1"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "Angles_DB1.dat" with pm3d
In [ ]:
# This loads the magics for gnuplot
%reload_ext gnuplot_kernel
#Configurando la salida para GNUplot
%gnuplot inline pngcairo transparent enhanced font "arial,20" fontscale 1.0 size 1280,960; set zeroaxis;;
In [ ]:
%%gnuplot
set output "db2_a1_a2.png"
set palette model RGB
set palette defined ( 0 '#000090',\
1 '#000fff',\
2 '#0090ff',\
3 '#0fffee',\
4 '#90ff70',\
5 '#ffee00',\
6 '#ff7000',\
7 '#ee0000',\
8 '#7f0000')
set view map
set dgrid3d
set pm3d interpolate 0,0
set xlabel "{/=30 C@^1_{/Symbol a}}-{/=30 C@^1_{/Symbol b}}-{/=30 S@^1_{/Symbol g}}"
set ylabel "{/=30 C@^2_{/Symbol a}}-{/=30 C@^2_{/Symbol b}}-{/=30 S@^2_{/Symbol g}}"
set title "Free Energy Surface Angles DB2"
##Descomentar la siguiente línea de código en caso de que la escala comience con valor de 1 y ejecutar nuevamente
#set cbrange[8:10]
splot "Angles_DB2.dat" with pm3d
In [ ]:
bonds_escale_y=[]
#Cargando valores del DB1
data_h_db1_left=np.loadtxt('bond_DB1_left.dat',comments=['#', '@'])
data_h_db1_right=np.loadtxt('bond_DB1_right.dat',comments=['#', '@'])
#Cargando valores del DB2
data_h_db2_left=np.loadtxt('bond_DB2_left.dat',comments=['#', '@'])
data_h_db2_right=np.loadtxt('bond_DB2_right.dat',comments=['#', '@'])
#Engrosar marco
figb=pl.figure(figsize=(12, 10), dpi=100, linewidth=3.0)
figb.subplots_adjust(hspace=.5)
ax = figb.add_subplot(221)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(222)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(223)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
ax = figb.add_subplot(224)
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(4)
#Formateando los valores de los ejes
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
bond1 = figb.add_subplot(221) # instead of plt.subplot(2, 2, 1)
#bond1.set_title('CA1 - CB1 - SY1')
bond1.set_xlabel('Angle (Degree)')
bond1.set_ylabel('P (Angle)')
n, bins, rectangles = bond1.hist(data_h_db1_left,100, normed=True, color='black',histtype='step',linewidth=3)
x1,x2,y1,y2=bond1.axis()
bonds_escale_y.append(y2)
bond2 = figb.add_subplot(222) # instead of plt.subplot(2, 2, 1)
#bond2.set_title('CA2 - CB2 - SY2')
bond2.set_xlabel('Angle (Degree)')
bond2.set_ylabel('P (Angle)')
n, bins, rectangles = bond2.hist(data_h_db1_right,100, normed=True, color='black',histtype='step', linewidth=3)
x1,x2,y1,y2=bond2.axis()
bonds_escale_y.append(y2)
bond3 = figb.add_subplot(223) # instead of plt.subplot(2, 2, 1)
#bond3.set_title('CA1 - CB1 - SY1')
bond3.set_xlabel('Angle (Degree)')
bond3.set_ylabel('P (Angle)')
n, bins, rectangles = bond3.hist(data_h_db2_left,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=bond3.axis()
bonds_escale_y.append(y2)
bond4 = figb.add_subplot(224) # instead of plt.subplot(2, 2, 1)
#bond4.set_title('CA2 - CB2 - SY2')
bond4.set_xlabel('Angle (Degree)')
bond4.set_ylabel('P (Angle)')
n, bins, rectangles = bond4.hist(data_h_db2_right,100, normed=True, color='red',histtype='step', linewidth=3)
x1,x2,y1,y2=bond4.axis()
bonds_escale_y.append(y2)
#escale_y
bonds_escale_y.sort(reverse=True)
bonds_escale_y
##Cambiando los ejes de las y
sub1.axis((x1,x2,y1,bonds_escale_y[0]))
sub2.axis((x1,x2,y1,bonds_escale_y[0]))
sub3.axis((x1,x2,y1,bonds_escale_y[0]))
sub4.axis((x1,x2,y1,bonds_escale_y[0]))
In [ ]:
### Creando el directorio para el análisis de los puentes
ruta_clusters = nuevaruta+'/clusters'
print ( ruta_clusters )
if not os.path.exists(ruta_clusters):
os.makedirs(ruta_clusters)
print ('Se ha creado la ruta ===>',ruta_clusters)
else:
print ("La ruta "+ruta_clusters+" existe..!!!")
print ( 'Nos vamos a ....', ruta_clusters)
os.chdir( ruta_clusters )
In [ ]:
!echo 1 1 | g_cluster -f ../output.xtc -s ../ionized.pdb -method gromos -cl out.pdb -g out.log -cutoff 0.2
In [ ]:
!vmd out.pdb
In [ ]:
### Creando el directorio para el análisis de colorByRMSF
ruta_colorByRMSF = nuevaruta+'/colorByRMSF'
print ( ruta_colorByRMSF )
if not os.path.exists(ruta_colorByRMSF):
os.makedirs(ruta_colorByRMSF)
print ('Se ha creado la ruta ===>',ruta_colorByRMSF)
else:
print ("La ruta "+ruta_colorByRMSF+" existe..!!!")
print ( 'Nos vamos a ....', ruta_colorByRMSF)
os.chdir( ruta_colorByRMSF )
In [ ]:
print ('\nCopiando el archivo colorByRMSF.vmd a '+ruta_colorByRMSF)
source_file=ruta_scripts+'/colorByRMSF/colorByRMSF.vmd'
dest_file=ruta_colorByRMSF+'/colorByRMSF.vmd'
shutil.copy(source_file,dest_file)
In [ ]:
print ('Ejecutando el análisis de rmsf...')
!echo 1 | g_rmsf -f ../output.xtc -s ../ionized.pdb -oq bfac.pdb -o rmsf.xvg
In [ ]:
#Calculando el mínimo y máximo del rmsf
#Cargando valores del RMSF
data_rmsf_gcolor=np.loadtxt('rmsf.xvg',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del RMSF
min_rmsf_gcolor=np.amin(data_rmsf_gcolor[:,1])
max_rmsf_gcolor=np.amax(data_rmsf_gcolor[:,1])
print ('Minimo_RMSF=>',min_rmsf_gcolor)
print ('Máximo_RMSF=>',max_rmsf_gcolor)
Arrancar VMD, dirigirse al menú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos en el cual pondremos los valores del Mínimo_RMSF y Máximo_RMSF calculado en la celda anterior:
source colorByRMSF.vmd
colorByRMSF top rmsf.xvg Mínimo_RMSF Máximo_RMSF
ESCALA DE COLOR
Dirigirse al menú Extensions -> Visualization -> Color Scale Bar y cambiar los valores de los siguientes campos:
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
REMOVER EJE X,Y,Z
Dirigirse al menú Display -> Axes -> Off, con el cual eliminaremos el eje de X,Y,Z.
In [ ]:
# Cargando el pdb con VMD
!vmd ../ionized.pdb
In [ ]:
print ( 'Nos vamos a ....', ruta_colorByRMSF )
os.chdir( ruta_colorByRMSF )
In [ ]:
#Inicializando vector
rmsf=[]
rmsf_x=[]
rmsf_y=[]
try:
file_Bfactor = open( 'bfac.pdb' )
new_bfactor=open('bfac_new.pdb','w')
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_Bfactor.readlines():
fila = linea.strip()
sl = fila.split()
cadena=sl[0]
if (cadena=='ATOM'):
if (len(sl)==12):
new_bfactor.write(linea)
else:
x=linea[0:60]
tempFactor=linea[60:66]
#print (x)
#print(tempFactor)
y=fila[67:]
#print (y)
enviar=x+' '+tempFactor+y
new_bfactor.write(enviar+'\n')
#print(enviar)
else:
#print (linea)
new_bfactor.write(linea)
new_bfactor.close()
In [ ]:
!gedit bfac_new.pdb
In [ ]:
#Inicializando vector
bfactors_color=[]
try:
file_bfactor_color = open( 'bfac_new.pdb' )
except IOError:
print ('No se pudo abrir el archivo o no existe·..')
i=0
for linea in file_bfactor_color.readlines():
fila = linea.strip()
sl = fila.split()
if (sl[0]=='ATOM'):
#print (sl[0])
idresidue=fila[23:26]
bfactor=fila[60:66]
#print (idresidue + '\t'+bfactor)
bfactors_color.append(idresidue+'\t'+bfactor+'\n')
#i=i+1
#Escribiendo el archivo BFACTOR.dat
f = open('protein_bfactor.dat', 'w')
#f.write('@ title "B-factors" \n')
f.write('@ xaxis label " Residue" \n')
f.write('@ xaxis label char size 1.480000\n')
f.write('@ xaxis bar linewidth 5.0\n')
f.write('@ xaxis ticklabel char size 1.480000\n')
f.write('@ yaxis label "B-factors (' +"\\"+'cE'+"\\"+'C)"\n')
f.write('@ yaxis label char size 1.480000\n')
f.write('@ yaxis bar linewidth 5.0\n')
f.write('@ yaxis ticklabel char size 1.480000\n')
f.write('@ s0 line linewidth 7\n')
f.write('@ s0 symbol 1\n')
f.write('@ s0 symbol size 1.000000\n')
f.write('@ s0 symbol color 1\n')
f.write('@ s0 symbol pattern 1\n')
f.write('@ s0 symbol fill color 2\n')
f.write('@ s0 symbol fill pattern 1\n')
f.write('@ s0 symbol linewidth 1.0\n')
f.write('@TYPE xy \n')
f.write("".join(bfactors_color))
f.close()
In [ ]:
!xmgrace protein_bfactor.dat
In [ ]:
#Cargando la imagen generada en xmgrace
Image(filename='protein_bfactor.png')
In [ ]:
#Calculando el mínimo y máximo del rmsf
#Cargando valores del RMSF
data_bfactor_color=np.loadtxt('protein_bfactor.dat',comments=['#', '@'])
#Obteniendo los valores máximo y mínimo del RMSF
min_bfactor_color=np.amin(data_bfactor_color[:,1])
max_bfactor_color=np.amax(data_bfactor_color[:,1])
print ('Minimo_B-Factor=>',min_bfactor_color)
print ('Máximo_B-Factor=>',max_bfactor_color)
In [ ]:
!chimera bfac_new.pdb
ESTABLECER EL MODO DE VISUALIZACIÓN
COLOREAR LOS B-FACTORS
Seleccionar Tools -> Depiction -> Render by Attribute.
Nos desplegará una ventana Render/Select by Attribute.
FONDO BLANCO
Para aplicar el fondo blanco:
SALVAR LA POSICIÓN DE LA IMAGEN
Una vez que se ha obtenido la imagen coloreada, ajustar la visualización rotando la imagen, con la finalidad de dejar los espacios adecuados para la inclusión de las etiquetas y la barra de color.
Para salvar la posición final de la imagen:
Si por alguna razón movemos la posición, para restaurarla hacer lo siguiente:
TITULO Y BARRA DE COLOR
Seleccionar del menú principal Tools -> Utilities -> Color Key. El cual desplegará la ventana 2D Labels/Color Key.
Para desplegar la barra de color:
Para desplegar el título de la barra:
Para desplegar el título de la imagen:
Notas:
SALVAR LA IMAGEN
Seleccionar del menú principal File -> Save Image. El cual desplegará la ventana Save image, en el cual en el campo File name dar el nombre de image.png.
SALVAR LA SESIÓN DE QUIMERA
Seleccionar del menú principal File -> Save Session as. El cual desplegará la ventana Choose Session Save File, en el cual en el campo File name colocar el nombre con la extensión .py.
In [ ]:
##Cargando la imagen generada
print ('Cargando el archivo...')
Image(filename='image.png')
In [ ]:
### Creando el directorio para el análisis del SASA en el directorio de VMD
print ('Nos vamos a ', ruta)
os.chdir( ruta )
output_find=!find /usr/local -maxdepth 2 -type d -name vmd
print (output_find)
ruta_vmd=output_find[0]
print (ruta_vmd)
ruta_vmd_sasa = ruta_vmd+'/plugins/noarch/tcl/iceVMD1.0'
print ( ruta_vmd_sasa )
if not os.path.exists(ruta_vmd_sasa):
os.makedirs(ruta_vmd_sasa)
print ('Se ha creado la ruta ===>',ruta_vmd_sasa)
else:
print ("La ruta "+ruta_vmd_sasa+" existe..!!!")
print ( 'Nos vamos a ....', ruta_vmd_sasa )
os.chdir( ruta_vmd_sasa )
In [ ]:
#Copiando los archivos generados a la carpeta plugins de VMD
print ('\nCopiando los archivos generados a '+ruta_vmd_sasa)
source_file=ruta_scripts+'/iceVMD1.0/colorplot.tcl'
dest_file=ruta_vmd_sasa+'/colorplot.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/multiplot.tcl'
dest_file=ruta_vmd_sasa+'/multiplot.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/pkgIndex.tcl'
dest_file=ruta_vmd_sasa+'/pkgIndex.tcl'
shutil.copy(source_file,dest_file)
source_file=ruta_scripts+'/iceVMD1.0/vmdICE.tcl'
dest_file=ruta_vmd_sasa+'/vmdICE.tcl'
shutil.copy(source_file,dest_file)
print('\nArchivos copiados.. Regresando a... '+nuevaruta)
os.chdir( nuevaruta )
In [ ]:
### Creando el directorio para la graficación del sasa
ruta_sasaColor = nuevaruta+'/sasaColor'
print ( ruta_sasaColor )
if not os.path.exists(ruta_sasaColor):
os.makedirs(ruta_sasaColor)
print ('Se ha creado la ruta ===>',ruta_sasaColor)
else:
print ("La ruta "+ruta_sasaColor+" existe..!!!")
print ( 'Nos vamos a ....', ruta_sasaColor )
os.chdir( ruta_sasaColor )
In [ ]:
print ('\nCopiando el archivo de configuracion a '+ruta_sasaColor)
source_file=ruta_scripts+'/iceVMD1.0/vmdrc'
dest_file=ruta_sasaColor+'/.vmdrc'
shutil.copy(source_file,dest_file)
Arrancar VMD.
Ventana vmdICE
Dirigirse al menú Extensions -> Analysis -> vmdICE, se presentará una ventana y se deberán cambiar los valores de los siguientes campos:
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
CAMBIAR RESOLUCIÓN DE ESFERAS
Dirigirse al menú Graphics - Representations, y en el campo Sphere Resolution cambiamos al valor de 50.
ROTAR LA IMAGEN PARA PRESENTAR UNA MEJOR VISTA Y GUARDARLA.
In [ ]:
!vmd ../ionized.psf ../output.xtc
In [ ]:
#Borrando los archivos del vmd
!rm -r $ruta_vmd_sasa
In [ ]:
### Creando el directorio para la graficación del rgyro
ruta_gyroColor = nuevaruta+'/color_rgyro'
print ( ruta_gyroColor )
if not os.path.exists(ruta_gyroColor):
os.makedirs(ruta_gyroColor)
print ('Se ha creado la ruta ===>',ruta_gyroColor)
else:
print ("La ruta "+ruta_gyroColor+" existe..!!!")
print ( 'Nos vamos a ....', ruta_gyroColor )
os.chdir( ruta_gyroColor )
In [ ]:
print ('\nCopiando el script colorRgyro.tcl a '+ruta_gyroColor)
source_file=ruta_scripts+'/colorRgyro/colorRgyro.tcl'
dest_file=ruta_gyroColor+'/colorRgyro.tcl'
shutil.copy(source_file,dest_file)
Arrancar VMD, dirigirse al manú Extensions -> Tk Console, copiar y ejecutar la siguiente secuencia de comandos:
source colorRgyro.tcl
CAMBIAR EL COLOR DE FONDO
Dirigirse al menú Graphics -> Colors , y realizar las siguientes selecciones:
ROTAR LA IMAGEN PARA PRESENTAR UNA MEJOR VISTA Y GUARDARLA.
In [ ]:
!vmd ../ionized.psf ../output.xtc
In [ ]: