saxs restraints
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 22 Dec 2017 00:21:02 +0000 (01:21 +0100)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 22 Dec 2017 00:21:02 +0000 (01:21 +0100)
django_simple/todo/forms.py
django_simple/todo/jobfiles.py
django_simple/todo/models.py
django_simple/todo/templates/details.html
django_simple/todo/templates/details1.html
django_simple/todo/templates/edit.html
django_simple/todo/views.py
files/pbs8.csh
files/pbs8_new.csh
files/plot_saxs.py [new file with mode: 0755]
files/saxs_dist.awk [new file with mode: 0644]

index 0a8f597..5eb67fb 100644 (file)
@@ -411,6 +411,16 @@ class TaskForm_remd_a(forms.Form):
                        help_text='box z dimension')
 
 
+     wsaxs = forms.FloatField(label='SAXS weight',initial=100.0,
+                            help_text='weight for SAXS restraint term')
+     scal_rad = forms.FloatField(label='Scal_rad (SAXS)',initial=1.0,
+                            help_text='downscaling factor of residue radii used in SAXS restraints')
+     saxs_data = forms.CharField(label='P(r) SAXS data',
+                     help_text='distance distribution from SAXS, two columns: r and P(r)',
+                     required=False,
+                     widget=forms.Textarea(attrs={'cols': 25, 'rows': 20}))
+
+
      def clean(self):
              cleaned_data = super(TaskForm_remd_a, self).clean()
 
index d09ab98..fd7fd7c 100644 (file)
@@ -48,15 +48,13 @@ def write_on_task_save(sender, instance, **kwargs):
 WSCLOC=0.10554 WTOR=1.34316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000        &
 WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000    &
 WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.25                                        &
