""" refmac.py: CCP4 GUI Project Copyright (C) 2010 University of York 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,S 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. """ from CCP4PluginScript import CPluginScript import CCP4ErrorHandling from CCP4ErrorHandling import * import CCP4Modules from lxml import etree class refmac_i2(CPluginScript): TASKMODULE = 'wrappers' TASKTITLE = 'Refinement (Refmac5)' TASKNAME = 'refmac' TASKCOMMAND = 'refmac5' TASKVERSION= 0.0 WHATNEXT = ['prosmart_refmac','buccaneer_build_refine_mr'] ASYNCHRONOUS = False ERROR_CODES = { 201 : {'description' : 'Refmac returned with non zero status' }, 202: {'description': 'New library created but strictly required' }, 203: {'description': 'New library created', 'severity':CCP4ErrorHandling.SEVERITY_WARNING}, 204: {'description': 'Program completed without generating XMLOUT.' }, } def __init__(self,*args, **kwargs): super(refmac_i2, self).__init__(*args, **kwargs) self._readyReadStandardOutputHandler = self.handleReadyReadStandardOutput self.xmlroot = etree.Element('REFMAC') from refmacLogScraper import logScraper self.logScraper = logScraper(xmlroot=self.xmlroot, flushXML=self.flushXML) self.xmlLength = 0 def handleReadyReadStandardOutput(self): if not hasattr(self,'logFileHandle'): self.logFileHandle = open(self.makeFileName('LOG'),'w') if not hasattr(self,'logFileBuffer'): self.logFileBuffer = '' pid = self.getProcessId() qprocess = CCP4Modules.PROCESSMANAGER().getJobData(pid,attribute='qprocess') availableStdout = qprocess.readAllStandardOutput() self.logFileHandle.write(availableStdout) self.logFileHandle.flush() self.logScraper.processLogChunk(str(availableStdout)) def flushXML(self): newXml = etree.tostring(self.xmlroot,pretty_print=True) if len(newXml)>self.xmlLength: self.xmlLength = len(newXml) with open (self.makeFileName('PROGRAMXML')+'_tmp','w') as programXmlFile: programXmlFile.write(newXml) import shutil shutil.move(self.makeFileName('PROGRAMXML')+'_tmp', self.makeFileName('PROGRAMXML')) def processInputFiles(self): import CCP4XtalData error = None self.hklin = None dataObjects = [] #print '\n\n\n***contentFlag',self.container.inputData.F_SIGF.contentFlag #Append Observation with representation dependent on whether we are detwining on Is or not obsTypeRoot = 'CONTENT_FLAG_F' if self.container.controlParameters.USE_TWIN and self.container.controlParameters.TWIN_TYPE.isSet() and self.container.controlParameters.TWIN_TYPE=="I": obsTypeRoot = 'CONTENT_FLAG_I' obsPairOrMean = 'MEAN' if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__() != 'NOTHING': obsPairOrMean = 'PAIR' obsType = getattr(CCP4XtalData.CObsDataFile, obsTypeRoot+obsPairOrMean) dataObjects += [['F_SIGF',obsType]] #Include phase estimates if called for if self.container.inputData.ABCD.isSet(): dataObjects += ['ABCD'] #Apply coordinate selection if set import os self.inputCoordPath = os.path.normpath(self.container.inputData.XYZIN.fullPath.__str__()) if self.container.inputData.XYZIN.isSelectionSet(): self.inputCoordPath = os.path.normpath(os.path.join(self.getWorkDirectory(),'selected.pdb')) self.container.inputData.XYZIN.getSelectedAtomsPdbFile(self.inputCoordPath) #Include FreeRflag if called for if self.container.inputData.FREERFLAG.isSet(): dataObjects += ['FREERFLAG'] self.hklin,error = self.makeHklin(dataObjects) if error.maxSeverity()>CCP4ErrorHandling.SEVERITY_WARNING: return CPluginScript.FAILED else: return CPluginScript.SUCCEEDED def processOutputFiles(self): if hasattr(self,'logFileHandle'): self.logFileHandle.close() else: self.xmlroot.clear() self.logScraper.scrapeFile( self.makeFileName('LOG') ) #First up check for exit status of the program from CCP4Modules import PROCESSMANAGER from lxml import etree exitStatus = 0 exitCode=0 try: exitStatus = PROCESSMANAGER().getJobData(pid=self.getProcessId(), attribute='exitStatus') except Exception as e: print e self.appendErrorReport(201,'Exit status: Unable to recover exitStatus') return CPluginScript.FAILED if exitStatus != 0: self.appendErrorReport(201,'Exit status: '+str(exitStatus)) return CPluginScript.FAILED #Now the exit codes...I think that non zero means Refmac identified an issue try: exitCode = PROCESSMANAGER().getJobData(pid=self.getProcessId(), attribute='exitCode') except: self.appendErrorReport(201,'Exit code: Unable to recover exitCode') return CPluginScript.FAILED if exitCode != 0: import os try: logFileText = open(self.makeFileName('LOG')).read() if 'Your coordinate file has a ligand which has either minimum or no description in the library' in logFileText and self.container.controlParameters.MAKE_NEW_LIGAND_EXIT.isSet() and self.container.controlParameters.MAKE_NEW_LIGAND_EXIT: self.appendErrorReport(201,'You did not supply a full ligand geometry file: either make and supply one (Make Ligand task), or set the appropriate flag in the advanced options') import re #Example line: * Plotfile: /tmp/martin/refmac5_temp1.64630_new_TM7.ps plotFiles = re.findall(r'^.*\* Plotfile:.*$',logFileText,re.MULTILINE) print plotFiles for plotFile in plotFiles: psfileName = plotFile.split()[-1] import shutil shutil.copyfile(psfileName, self.container.outputData.PSOUT.__str__()) self.container.outputData.PSOUT.annotation.set(psfileName+' from REFMAC') return CPluginScript.UNSATISFACTORY else: self.appendErrorReport(201,'Exit code: '+str(exitCode)) except: self.appendErrorReport(201,'Exit code: '+str(exitCode)) return CPluginScript.FAILED import CCP4XtalData import CCP4File import os # Need to set the expected content flag for phases data self.container.outputData.XYZOUT.annotation = 'Model from refinement' self.container.outputData.FPHIOUT.annotation = 'Weighted map from refinement' self.container.outputData.DIFFPHIOUT.annotation = 'Weighted difference map from refinement' self.container.outputData.ABCDOUT.annotation = 'Calculated phases from refinement' self.container.outputData.ABCDOUT.contentFlag = CCP4XtalData.CPhsDataFile.CONTENT_FLAG_HL self.container.outputData.TLSOUT.annotation = 'TLS parameters from refinement' self.container.outputData.LIBOUT.annotation = 'Generated dictionary from refinement' self.container.outputData.ANOMFPHIOUT.annotation = 'Anomalous map coefficients' self.container.outputData.DIFANOMFPHIOUT.annotation = 'Difference anomalous map coefficients (LLG map)' # Split out data objects that have been generated. Do this after applying the annotation, and flagging # above, since splitHklout needs to know the ABCDOUT contentFlag outputFiles = ['FPHIOUT','DIFFPHIOUT'] outputColumns = ['FWT,PHWT','DELFWT,PHDELWT'] if self.container.controlParameters.PHOUT: outputFiles+=['ABCDOUT'] outputColumns+=['HLACOMB,HLBCOMB,HLCCOMB,HLDCOMB'] if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__() != 'NOTHING': outputFiles += ['ANOMFPHIOUT'] outputColumns += ['FAN,PHAN'] import CCP4XtalData import os hkloutFile=CCP4XtalData.CMtzDataFile(os.path.join(self.getWorkDirectory(), "hklout.mtz")) hkloutFile.loadFile() columnLabelsInFile = [column.columnLabel.__str__() for column in hkloutFile.fileContent.listOfColumns] print 'columnLabelsInFile',columnLabelsInFile if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__() != 'NOTHING' and 'DELFAN' in columnLabelsInFile and 'PHDELAN' in columnLabelsInFile: outputFiles += ['DIFANOMFPHIOUT'] outputColumns += ['DELFAN,PHDELAN'] error = self.splitHklout(outputFiles,outputColumns) if error.maxSeverity()>CCP4ErrorHandling.SEVERITY_WARNING: return CPluginScript.FAILED #Use Refmacs XMLOUT as the basis for output XML. If not existent (probably due to failure), then create a new one import CCP4Utils rxml = None try: rxml = CCP4Utils.openFileToEtree(os.path.normpath(os.path.join(self.getWorkDirectory(),'XMLOUT.xml'))) except: rxml = etree.Element('REFMAC') self.appendErrorReport(204,self.makeFileName('PROGRAMXML')) return CPluginScript.FAILED #Copy across the segments of the scraped log file into this new xml root for childElement in self.xmlroot: rxml.append(childElement) #Extract performanceindictors from XML try: self.container.outputData.PERFORMANCEINDICATOR.RFactor = rxml.xpath("//new_cycle[last()]/r_factor")[-1].text self.container.outputData.PERFORMANCEINDICATOR.RFree = rxml.xpath("//new_cycle[last()]/r_free")[-1].text except: pass #Perform analysis of output coordinate file composition if os.path.isfile(str(self.container.outputData.XYZOUT.fullPath)): self.container.outputData.XYZOUT.fileContent.loadFile(self.container.outputData.XYZOUT.fullPath) modelCompositionNode = etree.SubElement(rxml,'ModelComposition') for chain in self.container.outputData.XYZOUT.fileContent.composition.chains: chainNode = etree.SubElement(modelCompositionNode,'Chain',id=chain) for monomer in self.container.outputData.XYZOUT.fileContent.composition.monomers: monomerNode = etree.SubElement(modelCompositionNode,'Monomer',id=monomer) #Skim smartie graphs from the log file smartieNode = etree.SubElement(rxml,'SmartieGraphs') self.scrapeSmartieGraphs(smartieNode) et = etree.ElementTree(rxml) #And write out the XML et.write(self.makeFileName('PROGRAMXML'), pretty_print=True) with open(self.container.outputData.COOTSCRIPTOUT.fullPath.__str__(),"w") as cootscript: #Write a GUI to regions that Refmac has identified as containing duffers badStretches = self.listOfTransgressingSegments(rxml) if len(badStretches) > 0: interestingBitsDef = 'interestingBits = [' for badStretch in badStretches: interestingBitsDef += ('{"chain":"%s","firstResidue":%s,"lastResidue":%s},'%(badStretch['chain'],badStretch['firstResidue'],badStretch['lastResidue'])) interestingBitsDef += ']\n' cootscript.write(interestingBitsDef) cootscript.write('ccp4i2Interface.addInterestingBitsMenu(title="Refmac-identified outliers", interestingBits=interestingBits)\n') return CPluginScript.SUCCEEDED def listOfTransgressingSegments(self, rxml): from sets import Set badResidueSet = Set() badStretches = [] outlierNodes = rxml.xpath('//OutliersByCriteria') if len(outlierNodes) == 0: return badStretches for outlierTypeNode in outlierNodes[0]: key = outlierTypeNode.tag for outlier in outlierTypeNode.xpath('Outlier'): res1Tuple = outlier.get('chainId1'),outlier.get('resId1') res2Tuple = outlier.get('chainId2'),outlier.get('resId2') badResidueSet.add(res1Tuple) badResidueSet.add(res2Tuple) badResidueTuple = [tuple for tuple in badResidueSet] #Sort on the basis of a string formed by adding the chain and residue elements of the tuple orderedBadResidues = sorted(badResidueTuple, key=lambda residue: residue[0]+residue[1]) for orderedBadResidue in orderedBadResidues: if len(badStretches) == 0 or badStretches[-1]['chain'] != orderedBadResidue[0] or int(orderedBadResidue[1])-int(badStretches[-1]['lastResidue']) > 1: badStretch = {"chain":orderedBadResidue[0], "firstResidue":orderedBadResidue[1], "lastResidue":orderedBadResidue[1]} badStretches.append(badStretch) else: badStretch["lastResidue"] = orderedBadResidue[1] return badStretches def scrapeSmartieGraphs(self, smartieNode): import sys, os import CCP4Utils smartiePath = os.path.join(CCP4Utils.getCCP4I2Dir(),'smartie') sys.path.append(smartiePath) import smartie from lxml import etree logfile = smartie.parselog(self.makeFileName('LOG')) for smartieTable in logfile.tables(): if smartieTable.ngraphs() > 0: tableelement = self.xmlForSmartieTable(smartieTable, smartieNode) return def xmlForSmartieTable(self, table, parent): import MGQTmatplotlib tableetree = MGQTmatplotlib.CCP4LogToEtree(table.rawtable()) parent.append(tableetree) return tableetree def makeCommandAndScript(self): import os self.hklout = os.path.join(self.workDirectory,"hklout.mtz") # make refmac command script self.appendCommandLine(['XYZIN',self.inputCoordPath]) self.appendCommandLine(['HKLIN',self.hklin]) if self.container.inputData.DICT.isSet(): self.appendCommandLine(['LIBIN',self.container.inputData.DICT.fullPath]) if self.container.inputData.TLSIN.isSet(): self.appendCommandLine(['TLSIN',self.container.inputData.TLSIN.fullPath]) self.appendCommandLine(['XYZOUT',self.container.outputData.XYZOUT.fullPath]) self.appendCommandLine(['HKLOUT',self.hklout]) #self.appendCommandLine(['LIBOUT',self.libout]) self.appendCommandLine(['LIBOUT',self.container.outputData.LIBOUT.fullPath]) self.appendCommandLine(['XMLOUT',os.path.normpath(os.path.join(self.getWorkDirectory(),'XMLOUT.xml'))]) if self.container.controlParameters.TITLE.isSet(): self.appendCommandScript("TITLE %s"%(str(self.container.controlParameters.TITLE))) # Main refinement options if self.container.controlParameters.NCYCLES.isSet(): self.appendCommandScript("NCYCLES %s"%(str(self.container.controlParameters.NCYCLES))) if str(self.container.controlParameters.WEIGHT_OPT) == 'AUTO': self.appendCommandScript("WEIGHT AUTO") elif self.container.controlParameters.WEIGHT.isSet(): self.appendCommandScript("WEIGHT MATRIX %s"%(str(self.container.controlParameters.WEIGHT))) if self.container.controlParameters.USE_TWIN: self.appendCommandScript("TWIN") if self.container.controlParameters.HYDR_USE: if str(self.container.controlParameters.HYDR_ALL) == 'ALL': self.appendCommandScript("MAKE HYDR ALL") else: self.appendCommandScript("MAKE HYDR YES") else: self.appendCommandScript("MAKE HYDR %s"%(str(self.container.controlParameters.HYDROGENS))) # Parameters if self.container.controlParameters.B_REFINEMENT_MODE.isSet(): self.appendCommandScript("REFI BREF %s"%(str(self.container.controlParameters.B_REFINEMENT_MODE))) if self.container.inputData.TLSIN.isSet() and self.container.controlParameters.NTLSCYCLES.isSet(): self.appendCommandScript("REFI TLSC %s"%(str(self.container.controlParameters.NTLSCYCLES))) elif self.container.controlParameters.AUTOTLS and self.container.controlParameters.NTLSCYCLES_AUTO.isSet(): self.appendCommandScript("REFI TLSC %s"%(str(self.container.controlParameters.NTLSCYCLES_AUTO))) if self.container.inputData.TLSIN.isSet() or self.container.controlParameters.AUTOTLS: self.appendCommandLine(["TLSOUT",self.container.outputData.TLSOUT.fullPath]) if self.container.controlParameters.TLSOUT_ADDU: self.appendCommandScript("TLSOUT ADDU") if self.container.controlParameters.TLSBFACSETUSE and self.container.controlParameters.TLSBFACSET.isSet(): self.appendCommandScript("BFAC SET "+str(self.container.controlParameters.TLSBFACSET)) elif self.container.controlParameters.BFACSETUSE and self.container.controlParameters.BFACSET.isSet(): self.appendCommandScript("BFAC SET "+str(self.container.controlParameters.BFACSET)) # Restraints if self.container.controlParameters.USE_NCS: if str(self.container.controlParameters.NCS_TYPE) == 'L': self.appendCommandScript("NCSR LOCAL") else: self.appendCommandScript("NCSR GLOBAL") elif self.container.controlParameters.USE_LOCAL_SYMMETRY: # USE_LOCAL_SYMMETRY is not used by the refinement pipeline, but is used by the bucaneer pipeline self.appendCommandScript("NCSR LOCAL") if self.container.controlParameters.USE_JELLY: if self.container.controlParameters.JELLY_SIGMA.isSet(): self.appendCommandScript("RIDG DIST SIGM %s"%(str(self.container.controlParameters.JELLY_SIGMA))) else: self.appendCommandScript("RIDG DIST SIGM 0.01") if self.container.controlParameters.JELLY_DIST.isSet(): self.appendCommandScript("RIDG DIST DMAX %s"%(str(self.container.controlParameters.JELLY_DIST))) # Output options if self.container.controlParameters.MAP_SHARP: if self.container.controlParameters.MAP_SHARP_CUSTOM and self.container.controlParameters.BSHARP.isSet(): self.appendCommandScript("MAPC SHAR "+str(self.container.controlParameters.BSHARP)) else: self.appendCommandScript("MAPC SHAR") # Advanced options if self.container.controlParameters.RES_CUSTOM: if self.container.controlParameters.RES_MIN.isSet() and self.container.controlParameters.RES_MAX.isSet(): self.appendCommandScript("REFI RESO %s %s"%(str(self.container.controlParameters.RES_MIN),str(self.container.controlParameters.RES_MAX))) elif self.container.controlParameters.RES_MIN.isSet() and not self.container.controlParameters.RES_MAX.isSet(): self.appendCommandScript("REFI RESO %s 999"%(str(self.container.controlParameters.RES_MIN))) elif not self.container.controlParameters.RES_MIN.isSet() and self.container.controlParameters.RES_MAX.isSet(): self.appendCommandScript("REFI RESO 0 %s"%(str(self.container.controlParameters.RES_MAX))) elif self.container.controlParameters.RESOLUTION.isSet(): # RESOLUTION is not used by the refinement pipeline, but is used by the bucaneer pipeline self.appendCommandScript("REFI RESO %s"%(str(self.container.controlParameters.RESOLUTION))) if self.container.controlParameters.EXTRAREFMACKEYWORDS.isSet(): for kwLine in str(self.container.controlParameters.EXTRAREFMACKEYWORDS).split('\n'): #print 'KwLine','['+str(kwLine)+']' self.appendCommandScript(kwLine.rstrip() + '\n') if self.container.controlParameters.MAKE_NEW_LIGAND_EXIT: self.appendCommandScript("MAKE NEWLIGAND EXIT") else: self.appendCommandScript("MAKE NEWLIGAND NOEXIT") if self.container.controlParameters.SCALETYPE.isSet(): if self.container.controlParameters.SCALETYPE.__str__() == 'REFMACDEFAULT': pass if self.container.controlParameters.SCALETYPE.__str__() == 'BABINET': self.appendCommandScript("SCALE TYPE BULK") self.appendCommandScript("SOLVENT NO") if self.container.controlParameters.SCALETYPE.__str__() == 'SIMPLE': self.appendCommandScript("SCALE TYPE SIMPLE") self.appendCommandScript("SOLVENT YES") if self.container.controlParameters.SCALETYPE.__str__() == 'WILSON': self.appendCommandScript("SCALE TYPE SIMPLE") self.appendCommandScript("SOLVENT NO") if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__ is not 'NONE': if self.container.controlParameters.USEANOMALOUSFOR.__str__() == 'OUTPUTMAPS': self.appendCommandScript("ANOMALOUS MAPONLY") elif self.container.controlParameters.USEANOMALOUSFOR.__str__() == 'SADREFINEMENT': self.appendCommandScript("REFINE OREFINE ANOMALOUS") if self.container.controlParameters.WAVELENGTH.isSet(): self.appendCommandScript("ANOMALOUS WAVELENGTH %10.5f"%self.container.controlParameters.WAVELENGTH) else: self.appendCommandScript("ANOMALOUS WAVELENGTH %10.5f"%self.container.inputData.F_SIGF.fileContent.getListOfWavelengths()[-1]) # Additional input/output if self.container.controlParameters.PHOUT: self.appendCommandScript("PHOUT") if self.container.inputData.EXTERNALRESTRAINTS.isSet(): self.appendCommandScript("@%s"%(str(self.container.inputData.EXTERNALRESTRAINTS.fullPath))) self.appendCommandScript("MONI DIST 1000000") labin = None if self.container.controlParameters.USE_TWIN and self.container.controlParameters.TWIN_TYPE.isSet() and self.container.controlParameters.TWIN_TYPE=="I": if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__() != 'NOTHING': labin = "LABIN I+=Iplus SIGI+=SIGIplus I-=Iminus SIGI-=SIGIminus" else: labin = "LABIN IP=I SIGIP=SIGI" else: if self.container.controlParameters.USEANOMALOUSFOR.isSet() and self.container.controlParameters.USEANOMALOUSFOR.__str__() != 'NOTHING': labin = "LABIN F+=Fplus SIGF+=SIGFplus F-=Fminus SIGF-=SIGFminus" else: labin = "LABIN FP=F SIGFP=SIGF" if self.container.inputData.ABCD.isSet(): import CCP4XtalData if self.container.inputData.ABCD.contentFlag == CCP4XtalData.CPhsDataFile.CONTENT_FLAG_HL: labin += " HLA=HLA HLB=HLB HLC=HLC HLD=HLD" else: labin += " PHIB=PHI FOM=FOM" if self.container.inputData.FREERFLAG.isSet(): labin += " FREE=FREER" print 'FreeR flag set' self.appendCommandScript(labin) self.appendCommandScript('END') return CPluginScript.SUCCEEDED def setProgramVersion(self): print 'refmac.getProgramVersion' return CPluginScript.setProgramVersion(self,'Refmac_5')