CODE set show confining sphere optionally in MCsim

analysis/MCsim.py

@@ -88,19 +88,19 @@ def __init__(self, path_to_data = "%s/data/" %(default_dir), trials = [''], traj
         
         Parameters
         ----------
-        path_to_data : string
+        path_to_data : string, optional
             path to either the wlcsim/data directory or the top level directory above the nested `trials`
         trials : list of strings
             list of the subdirectories in `path_to_data` that specify the different trials of the simulations 
-        trajectories : list 
+        trajectories : list, optional
             list of the "channel"/thread values of the simulation replicas, i.e. [1,2,3,4]
         time_min : int, optional
             default : 0
             minimum "time point" that you want to start reading in the data from
         time_max : int, optional
             default : 110
             maximum "time point" that you want to end reading in the data from. 
-        channel : int
+        channel : int, optional
             "channel"/thread value of a specific replica you want to read in  
         """
         # path to data directory
@@ -220,11 +220,9 @@ def __init__(self, path_to_data, time_min, time_max, trajectories, channel):
             path to either the wlcsim/data directory or the top level directory above the nested `trials`
         trajectories : list 
             list of the "channel"/thread values of the simulation replicas, i.e. [1,2,3,4]
-        time_min : int, optional
-            default : 0
+        time_min : int
             minimum "time point" that you want to start reading in the data from
-        time_max : int, optional
-            default : 110
+        time_max : int
             maximum "time point" that you want to end reading in the data from. 
         channel : int
             "channel"/thread value of a specific replica throughout the simulation "timecouse"
@@ -274,11 +272,9 @@ def __init__(self, path_to_data, time_min, time_max, channel, temperature = None
         ----------
         path_to_data : string
             path to either the wlcsim/data directory or the top level directory above the nested `trials`
-        time_min : int, optional
-            default : 0
+        time_min : int
             minimum "time point" that you want to start reading in the data from
-        time_max : int, optional
-            default : 110
+        time_max : int
             maximum "time point" that you want to end reading in the data from. 
         channel : int
             "channel"/thread value of a specific replica throughout the simulation "timecouse"
@@ -357,7 +353,7 @@ def playFineMovie(self,path=default_dir+'/analysis/pdb/',topo='linear',pymol='py
         os.system(pymol + " -r "+default_dir+"/analysis/movieFine.py -- " 
                     + str(self.time_max-self.time_min) + " " + path)
 
-    def playCoarseMovie(self,path=default_dir+'/analysis/pdb/',topo='linear',pymol='pymol'):
+    def playCoarseMovie(self, path = default_dir+'/analysis/pdb/', topo = 'linear', pymol = 'pymol', sphere_radius = 0, show_hull = True):
         """Play PyMol movie of the polymer throughout the simulation "timecourse" visualizing the excluded volume of the chain. 
         See the saveCoarseGrainedPDB method in the `Snapshot` class for more informtion on the calculation.
 
@@ -374,17 +370,23 @@ def playCoarseMovie(self,path=default_dir+'/analysis/pdb/',topo='linear',pymol='
         pymol : string, optional
             default : 'pymol'
             exectable command to initiaite PyMol, i.e. "~/Applications/pymol/pymol" 
+        sphere_radius : float, optional
+            default : 0
+            set what size radius to visualize a confining sphere, where 0 equates to no confinement
+        show_hull : boolean, optional
+            default : True
+            whether to construct the hulls of the excluded volume of the fiber or not
         """
         for time in range(self.time_min,self.time_max):
             self.snapshots[time].saveCoarseGrainedPDB(path=path,topo=topo)
         if (self.temperature != None):
             os.system(pymol + " -r "+default_dir+"/analysis/movieCoarse.py -- " 
-                        + str(self.time_max-self.time_min) + " PT " + str(self.temperature) + " " + path + " " + self.path_to_data)
+                        + str(self.time_max-self.time_min) + " PT " + str(self.temperature) + " " 
+                        + path + " " + self.path_to_data + " " + str(show_hull) + " " + str(sphere_radius))
         else:
             os.system(pymol + " -r "+default_dir+"/analysis/movieCoarse.py -- " 
                         + str(self.time_max-self.time_min) + " " + str(self.channel[-1]) + " 1 "  
-                        + path + " " + self.path_to_data)
-
+                        + path + " " + self.path_to_data + " " + str(show_hull) + " " + str(sphere_radius))
 class Snapshot:
     """
     `Snapshot` object to store all the positions and orientations of the computational beads for a given snapshot. 
@@ -423,7 +425,7 @@ def __init__(self,path_to_data,time,channel):
         ----------
         path_to_data : string
             path to either the wlcsim/data directory or the top level directory above the nested `trials`
-        time : int, optional
+        time : int
             "time point" of the current wlcsim output structure
         channel : int
             "channel"/thread value of the current wlcsim output structure

analysis/README

@@ -1,3 +1,14 @@
 this folder contains parser/helper files for wlcsim analysis
 
 it was written for the use of the risca lab, but anyone else is free to use/edit it!
+
+main code is MCsim which is a wlcsim output parser with some analysis methods, such as:
+    * visualize trajectory in pymol
+    * end-to-end distance
+    * pairwise nucleosome distances
+    * interpolate basepair resolution DNA structure between coarse-graied wlcsim beads
+    * determine ricc-seq break patterns from interpolated fiber structure
+
+some minimal helper functions also exist in utility, such as:
+    * apply rotation/translation of linker DNA into and out of nucleosome
+    * parametrize DNA as two offset helices to build a phosphate-base-phosphate representation of B-DNA

analysis/movieCoarse.py

@@ -3,6 +3,7 @@
 import numpy as np
 import colorsys
 import sys
+from pymol.cgo import *
 sys.path.append('../analysis/')
 from analysis.utility import *
 
@@ -13,6 +14,8 @@
 temperature = int(argv[3])
 pathPDB = argv[4]
 pathData = argv[5]
+hull = argv[6]
+rad = float(argv[7])
 if (channel == 'PT'):
     nodes = np.loadtxt(pathData+'nodeNumber')
     nodes = np.vstack((np.linspace(0, np.shape(nodes)[1]-1, np.shape(nodes)[1]), nodes))
@@ -35,70 +38,82 @@
     # load discretization data
     disc = np.loadtxt('/%sd%sv%s' %(input_folder,ind,channel[ind]))
     wrap = disc[0]; bps = disc[1]
-    for i in range(len(r)):
-        if bps[i] != 0:
-            poly = np.zeros(side*3).reshape([side,3])
-            uin = np.asarray(u[i,0:3]); vin = np.asarray(u[i,3:6]); cross = np.cross(uin, vin)
-            mat = np.matrix([vin, cross, uin]).reshape([3,3]).T
-            if (wrap[i]>1):
-                space = 0
-                height = 5.5
-                radius = 5.2
-                center = np.asarray([4.8455, -2.4445, 0.6694])
-                for j in range(int(side/2.0)):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.cos(space), -height/2.0, radius*np.sin(space)])
-                    poly[j,:] = r[i,:] + np.matmul(mat, center+vec)
-                    space = space + incr
-                space = 0
-                for j in range(int(side/2.0),side):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.cos(space), height/2.0, radius*np.sin(space)])
-                    poly[j,:] = r[i,:] + np.matmul(mat, center+vec)
-                    space = space + incr
+    if (hull == 'True'):
+        for i in range(len(r)):
+            if bps[i] != 0:
+                poly = np.zeros(side*3).reshape([side,3])
+                uin = np.asarray(u[i,0:3]); vin = np.asarray(u[i,3:6]); cross = np.cross(uin, vin)
+                mat = np.matrix([vin, cross, uin]).reshape([3,3]).T
+                if (wrap[i]>1):
+                    space = 0
+                    height = 5.5
+                    radius = 5.2
+                    center = np.asarray([4.8455, -2.4445, 0.6694])
+                    for j in range(int(side/2.0)):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.cos(space), -height/2.0, radius*np.sin(space)])
+                        poly[j,:] = r[i,:] + np.matmul(mat, center+vec)
+                        space = space + incr
+                    space = 0
+                    for j in range(int(side/2.0),side):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.cos(space), height/2.0, radius*np.sin(space)])
+                        poly[j,:] = r[i,:] + np.matmul(mat, center+vec)
+                        space = space + incr
+                    # make cylinders
+                    x1,y1,z1 = np.mean(poly[:int(side/2),0]), np.mean(poly[:int(side/2),1]), np.mean(poly[:int(side/2),2])
+                    (re, g, b) = colorsys.hsv_to_rgb(float(i)/(len(r)-1), 1.0, 1.0)
+                    x2,y2,z2 = np.mean(poly[int(side/2):,0]), np.mean(poly[int(side/2):,1]), np.mean(poly[int(side/2):,2])
+                    cmd.load_cgo( [ 25.0, 0.25, 9.0, x1, y1, z1, x2, y2, z2, radius, re, g, b, re, g, b ], "seg"+str(i+1)+'nuc')
+                    # add extruding linker
+                    space = 2*np.pi/side
+                    tempU, tempV, tempR = rotate_bead(uin,vin,r[i,:],int(bps[i]),int(wrap[i]))
+                    height = np.sqrt(np.dot(r[i+1,:]-tempR, r[i+1,:]-tempR))
+                    radius = 1.0
+                    center = np.zeros(3)
+                    tempMat = np.matrix([tempV, np.cross(tempU, tempV), tempU]).reshape([3,3]).T
+                    for j in range(int(side/2.0)):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.sin(space), radius*np.cos(space), -height/2.0])
+                        poly[j,:] = (tempR+r[i+1,:])/2.0 + np.matmul(tempMat, center+vec)
+                        space = space + incr
+                    space = 2*np.pi/side
+                    for j in range(int(side/2.0),side):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.sin(space), radius*np.cos(space), height/2.0])
+                        poly[j,:] = (tempR+r[i+1,:])/2.0 + np.matmul(tempMat, center+vec)
+                        space = space + incr
+                else:
+                    space = 2*np.pi/side
+                    height = np.sqrt(np.dot(r[i+1,:]-r[i,:], r[i+1,:]-r[i,:]))
+                    radius = 1.0
+                    center = np.zeros(3)
+                    for j in range(int(side/2.0)):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.sin(space), radius*np.cos(space), -height/2.0])
+                        poly[j,:] = (r[i,:]+r[i+1,:])/2.0 + np.matmul(mat, center+vec)
+                        space = space + incr
+                    space = 2*np.pi/side
+                    for j in range(int(side/2.0),side):
+                        # rotate into material frame 
+                        vec = np.asarray([radius*np.sin(space), radius*np.cos(space), height/2.0])
+                        poly[j,:] = (r[i,:]+r[i+1,:])/2.0 + np.matmul(mat, center+vec)
+                        space = space + incr
                 # make cylinders
                 x1,y1,z1 = np.mean(poly[:int(side/2),0]), np.mean(poly[:int(side/2),1]), np.mean(poly[:int(side/2),2])
                 (re, g, b) = colorsys.hsv_to_rgb(float(i)/(len(r)-1), 1.0, 1.0)
                 x2,y2,z2 = np.mean(poly[int(side/2):,0]), np.mean(poly[int(side/2):,1]), np.mean(poly[int(side/2):,2])