-CUTOFF=7.00000 WCORR4=0.00000
-"""
+CUTOFF=7.00000 WCORR4=0.00000"""
      else:
       w="""WSC=0.82686 WSCP=0.96947 WELEC=0.79373 WBOND=1.00000 WANG=0.46542              &
 WSCLOC=0.07969 WTOR=0.81684 WTORD=0.67806 WCORRH=0.00000 WCORR5=0.00000        &
 WCORR6=0.00000 WEL_LOC=0.71100 WTURN3=2.30298 WTURN4=0.86517 WTURN6=0.00000    &
 WSCCOR=0.14577 WVDWPP=0.17781 WHPB=1.00000 WSCP14=0.00000                      &
-CUTOFF=7.00000 WCORR4=0.00000
-"""
+CUTOFF=7.00000 WCORR4=0.00000"""
      
 
      with open(instance.jobdirname+'/file.inp','w') as f:
@@ -134,7 +132,16 @@ CUTOFF=7.00000 WCORR4=0.00000
            
           if instance.md_pdbref:
              control_line = control_line+'pdbref '
-           
+          
+          cntrl_saxs=''
+          if instance.saxs_data != '':
+             lines=instance.saxs_data.split('\n')
+             nsaxs=0
+             for line in lines:
+              if len(line.split())==2:
+               nsaxs+=1 
+             cntrl_saxs=' nsaxs='+str(nsaxs)+' scal_rad='+str(instance.scal_rad) 
+             control_line += cntrl_saxs             
             
           type_line = 'reset_vel='+ str(instance.remd_nstex)\
              +' nstep='+str(instance.md_nstep)\
@@ -195,8 +202,10 @@ CUTOFF=7.00000 WCORR4=0.00000
           for element in word_list[:-1]:
               f.write('{:79}'.format(element)+'&\n')
           f.write(word_list[-1]+'\n')
-
-       f.write(w)
+       
+       if instance.saxs_data != '':
+         w+=' wsaxs='+str(instance.wsaxs)
+       f.write(w+'\n')
 
 
        if instance.type == 'min' or instance.md_start == 'pdbstart':
@@ -223,7 +232,12 @@ CUTOFF=7.00000 WCORR4=0.00000
              f.write(seq[i:i+80]+'\n')
           f.write('0\n0\n')
 
-
+       if instance.saxs_data != '':
+             lines=instance.saxs_data.split('\n')
+             for line in lines:
+              tmp=line.split()
+              if len(tmp)==2:
+               f.write(tmp[0]+' '+tmp[1]+'\n') 
        
         
      if instance.type == 'min':
@@ -309,20 +323,29 @@ CUTOFF=7.00000 WCORR4=0.00000
             f.write('{:79}'.format('SEED='+str(instance.md_seed)+' isampl='+str(isampl)+
                ' einicheck=1 rescale=2 delta=0.02 cxfile classify')+'&\n')
             f.write('BOXX='+str(instance.boxx)+' BOXY='+str(instance.boxy)+
-                    ' BOXZ='+str(instance.boxz) +'\n')
+                    ' BOXZ='+str(instance.boxz)+cntrl_saxs +'\n')
           else:    
             f.write('{:79}'.format('SEED='+str(instance.md_seed)+' isampl='+str(isampl)+
               ' einicheck=1 rescale=2 delta=0.02 cxfile')+'&\n')
             f.write('BOXX='+str(instance.boxx)+' BOXY='+str(instance.boxy)+
-                    ' BOXZ='+str(instance.boxz) +'\n')
+                    ' BOXZ='+str(instance.boxz)+cntrl_saxs +'\n')
           f.write('nres='+str(len(seq))+' one_letter\n')
           for i in range(0,len(seq),80):
              f.write(seq[i:i+80]+'\n')
 
           f.write(write_ssbond(instance.ssbond))
 
-
-          f.write(w)
+          if instance.saxs_data != '':
+             fsaxs = open(instance.jobdirname+'/saxs.data', 'w')
+             lines=instance.saxs_data.split('\n')
+             for line in lines:
+              tmp=line.split()
+              if len(tmp)==2:
+               f.write(tmp[0]+' '+tmp[1]+'\n')
+               fsaxs.write(tmp[0]+' '+tmp[1]+'\n') 
+             fsaxs.close()
+
+          f.write(w+'\n')
           f.write('\n')
           f.write('nt='+str(instance.remd_nrep)+' replica read_iset\n')
           tmp1=json.loads(instance.remd_multi_t)
@@ -344,11 +367,17 @@ CUTOFF=7.00000 WCORR4=0.00000
              +' one_letter rescale=2 PRINT_CART PDBOUT=1 iopt=1'
              +' temper='+str(instance.remd_cluter_temp))+'&\n')
           f.write('BOXX='+str(instance.boxx)+' BOXY='+str(instance.boxy)+
-                  ' BOXZ='+str(instance.boxz) +'\n')
-          f.write(w)
+                  ' BOXZ='+str(instance.boxz)+cntrl_saxs +'\n')
+          f.write(w+'\n')
           for i in range(0,len(seq),80):
              f.write(seq[i:i+80]+'\n')
           f.write(write_ssbond(instance.ssbond))
+          if instance.saxs_data != '':
+             lines=instance.saxs_data.split('\n')
+             for line in lines:
+              tmp=line.split()
+              if len(tmp)==2:
+               f.write(tmp[0]+' '+tmp[1]+'\n') 
       
       os.chdir(instance.jobdirname)  
       ret_code = subprocess.Popen(' /opt/torque/bin/qsub pbs8.csh', shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
index e74eca0..da15522 100644 (file)
@@ -79,6 +79,12 @@ class Task(models.Model):
     remd_model3 = models.TextField(default='')
     remd_model4 = models.TextField(default='')
     remd_model5 = models.TextField(default='')
+
+#saxs
+    scal_rad = models.FloatField(default=1.0)
+    wsaxs = models.FloatField(default=100.0)
+    saxs_data = models.TextField(max_length=10000,default='')
+        
     
 # system    
     ready = models.BooleanField(default=False)
index 5c3e43c..c6ad8e9 100644 (file)
@@ -280,7 +280,20 @@ Created {{ task.created_date  }}
                            <div class="col-xs-10"> temperature for clustering</div>
                            <div class="col-xs-10">{{ task.remd_cluter_temp}}</div>
                        </li>           
-
+                       {% if task.saxs_data != "" %}                   
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> wsaxs </div>
+                           <div class="col-xs-10">{{ task.wsaxs }}</div>
+                       </li>           
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> scal_rad (saxs) </div>
+                           <div class="col-xs-10">{{ task.scal_rad }}</div>
+                       </li>           
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> saxs distribution </div>
+                           <div class="col-xs-10"><pre> {{ task.saxs_data }}</pre></div>
+                       </li>           
+                        {% endif %}
 
         {% endif %}
 
@@ -479,6 +492,16 @@ Created {{ task.created_date  }}
                          width="500"> </div>
                        </li>
                        
+                          {% if task.saxs_data != "" %}
+                       <li class="list-group-item task-item">
+                         <div class="col-xs-10"> P(r) of input SAXS data and
+                         calculated for 5 models </div>
+                         <div class="col-xs-10"> <img 
+                         src="/myfiles/download-file/{{task.jobdirname}}/saxs.png"
+                         width="500"> </div>
+                       </li>
+                          {% endif %}
+                       
        <li class="list-group-item task-item">
        <div class="col-xs-10">
         <button onclick="plusDivs(-1)">&#10094;</button> 
index 519ab0d..eaf6ff3 100644 (file)
@@ -280,7 +280,21 @@ Created {{ task.created_date  }}
                            <div class="col-xs-10"> temperature for clustering</div>
                            <div class="col-xs-10">{{ task.remd_cluter_temp}}</div>
                        </li>           
-
+                       
+                       {% if task.saxs_data != "" %}                   
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> wsaxs </div>
+                           <div class="col-xs-10">{{ task.wsaxs }}</div>
+                       </li>           
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> scal_rad (saxs) </div>
+                           <div class="col-xs-10">{{ task.scal_rad }}</div>
+                       </li>           
+                       <li class="list-group-item task-item">
+                           <div class="col-xs-10"> saxs distribution </div>
+                           <div class="col-xs-10"><pre> {{ task.saxs_data }}</pre></div>
+                       </li>           
+                        {% endif %}
 
         {% endif %}
 
@@ -478,6 +492,16 @@ Created {{ task.created_date  }}
                          src="/myfiles/download-file/{{task.jobdirname}}/remd_ex.png"
                          width="500"> </div>
                        </li>
+                          {% if task.saxs_data != "" %}
+                       <li class="list-group-item task-item">
+                         <div class="col-xs-10"> P(r) of input SAXS data and
+                         calculated for 5 models </div>
+                         <div class="col-xs-10"> <img 
+                         src="/myfiles/download-file/{{task.jobdirname}}/saxs.png"
+                         width="500"> </div>
+                       </li>
+                          {% endif %}
+                       
                        
        <li class="list-group-item task-item">
        <div class="col-xs-10">
index 90cec38..940ab03 100644 (file)
@@ -77,6 +77,10 @@ Calculation type {{ p_type }}
                </table>
                <button type="submit" class="save btn btn-default">Save & submit</button>
                <button type="submit" value="example" name="_example" class="save btn btn-default">Load example data</button>           
+               {% if p_type == "replica exchange molecular dynamics - advanced options" %}
+               <button type="submit" value="example_saxs"
+               name="_example_saxs" class="save btn btn-default">Load example SAXS data</button>               
+               {% endif %}
        </form>
 
        
index 61d5b61..c513fbf 100644 (file)
@@ -413,6 +413,7 @@ def add_remd(request,task_id):
 
 @login_required
 def add_remd_a(request,task_id):
+    from django.core.files.uploadedfile import UploadedFile
     task = get_object_or_404(Task, id=task_id)
     if request.method == 'POST':
      if '_example' in request.POST:
@@ -420,6 +421,49 @@ def add_remd_a(request,task_id):
          'md_nstep':500000,'md_lang':'berendsen','unres_ff':'opt-wtfsa-2',
          'remd_cluter_temp':280}
         form = TaskForm_remd_a(initial=data)     
+     elif '_example_saxs' in request.POST:
+        data= {'name':task.name,'pdbid':'5UJQ','md_pdbref':True,
+         'md_nstep':200000,'md_lang':'langevin','unres_ff':'E0LL2Y',
+         'remd_cluter_temp':270, 
+         'scal_rad':4.0,'wsaxs':100.0,'saxs_data':
+""" 0.5     1.33868e-02
+ 1.5     1.95880e-02
+ 2.5     2.68896e-02
+ 3.5     3.43737e-02
+ 4.5     4.07099e-02
+ 5.5     4.47875e-02
+ 6.5     4.63486e-02
+ 7.5     4.60514e-02
+ 8.5     4.49130e-02
+ 9.5     4.36744e-02
+10.5     4.26085e-02
+11.5     4.17464e-02
+12.5     4.11217e-02
+13.5     4.07835e-02
+14.5     4.06776e-02
+15.5     4.06060e-02
+16.5     4.03241e-02
+17.5     3.96655e-02
+18.5     3.85756e-02
+19.5     3.70537e-02
+20.5     3.50982e-02
+21.5     3.27236e-02
+22.5     3.00046e-02
+23.5     2.70643e-02
+24.5     2.40044e-02
+25.5     2.08595e-02
+26.5     1.76342e-02
+27.5     1.43802e-02
+28.5     1.12281e-02
+29.5     8.34574e-03
+30.5     5.87354e-03
+31.5     3.88732e-03
+32.5     2.39755e-03
+33.5     1.36323e-03
+34.5     7.06686e-04
+35.5     3.30592e-04
+36.5     1.38359e-04"""}
+        form = TaskForm_remd_a(initial=data)     
      else:
         form = TaskForm_remd_a(request.POST,request.FILES)
         if form.is_valid():
@@ -472,6 +516,10 @@ def add_remd_a(request,task_id):
              task.remd_cluter_temp=form.cleaned_data["remd_cluter_temp"]
              task.unres_ff=form.cleaned_data["unres_ff"]
              
+             task.scal_rad = form.cleaned_data["scal_rad"]
+             task.saxs_data = form.cleaned_data["saxs_data"]
+             task.wsaxs = form.cleaned_data["wsaxs"]
+             
              task.save()
              return redirect('addmlist',task_id=task.id)
     else:
index 0e3867d..82e7143 100755 (executable)
@@ -89,5 +89,14 @@ rm tmp.pdb
 /users2/local/bin/tmscore MODEL5.pdb plik1.pdb > tmscore5.out
 rm plik1.pdb
 
+if (-e saxs.data) then
+awk -f ../files/saxs_dist.awk MODEL1.pdb > MODEL1_saxs.data
+awk -f ../files/saxs_dist.awk MODEL2.pdb > MODEL2_saxs.data
+awk -f ../files/saxs_dist.awk MODEL3.pdb > MODEL3_saxs.data
+awk -f ../files/saxs_dist.awk MODEL4.pdb > MODEL4_saxs.data
+awk -f ../files/saxs_dist.awk MODEL5.pdb > MODEL5_saxs.data
+../files/plot_saxs.py
+endif
+
 #END
 touch finished
index ca528bb..67dcc8c 100755 (executable)
@@ -89,5 +89,15 @@ rm tmp.pdb
 /users2/local/bin/tmscore MODEL5.pdb plik1.pdb > tmscore5.out
 rm plik1.pdb
 
+if (-e saxs.data) then
+awk -f ../files/saxs_dist.awk MODEL1.pdb > MODEL1_saxs.data
+awk -f ../files/saxs_dist.awk MODEL2.pdb > MODEL2_saxs.data
+awk -f ../files/saxs_dist.awk MODEL3.pdb > MODEL3_saxs.data
+awk -f ../files/saxs_dist.awk MODEL4.pdb > MODEL4_saxs.data
+awk -f ../files/saxs_dist.awk MODEL5.pdb > MODEL5_saxs.data
+../files/plot_saxs.py
+endif
+
+
 #END
 touch finished
diff --git a/files/plot_saxs.py b/files/plot_saxs.py
new file mode 100755 (executable)
index 0000000..292550b
--- /dev/null
@@ -0,0 +1,30 @@
+#! /usr/bin/env python
+
+import matplotlib
+#matplotlib.use('GTK')
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import numpy as np
+import sys
+
+x,y= np.loadtxt('saxs.data',usecols=(0,1),unpack=True)
+x1,y1= np.loadtxt('MODEL1_saxs.data',usecols=(0,2),unpack=True)
+x2,y2= np.loadtxt('MODEL2_saxs.data',usecols=(0,2),unpack=True)
+x3,y3= np.loadtxt('MODEL3_saxs.data',usecols=(0,2),unpack=True)
+x4,y4= np.loadtxt('MODEL4_saxs.data',usecols=(0,2),unpack=True)
+x5,y5= np.loadtxt('MODEL5_saxs.data',usecols=(0,2),unpack=True)
+
+
+plt.clf()
+plt.xlabel('r')
+plt.ylabel('P(r)')
+plt.xlim(x[0], x[-1])
+plt.plot(x,y,'-',color='black') 
+plt.plot(x1,y1,'-',color='red')
+plt.plot(x2,y2,'-',color='lime')
+plt.plot(x3,y3,'-',color='blue')
+plt.plot(x4,y4,'-',color='magenta')
+plt.plot(x5,y5,'-',color='cyan')
+plt.legend(['reference','MODEL1','MODEL2','MODEL3','MODEL4','MODEL5'])
+plt.savefig('saxs.png')
diff --git a/files/saxs_dist.awk b/files/saxs_dist.awk
new file mode 100644 (file)
index 0000000..347dfde
--- /dev/null
@@ -0,0 +1,33 @@
+{
+  if ($1 == "ATOM" && $3 == "CA") {
+    nres = nres + 1
+    x[nres] = substr($0,31,8)
+    y[nres] = substr($0,39,8)
+    z[nres] = substr($0,47,8)
+            
+#    print "nres",nres," x",x[nres],y[nres],z[nres]
+  }
+}END{
+#  print "nres",nres
+  sigma = 2.50
+  for (i=1;i<=nres;i++) {
+    for (j=1;j<i;j++) {
+       dij = sqrt ((x[i]-x[j])^2+(y[i]-y[j])^2+(z[i]-z[j])^2)
+#       print i,j,dij
+       ibin = int(dij)
+       h[ibin] = h[ibin] + 1
+       for (k=0;k<=100;k++) {
+          dd = k+0.5
+#  print k,dd,dij,dd-dij,(dd-dij)^2/(2*sigma^2)
+          if ((dd-dij)^2/(2*sigma^2)<100) hh[k] = hh[k] + exp(-(dd-dij)^2/(2*sigma^2)) 
+       } 
+     }
+   }
+  norm = 0;
+  nnorm = 0;
+  for (k=0;k<=100;k++) {
+    norm = norm + h[k];
+    nnorm = nnorm + hh[k];
+  }
+  for (k=0;k<=100;k++) printf("%10.5f%15.5e%15.5e\n",k+0.5,h[k]/norm,hh[k]/nnorm)
+}