files/fluct_plot.py

   1 #! /usr/bin/env python
   2 import sys
   3 import math
   4
   5 import matplotlib
   6 #matplotlib.use('GTK')
   7 matplotlib.use('Agg')
   8 import matplotlib.pyplot as plt
   9 import os.path
  10 from itertools import cycle
  11
  12 iatom=0
  13 nmol=0
  14 sx={}
  15 sy={}
  16 sz={}
  17 s2x={}
  18 s2y={}
  19 s2z={}
  20 with open(sys.argv[1]) as f:
  21    for line in f:
  22       if line[0:4]=='ATOM' and line[13:15]=='CA':
  23 #         print line
  24          x=float(line[30:38])
  25          y=float(line[38:46])
  26          z=float(line[46:54])
  27 #         print x,y,z
  28          iatom=iatom+1
  29          sx[iatom]=sx.get(iatom,0)+x
  30          sy[iatom]=sy.get(iatom,0)+y
  31          sz[iatom]=sz.get(iatom,0)+z
  32          s2x[iatom]=s2x.get(iatom,0)+x*x
  33          s2y[iatom]=s2y.get(iatom,0)+y*y
  34          s2z[iatom]=s2z.get(iatom,0)+z*z
  35
  36
  37       if line[0:6]=='ENDMDL':
  38          natom=iatom
  39          iatom=0
  40          nmol=nmol+1
  41
  42 x=[]
  43 y=[]
  44 for i in range(1,natom+1):
  45    fluct=math.sqrt((s2x.get(i)-sx.get(i)*sx.get(i)/nmol\
  46                    +s2y.get(i)-sy.get(i)*sy.get(i)/nmol\
  47                    +s2z.get(i)-sz.get(i)*sz.get(i)/nmol  )/nmol)
  48    print fluct
  49    x.append(i)
  50    y.append(fluct)
  51
  52 b=[]
  53 newchain=True
  54 if os.path.exists('plik.pdb'):
  55  with open('plik.pdb') as f:
  56    for line in f:
  57       if line[0:4]=='ATOM' and line[13:15]=='CA':
  58 #         print line
  59         if newchain or int(line[22:26])>ires:
  60          b.append(math.sqrt(float(line[60:66])*3/8/math.pi/math.pi))
  61          ires=int(line[22:26])
  62          newchain=False
  63       if line[0:3]=='TER':
  64         newchain=True
  65       if line[0:3]=='END':
  66         break
  67
  68
  69
  70 plt.xlabel('residue')
  71 plt.ylabel('fluctuations')
  72 plt.xlim(0,natom+1)
  73 plt.plot(x,y,'-')
  74 if len(b)!=0:
  75  plt.plot(x,b,'-',c='red')
  76  plt.legend(['fluctuations','sqrt(3*bfactor/8*pi^2)'])
  77 plt.savefig('fluct_plot.png')
  78
  79 ycycle=cycle(y)
  80 with open('plik.pdb') as f, open ('plik_bf.pdb','w') as fw:
  81   prev_ires=None
  82   for line in f:
  83      if line[0:4]=='ATOM':
  84           if line[22:26] != prev_ires:
  85              prev_ires=line[22:26]
  86              bf=next(ycycle)
  87           fw.write(line[0:60]+'{:6.2f}'.format(bf)+line[66:])
  88      elif line[0:3] == 'TER' or line[0:6] == 'CONECT':
  89           fw.write(line)
  90