IMP  2.2.1
The Integrative Modeling Platform
atomicDomino.py
1 import IMP
2 import subprocess
3 import random
4 import IMP.domino
5 import IMP.core
6 import IMP.rmf
7 import RMF
8 import time
9 import IMP.algebra
10 import types
11 import re
12 import sys
13 import operator
14 import os
15 import resource
16 import atomicDominoUtilities
17 
18 
19 class AtomicDomino:
20 
21  def __init__(self, model, protein, parameterFileName):
22  self.model = model
23  self.protein = protein
24  self.namesToParticles = atomicDominoUtilities.makeNamesToParticles(
25  protein)
26  self.readParameters(parameterFileName)
27  self.wroteNativeProtein = 0
28  self.maxMem = 0
29 
30  def readParameters(self, parameterFileName):
31  self.parameters = atomicDominoUtilities.readParameters(
32  parameterFileName)
33 
34  def logMemory(self):
35  currentMem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
36  currentMem *= 10 ** -6
37  print "log memory: %s" % currentMem
38  if (currentMem > self.maxMem):
39  self.maxMem = currentMem
40 
41  # Get parameter value. Throws dictionary exception if not found
42  def getParam(self, paramName):
43  paramValue = self.parameters[paramName]
44  return paramValue
45 
46  def loadDominoHelpers(self):
47  self.subsetStateScores = {}
48  self.subsetRestraintScores = {}
49  self.subsetStateFailures = {}
50  self.nodesToAssignmentCount = {}
51 
52  self.totalAssignments = 0
53  outputDir = self.getParam("output_directory")
54  mtreeCytoscapeAssignmentFile = self.getParam(
55  "mtree_cytoscape_assignment_file")
56  self.mtreeCytoscapeAssignmentFh = open(
57  os.path.join(outputDir, mtreeCytoscapeAssignmentFile), 'w')
58  self.mtreeCytoscapeAssignmentFh.write("Assignments\n")
59 
60  self.mtNamesToIndices = {}
61 
62  self.subsetNamesToAssignmentFiles = {}
63 
64  # Simple way to calculate run-time for certain methods and procedures.
65  # When reset is 0, compares the previous saved time to the current one
66  # and returns the difference (saving the current time too).
67  # When reset is 1, just saves the current time.
68  def logTimePoint(self, reset):
69  newTime = time.time()
70  if (reset == 0):
71  timeDifference = newTime - self.timePoint
72  timeDifference = round(timeDifference, 0)
73  self.timePoint = newTime
74  return timeDifference
75  else:
76 
77  self.timePoint = newTime
78  return newTime
79 
80  def createParticleStatesTable(self):
81 
82  for particleName in self.particleInfo.keys():
83 
84  statesToCenters = {}
85  dataArray = self.particleInfo[particleName]
86  for i in range(len(dataArray)):
87  [state, center] = dataArray[i]
88  statesToCenters[state] = center
89 
90  statesList = []
91  for stateIndex in sorted(statesToCenters.keys()):
92  vector3d = IMP.algebra.Vector3D(statesToCenters[stateIndex])
93  statesList.append(vector3d)
94  # print "appending particle %s state %s center %s" %
95  # (particleName, stateIndex, vector3d)
96  xyzStates = IMP.domino.XYZStates(statesList)
97  self.dominoPst.set_particle_states(
98  self.namesToParticles[particleName],
99  xyzStates)
100 
101  def quickParticleName(self, particle):
102  return atomicDominoUtilities.quickParticleName(particle)
103 
104  def filterAssignments(self, assignments, subset, nodeIndex, rssft):
105 
106  filteredSubsets = []
107  restraintList = []
108 
109  # make dependency graph for stats
110  for r in IMP.get_restraints([self.model.get_root_restraint_set()]):
111  restraintList.append(r)
112  dg = IMP.get_dependency_graph(restraintList)
113 
114  stateCounter = 0
115  passedCounter = 0
116 
117  # create hdf5AssignmentContainer
118 
119  sFilter = rssft.get_subset_filter(subset, filteredSubsets)
120  # check each unique state to see if passes filter
121  filteredAssignments = []
122  for assignment in assignments:
123 
124  if (sFilter is None or sFilter.get_is_ok(assignment)):
125  # add to assignment container if it passes
126  passedCounter += 1
127  filteredAssignments.append(assignment)
128 
129  stateCounter += 1
130  fraction = (passedCounter * 1.0) / (stateCounter * 1.0)
131  print "%s states passed out of %s total for this subset (fraction %s)" % (passedCounter, stateCounter, fraction)
132  if (passedCounter == 0):
133  print "subset %s had 0 assignments (out of %s) pass. Exiting..." % (subset, stateCounter)
134  sys.exit()
135 
136  return filteredAssignments
137 
138  # Convert the name of the subset to something more readable. Currently
139  # returning name as sorted list of atoms
140  def quickSubsetName(self, subset):
141  cleanName = self.cleanVertexName(subset)
142  atomNameList = cleanName.split(" ")
143  sortedList = sorted(atomNameList)
144 
145  name = " ".join(sortedList)
146  return name
147 
148  def getNodeIndexList(self):
149  return self.mt.get_vertices()
150 
151  def getLeafNodeIndexList(self):
152  leafNodeIndexList = []
153  for subset in self.subsets:
154  index = self.getMtIndexForSubset(self.quickSubsetName(subset))
155  leafNodeIndexList.append(index)
156  return leafNodeIndexList
157 
158  # Create Domino sampler and give it all the other domino objects it needs
159  def createSampler(self):
160  s = IMP.domino.DominoSampler(self.model, self.dominoPst)
161 
162  s.set_merge_tree(self.mt)
163  filterTables = []
164  filterTables.append(self.rssft)
165  s.set_subset_filter_tables(filterTables)
166 
167  s.set_assignments_table(self.lat)
168 
169  if (self.getParam("cross_subset_filtering") == 1):
170  s.set_use_cross_subset_filtering(1)
171 
172  self.sampler = s
173 
174  # Use the model restraint set to get the interaction graph, junction tree, and merge tree, and also get subsets
175  # from the junction tree and return them
176  def createSubsets(self):
177  self.initializeParticleStatesTable()
179  [self.model.get_root_restraint_set()],
180  self.dominoPst)
181  print "interaction graph:"
182  ig.show()
184  print "junction tree:"
185  jt.show()
186  self.subsetNamesToSubsets = {}
188 
189  # make map of vertex indices to atoms in subsets
190  for index in mt.get_vertices():
191  subset = mt.get_vertex_name(index)
192  subsetName = self.cleanVertexName(subset)
193  self.mtNamesToIndices[subsetName] = index
194 
195  subsets = IMP.domino.get_subsets(jt)
196  for subset in subsets:
197  print "created subset %s" % subset
198  subsetName = self.quickSubsetName(subset)
199  self.subsetNamesToSubsets[subsetName] = subset
200 
201  print "merge tree:"
202  mt.show()
203 
204  self.ig = ig
205  self.jt = jt
206  self.mt = mt
207  self.subsets = subsets
208 
209  self.parentSiblingMap = {}
210  self.parentSiblingMap[self.mt.get_vertices()[-1]] = {}
211  self.createSiblingMap(self.mt.get_vertices()[-1])
212 
213  def getMtIndexForSubset(self, subsetName):
214  for index in self.mt.get_vertices():
215  mtNode = self.mt.get_vertex_name(index)
216  mtName = self.cleanVertexName(mtNode)
217  mtName = mtName.rstrip()
218  mtName = mtName.lstrip()
219 
220  mtNameList = mtName.split(" ")
221  subsetNameList = subsetName.split(" ")
222  mtNameListSorted = sorted(mtNameList)
223  subsetNameListSorted = sorted(subsetNameList)
224  if (" ".join(subsetNameListSorted) == " ".join(mtNameListSorted)):
225 
226  return index
227  print "did not find merge tree index for subset name %s" % subsetName
228  sys.exit()
229 
230  def getDominoParticleNames(self):
231  particles = self.dominoPst.get_particles()
232  particleNameList = []
233  for particle in particles:
234  pName = self.quickParticleName(particle)
235  particleNameList.append(pName)
236 
237  return particleNameList
238 
239  def getSubsetNameList(self):
240  subsetNameList = []
241  for subset in self.subsets:
242  name = self.quickSubsetName(subset)
243  subsetNameList.append(name)
244  return subsetNameList
245 
246  def getMtIndexToParticles(self):
247 
248  mtIndexToParticles = {}
249  allVertices = self.mt.get_vertices()
250  for nodeIndex in allVertices:
251  subset = self.mt.get_vertex_name(nodeIndex)
252  particleList = self.quickSubsetName(subset)
253  mtIndexToParticles[nodeIndex] = particleList
254  print "createtMtIndexToParticles: index %s is particleList %s " % (nodeIndex, particleList)
255 
256  return mtIndexToParticles
257 
258  def readTrajectoryFile(
259  self,
260  atomList,
261  rh,
262  frames,
263  scoreOutputFile,
264  skipDomino,
265  flexibleAtoms):
266 
267  bestScore = 100000
268  bestRmsd = 10000
269  bestScoreFrame = 0
270  bestRmsdFrame = 0
271  scoreOutputFh = open(scoreOutputFile, 'w')
272  for i in frames:
273  try:
274  IMP.rmf.load_frame(rh, i, self.protein)
275  except Exception:
276  print "execption in loading frame %s" % i
277  continue
278  score = self.model.evaluate(False)
279  leaves = IMP.atom.get_leaves(self.protein)
280  # for leaf in leaves:
281  # xyzD = IMP.core.XYZ.decorate_particle(leaf)
282  # print "read trajectory file: coordinates for frame %s particle %s
283  # are %s" % (i, self.quickParticleName(leaf),
284  # xyzD.get_coordinates())
285 
286  rmsd = self.calculateNativeRmsd(flexibleAtoms)
287  scoreOutputFh.write("%s\t%s\t%s\n" % (i, score, rmsd))
288  if (score < bestScore):
289  bestScore = score
290  bestScoreFrame = i
291 
292  if (rmsd < bestRmsd):
293  bestRmsd = rmsd
294  bestRmsdFrame = i
295  if (skipDomino == 0):
296  for atomName in atomList:
297 
298  particle = self.namesToParticles[atomName]
299 
300  # get grid index and coordinates
301  gridIndex = self.snapToGrid(particle)
302  center = self.grid.get_center(gridIndex)
303  pythonCenter = []
304  for coordinate in center:
305  pythonCenter.append(coordinate)
306 
307  # grid indices are mapped to coordinate states. Check if
308  # we've seen this grid index
309  if ((gridIndex in self.particleStatesSeen[atomName]) == 0):
310  # set this particle state index to size of the
311  # dictionary, which effectively increments the index
312  # with each new state
313  currentSize = len(
314  self.particleStatesSeen[atomName].keys())
315  self.particleStatesSeen[
316  atomName][
317  gridIndex] = currentSize
318  state = self.particleStatesSeen[atomName][gridIndex]
319  self.particleInfo[atomName].append([state, pythonCenter])
320  else:
321  print "didn't add domino states due to skip parameters"
322  return [bestScore, bestScoreFrame, bestRmsd, bestRmsdFrame]
323 
324  # Get the grid index for the given particle. Returns an integer that can be used to get the center
325  # of the grid space for discretizing the particle
326  def getParticleGridIndex(self, leaf):
327  xyzDecorator = IMP.core.XYZ.decorate_particle(leaf)
328  coordinates = xyzDecorator.get_coordinates()
329  extendedIndex = 0
330  extendedIndex = self.grid.get_extended_index(coordinates)
331  if (self.grid.get_has_index(extendedIndex) == 0):
332  # mostly working with sparse grids now
333  if (self.getParam("grid_type") == "sparse"):
334  self.grid.add_voxel(extendedIndex, 1)
335  else:
336  self.grid.add_voxel(extendedIndex)
337 
338  index = self.grid.get_index(extendedIndex)
339  return index
340 
341  # Set coordinates of this atom to be the center of the grid space
342  # containing it (effectiely discretizing the system)
343  def snapToGrid(self, particle):
344 
345  index = self.getParticleGridIndex(particle)
346  center = self.grid.get_center(index)
347  xyzDecorator = IMP.core.XYZ.decorate_particle(particle)
348  xyzDecorator.set_coordinates(center)
349 
350  return index
351 
352  # Take initial protein and snap every atom to the center of its gridpoint
353  def discretizeNativeProtein(self):
354 
355  outputDir = self.getParam("output_directory")
356  nativeSnappedFile = self.getParam("native_protein_snapped_output_file")
357 
358  leaves = IMP.atom.get_leaves(self.protein)
359  for leaf in leaves:
360  self.snapToGrid(leaf)
361 
363  self.protein,
364  os.path.join(
365  outputDir,
366  nativeSnappedFile))
367 
368  # Get particle representing the alpha carbon at the center of the peptide
369  def getPeptideCa(self):
370 
371  # get all residue indices in the peptide
372  residues = IMP.atom.get_by_type(self.protein, IMP.atom.RESIDUE_TYPE)
373  peptideIndicesToResidues = {}
374  for residueH in residues:
375  chain = IMP.atom.get_chain(residueH)
376  chainId = chain.get_id()
377  residue = residueH.get_as_residue()
378  if (chainId == self.getParam("peptide_chain")):
379  peptideIndicesToResidues[residue.get_index()] = residue
380 
381  # use the min and max residue indices to get the residue in the middle
382  # (rounding down)
383  minPeptide = min(sorted(peptideIndicesToResidues.keys()))
384  maxPeptide = max(sorted(peptideIndicesToResidues.keys()))
385  centerPeptide = round(((maxPeptide - minPeptide) / 2 + minPeptide), 0)
386 
387  # get the particle corresponding to the ca atom at the middle residue
388  # and return it,
389  centerName = atomicDominoUtilities.makeParticleName(
390  self.getParam("peptide_chain"),
391  int(centerPeptide),
392  "CA")
393  centerParticle = self.namesToParticles[centerName]
394  return centerParticle
395 
396  # Create grid object that will be used to create discrete states for each
397  # particle
398  def createGrid(self):
399 
400  protBb = IMP.atom.get_bounding_box(self.protein)
401 
402  gridSpacing = float(self.getParam("grid_spacing"))
403  bufferSpace = float(self.getParam("grid_buffer_space"))
404 
405  protBb += bufferSpace # add buffer around grid
406  g = 0
407  if (self.getParam("grid_type") == "sparse"):
408  ca = self.getPeptideCa()
411  gridSpacing, xyzCa.get_coordinates())
412  else:
413  g = IMP.algebra.DenseDoubleGrid3D(gridSpacing, protBb)
414 
415  self.grid = g
416 
417  # Create Particle States Table and for each particle in the system, add XYZStates with that particle's
418  # initial location
419  def initializeParticleStatesTable(self):
420 
421  dominoPst = IMP.domino.ParticleStatesTable()
422  restrainedParticles = atomicDominoUtilities.getRestrainedParticles(
423  self.protein, self.model, self.namesToParticles)
424 
425  for p in restrainedParticles:
426 
428  xyz = IMP.core.XYZ(p).get_coordinates()
429  xyzStates = IMP.domino.XYZStates([xyz])
430  dominoPst.set_particle_states(p, xyzStates)
431 
432  self.dominoPst = dominoPst
433 
434  def createUniqueLeafAssignments(self, particleNameList, particleInfo):
435 
436  size = len(particleInfo[particleNameList[0]])
437  allAssignments = []
438 
439  for i in range(size):
440  nextAssignment = []
441  for particleName in particleNameList:
442  dataArray = particleInfo[particleName]
443  [state, center] = dataArray[i]
444  nextAssignment.append(state)
445  allAssignments.append(nextAssignment)
446 
447  # make unique assignments to avoid duplicates
448  uniqueAssignments = {}
449  for assignment in allAssignments:
450 
451  stateString = ""
452  for index in assignment:
453  stateString = stateString + str(index) + "_"
454  uniqueAssignments[stateString] = assignment
455  finalAssignments = []
456 
457  # add all assignments to final list
458  for stateString in uniqueAssignments.keys():
459  assignmentList = uniqueAssignments[stateString]
460  assignment = IMP.domino.Assignment(assignmentList)
461  finalAssignments.append(assignment)
462  return finalAssignments
463 
464  # Read MD trajectory; for each particle, save all unique states, and for
465  # each subset, save all assignments
466  def readMdTrajectory(self, atomList, flexibleAtoms):
467 
468  # open trajectory file
469  outputDir = self.getParam("output_directory")
470  trajectoryFile = self.getParam("md_trajectory_output_file")
471  fullFile = os.path.join(outputDir, trajectoryFile)
472  rh = RMF.open_rmf_file(fullFile)
473  IMP.rmf.set_hierarchies(rh, [self.protein])
474  framesToRead = atomicDominoUtilities.getMdIntervalFrames(
475  rh, int(self.getParam("md_interval")), self.protein)
476  print "preparing to read md frames %s" % framesToRead
477  # prepare data structures for tracking
478  particleInfo = {}
479  # for each particle, list where each entry corresponds to an md step, and its value [domino state, coordinates]
480  # for each particle, dictionary where the key is the grid index and
481  # value is domino state
482  particleStatesSeen = {}
483  for atomName in atomList:
484  particle = self.namesToParticles[atomName]
485  particleInfo[atomName] = []
486  particleStatesSeen[atomName] = {}
487 
488  self.particleStatesSeen = particleStatesSeen
489  self.particleInfo = particleInfo
490 
491  # read trajectory file
492  mdScoreOutputFile = os.path.join(
493  outputDir, "%s" %
494  self.getParam("md_score_output_file"))
495  [bestMdScore,
496  bestScoreFrame,
497  bestRmsd,
498  bestRmsdFrame] = self.readTrajectoryFile(
499  atomList,
500  rh,
501  framesToRead,
502  mdScoreOutputFile,
503  0,
504  flexibleAtoms)
505 
506  # output best score information
507  # print "try loading bad frame"
508  self.singlePdbResults(
509  fullFile,
510  bestScoreFrame,
511  self.getParam("best_md_score_output_file"))
512  #self.singlePdbResults(fullFile, 10000, self.getParam("best_md_score_output_file"))
513 
514  self.singlePdbResults(
515  fullFile,
516  bestRmsdFrame,
517  self.getParam("best_md_rmsd_output_file"))
518  self.singlePdbResults(
519  fullFile,
520  -1,
521  self.getParam("final_md_frame_output_file"))
522 
523  self.bestMdScore = round(bestMdScore, 2)
524  self.bestMdRmsd = round(bestRmsd, 2)
525  self.bestMdScoreFrame = bestScoreFrame
526  self.bestMdRmsdFrame = bestRmsdFrame
527 
528  def readCgTrajectories(self, atomList, flexibleAtoms):
529 
530  cgFileName = self.getParam("cg_output_file")
531  bestCgScore = 10000000
532  bestCgScoreFile = ""
533  bestCgRmsd = 10000000
534  bestCgRmsdFile = ""
535 
536  outputDir = self.getParam("output_directory")
537  trajectoryFile = self.getParam("md_trajectory_output_file")
538  fullFile = os.path.join(outputDir, trajectoryFile)
539  rh = RMF.open_rmf_file(fullFile)
540  IMP.rmf.set_hierarchies(rh, [self.protein])
541  framesToRead = atomicDominoUtilities.getMdIntervalFrames(
542  rh, int(self.getParam("cg_interval")), self.protein)
543 
544  skipCgDomino = int(self.getParam("skip_cg_domino"))
545 
546  if (len(framesToRead) > 0):
547  for cgNumber in framesToRead:
548  # Open next cg trajectory
549  outputDir = self.getParam("output_directory")
550  fullCgFileName = os.path.join(
551  outputDir, "%s%s" %
552  (cgFileName, cgNumber))
553  rh = RMF.open_rmf_file(fullCgFileName)
554  IMP.rmf.set_hierarchies(rh, [self.protein])
555 
556  # Only look at the bottom 20 frames
557  frameCount = IMP.rmf.get_number_of_frames(rh, self.protein)
558  cgFrames = []
559  startFrameCount = 0
560  if (frameCount > 20):
561  startFrameCount = frameCount - 20
562 
563  for i in range(startFrameCount, frameCount):
564  cgFrames.append(i)
565 
566  # Process trajectory
567  cgScoreOutputFile = os.path.join(
568  outputDir, "%s%s" %
569  (self.getParam("cg_score_output_file"), cgNumber))
570  [cgScore, cgScoreFrame, cgRmsd, cgRmsdFrame] = self.readTrajectoryFile(
571  atomList, rh, cgFrames, cgScoreOutputFile, skipCgDomino, flexibleAtoms)
572  print "cg number %s rmsd %s score %s" % (cgNumber, cgRmsd, cgScore)
573  # Update best score
574  if (cgScore < bestCgScore):
575  bestCgScore = cgScore
576  bestCgScoreFile = fullCgFileName
577  if (cgRmsd < bestCgRmsd):
578  bestCgRmsd = cgRmsd
579  bestCgRmsdFile = fullCgFileName
580 
581  # output best score information
582  self.singlePdbResults(
583  bestCgScoreFile,
584  -1,
585  self.getParam("best_cg_score_output_file"))
586  self.singlePdbResults(
587  bestCgRmsdFile,
588  -1,
589  self.getParam("best_cg_rmsd_output_file"))
590  self.singlePdbResults(
591  "%s%s" %
592  (cgFileName, framesToRead[-1]), -1, self.getParam("final_cg_frame_output_file"))
593  finalCgRmsd = self.calculateNativeRmsd(flexibleAtoms)
594  print "final cg rmsd is %s " % finalCgRmsd
595  self.bestCgScore = round(bestCgScore, 2)
596  self.bestCgRmsd = round(bestCgRmsd, 2)
597  self.bestCgScoreFile = bestCgScoreFile
598  self.bestCgRmsdFile = bestCgRmsdFile
599 
600  def singlePdbResults(self, trajectoryFile, frame, outputPdbFile):
601 
602  fullTrajectoryFile = os.path.join(
603  self.getParam("output_directory"),
604  trajectoryFile)
605  fullOutputFile = os.path.join(
606  self.getParam("output_directory"),
607  outputPdbFile)
608  rh = RMF.open_rmf_file(fullTrajectoryFile)
609  IMP.rmf.set_hierarchies(rh, [self.protein])
610  if (frame == -1):
611  frame = IMP.rmf.get_number_of_frames(rh, self.protein) - 1
612  IMP.rmf.load_frame(rh, frame, self.protein)
613  IMP.atom.write_pdb(self.protein, fullOutputFile)
614 
615  def calculateRmsd(self, otherProtein, flexibleAtoms):
616  otherNamesToParticles = atomicDominoUtilities.makeNamesToParticles(
617  otherProtein)
618  otherVector = []
619  modelVector = []
620  for pName in otherNamesToParticles.keys():
621  if ((pName in flexibleAtoms) == 0):
622  continue
623  otherParticle = otherNamesToParticles[pName]
624  modelParticle = self.namesToParticles[pName]
625  otherVector.append(
627  otherParticle).get_coordinates(
628  ))
629  modelVector.append(
631  modelParticle).get_coordinates(
632  ))
633  rmsd = IMP.atom.get_rmsd(otherVector, modelVector)
634  return rmsd
635 
636  def calculateNativeRmsd(self, flexibleAtoms):
637 
638  if (self.wroteNativeProtein == 0):
639  pdbName = self.getParam("native_pdb_input_file")
640  self.nativeModel = IMP.kernel.Model()
641  self.nativeProtein = IMP.atom.read_pdb(
642  pdbName,
643  self.nativeModel,
645  self.wroteNativeProtein = 1
646 
647  return self.calculateRmsd(self.nativeProtein, flexibleAtoms)
648 
649  def calculateTrajectoryRmsd(
650  self,
651  trajectoryFile,
652  trajectoryFrame,
653  flexibleAtoms):
654  pdbName = self.getParam("native_pdb_input_file")
655  otherModel = IMP.kernel.Model()
656  otherProtein = IMP.atom.read_pdb(
657  pdbName,
658  self.nativeModel,
660  outputDir = self.getParam("output_directory")
661  fullFile = os.path.join(outputDir, trajectoryFile)
662  print "open calculate traj rmf %s" % fullFile
663  rh = RMF.open_rmf_file(fullFile)
664  IMP.rmf.set_hierarchies(rh, [otherProtein])
665  if (trajectoryFrame == -1):
666  trajectoryFrame = IMP.rmf.get_number_of_frames(
667  rh, otherProtein) - 1
668  IMP.rmf.load_frame(rh, trajectoryFrame, otherProtein)
669  return self.calculateRmsd(otherProtein, flexibleAtoms)
670 
671  def createAllSubsetAssignments(self):
672 
675  self.model.get_root_restraint_set(),
676  self.dominoPst)
677 
678  leafNodeIndexList = self.getLeafNodeIndexList()
679 
680  for nodeIndex in leafNodeIndexList:
681  # get subset for this leaf
682  subset = self.mt.get_vertex_name(nodeIndex)
683  particleNameList = []
684  for particle in subset:
685  particleNameList.append(self.quickParticleName(particle))
686  print "creating initial assignments for leaf %s" % nodeIndex
687  # use particleInfo to create assignments and filter them
688  assignments = self.createUniqueLeafAssignments(
689  particleNameList, self.particleInfo)
690  filteredAssignments = self.filterAssignments(
691  assignments, subset, nodeIndex, rssft)
692 
693  # add assignemtns to container and listAssignmentTable
694  packedAssignmentContainer = IMP.domino.PackedAssignmentContainer()
695  for assignment in filteredAssignments:
696  packedAssignmentContainer.add_assignment(assignment)
697  lat.set_assignments(subset, packedAssignmentContainer)
698 
699  self.lat = lat
700  self.rssft = rssft
701 
702  def runDomino(self):
703  root = self.mt.get_vertices()[-1]
704  completeAc = self.loadAssignments(root)
705  self.completeAc = completeAc
706 
707  def loadAssignments(self, nodeIndex):
708 
709  children = self.mt.get_out_neighbors(nodeIndex)
710  subset = self.mt.get_vertex_name(nodeIndex)
711  heapCount = int(self.getParam("heap_count"))
713  heapCount, self.rssft.get_subset_filter(subset, []))
714  if len(children) == 0:
715  print "computing assignments for leaf %s" % nodeIndex
716 
717  self.sampler.load_vertex_assignments(nodeIndex, mine)
718  print "leaf node %s has %s leaf assignments" % (nodeIndex, mine.get_number_of_assignments())
719  else:
720  if (children[0] > children[1]):
721  children = [children[1], children[0]]
722  # recurse on the two children
723  firstAc = self.loadAssignments(children[0])
724  secondAc = self.loadAssignments(children[1])
725  self.logTimePoint(1)
726  self.sampler.load_vertex_assignments(
727  nodeIndex, firstAc, secondAc, mine)
728 
729  timeDifference = self.logTimePoint(0)
730  print "Done Parent %s Assignments %s first child %s second child %s time %s" % (nodeIndex, mine.get_number_of_assignments(), firstAc.get_number_of_assignments(),
731  secondAc.get_number_of_assignments(), timeDifference)
732  self.totalAssignments += mine.get_number_of_assignments()
733  self.logMemory()
734  return mine
735 
736  def writeOutput(self, flexibleAtoms, startTime):
737  bestAssignment = -1
738  bestDominoScore = 100000
739  bestAssignment = 0
740  finalAssignments = self.completeAc.get_assignments()
741  for assignment in finalAssignments:
743  self.dominoPst.get_subset(),
744  assignment,
745  self.dominoPst)
746  score = self.model.evaluate(False)
747  if (score < bestDominoScore):
748  bestAssignment = assignment
749  bestDominoScore = round(score, 2)
750  print "best domino score is %s " % bestDominoScore
751  print "best md score is %s" % self.bestMdScore
752  print "best md rmsd is %s" % self.bestMdRmsd
753  print "best cg score is %s" % self.bestCgScore
754  print "best cg rmsd is %s" % self.bestCgRmsd
755  print "merge tree contained %s total assignments" % self.totalAssignments
756 
758  self.dominoPst.get_subset(),
759  bestAssignment,
760  self.dominoPst)
761  dominoVsMdRmsd = round(
762  self.calculateTrajectoryRmsd(
763  self.getParam(
764  "md_trajectory_output_file"),
765  self.bestMdScoreFrame,
766  flexibleAtoms),
767  2)
768  cg = IMP.core.ConjugateGradients(self.model)
769  cg.optimize(100)
771  self.protein,
772  os.path.join(self.getParam("output_directory"),
773  self.getParam("minimum_domino_score_pdb")))
774 
775  dominoVsCgRmsd = round(
776  self.calculateTrajectoryRmsd(
777  self.bestCgScoreFile,
778  -1,
779  flexibleAtoms),
780  2)
781  dominoMinimizedScore = round(self.model.evaluate(False), 2)
782  dominoRmsd = round(self.calculateNativeRmsd(flexibleAtoms), 2)
783 
784  runTime = round(time.time() - startTime, 2)
785  print "final domino score (after cg): %s" % dominoMinimizedScore
786  print "final domino rmsd: %s" % dominoRmsd
787  print "best domino rmsd with best md score: %s" % dominoVsMdRmsd
788  print "domino rmsd with best cg score: %s" % dominoVsCgRmsd
789  print "Final Results\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.bestMdScore, self.bestMdRmsd, self.bestCgScore, self.bestCgRmsd, bestDominoScore, dominoRmsd, dominoMinimizedScore, dominoVsCgRmsd, self.totalAssignments, self.maxMem, runTime)
790 
791  #
792  # Parallel methods
793  #
794  def createSubsetFromParticles(self, particleNames):
795  particleNameList = particleNames.split(" ")
796  particleList = []
797  for particleName in particleNameList:
798  particle = self.namesToParticles[particleName]
799  particleList.append(particle)
800 
801  subset = IMP.domino.Subset(particleList)
802  return [subset, particleList]
803 
804  def createHdf5AssignmentContainer(self, index, particleNames, read):
805  root = self.getAssignmentContainerRoot(index, read)
806  print "got root for index %s" % index
807  dataset = 0
808  if (read == 1):
809  dataset = root.get_child_index_data_set_2d(str(index))
810  else:
811  dataset = root.add_child_index_data_set_2d(str(index))
812  print "added child index dataset"
813  # firstDataset.set_size([
814  [subset, particleOrder] = self.createSubsetFromParticles(particleNames)
815  print "created subset for index %s" % index
816  hdf5Ac = IMP.domino.create_assignments_container(
817  dataset, subset, particleOrder)
818  print "returning from create"
819  order = IMP.domino.get_order(subset, particleOrder)
820  for nextInt in order:
821  print "next int is %s" % nextInt
822  return [subset, hdf5Ac]
823 
824  def loadAssignmentsParallel(
825  self,
826  nodeIndex,
827  particleInfo,
828  mtIndexToNodeInfo,
829  mtIndexToSubsetOrder,
830  mtIndexToParticles):
831  IMP.base.set_log_level(IMP.WARNING)
832 
833  if (("firstChild" in mtIndexToNodeInfo[nodeIndex]) == 0):
834  print "writing file for leaf index %s" % nodeIndex
835  return (
836  self.createAssignmentsParallel(
837  particleInfo,
838  nodeIndex,
839  mtIndexToParticles)
840  )
841 
842  else:
843  beginTime = self.logTimePoint(1)
844  firstChildIndex = mtIndexToNodeInfo[nodeIndex]["firstChild"]
845  [firstSubset, firstAc] = self.createHdf5AssignmentContainer(
846  firstChildIndex, mtIndexToSubsetOrder[firstChildIndex], 1)
847  firstAcCreateTime = self.logTimePoint(0)
848 
849  secondChildIndex = mtIndexToNodeInfo[nodeIndex]["secondChild"]
850  [secondSubset, secondAc] = self.createHdf5AssignmentContainer(
851  secondChildIndex, mtIndexToSubsetOrder[secondChildIndex], 1)
852  secondAcCreateTime = self.logTimePoint(0)
853 
854  print "getting assignments for nodeIndex %s first child %s second child %s" % (nodeIndex, firstChildIndex, secondChildIndex)
855  firstChildParticles = mtIndexToSubsetOrder[firstChildIndex]
856  secondChildParticles = mtIndexToSubsetOrder[secondChildIndex]
857  myParticles = mtIndexToParticles[nodeIndex]
858 
859  print "first child particles %s\nsecond child particles %s\nmy particles; %s\n" % (firstChildParticles, secondChildParticles, myParticles)
860  for p in firstSubset:
861  print "next particle in first subset: %s" % self.quickParticleName(p)
862 
863  for p in secondSubset:
864  print "next particle in second subset: %s" % self.quickParticleName(p)
865 
866  for assignment in firstAc.get_assignments():
867  print "next assignment for first child %s: %s" % (firstChildIndex, assignment)
868 
869  for assignment in secondAc.get_assignments():
870  print "next assignment for second child %s: %s" % (secondChildIndex, assignment)
871 
872  [mySubset, myAc] = self.createHdf5AssignmentContainer(
873  nodeIndex, mtIndexToParticles[nodeIndex], 0)
874  print "done creating hdf5"
875  prepTime = self.logTimePoint(0)
877  self.model.get_root_restraint_set(),
878  self.dominoPst)
879  rssf = rssft.get_subset_filter(mySubset, [])
880 
881  #heapAc = IMP.domino.HeapAssignmentContainer(1000, rssf)
882 
884  firstSubset,
885  firstAc,
886  secondSubset,
887  secondAc,
888  [rssft],
889  myAc)
890 
891  heapTime = self.logTimePoint(0)
892  # myAc.add_assignments(heapAc.get_assignments())
893  addTime = self.logTimePoint(0)
894  for assignment in myAc.get_assignments():
895  print "loadAssignmentsParallel next assignment for %s: %s" % (nodeIndex, assignment)
896  doneTime = self.logTimePoint(0)
897  firstChildCount = firstAc.get_number_of_assignments()
898  secondChildCount = secondAc.get_number_of_assignments()
899  print "first count: %s second count: %s begin: %s firstAc: %s secondAc: %s prep: %s heap: %s add: %s done: %s" % (firstChildCount, secondChildCount, beginTime, firstAcCreateTime, secondAcCreateTime, prepTime, heapTime, addTime, doneTime)
900  subsetOrder = self.getSubsetOrderList(mySubset)
901  return subsetOrder
902 
903  def createAssignmentsParallel(
904  self,
905  particleInfo,
906  nodeIndex,
907  mtIndexToParticles):
908 
909  subsetName = mtIndexToParticles[nodeIndex]
910  print "starting assignments parallel leaf index %s subset name %s" % (nodeIndex, subsetName)
911  [subset, particleList] = self.createSubsetFromParticles(subsetName)
912 
913  particleNameList = subsetName.split(" ")
914 
915  # create assignment by reading states
916  finalAssignments = self.createUniqueLeafAssignments(
917  particleNameList, particleInfo)
918 
919  #lat = IMP.domino.ListAssignmentsTable()
920  #lat.set_assignments(subset, finalAssignmentContainer)
921  # self.sampler.set_assignments_table(lat)
922  #hdf5AssignmentContainer = IMP.domino.HDF5AssignmentContainer(dataset, subset, self.dominoPst.get_particles(), subsetName)
924  self.model.get_root_restraint_set(),
925  self.dominoPst)
926  filteredAssignments = self.filterAssignments(
927  finalAssignments, subset, nodeIndex, rssft)
928  root = self.getAssignmentContainerRoot(nodeIndex, 0)
929  dataset = root.add_child_index_data_set_2d(str(nodeIndex))
930  dataset.set_size([0, len(subset)])
931 
932  hdf5AssignmentContainer = IMP.domino.create_assignments_container(
933  dataset, subset, particleOrder)
934  for assignment in filteredAssignments:
935  hdf5AssignmentContainer.add_assignment(assignment)
936  for assignment in hdf5AssignmentContainer.get_assignments():
937  print "hdf5 assignment container node %s next assignmet %s" % (nodeIndex, assignment)
938 
939  #self.checkAssignments(subset, nodeIndex, particleOrder)
940  subsetOrder = self.getSubsetOrderList(subset)
941  print "leaf node returning with order %s" % subsetOrder
942  return subsetOrder
943 
944  # Write pymol session files for the interactions across atoms and all
945  # subsets
946  def writePymolData(self):
947 
948  outputDir = self.getParam("output_directory")
950  pymolInteractions = self.getParam("pymol_interactions_file")
951  w = IMP.display.PymolWriter(os.path.join(outputDir, pymolInteractions))
952 
953  for gg in geometry:
954  w.add_geometry(gg)
955 
956  pymolSubsets = self.getParam("pymol_subsets_file")
957  geometry = IMP.domino.get_subset_graph_geometry(self.jt)
958  w = IMP.display.PymolWriter(os.path.join(outputDir, pymolSubsets))
959  for gg in geometry:
960  w.add_geometry(gg)
961 
962  # Clean the default name of the vertex (in brackets and with each atom
963  # contained in quotes) and return a string where [] and " are removed
964  def cleanVertexName(self, vertexName):
965 
966  # not sure if any vertices still have a Subset prefix but keeping
967  # anyway
968  nodeRe = re.compile('Subset\(\[(.*?)\s*\]')
969  secondNodeRe = re.compile('\[(.*?)\s*\]') # atom name
970  node = nodeRe.search(str(vertexName))
971  secondNode = secondNodeRe.search(str(vertexName))
972  vertexNameFinal = ""
973  foundName = 0
974  if node:
975  foundName = node.group(1)
976  if secondNode:
977  foundName = secondNode.group(1)
978  vertexNameFinal = foundName.replace('"', '')
979  return vertexNameFinal
980 
981  def getSubsetOrderList(self, subset):
982  subsetOrderList = []
983  for particle in subset:
984  name = self.quickParticleName(particle)
985  subsetOrderList.append(name)
986  subsetOrder = " ".join(subsetOrderList)
987  return subsetOrder
988 
989  def checkAssignments(self, subset, nodeIndex, particleOrder):
990  print "reading back in assignments for leaf index %s" % nodeIndex
991  root = self.getAssignmentContainerRoot(nodeIndex, 1)
992  dataset = root.get_child_index_data_set_2d(str(nodeIndex))
993  hdf5 = IMP.domino.create_assignments_container(
994  dataset, subset, particleOrder)
995  for assignment in hdf5.get_assignments():
996  print "leaf index %s read back in %s" % (nodeIndex, assignment)
997 
998  def createSamplerLite(self):
999  s = IMP.domino.DominoSampler(self.model, self.dominoPst)
1000 
1001  if (self.getParam("cross_subset_filtering") == 1):
1002  s.set_use_cross_subset_filtering(1)
1003 
1004  self.sampler = s
1005 
1006  def createSiblingMap(self, parentIndex):
1007 
1008  children = self.mt.get_out_neighbors(parentIndex)
1009  if (len(children) > 0):
1010  firstChild = children[0]
1011  secondChild = children[1]
1012  self.parentSiblingMap[firstChild] = {}
1013  self.parentSiblingMap[firstChild]["sibling"] = secondChild
1014  self.parentSiblingMap[firstChild]["parent"] = parentIndex
1015 
1016  self.parentSiblingMap[secondChild] = {}
1017  self.parentSiblingMap[secondChild]["sibling"] = firstChild
1018  self.parentSiblingMap[secondChild]["parent"] = parentIndex
1019  print "created map for parent %s first child %s second child %s" % (parentIndex, firstChild, secondChild)
1020 
1021  self.parentSiblingMap[parentIndex]["firstChild"] = firstChild
1022  self.parentSiblingMap[parentIndex]["secondChild"] = secondChild
1023 
1024  self.createSiblingMap(firstChild)
1025  self.createSiblingMap(secondChild)
1026 
1027  def getMtIndexToNodeInfo(self):
1028  return self.parentSiblingMap
1029 
1030  def getLeafParentNodeIndexList(self):
1031  leafParentNodeIndexList = {}
1032  leafNodeIndexList = self.getLeafNodeIndexList()
1033  for leafIndex in leafNodeIndexList:
1034  parent = self.parentSiblingMap[leafIndex]["parent"]
1035  leafParentNodeIndexList[parent] = 1
1036  return leafParentNodeIndexList
1037 
1038  def getMtIndexToNameList(selt):
1039  mtIndexToNames = {}
1040  for index in self.mt.get_vertices():
1041  name = self.mt.get_vertex_name(index)
1042  mtIndexToNames[index] = name
1043  return mtIndexToNames
1044 
1045  def getAssignmentContainerRoot(self, subsetIndex, read):
1046  outputDir = self.getParam("output_directory")
1047  filePrefix = self.getParam("subset_assignment_output_file")
1048  assignmentFileName = os.path.join(
1049  outputDir, "%s%s" %
1050  (filePrefix, subsetIndex))
1051  print "creating hdf5 file with name %s" % assignmentFileName
1052  root = 0
1053  if (read == 1):
1054  root = RMF.open_hdf5_file(assignmentFileName)
1055  else:
1056  root = RMF.create_hdf5_file(assignmentFileName)
1057  return root
1058 
1059  #
1060  # Begin Cytoscape Methods
1061  #
1062 
1063  def writeVisualization(self):
1064 
1065  self.writeCytoscapeIgInput()
1066  self.writeCytoscapeJtInput()
1067  self.writeCytoscapeMtInput()
1068  self.writeCytoscapeScripts()
1069  self.writePymolData()
1070 
1071  def writeCytoscapeScripts(self):
1072  outputDir = self.getParam("output_directory")
1073  mTreeCytoscapeInput = self.getParam("mtree_cytoscape_input_file")
1074  mTreeCytoscapeAssignments = self.getParam(
1075  "mtree_cytoscape_assignment_file")
1076  mTreeCytoscapeAtomChains = self.getParam(
1077  "mtree_cytoscape_atom_chain_file")
1078  mTreeCytoscapeAtomSummary = self.getParam(
1079  "mtree_cytoscape_atom_summary_file")
1080 
1081  mTreeCytoscapeScript = self.getParam("mtree_cytoscape_script")
1082  mTreeFh = open(os.path.join(outputDir, mTreeCytoscapeScript), 'w')
1083  mTreeFh.write(
1084  "network import file=%s\n" %
1085  os.path.join(outputDir, mTreeCytoscapeInput))
1086  mTreeFh.write(
1087  "node import attributes file=\"%s\"\n" %
1088  os.path.join(outputDir, mTreeCytoscapeAssignments))
1089  mTreeFh.write(
1090  "node import attributes file=\"%s\"\n" %
1091  os.path.join(outputDir, mTreeCytoscapeAtomSummary))
1092  mTreeFh.write(
1093  "node import attributes file=\"%s\"\n" %
1094  os.path.join(outputDir, mTreeCytoscapeAtomChains))
1095  mTreeFh.write("layout jgraph-tree\n")
1096 
1097  def getGraphStructure(self, graph, fileName, separator):
1098  # write junction tree to file
1099  graphLogWrite = open(fileName, 'w')
1100  graph.show(graphLogWrite)
1101  graphLogWrite.close()
1102 
1103  # read file
1104  graphLogRead = open(fileName, 'r')
1105 
1106  nodeRe = re.compile('^(\d+)\[label\=\"\[*(.*?)\s*\]*\"') # atom name
1107  separatorEscape = "\\" + separator
1108  edgeString = "^(\d+)\-%s(\d+)" % separatorEscape
1109  edgeRe = re.compile(edgeString)
1110 
1111  nodesToNodes = {}
1112  # keys: source node, value: target node (track edges)
1113  nodesToNames = {} # keys: node number, value; string parsed from file
1114 
1115  for line in graphLogRead:
1116 
1117  # search for nodes
1118  node = nodeRe.search(line)
1119  if node:
1120  nodeNumber = node.group(1)
1121  atomString = node.group(2)
1122  nodesToNames[nodeNumber] = atomString
1123  continue
1124 
1125  # search for edges
1126  edge = edgeRe.search(line)
1127  if edge:
1128  firstNode = edge.group(1)
1129  secondNode = edge.group(2)
1130  firstNodeDict = {}
1131  if (firstNode in nodesToNodes):
1132  firstNodeDict = nodesToNodes[firstNode]
1133  firstNodeDict[secondNode] = 1
1134  nodesToNodes[firstNode] = firstNodeDict
1135 
1136  return [nodesToNames, nodesToNodes]
1137 
1138  def writeEdgeFile(self, nodesToNodes, edgeFileName):
1139  # write edge file
1140  outputDir = self.getParam("output_directory")
1141  graphInputFile = open(os.path.join(outputDir, edgeFileName), 'w')
1142  for firstNode in nodesToNodes.keys():
1143  nodeDict = nodesToNodes[firstNode]
1144  for secondNode in nodeDict.keys():
1145  graphInputFile.write("%s ttt %s\n" % (firstNode, secondNode))
1146  graphInputFile.close()
1147 
1148  def writeCytoscapeIgInput(self):
1149  outputDir = self.getParam("output_directory")
1150  igOutputFile = self.getParam("ig_output_file")
1151  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1152  self.ig, os.path.join(outputDir, igOutputFile), "-")
1153 
1154  self.writeEdgeFile(
1155  nodesToNodes,
1156  self.getParam("ig_cytoscape_input_file"))
1157 
1158  # write residue numbers for each node
1159  igResiduesFile = self.getParam("ig_cytoscape_residues_file")
1160  peptideChain = self.getParam("peptide_chain")
1161 
1162  subsetResidueFile = open(os.path.join(outputDir, igResiduesFile), 'w')
1163  subsetResidueFile.write("ResidueNumber\n")
1164  for nodeNumber in nodesToNames.keys():
1165  nodeName = nodesToNames[nodeNumber]
1166 
1167  [nodeChain, residueNumber,
1168  nodeAtom] = atomicDominoUtilities.getAtomInfoFromName(nodeName)
1169  if (nodeChain == peptideChain): # peptideAtom
1170  subsetResidueFile.write(
1171  "%s = %s\n" %
1172  (nodeNumber, residueNumber))
1173  else:
1174  # for now just write arbitrary number to designate as protein
1175  # atom
1176  subsetResidueFile.write("%s = 100\n" % nodeNumber)
1177 
1178  def writeCytoscapeMtInput(self):
1179  outputDir = self.getParam("output_directory")
1180  mTreeOutputFile = self.getParam("mtree_output_file")
1181  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1182  self.mt, os.path.join(outputDir, mTreeOutputFile), ">")
1183 
1184  self.writeEdgeFile(
1185  nodesToNodes,
1186  self.getParam("mtree_cytoscape_input_file"))
1187  self.writeNodeNameAttributes(
1188  nodesToNames, self.getParam("mtree_cytoscape_atom_name_file"), self.getParam(
1189  "mtree_cytoscape_atom_summary_file"),
1190  self.getParam("mtree_cytoscape_atom_chain_file"))
1191 
1192  def writeCytoscapeJtInput(self):
1193 
1194  outputDir = self.getParam("output_directory")
1195  jTreeOutputFile = self.getParam("jtree_output_file")
1196  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1197  self.jt, os.path.join(outputDir, jTreeOutputFile), "-")
1198 
1199  self.writeEdgeFile(
1200  nodesToNodes,
1201  self.getParam("jtree_cytoscape_input_file"))
1202 
1203  self.writeNodeNameAttributes(
1204  nodesToNames, self.getParam("jtree_cytoscape_atom_name_file"), self.getParam(
1205  "jtree_cytoscape_atom_summary_file"),
1206  self.getParam("jtree_cytoscape_atom_chain_file"))
1207 
1208  # write edge weight file -- weights are number of shared particles
1209  # across nodes
1210  jtreeCytoscapeEdgeFile = self.getParam("jtree_cytoscape_edge_file")
1211  edgeWeightFile = open(
1212  os.path.join(outputDir,
1213  jtreeCytoscapeEdgeFile),
1214  'w')
1215  edgeWeightFile.write("SubsetOverlap (class=Integer)\n")
1216  for firstNode in nodesToNodes.keys():
1217  nodeDict = nodesToNodes[firstNode]
1218  for secondNode in nodeDict.keys():
1219  firstNodeAtoms = nodesToNames[firstNode].split(" ")
1220  secondNodeAtoms = nodesToNames[secondNode].split(" ")
1221  intersection = [
1222  val for val in firstNodeAtoms if val in secondNodeAtoms]
1223  edgeWeightFile.write(
1224  "%s (pp) %s = %s\n" %
1225  (firstNode, secondNode, len(intersection)))
1226  edgeWeightFile.close()
1227 
1228  def getAtomTypeCounts(self, atomNames):
1229 
1230  atomNames = atomNames.lstrip('[')
1231  atomNames = atomNames.rstrip(']')
1232  atomNames = atomNames.rstrip()
1233 
1234  atomList = atomNames.split(" ") # atom names
1235  peptideAtomCount = 0
1236  proteinAtomCount = 0
1237  for atom in atomList:
1238  [chain, residue, atom] = atomicDominoUtilities.getAtomInfoFromName(
1239  atom)
1240  if (chain == self.getParam("peptide_chain")):
1241  peptideAtomCount += 1
1242  else:
1243  proteinAtomCount += 1
1244  return [peptideAtomCount, proteinAtomCount]
1245 
1246  def writeNodeNameAttributes(
1247  self,
1248  nodesToNames,
1249  atomNameFile,
1250  atomSummaryFile,
1251  atomChainFile):
1252  # write attribute file (atom names for each node)
1253  outputDir = self.getParam("output_directory")
1254  subsetAtomNameFile = open(os.path.join(outputDir, atomNameFile), 'w')
1255  subsetAtomSummaryFile = open(
1256  os.path.join(outputDir, atomSummaryFile), 'w')
1257  subsetAtomChainFile = open(os.path.join(outputDir, atomChainFile), 'w')
1258 
1259  subsetAtomNameFile.write("Atom names\n")
1260  subsetAtomSummaryFile.write("Atom Summary\n")
1261  subsetAtomChainFile.write("Atom chain\n")
1262  for node in nodesToNames.keys():
1263  atomNames = nodesToNames[node]
1264  subsetAtomNameFile.write("%s = %s\n" % (node, atomNames))
1265  [peptideAtomCount, proteinAtomCount] = self.getAtomTypeCounts(
1266  atomNames)
1267  # Number of protein and peptide atoms in each subset
1268  subsetAtomSummaryFile.write(
1269  "%s = %sp %sl\n" %
1270  (node, proteinAtomCount, peptideAtomCount))
1271 
1272  # whether each subset has protein, peptide, or a mix of atoms
1273  if (proteinAtomCount == 0):
1274  subsetAtomChainFile.write("%s = 1\n" % node)
1275  elif (peptideAtomCount == 0):
1276  subsetAtomChainFile.write("%s = 2\n" % node)
1277  else:
1278  subsetAtomChainFile.write("%s = 3\n" % node)
1279 
1280  subsetAtomChainFile.close()
1281  subsetAtomSummaryFile.close()
1282  subsetAtomNameFile.close()
1283 
1284 
1285  #
1286  # End Cytoscape Methods
1287  #
Chain get_chain(Hierarchy h)
InteractionGraph get_interaction_graph(ScoringFunctionAdaptor rs, const kernel::ParticlesTemp &pst)
void write_pdb(const Selection &mhd, base::TextOutput out, unsigned int model=1)
Grid3D< int, SparseGridStorage3D< int, UnboundedGridStorage3D > > SparseUnboundedIntGrid3D
kernel::ParticlesTemp get_order(const Subset &s, const SubsetFilterTables &sft)
void set_log_level(LogLevel l)
Set the current global log level.
SubsetGraph get_junction_tree(const InteractionGraph &ig)
Simple conjugate gradients optimizer.
Sample best solutions using Domino.
Definition: DominoSampler.h:32
algebra::BoundingBoxD< 3 > get_bounding_box(const Hierarchy &h)
Get a bounding box for the Hierarchy.
Subsets get_subsets(const SubsetGraph &g)
Gets all of the Subsets of a SubsetGraph.
Filter a configuration of the subset using the kernel::Model thresholds.
Represent a subset of the particles being optimized.
Definition: Subset.h:33
void load_frame(RMF::FileConstHandle file, unsigned int frame)
Definition: frames.h:27
Grid3D< double, DenseGridStorage3D< double > > DenseDoubleGrid3D
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
MergeTree get_balanced_merge_tree(const SubsetGraph &junction_tree)
display::Geometries get_subset_graph_geometry(const SubsetGraph &ig)
Select all non-alternative ATOM records.
Definition: pdb.h:63
Hierarchies get_by_type(Hierarchy mhd, GetByType t)
void load_merged_assignments(const Subset &first_subset, AssignmentContainer *first, const Subset &second_subset, AssignmentContainer *second, const SubsetFilterTablesTemp &filters, AssignmentContainer *ret)
Fill in assignments for an internal node.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
See IMP.core for more information.
See IMP.algebra for more information.
Store a configuration of a subset.
Definition: Assignment.h:32
VectorD< 3 > Vector3D
Definition: VectorD.h:395
Write a CGO file with the geometry.
Definition: PymolWriter.h:34
static XYZ decorate_particle(::IMP::kernel::Particle *p)
Definition: XYZ.h:49
display::Geometries get_interaction_graph_geometry(const InteractionGraph &ig)
void load_particle_states(const Subset &s, const Assignment &ss, const ParticleStatesTable *pst)
void read_pdb(base::TextInput input, int model, Hierarchy h)
Hierarchies get_leaves(const Selection &h)
See IMP.rmf for more information.
Definition: associations.h:20
See IMP.domino for more information.
Definition: analysis.h:15
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:72