16 import atomicDominoUtilities
21 def __init__(self, model, protein, parameterFileName):
23 self.protein = protein
24 self.namesToParticles = atomicDominoUtilities.makeNamesToParticles(
26 self.readParameters(parameterFileName)
27 self.wroteNativeProtein = 0
30 def readParameters(self, parameterFileName):
31 self.parameters = atomicDominoUtilities.readParameters(
35 currentMem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
36 currentMem *= 10 ** -6
37 print "log memory: %s" % currentMem
38 if (currentMem > self.maxMem):
39 self.maxMem = currentMem
42 def getParam(self, paramName):
43 paramValue = self.parameters[paramName]
46 def loadDominoHelpers(self):
47 self.subsetStateScores = {}
48 self.subsetRestraintScores = {}
49 self.subsetStateFailures = {}
50 self.nodesToAssignmentCount = {}
52 self.totalAssignments = 0
53 outputDir = self.getParam(
"output_directory")
54 mtreeCytoscapeAssignmentFile = self.getParam(
55 "mtree_cytoscape_assignment_file")
56 self.mtreeCytoscapeAssignmentFh = open(
57 os.path.join(outputDir, mtreeCytoscapeAssignmentFile),
'w')
58 self.mtreeCytoscapeAssignmentFh.write(
"Assignments\n")
60 self.mtNamesToIndices = {}
62 self.subsetNamesToAssignmentFiles = {}
68 def logTimePoint(self, reset):
71 timeDifference = newTime - self.timePoint
72 timeDifference = round(timeDifference, 0)
73 self.timePoint = newTime
77 self.timePoint = newTime
80 def createParticleStatesTable(self):
82 for particleName
in self.particleInfo.keys():
85 dataArray = self.particleInfo[particleName]
86 for i
in range(len(dataArray)):
87 [state, center] = dataArray[i]
88 statesToCenters[state] = center
91 for stateIndex
in sorted(statesToCenters.keys()):
93 statesList.append(vector3d)
97 self.dominoPst.set_particle_states(
98 self.namesToParticles[particleName],
101 def quickParticleName(self, particle):
102 return atomicDominoUtilities.quickParticleName(particle)
104 def filterAssignments(self, assignments, subset, nodeIndex, rssft):
110 for r
in IMP.get_restraints([self.model.get_root_restraint_set()]):
111 restraintList.append(r)
112 dg = IMP.get_dependency_graph(restraintList)
119 sFilter = rssft.get_subset_filter(subset, filteredSubsets)
121 filteredAssignments = []
122 for assignment
in assignments:
124 if (sFilter
is None or sFilter.get_is_ok(assignment)):
127 filteredAssignments.append(assignment)
130 fraction = (passedCounter * 1.0) / (stateCounter * 1.0)
131 print "%s states passed out of %s total for this subset (fraction %s)" % (passedCounter, stateCounter, fraction)
132 if (passedCounter == 0):
133 print "subset %s had 0 assignments (out of %s) pass. Exiting..." % (subset, stateCounter)
136 return filteredAssignments
140 def quickSubsetName(self, subset):
141 cleanName = self.cleanVertexName(subset)
142 atomNameList = cleanName.split(
" ")
143 sortedList = sorted(atomNameList)
145 name =
" ".join(sortedList)
148 def getNodeIndexList(self):
149 return self.mt.get_vertices()
151 def getLeafNodeIndexList(self):
152 leafNodeIndexList = []
153 for subset
in self.subsets:
154 index = self.getMtIndexForSubset(self.quickSubsetName(subset))
155 leafNodeIndexList.append(index)
156 return leafNodeIndexList
159 def createSampler(self):
162 s.set_merge_tree(self.mt)
164 filterTables.append(self.rssft)
165 s.set_subset_filter_tables(filterTables)
167 s.set_assignments_table(self.lat)
169 if (self.getParam(
"cross_subset_filtering") == 1):
170 s.set_use_cross_subset_filtering(1)
176 def createSubsets(self):
177 self.initializeParticleStatesTable()
179 [self.model.get_root_restraint_set()],
181 print "interaction graph:"
184 print "junction tree:"
186 self.subsetNamesToSubsets = {}
190 for index
in mt.get_vertices():
191 subset = mt.get_vertex_name(index)
192 subsetName = self.cleanVertexName(subset)
193 self.mtNamesToIndices[subsetName] = index
196 for subset
in subsets:
197 print "created subset %s" % subset
198 subsetName = self.quickSubsetName(subset)
199 self.subsetNamesToSubsets[subsetName] = subset
207 self.subsets = subsets
209 self.parentSiblingMap = {}
210 self.parentSiblingMap[self.mt.get_vertices()[-1]] = {}
211 self.createSiblingMap(self.mt.get_vertices()[-1])
213 def getMtIndexForSubset(self, subsetName):
214 for index
in self.mt.get_vertices():
215 mtNode = self.mt.get_vertex_name(index)
216 mtName = self.cleanVertexName(mtNode)
217 mtName = mtName.rstrip()
218 mtName = mtName.lstrip()
220 mtNameList = mtName.split(
" ")
221 subsetNameList = subsetName.split(
" ")
222 mtNameListSorted = sorted(mtNameList)
223 subsetNameListSorted = sorted(subsetNameList)
224 if (
" ".join(subsetNameListSorted) ==
" ".join(mtNameListSorted)):
227 print "did not find merge tree index for subset name %s" % subsetName
230 def getDominoParticleNames(self):
231 particles = self.dominoPst.get_particles()
232 particleNameList = []
233 for particle
in particles:
234 pName = self.quickParticleName(particle)
235 particleNameList.append(pName)
237 return particleNameList
239 def getSubsetNameList(self):
241 for subset
in self.subsets:
242 name = self.quickSubsetName(subset)
243 subsetNameList.append(name)
244 return subsetNameList
246 def getMtIndexToParticles(self):
248 mtIndexToParticles = {}
249 allVertices = self.mt.get_vertices()
250 for nodeIndex
in allVertices:
251 subset = self.mt.get_vertex_name(nodeIndex)
252 particleList = self.quickSubsetName(subset)
253 mtIndexToParticles[nodeIndex] = particleList
254 print "createtMtIndexToParticles: index %s is particleList %s " % (nodeIndex, particleList)
256 return mtIndexToParticles
258 def readTrajectoryFile(
271 scoreOutputFh = open(scoreOutputFile,
'w')
276 print "execption in loading frame %s" % i
278 score = self.model.evaluate(
False)
286 rmsd = self.calculateNativeRmsd(flexibleAtoms)
287 scoreOutputFh.write(
"%s\t%s\t%s\n" % (i, score, rmsd))
288 if (score < bestScore):
292 if (rmsd < bestRmsd):
295 if (skipDomino == 0):
296 for atomName
in atomList:
298 particle = self.namesToParticles[atomName]
301 gridIndex = self.snapToGrid(particle)
302 center = self.grid.get_center(gridIndex)
304 for coordinate
in center:
305 pythonCenter.append(coordinate)
309 if ((gridIndex
in self.particleStatesSeen[atomName]) == 0):
314 self.particleStatesSeen[atomName].keys())
315 self.particleStatesSeen[
317 gridIndex] = currentSize
318 state = self.particleStatesSeen[atomName][gridIndex]
319 self.particleInfo[atomName].append([state, pythonCenter])
321 print "didn't add domino states due to skip parameters"
322 return [bestScore, bestScoreFrame, bestRmsd, bestRmsdFrame]
326 def getParticleGridIndex(self, leaf):
328 coordinates = xyzDecorator.get_coordinates()
330 extendedIndex = self.grid.get_extended_index(coordinates)
331 if (self.grid.get_has_index(extendedIndex) == 0):
333 if (self.getParam(
"grid_type") ==
"sparse"):
334 self.grid.add_voxel(extendedIndex, 1)
336 self.grid.add_voxel(extendedIndex)
338 index = self.grid.get_index(extendedIndex)
343 def snapToGrid(self, particle):
345 index = self.getParticleGridIndex(particle)
346 center = self.grid.get_center(index)
348 xyzDecorator.set_coordinates(center)
353 def discretizeNativeProtein(self):
355 outputDir = self.getParam(
"output_directory")
356 nativeSnappedFile = self.getParam(
"native_protein_snapped_output_file")
360 self.snapToGrid(leaf)
369 def getPeptideCa(self):
373 peptideIndicesToResidues = {}
374 for residueH
in residues:
376 chainId = chain.get_id()
377 residue = residueH.get_as_residue()
378 if (chainId == self.getParam(
"peptide_chain")):
379 peptideIndicesToResidues[residue.get_index()] = residue
383 minPeptide = min(sorted(peptideIndicesToResidues.keys()))
384 maxPeptide = max(sorted(peptideIndicesToResidues.keys()))
385 centerPeptide = round(((maxPeptide - minPeptide) / 2 + minPeptide), 0)
389 centerName = atomicDominoUtilities.makeParticleName(
390 self.getParam(
"peptide_chain"),
393 centerParticle = self.namesToParticles[centerName]
394 return centerParticle
398 def createGrid(self):
402 gridSpacing = float(self.getParam(
"grid_spacing"))
403 bufferSpace = float(self.getParam(
"grid_buffer_space"))
405 protBb += bufferSpace
407 if (self.getParam(
"grid_type") ==
"sparse"):
408 ca = self.getPeptideCa()
411 gridSpacing, xyzCa.get_coordinates())
419 def initializeParticleStatesTable(self):
422 restrainedParticles = atomicDominoUtilities.getRestrainedParticles(
423 self.protein, self.model, self.namesToParticles)
425 for p
in restrainedParticles:
430 dominoPst.set_particle_states(p, xyzStates)
432 self.dominoPst = dominoPst
434 def createUniqueLeafAssignments(self, particleNameList, particleInfo):
436 size = len(particleInfo[particleNameList[0]])
439 for i
in range(size):
441 for particleName
in particleNameList:
442 dataArray = particleInfo[particleName]
443 [state, center] = dataArray[i]
444 nextAssignment.append(state)
445 allAssignments.append(nextAssignment)
448 uniqueAssignments = {}
449 for assignment
in allAssignments:
452 for index
in assignment:
453 stateString = stateString + str(index) +
"_"
454 uniqueAssignments[stateString] = assignment
455 finalAssignments = []
458 for stateString
in uniqueAssignments.keys():
459 assignmentList = uniqueAssignments[stateString]
461 finalAssignments.append(assignment)
462 return finalAssignments
466 def readMdTrajectory(self, atomList, flexibleAtoms):
469 outputDir = self.getParam(
"output_directory")
470 trajectoryFile = self.getParam(
"md_trajectory_output_file")
471 fullFile = os.path.join(outputDir, trajectoryFile)
472 rh = RMF.open_rmf_file(fullFile)
473 IMP.rmf.set_hierarchies(rh, [self.protein])
474 framesToRead = atomicDominoUtilities.getMdIntervalFrames(
475 rh, int(self.getParam(
"md_interval")), self.protein)
476 print "preparing to read md frames %s" % framesToRead
482 particleStatesSeen = {}
483 for atomName
in atomList:
484 particle = self.namesToParticles[atomName]
485 particleInfo[atomName] = []
486 particleStatesSeen[atomName] = {}
488 self.particleStatesSeen = particleStatesSeen
489 self.particleInfo = particleInfo
492 mdScoreOutputFile = os.path.join(
494 self.getParam(
"md_score_output_file"))
498 bestRmsdFrame] = self.readTrajectoryFile(
508 self.singlePdbResults(
511 self.getParam(
"best_md_score_output_file"))
514 self.singlePdbResults(
517 self.getParam(
"best_md_rmsd_output_file"))
518 self.singlePdbResults(
521 self.getParam(
"final_md_frame_output_file"))
523 self.bestMdScore = round(bestMdScore, 2)
524 self.bestMdRmsd = round(bestRmsd, 2)
525 self.bestMdScoreFrame = bestScoreFrame
526 self.bestMdRmsdFrame = bestRmsdFrame
528 def readCgTrajectories(self, atomList, flexibleAtoms):
530 cgFileName = self.getParam(
"cg_output_file")
531 bestCgScore = 10000000
533 bestCgRmsd = 10000000
536 outputDir = self.getParam(
"output_directory")
537 trajectoryFile = self.getParam(
"md_trajectory_output_file")
538 fullFile = os.path.join(outputDir, trajectoryFile)
539 rh = RMF.open_rmf_file(fullFile)
540 IMP.rmf.set_hierarchies(rh, [self.protein])
541 framesToRead = atomicDominoUtilities.getMdIntervalFrames(
542 rh, int(self.getParam(
"cg_interval")), self.protein)
544 skipCgDomino = int(self.getParam(
"skip_cg_domino"))
546 if (len(framesToRead) > 0):
547 for cgNumber
in framesToRead:
549 outputDir = self.getParam(
"output_directory")
550 fullCgFileName = os.path.join(
552 (cgFileName, cgNumber))
553 rh = RMF.open_rmf_file(fullCgFileName)
554 IMP.rmf.set_hierarchies(rh, [self.protein])
557 frameCount = IMP.rmf.get_number_of_frames(rh, self.protein)
560 if (frameCount > 20):
561 startFrameCount = frameCount - 20
563 for i
in range(startFrameCount, frameCount):
567 cgScoreOutputFile = os.path.join(
569 (self.getParam(
"cg_score_output_file"), cgNumber))
570 [cgScore, cgScoreFrame, cgRmsd, cgRmsdFrame] = self.readTrajectoryFile(
571 atomList, rh, cgFrames, cgScoreOutputFile, skipCgDomino, flexibleAtoms)
572 print "cg number %s rmsd %s score %s" % (cgNumber, cgRmsd, cgScore)
574 if (cgScore < bestCgScore):
575 bestCgScore = cgScore
576 bestCgScoreFile = fullCgFileName
577 if (cgRmsd < bestCgRmsd):
579 bestCgRmsdFile = fullCgFileName
582 self.singlePdbResults(
585 self.getParam(
"best_cg_score_output_file"))
586 self.singlePdbResults(
589 self.getParam(
"best_cg_rmsd_output_file"))
590 self.singlePdbResults(
592 (cgFileName, framesToRead[-1]), -1, self.getParam(
"final_cg_frame_output_file"))
593 finalCgRmsd = self.calculateNativeRmsd(flexibleAtoms)
594 print "final cg rmsd is %s " % finalCgRmsd
595 self.bestCgScore = round(bestCgScore, 2)
596 self.bestCgRmsd = round(bestCgRmsd, 2)
597 self.bestCgScoreFile = bestCgScoreFile
598 self.bestCgRmsdFile = bestCgRmsdFile
600 def singlePdbResults(self, trajectoryFile, frame, outputPdbFile):
602 fullTrajectoryFile = os.path.join(
603 self.getParam(
"output_directory"),
605 fullOutputFile = os.path.join(
606 self.getParam(
"output_directory"),
608 rh = RMF.open_rmf_file(fullTrajectoryFile)
609 IMP.rmf.set_hierarchies(rh, [self.protein])
611 frame = IMP.rmf.get_number_of_frames(rh, self.protein) - 1
615 def calculateRmsd(self, otherProtein, flexibleAtoms):
616 otherNamesToParticles = atomicDominoUtilities.makeNamesToParticles(
620 for pName
in otherNamesToParticles.keys():
621 if ((pName
in flexibleAtoms) == 0):
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())
630 def calculateNativeRmsd(self, flexibleAtoms):
632 if (self.wroteNativeProtein == 0):
633 pdbName = self.getParam(
"native_pdb_input_file")
639 self.wroteNativeProtein = 1
641 return self.calculateRmsd(self.nativeProtein, flexibleAtoms)
643 def calculateTrajectoryRmsd(
648 pdbName = self.getParam(
"native_pdb_input_file")
654 outputDir = self.getParam(
"output_directory")
655 fullFile = os.path.join(outputDir, trajectoryFile)
656 print "open calculate traj rmf %s" % fullFile
657 rh = RMF.open_rmf_file(fullFile)
658 IMP.rmf.set_hierarchies(rh, [otherProtein])
659 if (trajectoryFrame == -1):
660 trajectoryFrame = IMP.rmf.get_number_of_frames(
661 rh, otherProtein) - 1
663 return self.calculateRmsd(otherProtein, flexibleAtoms)
665 def createAllSubsetAssignments(self):
669 self.model.get_root_restraint_set(),
672 leafNodeIndexList = self.getLeafNodeIndexList()
674 for nodeIndex
in leafNodeIndexList:
676 subset = self.mt.get_vertex_name(nodeIndex)
677 particleNameList = []
678 for particle
in subset:
679 particleNameList.append(self.quickParticleName(particle))
680 print "creating initial assignments for leaf %s" % nodeIndex
682 assignments = self.createUniqueLeafAssignments(
683 particleNameList, self.particleInfo)
684 filteredAssignments = self.filterAssignments(
685 assignments, subset, nodeIndex, rssft)
689 for assignment
in filteredAssignments:
690 packedAssignmentContainer.add_assignment(assignment)
691 lat.set_assignments(subset, packedAssignmentContainer)
697 root = self.mt.get_vertices()[-1]
698 completeAc = self.loadAssignments(root)
699 self.completeAc = completeAc
701 def loadAssignments(self, nodeIndex):
703 children = self.mt.get_out_neighbors(nodeIndex)
704 subset = self.mt.get_vertex_name(nodeIndex)
705 heapCount = int(self.getParam(
"heap_count"))
707 heapCount, self.rssft.get_subset_filter(subset, []))
708 if len(children) == 0:
709 print "computing assignments for leaf %s" % nodeIndex
711 self.sampler.load_vertex_assignments(nodeIndex, mine)
712 print "leaf node %s has %s leaf assignments" % (nodeIndex, mine.get_number_of_assignments())
714 if (children[0] > children[1]):
715 children = [children[1], children[0]]
717 firstAc = self.loadAssignments(children[0])
718 secondAc = self.loadAssignments(children[1])
720 self.sampler.load_vertex_assignments(
721 nodeIndex, firstAc, secondAc, mine)
723 timeDifference = self.logTimePoint(0)
724 print "Done Parent %s Assignments %s first child %s second child %s time %s" % (nodeIndex, mine.get_number_of_assignments(), firstAc.get_number_of_assignments(),
725 secondAc.get_number_of_assignments(), timeDifference)
726 self.totalAssignments += mine.get_number_of_assignments()
730 def writeOutput(self, flexibleAtoms, startTime):
732 bestDominoScore = 100000
734 finalAssignments = self.completeAc.get_assignments()
735 for assignment
in finalAssignments:
737 self.dominoPst.get_subset(),
740 score = self.model.evaluate(
False)
741 if (score < bestDominoScore):
742 bestAssignment = assignment
743 bestDominoScore = round(score, 2)
744 print "best domino score is %s " % bestDominoScore
745 print "best md score is %s" % self.bestMdScore
746 print "best md rmsd is %s" % self.bestMdRmsd
747 print "best cg score is %s" % self.bestCgScore
748 print "best cg rmsd is %s" % self.bestCgRmsd
749 print "merge tree contained %s total assignments" % self.totalAssignments
752 self.dominoPst.get_subset(),
755 dominoVsMdRmsd = round(
756 self.calculateTrajectoryRmsd(
758 "md_trajectory_output_file"),
759 self.bestMdScoreFrame,
766 os.path.join(self.getParam(
"output_directory"),
767 self.getParam(
"minimum_domino_score_pdb")))
769 dominoVsCgRmsd = round(
770 self.calculateTrajectoryRmsd(
771 self.bestCgScoreFile,
775 dominoMinimizedScore = round(self.model.evaluate(
False), 2)
776 dominoRmsd = round(self.calculateNativeRmsd(flexibleAtoms), 2)
778 runTime = round(time.time() - startTime, 2)
779 print "final domino score (after cg): %s" % dominoMinimizedScore
780 print "final domino rmsd: %s" % dominoRmsd
781 print "best domino rmsd with best md score: %s" % dominoVsMdRmsd
782 print "domino rmsd with best cg score: %s" % dominoVsCgRmsd
783 print "Final Results\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.bestMdScore, self.bestMdRmsd, self.bestCgScore, self.bestCgRmsd, bestDominoScore, dominoRmsd, dominoMinimizedScore, dominoVsCgRmsd, self.totalAssignments, self.maxMem, runTime)
788 def createSubsetFromParticles(self, particleNames):
789 particleNameList = particleNames.split(
" ")
791 for particleName
in particleNameList:
792 particle = self.namesToParticles[particleName]
793 particleList.append(particle)
796 return [subset, particleList]
798 def createHdf5AssignmentContainer(self, index, particleNames, read):
799 root = self.getAssignmentContainerRoot(index, read)
800 print "got root for index %s" % index
803 dataset = root.get_child_index_data_set_2d(str(index))
805 dataset = root.add_child_index_data_set_2d(str(index))
806 print "added child index dataset"
808 [subset, particleOrder] = self.createSubsetFromParticles(particleNames)
809 print "created subset for index %s" % index
810 hdf5Ac = IMP.domino.create_assignments_container(
811 dataset, subset, particleOrder)
812 print "returning from create"
814 for nextInt
in order:
815 print "next int is %s" % nextInt
816 return [subset, hdf5Ac]
818 def loadAssignmentsParallel(
823 mtIndexToSubsetOrder,
827 if ((
"firstChild" in mtIndexToNodeInfo[nodeIndex]) == 0):
828 print "writing file for leaf index %s" % nodeIndex
830 self.createAssignmentsParallel(
837 beginTime = self.logTimePoint(1)
838 firstChildIndex = mtIndexToNodeInfo[nodeIndex][
"firstChild"]
839 [firstSubset, firstAc] = self.createHdf5AssignmentContainer(
840 firstChildIndex, mtIndexToSubsetOrder[firstChildIndex], 1)
841 firstAcCreateTime = self.logTimePoint(0)
843 secondChildIndex = mtIndexToNodeInfo[nodeIndex][
"secondChild"]
844 [secondSubset, secondAc] = self.createHdf5AssignmentContainer(
845 secondChildIndex, mtIndexToSubsetOrder[secondChildIndex], 1)
846 secondAcCreateTime = self.logTimePoint(0)
848 print "getting assignments for nodeIndex %s first child %s second child %s" % (nodeIndex, firstChildIndex, secondChildIndex)
849 firstChildParticles = mtIndexToSubsetOrder[firstChildIndex]
850 secondChildParticles = mtIndexToSubsetOrder[secondChildIndex]
851 myParticles = mtIndexToParticles[nodeIndex]
853 print "first child particles %s\nsecond child particles %s\nmy particles; %s\n" % (firstChildParticles, secondChildParticles, myParticles)
854 for p
in firstSubset:
855 print "next particle in first subset: %s" % self.quickParticleName(p)
857 for p
in secondSubset:
858 print "next particle in second subset: %s" % self.quickParticleName(p)
860 for assignment
in firstAc.get_assignments():
861 print "next assignment for first child %s: %s" % (firstChildIndex, assignment)
863 for assignment
in secondAc.get_assignments():
864 print "next assignment for second child %s: %s" % (secondChildIndex, assignment)
866 [mySubset, myAc] = self.createHdf5AssignmentContainer(
867 nodeIndex, mtIndexToParticles[nodeIndex], 0)
868 print "done creating hdf5"
869 prepTime = self.logTimePoint(0)
871 self.model.get_root_restraint_set(),
873 rssf = rssft.get_subset_filter(mySubset, [])
885 heapTime = self.logTimePoint(0)
887 addTime = self.logTimePoint(0)
888 for assignment
in myAc.get_assignments():
889 print "loadAssignmentsParallel next assignment for %s: %s" % (nodeIndex, assignment)
890 doneTime = self.logTimePoint(0)
891 firstChildCount = firstAc.get_number_of_assignments()
892 secondChildCount = secondAc.get_number_of_assignments()
893 print "first count: %s second count: %s begin: %s firstAc: %s secondAc: %s prep: %s heap: %s add: %s done: %s" % (firstChildCount, secondChildCount, beginTime, firstAcCreateTime, secondAcCreateTime, prepTime, heapTime, addTime, doneTime)
894 subsetOrder = self.getSubsetOrderList(mySubset)
897 def createAssignmentsParallel(
903 subsetName = mtIndexToParticles[nodeIndex]
904 print "starting assignments parallel leaf index %s subset name %s" % (nodeIndex, subsetName)
905 [subset, particleList] = self.createSubsetFromParticles(subsetName)
907 particleNameList = subsetName.split(
" ")
910 finalAssignments = self.createUniqueLeafAssignments(
911 particleNameList, particleInfo)
918 self.model.get_root_restraint_set(),
920 filteredAssignments = self.filterAssignments(
921 finalAssignments, subset, nodeIndex, rssft)
922 root = self.getAssignmentContainerRoot(nodeIndex, 0)
923 dataset = root.add_child_index_data_set_2d(str(nodeIndex))
924 dataset.set_size([0, len(subset)])
926 hdf5AssignmentContainer = IMP.domino.create_assignments_container(
927 dataset, subset, particleOrder)
928 for assignment
in filteredAssignments:
929 hdf5AssignmentContainer.add_assignment(assignment)
930 for assignment
in hdf5AssignmentContainer.get_assignments():
931 print "hdf5 assignment container node %s next assignmet %s" % (nodeIndex, assignment)
934 subsetOrder = self.getSubsetOrderList(subset)
935 print "leaf node returning with order %s" % subsetOrder
940 def writePymolData(self):
942 outputDir = self.getParam(
"output_directory")
944 pymolInteractions = self.getParam(
"pymol_interactions_file")
950 pymolSubsets = self.getParam(
"pymol_subsets_file")
958 def cleanVertexName(self, vertexName):
962 nodeRe = re.compile(
'Subset\(\[(.*?)\s*\]')
963 secondNodeRe = re.compile(
'\[(.*?)\s*\]')
964 node = nodeRe.search(str(vertexName))
965 secondNode = secondNodeRe.search(str(vertexName))
969 foundName = node.group(1)
971 foundName = secondNode.group(1)
972 vertexNameFinal = foundName.replace(
'"',
'')
973 return vertexNameFinal
975 def getSubsetOrderList(self, subset):
977 for particle
in subset:
978 name = self.quickParticleName(particle)
979 subsetOrderList.append(name)
980 subsetOrder =
" ".join(subsetOrderList)
983 def checkAssignments(self, subset, nodeIndex, particleOrder):
984 print "reading back in assignments for leaf index %s" % nodeIndex
985 root = self.getAssignmentContainerRoot(nodeIndex, 1)
986 dataset = root.get_child_index_data_set_2d(str(nodeIndex))
987 hdf5 = IMP.domino.create_assignments_container(
988 dataset, subset, particleOrder)
989 for assignment
in hdf5.get_assignments():
990 print "leaf index %s read back in %s" % (nodeIndex, assignment)
992 def createSamplerLite(self):
995 if (self.getParam(
"cross_subset_filtering") == 1):
996 s.set_use_cross_subset_filtering(1)
1000 def createSiblingMap(self, parentIndex):
1002 children = self.mt.get_out_neighbors(parentIndex)
1003 if (len(children) > 0):
1004 firstChild = children[0]
1005 secondChild = children[1]
1006 self.parentSiblingMap[firstChild] = {}
1007 self.parentSiblingMap[firstChild][
"sibling"] = secondChild
1008 self.parentSiblingMap[firstChild][
"parent"] = parentIndex
1010 self.parentSiblingMap[secondChild] = {}
1011 self.parentSiblingMap[secondChild][
"sibling"] = firstChild
1012 self.parentSiblingMap[secondChild][
"parent"] = parentIndex
1013 print "created map for parent %s first child %s second child %s" % (parentIndex, firstChild, secondChild)
1015 self.parentSiblingMap[parentIndex][
"firstChild"] = firstChild
1016 self.parentSiblingMap[parentIndex][
"secondChild"] = secondChild
1018 self.createSiblingMap(firstChild)
1019 self.createSiblingMap(secondChild)
1021 def getMtIndexToNodeInfo(self):
1022 return self.parentSiblingMap
1024 def getLeafParentNodeIndexList(self):
1025 leafParentNodeIndexList = {}
1026 leafNodeIndexList = self.getLeafNodeIndexList()
1027 for leafIndex
in leafNodeIndexList:
1028 parent = self.parentSiblingMap[leafIndex][
"parent"]
1029 leafParentNodeIndexList[parent] = 1
1030 return leafParentNodeIndexList
1032 def getMtIndexToNameList(selt):
1034 for index
in self.mt.get_vertices():
1035 name = self.mt.get_vertex_name(index)
1036 mtIndexToNames[index] = name
1037 return mtIndexToNames
1039 def getAssignmentContainerRoot(self, subsetIndex, read):
1040 outputDir = self.getParam(
"output_directory")
1041 filePrefix = self.getParam(
"subset_assignment_output_file")
1042 assignmentFileName = os.path.join(
1044 (filePrefix, subsetIndex))
1045 print "creating hdf5 file with name %s" % assignmentFileName
1048 root = RMF.open_hdf5_file(assignmentFileName)
1050 root = RMF.create_hdf5_file(assignmentFileName)
1057 def writeVisualization(self):
1059 self.writeCytoscapeIgInput()
1060 self.writeCytoscapeJtInput()
1061 self.writeCytoscapeMtInput()
1062 self.writeCytoscapeScripts()
1063 self.writePymolData()
1065 def writeCytoscapeScripts(self):
1066 outputDir = self.getParam(
"output_directory")
1067 mTreeCytoscapeInput = self.getParam(
"mtree_cytoscape_input_file")
1068 mTreeCytoscapeAssignments = self.getParam(
1069 "mtree_cytoscape_assignment_file")
1070 mTreeCytoscapeAtomChains = self.getParam(
1071 "mtree_cytoscape_atom_chain_file")
1072 mTreeCytoscapeAtomSummary = self.getParam(
1073 "mtree_cytoscape_atom_summary_file")
1075 mTreeCytoscapeScript = self.getParam(
"mtree_cytoscape_script")
1076 mTreeFh = open(os.path.join(outputDir, mTreeCytoscapeScript),
'w')
1078 "network import file=%s\n" %
1079 os.path.join(outputDir, mTreeCytoscapeInput))
1081 "node import attributes file=\"%s\"\n" %
1082 os.path.join(outputDir, mTreeCytoscapeAssignments))
1084 "node import attributes file=\"%s\"\n" %
1085 os.path.join(outputDir, mTreeCytoscapeAtomSummary))
1087 "node import attributes file=\"%s\"\n" %
1088 os.path.join(outputDir, mTreeCytoscapeAtomChains))
1089 mTreeFh.write(
"layout jgraph-tree\n")
1091 def getGraphStructure(self, graph, fileName, separator):
1093 graphLogWrite = open(fileName,
'w')
1094 graph.show(graphLogWrite)
1095 graphLogWrite.close()
1098 graphLogRead = open(fileName,
'r')
1100 nodeRe = re.compile('^(\d+)\[label\=\"\[*(.*?)\s*\]*\"')
1101 separatorEscape =
"\\" + separator
1102 edgeString =
"^(\d+)\-%s(\d+)" % separatorEscape
1103 edgeRe = re.compile(edgeString)
1109 for line
in graphLogRead:
1112 node = nodeRe.search(line)
1114 nodeNumber = node.group(1)
1115 atomString = node.group(2)
1116 nodesToNames[nodeNumber] = atomString
1120 edge = edgeRe.search(line)
1122 firstNode = edge.group(1)
1123 secondNode = edge.group(2)
1125 if (firstNode
in nodesToNodes):
1126 firstNodeDict = nodesToNodes[firstNode]
1127 firstNodeDict[secondNode] = 1
1128 nodesToNodes[firstNode] = firstNodeDict
1130 return [nodesToNames, nodesToNodes]
1132 def writeEdgeFile(self, nodesToNodes, edgeFileName):
1134 outputDir = self.getParam(
"output_directory")
1135 graphInputFile = open(os.path.join(outputDir, edgeFileName),
'w')
1136 for firstNode
in nodesToNodes.keys():
1137 nodeDict = nodesToNodes[firstNode]
1138 for secondNode
in nodeDict.keys():
1139 graphInputFile.write(
"%s ttt %s\n" % (firstNode, secondNode))
1140 graphInputFile.close()
1142 def writeCytoscapeIgInput(self):
1143 outputDir = self.getParam(
"output_directory")
1144 igOutputFile = self.getParam(
"ig_output_file")
1145 [nodesToNames, nodesToNodes] = self.getGraphStructure(
1146 self.ig, os.path.join(outputDir, igOutputFile),
"-")
1150 self.getParam(
"ig_cytoscape_input_file"))
1153 igResiduesFile = self.getParam(
"ig_cytoscape_residues_file")
1154 peptideChain = self.getParam(
"peptide_chain")
1156 subsetResidueFile = open(os.path.join(outputDir, igResiduesFile),
'w')
1157 subsetResidueFile.write(
"ResidueNumber\n")
1158 for nodeNumber
in nodesToNames.keys():
1159 nodeName = nodesToNames[nodeNumber]
1161 [nodeChain, residueNumber,
1162 nodeAtom] = atomicDominoUtilities.getAtomInfoFromName(nodeName)
1163 if (nodeChain == peptideChain):
1164 subsetResidueFile.write(
1166 (nodeNumber, residueNumber))
1170 subsetResidueFile.write(
"%s = 100\n" % nodeNumber)
1172 def writeCytoscapeMtInput(self):
1173 outputDir = self.getParam(
"output_directory")
1174 mTreeOutputFile = self.getParam(
"mtree_output_file")
1175 [nodesToNames, nodesToNodes] = self.getGraphStructure(
1176 self.mt, os.path.join(outputDir, mTreeOutputFile),
">")
1180 self.getParam(
"mtree_cytoscape_input_file"))
1181 self.writeNodeNameAttributes(
1182 nodesToNames, self.getParam(
"mtree_cytoscape_atom_name_file"), self.getParam(
1183 "mtree_cytoscape_atom_summary_file"),
1184 self.getParam(
"mtree_cytoscape_atom_chain_file"))
1186 def writeCytoscapeJtInput(self):
1188 outputDir = self.getParam(
"output_directory")
1189 jTreeOutputFile = self.getParam(
"jtree_output_file")
1190 [nodesToNames, nodesToNodes] = self.getGraphStructure(
1191 self.jt, os.path.join(outputDir, jTreeOutputFile),
"-")
1195 self.getParam(
"jtree_cytoscape_input_file"))
1197 self.writeNodeNameAttributes(
1198 nodesToNames, self.getParam(
"jtree_cytoscape_atom_name_file"), self.getParam(
1199 "jtree_cytoscape_atom_summary_file"),
1200 self.getParam(
"jtree_cytoscape_atom_chain_file"))
1204 jtreeCytoscapeEdgeFile = self.getParam(
"jtree_cytoscape_edge_file")
1205 edgeWeightFile = open(
1206 os.path.join(outputDir,
1207 jtreeCytoscapeEdgeFile),
1209 edgeWeightFile.write(
"SubsetOverlap (class=Integer)\n")
1210 for firstNode
in nodesToNodes.keys():
1211 nodeDict = nodesToNodes[firstNode]
1212 for secondNode
in nodeDict.keys():
1213 firstNodeAtoms = nodesToNames[firstNode].split(
" ")
1214 secondNodeAtoms = nodesToNames[secondNode].split(
" ")
1216 val
for val
in firstNodeAtoms
if val
in secondNodeAtoms]
1217 edgeWeightFile.write(
1218 "%s (pp) %s = %s\n" %
1219 (firstNode, secondNode, len(intersection)))
1220 edgeWeightFile.close()
1222 def getAtomTypeCounts(self, atomNames):
1224 atomNames = atomNames.lstrip(
'[')
1225 atomNames = atomNames.rstrip(
']')
1226 atomNames = atomNames.rstrip()
1228 atomList = atomNames.split(
" ")
1229 peptideAtomCount = 0
1230 proteinAtomCount = 0
1231 for atom
in atomList:
1232 [chain, residue, atom] = atomicDominoUtilities.getAtomInfoFromName(
1234 if (chain == self.getParam(
"peptide_chain")):
1235 peptideAtomCount += 1
1237 proteinAtomCount += 1
1238 return [peptideAtomCount, proteinAtomCount]
1240 def writeNodeNameAttributes(
1247 outputDir = self.getParam(
"output_directory")
1248 subsetAtomNameFile = open(os.path.join(outputDir, atomNameFile),
'w')
1249 subsetAtomSummaryFile = open(
1250 os.path.join(outputDir, atomSummaryFile),
'w')
1251 subsetAtomChainFile = open(os.path.join(outputDir, atomChainFile),
'w')
1253 subsetAtomNameFile.write(
"Atom names\n")
1254 subsetAtomSummaryFile.write(
"Atom Summary\n")
1255 subsetAtomChainFile.write(
"Atom chain\n")
1256 for node
in nodesToNames.keys():
1257 atomNames = nodesToNames[node]
1258 subsetAtomNameFile.write(
"%s = %s\n" % (node, atomNames))
1259 [peptideAtomCount, proteinAtomCount] = self.getAtomTypeCounts(
1262 subsetAtomSummaryFile.write(
1264 (node, proteinAtomCount, peptideAtomCount))
1267 if (proteinAtomCount == 0):
1268 subsetAtomChainFile.write(
"%s = 1\n" % node)
1269 elif (peptideAtomCount == 0):
1270 subsetAtomChainFile.write(
"%s = 2\n" % node)
1272 subsetAtomChainFile.write(
"%s = 3\n" % node)
1274 subsetAtomChainFile.close()
1275 subsetAtomSummaryFile.close()
1276 subsetAtomNameFile.close()
Chain get_chain(Hierarchy h)
InteractionGraph get_interaction_graph(ScoringFunctionAdaptor rs, const kernel::ParticlesTemp &pst)
void write_pdb(const Selection &mhd, base::TextOutput out, unsigned int model=1)
Grid3D< int, SparseGridStorage3D< int, UnboundedGridStorage3D > > SparseUnboundedIntGrid3D
kernel::ParticlesTemp get_order(const Subset &s, const SubsetFilterTables &sft)
void set_log_level(LogLevel l)
Set the current global log level.
SubsetGraph get_junction_tree(const InteractionGraph &ig)
Simple conjugate gradients optimizer.
Sample best solutions using Domino.
algebra::BoundingBoxD< 3 > get_bounding_box(const Hierarchy &h)
Get a bounding box for the Hierarchy.
Subsets get_subsets(const SubsetGraph &g)
Gets all of the Subsets of a SubsetGraph.
Filter a configuration of the subset using the kernel::Model thresholds.
Represent a subset of the particles being optimized.
void load_frame(RMF::FileConstHandle file, unsigned int frame)
Grid3D< double, DenseGridStorage3D< double > > DenseDoubleGrid3D
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
MergeTree get_balanced_merge_tree(const SubsetGraph &junction_tree)
display::Geometries get_subset_graph_geometry(const SubsetGraph &ig)
Select all non-alternative ATOM records.
Hierarchies get_by_type(Hierarchy mhd, GetByType t)
void load_merged_assignments(const Subset &first_subset, AssignmentContainer *first, const Subset &second_subset, AssignmentContainer *second, const SubsetFilterTablesTemp &filters, AssignmentContainer *ret)
Fill in assignments for an internal node.
A decorator for a particle with x,y,z coordinates.
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.
Write a CGO file with the geometry.
display::Geometries get_interaction_graph_geometry(const InteractionGraph &ig)
void load_particle_states(const Subset &s, const Assignment &ss, const ParticleStatesTable *pst)
void read_pdb(base::TextInput input, int model, Hierarchy h)
Hierarchies get_leaves(const Selection &h)
Support for the RMF file format for storing hierarchical molecular data and markup.
Divide-and-conquer inferential optimization in discrete space.
Class for storing model, its restraints, constraints, and particles.