-                cmd.load_cgo( [ 25.0, 0.25, 9.0, x1, y1, z1, x2, y2, z2, radius, re, g, b, re, g, b ], "seg"+str(i+1)+'nuc')
-                # add extruding linker
-                space = 2*np.pi/side
-                tempU, tempV, tempR = rotate_bead(uin,vin,r[i,:],int(bps[i]),int(wrap[i]))
-                height = np.sqrt(np.dot(r[i+1,:]-tempR, r[i+1,:]-tempR))
-                radius = 1.0
-                center = np.zeros(3)
-                tempMat = np.matrix([tempV, np.cross(tempU, tempV), tempU]).reshape([3,3]).T
-                for j in range(int(side/2.0)):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.sin(space), radius*np.cos(space), -height/2.0])
-                    poly[j,:] = (tempR+r[i+1,:])/2.0 + np.matmul(tempMat, center+vec)
-                    space = space + incr
-                space = 2*np.pi/side
-                for j in range(int(side/2.0),side):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.sin(space), radius*np.cos(space), height/2.0])
-                    poly[j,:] = (tempR+r[i+1,:])/2.0 + np.matmul(tempMat, center+vec)
-                    space = space + incr
-            else:
-                space = 2*np.pi/side
-                height = np.sqrt(np.dot(r[i+1,:]-r[i,:], r[i+1,:]-r[i,:]))
-                radius = 1.0
-                center = np.zeros(3)
-                for j in range(int(side/2.0)):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.sin(space), radius*np.cos(space), -height/2.0])
-                    poly[j,:] = (r[i,:]+r[i+1,:])/2.0 + np.matmul(mat, center+vec)
-                    space = space + incr
-                space = 2*np.pi/side
-                for j in range(int(side/2.0),side):
-                    # rotate into material frame 
-                    vec = np.asarray([radius*np.sin(space), radius*np.cos(space), height/2.0])
-                    poly[j,:] = (r[i,:]+r[i+1,:])/2.0 + np.matmul(mat, center+vec)
-                    space = space + incr
-            # make cylinders
-            x1,y1,z1 = np.mean(poly[:int(side/2),0]), np.mean(poly[:int(side/2),1]), np.mean(poly[:int(side/2),2])
-            (re, g, b) = colorsys.hsv_to_rgb(float(i)/(len(r)-1), 1.0, 1.0)
-            x2,y2,z2 = np.mean(poly[int(side/2):,0]), np.mean(poly[int(side/2):,1]), np.mean(poly[int(side/2):,2])
-            cmd.load_cgo( [ 25.0, 0.25, 9.0, x1, y1, z1, x2, y2, z2, radius, re, g, b, re, g, b ], "seg"+str(i+1))
+                cmd.load_cgo( [ 25.0, 0.25, 9.0, x1, y1, z1, x2, y2, z2, radius, re, g, b, re, g, b ], "seg"+str(i+1))
 #cmd.mplay()
 cmd.orient('snap')
+
+# add sphere volume if set 
+if (rad>0):
+   spherelist = [
+      COLOR,    1,   0,    0,
+      SPHERE,   rad, rad, rad, rad,
+   ]
+   cmd.load_cgo(spherelist, 'sphere',   1)
+   cmd.set('cgo_transparency', 0.6)
+   cmd.set('cgo_sphere_quality', 100)
+