IMP  2.3.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(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(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()
409  xyzCa = IMP.core.XYZ(ca)
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 
427  xyzD = IMP.core.XYZ(p)
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(IMP.core.XYZ(otherParticle).get_coordinates())
626  modelVector.append(IMP.core.XYZ(modelParticle).get_coordinates())
627  rmsd = IMP.atom.get_rmsd(otherVector, modelVector)
628  return rmsd
629 
630  def calculateNativeRmsd(self, flexibleAtoms):
631 
632  if (self.wroteNativeProtein == 0):
633  pdbName = self.getParam("native_pdb_input_file")
634  self.nativeModel = IMP.kernel.Model()
635  self.nativeProtein = IMP.atom.read_pdb(
636  pdbName,
637  self.nativeModel,
639  self.wroteNativeProtein = 1
640 
641  return self.calculateRmsd(self.nativeProtein, flexibleAtoms)
642 
643  def calculateTrajectoryRmsd(
644  self,
645  trajectoryFile,
646  trajectoryFrame,
647  flexibleAtoms):
648  pdbName = self.getParam("native_pdb_input_file")
649  otherModel = IMP.kernel.Model()
650  otherProtein = IMP.atom.read_pdb(
651  pdbName,
652  self.nativeModel,
654  outputDir = self.getParam("output_directory")
655  fullFile = os.path.join(outputDir, trajectoryFile)
656  print "open calculate traj rmf %s" % fullFile
657  rh = RMF.open_rmf_file(fullFile)
658  IMP.rmf.set_hierarchies(rh, [otherProtein])
659  if (trajectoryFrame == -1):
660  trajectoryFrame = IMP.rmf.get_number_of_frames(
661  rh, otherProtein) - 1
662  IMP.rmf.load_frame(rh, trajectoryFrame, otherProtein)
663  return self.calculateRmsd(otherProtein, flexibleAtoms)
664 
665  def createAllSubsetAssignments(self):
666 
669  self.model.get_root_restraint_set(),
670  self.dominoPst)
671 
672  leafNodeIndexList = self.getLeafNodeIndexList()
673 
674  for nodeIndex in leafNodeIndexList:
675  # get subset for this leaf
676  subset = self.mt.get_vertex_name(nodeIndex)
677  particleNameList = []
678  for particle in subset:
679  particleNameList.append(self.quickParticleName(particle))
680  print "creating initial assignments for leaf %s" % nodeIndex
681  # use particleInfo to create assignments and filter them
682  assignments = self.createUniqueLeafAssignments(
683  particleNameList, self.particleInfo)
684  filteredAssignments = self.filterAssignments(
685  assignments, subset, nodeIndex, rssft)
686 
687  # add assignemtns to container and listAssignmentTable
688  packedAssignmentContainer = IMP.domino.PackedAssignmentContainer()
689  for assignment in filteredAssignments:
690  packedAssignmentContainer.add_assignment(assignment)
691  lat.set_assignments(subset, packedAssignmentContainer)
692 
693  self.lat = lat
694  self.rssft = rssft
695 
696  def runDomino(self):
697  root = self.mt.get_vertices()[-1]
698  completeAc = self.loadAssignments(root)
699  self.completeAc = completeAc
700 
701  def loadAssignments(self, nodeIndex):
702 
703  children = self.mt.get_out_neighbors(nodeIndex)
704  subset = self.mt.get_vertex_name(nodeIndex)
705  heapCount = int(self.getParam("heap_count"))
707  heapCount, self.rssft.get_subset_filter(subset, []))
708  if len(children) == 0:
709  print "computing assignments for leaf %s" % nodeIndex
710 
711  self.sampler.load_vertex_assignments(nodeIndex, mine)
712  print "leaf node %s has %s leaf assignments" % (nodeIndex, mine.get_number_of_assignments())
713  else:
714  if (children[0] > children[1]):
715  children = [children[1], children[0]]
716  # recurse on the two children
717  firstAc = self.loadAssignments(children[0])
718  secondAc = self.loadAssignments(children[1])
719  self.logTimePoint(1)
720  self.sampler.load_vertex_assignments(
721  nodeIndex, firstAc, secondAc, mine)
722 
723  timeDifference = self.logTimePoint(0)
724  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(),
725  secondAc.get_number_of_assignments(), timeDifference)
726  self.totalAssignments += mine.get_number_of_assignments()
727  self.logMemory()
728  return mine
729 
730  def writeOutput(self, flexibleAtoms, startTime):
731  bestAssignment = -1
732  bestDominoScore = 100000
733  bestAssignment = 0
734  finalAssignments = self.completeAc.get_assignments()
735  for assignment in finalAssignments:
737  self.dominoPst.get_subset(),
738  assignment,
739  self.dominoPst)
740  score = self.model.evaluate(False)
741  if (score < bestDominoScore):
742  bestAssignment = assignment
743  bestDominoScore = round(score, 2)
744  print "best domino score is %s " % bestDominoScore
745  print "best md score is %s" % self.bestMdScore
746  print "best md rmsd is %s" % self.bestMdRmsd
747  print "best cg score is %s" % self.bestCgScore
748  print "best cg rmsd is %s" % self.bestCgRmsd
749  print "merge tree contained %s total assignments" % self.totalAssignments
750 
752  self.dominoPst.get_subset(),
753  bestAssignment,
754  self.dominoPst)
755  dominoVsMdRmsd = round(
756  self.calculateTrajectoryRmsd(
757  self.getParam(
758  "md_trajectory_output_file"),
759  self.bestMdScoreFrame,
760  flexibleAtoms),
761  2)
762  cg = IMP.core.ConjugateGradients(self.model)
763  cg.optimize(100)
765  self.protein,
766  os.path.join(self.getParam("output_directory"),
767  self.getParam("minimum_domino_score_pdb")))
768 
769  dominoVsCgRmsd = round(
770  self.calculateTrajectoryRmsd(
771  self.bestCgScoreFile,
772  -1,
773  flexibleAtoms),
774  2)
775  dominoMinimizedScore = round(self.model.evaluate(False), 2)
776  dominoRmsd = round(self.calculateNativeRmsd(flexibleAtoms), 2)
777 
778  runTime = round(time.time() - startTime, 2)
779  print "final domino score (after cg): %s" % dominoMinimizedScore
780  print "final domino rmsd: %s" % dominoRmsd
781  print "best domino rmsd with best md score: %s" % dominoVsMdRmsd
782  print "domino rmsd with best cg score: %s" % dominoVsCgRmsd
783  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)
784 
785  #
786  # Parallel methods
787  #
788  def createSubsetFromParticles(self, particleNames):
789  particleNameList = particleNames.split(" ")
790  particleList = []
791  for particleName in particleNameList:
792  particle = self.namesToParticles[particleName]
793  particleList.append(particle)
794 
795  subset = IMP.domino.Subset(particleList)
796  return [subset, particleList]
797 
798  def createHdf5AssignmentContainer(self, index, particleNames, read):
799  root = self.getAssignmentContainerRoot(index, read)
800  print "got root for index %s" % index
801  dataset = 0
802  if (read == 1):
803  dataset = root.get_child_index_data_set_2d(str(index))
804  else:
805  dataset = root.add_child_index_data_set_2d(str(index))
806  print "added child index dataset"
807  # firstDataset.set_size([
808  [subset, particleOrder] = self.createSubsetFromParticles(particleNames)
809  print "created subset for index %s" % index
810  hdf5Ac = IMP.domino.create_assignments_container(
811  dataset, subset, particleOrder)
812  print "returning from create"
813  order = IMP.domino.get_order(subset, particleOrder)
814  for nextInt in order:
815  print "next int is %s" % nextInt
816  return [subset, hdf5Ac]
817 
818  def loadAssignmentsParallel(
819  self,
820  nodeIndex,
821  particleInfo,
822  mtIndexToNodeInfo,
823  mtIndexToSubsetOrder,
824  mtIndexToParticles):
825  IMP.base.set_log_level(IMP.WARNING)
826 
827  if (("firstChild" in mtIndexToNodeInfo[nodeIndex]) == 0):
828  print "writing file for leaf index %s" % nodeIndex
829  return (
830  self.createAssignmentsParallel(
831  particleInfo,
832  nodeIndex,
833  mtIndexToParticles)
834  )
835 
836  else:
837  beginTime = self.logTimePoint(1)
838  firstChildIndex = mtIndexToNodeInfo[nodeIndex]["firstChild"]
839  [firstSubset, firstAc] = self.createHdf5AssignmentContainer(
840  firstChildIndex, mtIndexToSubsetOrder[firstChildIndex], 1)
841  firstAcCreateTime = self.logTimePoint(0)
842 
843  secondChildIndex = mtIndexToNodeInfo[nodeIndex]["secondChild"]
844  [secondSubset, secondAc] = self.createHdf5AssignmentContainer(
845  secondChildIndex, mtIndexToSubsetOrder[secondChildIndex], 1)
846  secondAcCreateTime = self.logTimePoint(0)
847 
848  print "getting assignments for nodeIndex %s first child %s second child %s" % (nodeIndex, firstChildIndex, secondChildIndex)
849  firstChildParticles = mtIndexToSubsetOrder[firstChildIndex]
850  secondChildParticles = mtIndexToSubsetOrder[secondChildIndex]
851  myParticles = mtIndexToParticles[nodeIndex]
852 
853  print "first child particles %s\nsecond child particles %s\nmy particles; %s\n" % (firstChildParticles, secondChildParticles, myParticles)
854  for p in firstSubset:
855  print "next particle in first subset: %s" % self.quickParticleName(p)
856 
857  for p in secondSubset:
858  print "next particle in second subset: %s" % self.quickParticleName(p)
859 
860  for assignment in firstAc.get_assignments():
861  print "next assignment for first child %s: %s" % (firstChildIndex, assignment)
862 
863  for assignment in secondAc.get_assignments():
864  print "next assignment for second child %s: %s" % (secondChildIndex, assignment)
865 
866  [mySubset, myAc] = self.createHdf5AssignmentContainer(
867  nodeIndex, mtIndexToParticles[nodeIndex], 0)
868  print "done creating hdf5"
869  prepTime = self.logTimePoint(0)
871  self.model.get_root_restraint_set(),
872  self.dominoPst)
873  rssf = rssft.get_subset_filter(mySubset, [])
874 
875  #heapAc = IMP.domino.HeapAssignmentContainer(1000, rssf)
876 
878  firstSubset,
879  firstAc,
880  secondSubset,
881  secondAc,
882  [rssft],
883  myAc)
884 
885  heapTime = self.logTimePoint(0)
886  # myAc.add_assignments(heapAc.get_assignments())
887  addTime = self.logTimePoint(0)
888  for assignment in myAc.get_assignments():
889  print "loadAssignmentsParallel next assignment for %s: %s" % (nodeIndex, assignment)
890  doneTime = self.logTimePoint(0)
891  firstChildCount = firstAc.get_number_of_assignments()
892  secondChildCount = secondAc.get_number_of_assignments()
893  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)
894  subsetOrder = self.getSubsetOrderList(mySubset)
895  return subsetOrder
896 
897  def createAssignmentsParallel(
898  self,
899  particleInfo,
900  nodeIndex,
901  mtIndexToParticles):
902 
903  subsetName = mtIndexToParticles[nodeIndex]
904  print "starting assignments parallel leaf index %s subset name %s" % (nodeIndex, subsetName)
905  [subset, particleList] = self.createSubsetFromParticles(subsetName)
906 
907  particleNameList = subsetName.split(" ")
908 
909  # create assignment by reading states
910  finalAssignments = self.createUniqueLeafAssignments(
911  particleNameList, particleInfo)
912 
913  #lat = IMP.domino.ListAssignmentsTable()
914  #lat.set_assignments(subset, finalAssignmentContainer)
915  # self.sampler.set_assignments_table(lat)
916  #hdf5AssignmentContainer = IMP.domino.HDF5AssignmentContainer(dataset, subset, self.dominoPst.get_particles(), subsetName)
918  self.model.get_root_restraint_set(),
919  self.dominoPst)
920  filteredAssignments = self.filterAssignments(
921  finalAssignments, subset, nodeIndex, rssft)
922  root = self.getAssignmentContainerRoot(nodeIndex, 0)
923  dataset = root.add_child_index_data_set_2d(str(nodeIndex))
924  dataset.set_size([0, len(subset)])
925 
926  hdf5AssignmentContainer = IMP.domino.create_assignments_container(
927  dataset, subset, particleOrder)
928  for assignment in filteredAssignments:
929  hdf5AssignmentContainer.add_assignment(assignment)
930  for assignment in hdf5AssignmentContainer.get_assignments():
931  print "hdf5 assignment container node %s next assignmet %s" % (nodeIndex, assignment)
932 
933  #self.checkAssignments(subset, nodeIndex, particleOrder)
934  subsetOrder = self.getSubsetOrderList(subset)
935  print "leaf node returning with order %s" % subsetOrder
936  return subsetOrder
937 
938  # Write pymol session files for the interactions across atoms and all
939  # subsets
940  def writePymolData(self):
941 
942  outputDir = self.getParam("output_directory")
944  pymolInteractions = self.getParam("pymol_interactions_file")
945  w = IMP.display.PymolWriter(os.path.join(outputDir, pymolInteractions))
946 
947  for gg in geometry:
948  w.add_geometry(gg)
949 
950  pymolSubsets = self.getParam("pymol_subsets_file")
951  geometry = IMP.domino.get_subset_graph_geometry(self.jt)
952  w = IMP.display.PymolWriter(os.path.join(outputDir, pymolSubsets))
953  for gg in geometry:
954  w.add_geometry(gg)
955 
956  # Clean the default name of the vertex (in brackets and with each atom
957  # contained in quotes) and return a string where [] and " are removed
958  def cleanVertexName(self, vertexName):
959 
960  # not sure if any vertices still have a Subset prefix but keeping
961  # anyway
962  nodeRe = re.compile('Subset\(\[(.*?)\s*\]')
963  secondNodeRe = re.compile('\[(.*?)\s*\]') # atom name
964  node = nodeRe.search(str(vertexName))
965  secondNode = secondNodeRe.search(str(vertexName))
966  vertexNameFinal = ""
967  foundName = 0
968  if node:
969  foundName = node.group(1)
970  if secondNode:
971  foundName = secondNode.group(1)
972  vertexNameFinal = foundName.replace('"', '')
973  return vertexNameFinal
974 
975  def getSubsetOrderList(self, subset):
976  subsetOrderList = []
977  for particle in subset:
978  name = self.quickParticleName(particle)
979  subsetOrderList.append(name)
980  subsetOrder = " ".join(subsetOrderList)
981  return subsetOrder
982 
983  def checkAssignments(self, subset, nodeIndex, particleOrder):
984  print "reading back in assignments for leaf index %s" % nodeIndex
985  root = self.getAssignmentContainerRoot(nodeIndex, 1)
986  dataset = root.get_child_index_data_set_2d(str(nodeIndex))
987  hdf5 = IMP.domino.create_assignments_container(
988  dataset, subset, particleOrder)
989  for assignment in hdf5.get_assignments():
990  print "leaf index %s read back in %s" % (nodeIndex, assignment)
991 
992  def createSamplerLite(self):
993  s = IMP.domino.DominoSampler(self.model, self.dominoPst)
994 
995  if (self.getParam("cross_subset_filtering") == 1):
996  s.set_use_cross_subset_filtering(1)
997 
998  self.sampler = s
999 
1000  def createSiblingMap(self, parentIndex):
1001 
1002  children = self.mt.get_out_neighbors(parentIndex)
1003  if (len(children) > 0):
1004  firstChild = children[0]
1005  secondChild = children[1]
1006  self.parentSiblingMap[firstChild] = {}
1007  self.parentSiblingMap[firstChild]["sibling"] = secondChild
1008  self.parentSiblingMap[firstChild]["parent"] = parentIndex
1009 
1010  self.parentSiblingMap[secondChild] = {}
1011  self.parentSiblingMap[secondChild]["sibling"] = firstChild
1012  self.parentSiblingMap[secondChild]["parent"] = parentIndex
1013  print "created map for parent %s first child %s second child %s" % (parentIndex, firstChild, secondChild)
1014 
1015  self.parentSiblingMap[parentIndex]["firstChild"] = firstChild
1016  self.parentSiblingMap[parentIndex]["secondChild"] = secondChild
1017 
1018  self.createSiblingMap(firstChild)
1019  self.createSiblingMap(secondChild)
1020 
1021  def getMtIndexToNodeInfo(self):
1022  return self.parentSiblingMap
1023 
1024  def getLeafParentNodeIndexList(self):
1025  leafParentNodeIndexList = {}
1026  leafNodeIndexList = self.getLeafNodeIndexList()
1027  for leafIndex in leafNodeIndexList:
1028  parent = self.parentSiblingMap[leafIndex]["parent"]
1029  leafParentNodeIndexList[parent] = 1
1030  return leafParentNodeIndexList
1031 
1032  def getMtIndexToNameList(selt):
1033  mtIndexToNames = {}
1034  for index in self.mt.get_vertices():
1035  name = self.mt.get_vertex_name(index)
1036  mtIndexToNames[index] = name
1037  return mtIndexToNames
1038 
1039  def getAssignmentContainerRoot(self, subsetIndex, read):
1040  outputDir = self.getParam("output_directory")
1041  filePrefix = self.getParam("subset_assignment_output_file")
1042  assignmentFileName = os.path.join(
1043  outputDir, "%s%s" %
1044  (filePrefix, subsetIndex))
1045  print "creating hdf5 file with name %s" % assignmentFileName
1046  root = 0
1047  if (read == 1):
1048  root = RMF.open_hdf5_file(assignmentFileName)
1049  else:
1050  root = RMF.create_hdf5_file(assignmentFileName)
1051  return root
1052 
1053  #
1054  # Begin Cytoscape Methods
1055  #
1056 
1057  def writeVisualization(self):
1058 
1059  self.writeCytoscapeIgInput()
1060  self.writeCytoscapeJtInput()
1061  self.writeCytoscapeMtInput()
1062  self.writeCytoscapeScripts()
1063  self.writePymolData()
1064 
1065  def writeCytoscapeScripts(self):
1066  outputDir = self.getParam("output_directory")
1067  mTreeCytoscapeInput = self.getParam("mtree_cytoscape_input_file")
1068  mTreeCytoscapeAssignments = self.getParam(
1069  "mtree_cytoscape_assignment_file")
1070  mTreeCytoscapeAtomChains = self.getParam(
1071  "mtree_cytoscape_atom_chain_file")
1072  mTreeCytoscapeAtomSummary = self.getParam(
1073  "mtree_cytoscape_atom_summary_file")
1074 
1075  mTreeCytoscapeScript = self.getParam("mtree_cytoscape_script")
1076  mTreeFh = open(os.path.join(outputDir, mTreeCytoscapeScript), 'w')
1077  mTreeFh.write(
1078  "network import file=%s\n" %
1079  os.path.join(outputDir, mTreeCytoscapeInput))
1080  mTreeFh.write(
1081  "node import attributes file=\"%s\"\n" %
1082  os.path.join(outputDir, mTreeCytoscapeAssignments))
1083  mTreeFh.write(
1084  "node import attributes file=\"%s\"\n" %
1085  os.path.join(outputDir, mTreeCytoscapeAtomSummary))
1086  mTreeFh.write(
1087  "node import attributes file=\"%s\"\n" %
1088  os.path.join(outputDir, mTreeCytoscapeAtomChains))
1089  mTreeFh.write("layout jgraph-tree\n")
1090 
1091  def getGraphStructure(self, graph, fileName, separator):
1092  # write junction tree to file
1093  graphLogWrite = open(fileName, 'w')
1094  graph.show(graphLogWrite)
1095  graphLogWrite.close()
1096 
1097  # read file
1098  graphLogRead = open(fileName, 'r')
1099 
1100  nodeRe = re.compile('^(\d+)\[label\=\"\[*(.*?)\s*\]*\"') # atom name
1101  separatorEscape = "\\" + separator
1102  edgeString = "^(\d+)\-%s(\d+)" % separatorEscape
1103  edgeRe = re.compile(edgeString)
1104 
1105  nodesToNodes = {}
1106  # keys: source node, value: target node (track edges)
1107  nodesToNames = {} # keys: node number, value; string parsed from file
1108 
1109  for line in graphLogRead:
1110 
1111  # search for nodes
1112  node = nodeRe.search(line)
1113  if node:
1114  nodeNumber = node.group(1)
1115  atomString = node.group(2)
1116  nodesToNames[nodeNumber] = atomString
1117  continue
1118 
1119  # search for edges
1120  edge = edgeRe.search(line)
1121  if edge:
1122  firstNode = edge.group(1)
1123  secondNode = edge.group(2)
1124  firstNodeDict = {}
1125  if (firstNode in nodesToNodes):
1126  firstNodeDict = nodesToNodes[firstNode]
1127  firstNodeDict[secondNode] = 1
1128  nodesToNodes[firstNode] = firstNodeDict
1129 
1130  return [nodesToNames, nodesToNodes]
1131 
1132  def writeEdgeFile(self, nodesToNodes, edgeFileName):
1133  # write edge file
1134  outputDir = self.getParam("output_directory")
1135  graphInputFile = open(os.path.join(outputDir, edgeFileName), 'w')
1136  for firstNode in nodesToNodes.keys():
1137  nodeDict = nodesToNodes[firstNode]
1138  for secondNode in nodeDict.keys():
1139  graphInputFile.write("%s ttt %s\n" % (firstNode, secondNode))
1140  graphInputFile.close()
1141 
1142  def writeCytoscapeIgInput(self):
1143  outputDir = self.getParam("output_directory")
1144  igOutputFile = self.getParam("ig_output_file")
1145  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1146  self.ig, os.path.join(outputDir, igOutputFile), "-")
1147 
1148  self.writeEdgeFile(
1149  nodesToNodes,
1150  self.getParam("ig_cytoscape_input_file"))
1151 
1152  # write residue numbers for each node
1153  igResiduesFile = self.getParam("ig_cytoscape_residues_file")
1154  peptideChain = self.getParam("peptide_chain")
1155 
1156  subsetResidueFile = open(os.path.join(outputDir, igResiduesFile), 'w')
1157  subsetResidueFile.write("ResidueNumber\n")
1158  for nodeNumber in nodesToNames.keys():
1159  nodeName = nodesToNames[nodeNumber]
1160 
1161  [nodeChain, residueNumber,
1162  nodeAtom] = atomicDominoUtilities.getAtomInfoFromName(nodeName)
1163  if (nodeChain == peptideChain): # peptideAtom
1164  subsetResidueFile.write(
1165  "%s = %s\n" %
1166  (nodeNumber, residueNumber))
1167  else:
1168  # for now just write arbitrary number to designate as protein
1169  # atom
1170  subsetResidueFile.write("%s = 100\n" % nodeNumber)
1171 
1172  def writeCytoscapeMtInput(self):
1173  outputDir = self.getParam("output_directory")
1174  mTreeOutputFile = self.getParam("mtree_output_file")
1175  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1176  self.mt, os.path.join(outputDir, mTreeOutputFile), ">")
1177 
1178  self.writeEdgeFile(
1179  nodesToNodes,
1180  self.getParam("mtree_cytoscape_input_file"))
1181  self.writeNodeNameAttributes(
1182  nodesToNames, self.getParam("mtree_cytoscape_atom_name_file"), self.getParam(
1183  "mtree_cytoscape_atom_summary_file"),
1184  self.getParam("mtree_cytoscape_atom_chain_file"))
1185 
1186  def writeCytoscapeJtInput(self):
1187 
1188  outputDir = self.getParam("output_directory")
1189  jTreeOutputFile = self.getParam("jtree_output_file")
1190  [nodesToNames, nodesToNodes] = self.getGraphStructure(
1191  self.jt, os.path.join(outputDir, jTreeOutputFile), "-")
1192 
1193  self.writeEdgeFile(
1194  nodesToNodes,
1195  self.getParam("jtree_cytoscape_input_file"))
1196 
1197  self.writeNodeNameAttributes(
1198  nodesToNames, self.getParam("jtree_cytoscape_atom_name_file"), self.getParam(
1199  "jtree_cytoscape_atom_summary_file"),
1200  self.getParam("jtree_cytoscape_atom_chain_file"))
1201 
1202  # write edge weight file -- weights are number of shared particles
1203  # across nodes
1204  jtreeCytoscapeEdgeFile = self.getParam("jtree_cytoscape_edge_file")
1205  edgeWeightFile = open(
1206  os.path.join(outputDir,
1207  jtreeCytoscapeEdgeFile),
1208  'w')
1209  edgeWeightFile.write("SubsetOverlap (class=Integer)\n")
1210  for firstNode in nodesToNodes.keys():
1211  nodeDict = nodesToNodes[firstNode]
1212  for secondNode in nodeDict.keys():
1213  firstNodeAtoms = nodesToNames[firstNode].split(" ")
1214  secondNodeAtoms = nodesToNames[secondNode].split(" ")
1215  intersection = [
1216  val for val in firstNodeAtoms if val in secondNodeAtoms]
1217  edgeWeightFile.write(
1218  "%s (pp) %s = %s\n" %
1219  (firstNode, secondNode, len(intersection)))
1220  edgeWeightFile.close()
1221 
1222  def getAtomTypeCounts(self, atomNames):
1223 
1224  atomNames = atomNames.lstrip('[')
1225  atomNames = atomNames.rstrip(']')
1226  atomNames = atomNames.rstrip()
1227 
1228  atomList = atomNames.split(" ") # atom names
1229  peptideAtomCount = 0
1230  proteinAtomCount = 0
1231  for atom in atomList:
1232  [chain, residue, atom] = atomicDominoUtilities.getAtomInfoFromName(
1233  atom)
1234  if (chain == self.getParam("peptide_chain")):
1235  peptideAtomCount += 1
1236  else:
1237  proteinAtomCount += 1
1238  return [peptideAtomCount, proteinAtomCount]
1239 
1240  def writeNodeNameAttributes(
1241  self,
1242  nodesToNames,
1243  atomNameFile,
1244  atomSummaryFile,
1245  atomChainFile):
1246  # write attribute file (atom names for each node)
1247  outputDir = self.getParam("output_directory")
1248  subsetAtomNameFile = open(os.path.join(outputDir, atomNameFile), 'w')
1249  subsetAtomSummaryFile = open(
1250  os.path.join(outputDir, atomSummaryFile), 'w')
1251  subsetAtomChainFile = open(os.path.join(outputDir, atomChainFile), 'w')
1252 
1253  subsetAtomNameFile.write("Atom names\n")
1254  subsetAtomSummaryFile.write("Atom Summary\n")
1255  subsetAtomChainFile.write("Atom chain\n")
1256  for node in nodesToNames.keys():
1257  atomNames = nodesToNames[node]
1258  subsetAtomNameFile.write("%s = %s\n" % (node, atomNames))
1259  [peptideAtomCount, proteinAtomCount] = self.getAtomTypeCounts(
1260  atomNames)
1261  # Number of protein and peptide atoms in each subset
1262  subsetAtomSummaryFile.write(
1263  "%s = %sp %sl\n" %
1264  (node, proteinAtomCount, peptideAtomCount))
1265 
1266  # whether each subset has protein, peptide, or a mix of atoms
1267  if (proteinAtomCount == 0):
1268  subsetAtomChainFile.write("%s = 1\n" % node)
1269  elif (peptideAtomCount == 0):
1270  subsetAtomChainFile.write("%s = 2\n" % node)
1271  else:
1272  subsetAtomChainFile.write("%s = 3\n" % node)
1273 
1274  subsetAtomChainFile.close()
1275  subsetAtomSummaryFile.close()
1276  subsetAtomNameFile.close()
1277 
1278 
1279  #
1280  # End Cytoscape Methods
1281  #
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
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
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
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)
Support for the RMF file format for storing hierarchical molecular data and markup.
Divide-and-conquer inferential optimization in discrete space.
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:73