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