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