""" molrep_pipe.py: CCP4 GUI Project Copyright (C) 2015 STFC This library is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License version 3, modified in accordance with the provisions of the license to address the requirements of UK law. You should have received a copy of the modified GNU Lesser General Public License along with this library. If not, copies may be downloaded from http://www.ccp4.ac.uk/ccp4license.php This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. """ import os, shutil from CCP4PluginScript import CPluginScript import CCP4Utils from CCP4ErrorHandling import * from xml.etree import ElementTree as ET class molrep_pipe(CPluginScript): TASKNAME = 'molrep_pipe' MAINTAINER = 'liz.potterton@york.ac.uk' PERFORMANCECLASS = 'CRefinementPerformance' WHATNEXT = [] ERROR_CODES = { 301 : { 'description' : 'Error reading program xml output from first molrep run' }, 302 : { 'description' : 'No Laue results in program xml output from first molrep run' } } PURGESEARCHLIST = [ [ 'molrep_mr%*/align.pdb' , 1], [ 'molrep_mr%*/molrep_mtz.cif' , 1 ], [ 'molrep_mr%*/molrep.doc.txt' , 5 ] ] def process(self): if str(self.container.controlParameters.SG_OPTIONS) == 'specify': self.newspacegroup = str(self.container.controlParameters.SG) else: self.newspacegroup = str(self.container.inputData.F_SIGF.fileContent.spaceGroup) self.xmlroot = ET.Element('MolrepPipe') self.xmlroot.text = '\n' with open(str(self.makeFileName('PROGRAMXML')), 'w') as ostream: print >>ostream, ET.tostring(self.xmlroot) if str(self.container.controlParameters.SG_OPTIONS) == 'laue': self.run1() else: self.run2() return CPluginScript.SUCCEEDED def run1(self): self.molrep1 = self.makePluginObject('molrep_mr') self.molrep1.container.inputData.copyData(self.container.inputData) self.molrep1.container.controlParameters.copyData(self.container.controlParameters) self.molrep1.container.guiParameters.copyData(self.container.guiParameters) self.connectSignal(self.molrep1,'finished', self.fin1) self.molrep1.process() def readxml(self, xml): try: tree = ET.parse(xml).getroot() return tree except: self.appendErrorReport(301, xml) self.reportStatus(CPluginScript.FAILED) def unique_tags(self, tree): cou = 0 for e1 in tree.findall('RFpeaks'): e1.tag = 'RFpeaks' + str(cou) rf_list = list() for e2 in e1.findall('RFpeak'): rft = e2.findtext('RF') rf = float(rft) rf_list.append((rf, e2)) e1 = ET.SubElement(tree, 'RFsorted' + str(cou)) e1.text = '\n ' e1.tail = '\n ' for rf, e2 in sorted(rf_list): e1.append(e2) cou += 1 def appendxml(self, tree, tag): self.xmlroot.append(tree) tree.tail = '\n' tree.tag = tag with open(str(self.makeFileName('PROGRAMXML')), 'w') as ostream: print >>ostream, ET.tostring(self.xmlroot) def fin1(self, status): if status is not None and status.get('finishStatus') == CPluginScript.FAILED: self.reportStatus(status) return xml = self.molrep1.makeFileName('PROGRAMXML') tree = self.readxml(xml) if tree is None: return self.unique_tags(tree) laueE = tree.find('laue_group_alternatives') if laueE is None: self.appendErrorReport(302, xml) self.reportStatus(CPluginScript.FAILED) return laueData = [] testList = laueE.findall('test') eleNames = ['space_group', 'score', 'contrast'] for testE in testList: e1 = ET.SubElement(testE, 'selected') e1.text = '-' laueData.append({}) for name in eleNames: e = testE.find(name) if name == 'space_group': laueData[-1][name] = str(e.text) else: laueData[-1][name] = float(e.text) best = 0 for ind in range(1, len(laueData)): if (laueData[ind]['score'])>laueData[best]['score']: best = ind print 'Proceeding for best space group', laueData[best]['space_group'] self.newspacegroup = laueData[best]['space_group'] e1 = testList[best].find('selected') e1.text = 'yes' e0 = ET.SubElement(laueE, 'selected') e0.text = self.newspacegroup e0.tail = '\n' self.appendxml(tree, 'MolrepSpaceGroup') self.run2() def reindex(self, mtzin, mtzout): import subprocess as SP cmd = ('reindex', 'hklin', mtzin, 'hklout', mtzout) stdi = 'symm \'%s\'\nend\n' %self.newspacegroup sp = SP.Popen(cmd, stdin=SP.PIPE) sp.stdin.write(stdi) return sp.wait() def run2(self): if self.newspacegroup == str(self.container.inputData.F_SIGF.fileContent.spaceGroup): self.fobs = str(self.container.inputData.F_SIGF) self.free = str(self.container.inputData.FREERFLAG) else: self.fobs = str(self.container.outputData.F_SIGF) self.free = str(self.container.outputData.FREERFLAG) ret1 = self.reindex(str(self.container.inputData.F_SIGF), self.fobs) ret2 = self.reindex(str(self.container.inputData.FREERFLAG), self.free) if ret1 or ret2: self.reportStatus(CPluginScript.FAILED) return else: self.container.outputData.F_SIGF.setContentFlag(reset=True) self.container.outputData.F_SIGF.annotation = 'Observed data reindexed to'+str(self.newspacegroup) self.container.outputData.FREERFLAG.annotation = 'FreeR reindexed to'+str(self.newspacegroup) self.molrep2 = self.makePluginObject('molrep_mr') self.molrep2.container.inputData.copyData(self.container.inputData) self.molrep2.container.inputData.F_SIGF.set(self.fobs) self.molrep2.container.controlParameters.copyData(self.container.controlParameters) self.molrep2.container.controlParameters.SG_OPTIONS = 'no' self.molrep2.container.guiParameters.copyData(self.container.guiParameters) self.connectSignal(self.molrep2, 'finished', self.fin2) self.molrep2.process() def fin2(self, status): if status is not None and status.get('finishStatus') == CPluginScript.FAILED: self.reportStatus(status) return tree = self.readxml(self.molrep2.makeFileName('PROGRAMXML')) if tree is None: return self.unique_tags(tree) self.appendxml(tree, 'MolrepSearch') self.xyz = str(self.molrep2.container.outputData.XYZOUT) if os.path.isfile(self.xyz): self.run3() else: self.reportStatus(CPluginScript.SUCCEEDED) def run3(self): self.refmac = self.makePluginObject('refmac') self.refmac.container.inputData.XYZIN.set(self.xyz) self.refmac.container.inputData.F_SIGF.set(self.fobs) self.refmac.container.inputData.FREERFLAG.set(self.free) self.refmac.container.controlParameters.HYDROGENS = 'NO' self.refmac.container.controlParameters.NCYCLES = '5' self.refmac.container.controlParameters.PHOUT = False self.connectSignal(self.refmac, 'finished', self.fin3) self.refmac.process() def fin3(self, status): if status is not None and status.get('finishStatus') == CPluginScript.FAILED: self.reportStatus(status) return tree = self.readxml(self.refmac.makeFileName('PROGRAMXML')) if tree is None: return self.appendxml(tree, 'Refmac') try: shutil.copyfile(str(self.refmac.container.outputData.XYZOUT), str(self.container.outputData.XYZOUT)) shutil.copyfile(str(self.refmac.container.outputData.DIFFPHIOUT), str(self.container.outputData.DIFFPHIOUT)) shutil.copyfile(str(self.refmac.container.outputData.FPHIOUT), str(self.container.outputData.FPHIOUT)) cycles = list() for e3 in tree.findall('./Overall_stats/stats_vs_cycle/new_cycle'): cycles.append((int(e3.findtext('cycle')), float(e3.findtext('r_factor')), float(e3.findtext('r_free')))) cycle_no, Rcrist, Rfree = sorted(cycles)[-1] self.container.outputData.PERFORMANCE.RFactor = Rcrist self.container.outputData.PERFORMANCE.RFree = Rfree self.reportStatus(CPluginScript.SUCCEEDED) except: self.reportStatus(CPluginScript.FAILED)