diff --git a/analysis/MCsim.py b/analysis/MCsim.py index 98514685..5f981775 100644 --- a/analysis/MCsim.py +++ b/analysis/MCsim.py @@ -88,11 +88,11 @@ 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 @@ -100,7 +100,7 @@ def __init__(self, path_to_data = "%s/data/" %(default_dir), trials = [''], traj 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 diff --git a/analysis/README b/analysis/README index 1a168fbb..cbb82ec2 100644 --- a/analysis/README +++ b/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 diff --git a/analysis/drawSphere.py b/analysis/drawSphere.py deleted file mode 100644 index 5fc015f2..00000000 --- a/analysis/drawSphere.py +++ /dev/null @@ -1,12 +0,0 @@ -from pymol.cgo import * -from pymol import cmd - - -spherelist = [ - COLOR, 1, 0, 0, - SPHERE, 250/2, 250/2, 250/2, 250/2, - ] - -cmd.load_cgo(spherelist, 'sphere', 1) -cmd.set('cgo_transparency', 0.6) -cmd.set('cgo_sphere_quality', 100) diff --git a/analysis/movieCoarse.py b/analysis/movieCoarse.py index 87a31984..d5dba716 100644 --- a/analysis/movieCoarse.py +++ b/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) +