#!/usr/bin/python import os,sys,shutil,subprocess,copy from process import process import common import fileinput par=common.parameter class mbref(process): name="iterative model building and refinement" short_name="model tracing with refinement" supported_procs=["ref", "mb"] supported_progs=["arpwarp"] supported_params={} supported_params['target'] = par( desc='Experiment/refinement target', typ=str, cap=True, share=True ) supported_params['bigcyc'] = par( desc='Number of cycles of model building', typ=int ) supported_params['low_res_par'] = par( desc='Use low resolution parameters in refinement and model building.', typ=bool, share=True ) supported_params['arp_sad'] = par( desc='Use the original ARP/wARP SAD procedure', typ=bool ) def RunPreprocess(self, rundir=None, **kwargs): process.RunPreprocess(self,rundir,kwargs) # get the number of monomers for arpwarp if that may be needed if 'arpwarp' in self.programs and not self.inp.Get(has_monomers=True) and self.inp.Get('sequence'): matthews=self.AddProcess('matthews', propagate_out=False) matthews.Run() self.processes.remove(matthews) def Init(self): self.mb_res_all, self.ref_res_all = [], [] def TreatInOutPar(self, set_all_par=False): # defaults if no program and no subprocess is specified if not self.GetProg(supported=True) and not self.GetProcess(supported=True): # maybe we could use arpwarp by default if a) available, b) high res (<2.5) ? self.AddProcess('ref') self.AddProcess('mb') if self.GetProcess(supported=True): # assign the child processes to specific class attributes for convenience self.ph = self.GetOrAddProcess('ref') self.mb = self.GetOrAddProcess('mb') # check and set some defaults if not self.GetParam('bigcyc'): self.SetParam('bigcyc',5) self.ph.SetRunDir() self.ph.TreatInOutPar() self.mb.SetRunDir() self.mb.TreatInOutPar() # always merge for MLHL (to prevent the originally input HL's be rewritten by the new HL's) if self.ph.GetParam('target')=='MLHL': self.ph.GetProg().always_merge=True #if os.environ['CCP4_LIB'] and os.path.isfile(os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.pdb')): # self.mb.GetProg().SetKey('pdbin-ref',os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.pdb')) # self.mb.GetProg().SetKey('mtzin-ref',os.path.join(os.environ['CCP4_LIB'],'data','reference_structures','reference-1ajr.mtz')) if not self.mb.GetProg().GetKey('sequence-reliability'): self.mb.GetProg().SetKey('sequence-reliability',0.99) #if not self.mb.GetProg().GetKey('resolution'): # self.mb.GetProg().SetKey('resolution',2.0) if not self.ph.GetParam('cycles'): self.ph.SetParam('cycles',10) if not self.ph.inp.Get('model', typ='substr', has_atomtypes=True) and \ self.ph.GetParam('target',capitalized=True) in ('SAD','SIRAS'): common.Error('No substructure model or atomtype(s) inputted for {0}.'.format(self.name)) else: # the arpwarp specific part shall follow arp=self.GetProg('arpwarp') arp_auto=self.GetProg('arpwarp_auto') if arp and not arp_auto and not arp.inp.Get('datafile',typ='par'): arp_auto=self.AddProg('arpwarp_auto',ind=0) arp_auto.inp=copy.copy(arp.inp) arp_auto.SetArg('parfile','arp.par') process.TreatInOutPar(self,set_all_par) def RunBody(self,*args,**kwargs): if self.processes: for self.mb.cyc in range(1, self.GetParam('bigcyc')+1): self.out.ClearAll() # build model if self.ph.out.mapcoef: self.mb.inp.Set(self.ph.out.mapcoef) if not self.mb.GetProg(): common.Error('No program specified for {0}'.format(self.mb.name)) self.mb.GetProg().runname=self.mb.GetProg().name+'_cyc'+str(self.mb.cyc) self.mb.Run() self.UpdateResults('mb') # run refinement if self.ph.out.Get('model',typ='substr',filetype='pdb'): self.ph.inp.Set(self.ph.out.Get('model',typ='substr',filetype='pdb')) self.ph.inp.Add(self.mb.out.Get('model')) self.ph.Run() self.UpdateResults('ph') # separate substructure/partial if needed if self.ph.out.Get('model',typ='partial+substr') and not self.ph.out.Get('model',typ='substr'): separate_models=self.ph.GetOrAddProcess('sepsubstrprot') separate_models.inp.Set(self.ph.out.Get('model',typ='partial+substr')) separate_models.Run() # set correlation mode for Buccaneer if self.mb.GetProg('buccaneer'): self.mb.GetProg('buccaneer').inp.Add( self.out.Get('model',typ='partial') ) else: # arpwarp (or possibly other mb+ref programs in the future) arp_auto=self.GetProg('arpwarp_auto') arp=self.GetProg('arpwarp') if arp_auto: arp_auto.Run() if arp: arp.inp.Add(arp_auto.out.Get('datafile',typ='par')) arp.fsf=arp_auto.fsf arp.SetRunDir(arp_auto.rundir) for line in fileinput.input(arp.inp.Get(typ='par').GetFileName(), inplace=True): if line.startswith('set JOB_ID'): line='set JOB_ID = arp\n' sys.stdout.write(line) if arp: if self.GetParam('target') == 'SAD' and not self.GetParam('arp_sad'): parf=arp.inp.Get(typ='par').GetFileName() obj,N=self.inp.GetTargetObjects(self.GetParam('target'), modtyp='pdb', fsftyp='mtz') if obj is None: common.Error('Input data/model for ARP/wARP SAD could not be retrived.') with open(parf,'a') as f: f.write('set phaseref = SADH\n') f.write("set F_1 = '{0}'\n".format(obj['f+'].GetLabel('f'))) f.write("set SIGF_1 = '{0}'\n".format(obj['f+'].GetLabel('sigf'))) f.write("set F_2 = '{0}'\n".format(obj['f-'].GetLabel('f'))) f.write("set SIGF_2 = '{0}'\n".format(obj['f-'].GetLabel('sigf'))) fp,fpp = obj['mod'].atomtypes[obj['mod'].GetAtomType()] if fp is not None and fpp is not None: f.write("set ANOM = ' FORM {0} {1} {2}'\n".format(obj['mod'].GetAtomType(),fp,fpp)) f.write("set HEAVYIN = '{0}'\n".format(obj['mod'].GetFileName('pdb'))) ### this should be disabled and the various heavy atom exceptions may be removed from dummy/refmac #f.write("set heavyin = '{0}'\n".format(obj['mod'].GetFileName('pdb'))) dummy_env = os.path.join( os.getenv('CRANK2'), 'bin', 'dummy' ) inp_parf=os.path.join(os.path.dirname(parf),'input.par') shutil.copy(parf, inp_parf) shutil.copy(obj['mod'].GetFileName('pdb'), obj['mod'].GetFileName('pdb')+'.ref') arp.env['PATH'] = "{0}{2}{1}".format(dummy_env, os.getenv('PATH'), os.pathsep) arp.env['CCP4I_TOP'] = dummy_env arp.env['CBIN'] = dummy_env arp.Run() #with open('arp_warp.log','w') as g: # subprocess.call( [os.path.join(os.getenv('warpbin'),'warp_tracing.sh'),inp_parf,'1'], stdout=g ) # the following two loggraph procedures are basically the same as in comb_phdmmb # would be nice to merge them somehow? def UpdateResults(self, update): # defining variables just for the convenience mb_stat=self.mb.GetProg(supported=True).GetStat ph_stat=self.ph.GetProg(supported=True).GetStat if update=='mb': self.mb_res_all.append( (self.mb.cyc, mb_stat('res_built'), mb_stat('frag_built'), mb_stat('compl_chain'), mb_stat('compl_res')) ) elif update=='ph': self.ref_res_all.append( (ph_stat('rfact'), ph_stat('fom')) ) self.PrintActualLogGraph() def PrintActualLogGraph(self, finished=False): if self.opened_loggraph: self.GetLogGraphHandle().seek(0,0) self.LGInfo(self.GetCCP4Header()) self.LGInfo(self.GetLogGraphPrograms()) if self.ref_res_all: self.LGInfo('\n $TABLE : R-factor and FOM per cycle:') self.LGInfo('$GRAPHS : R-factor and FOM per cycle :A:1,2,3:\n$$') self.LGInfo('Cycle R-factor FOM $$ $$') for cyc,(R,fom) in enumerate(self.ref_res_all): self.LGInfo('{0:3} {1:8.3f} {2:8.3f}'.format(cyc+1,R,fom)) self.LGInfo('$$\n') if self.mb_res_all: self.LGInfo('\n $TABLE : Built per building cycle:') self.LGInfo('$GRAPHS : Residues and fragments per cycle :A:1,2,3:\n$$') self.LGInfo('Cycle Residues Fragments $$ $$') for cyc,res,fr,c1,c2 in self.mb_res_all: self.LGInfo('{0:3} {1:8} {2:8}'.format(cyc,res,fr)) self.LGInfo('$$\n') if finished: self.LGInfo('\n$TEXT:Result: $$ Final result $$') self.LGInfo("{0}\n$$\n".format(self.result_str))