#!/usr/bin/python import os,sys,copy from process import process from program import program import common,data class faest(process): name="FA estimation and/or other substructure detection preparations" short_name="FA estimation" supported_progs=["afro","shelxc","ecalc"] supported_params = {} supported_params['target'] = common.parameter( desc='Experiment/refinement target', typ=str, cap=True ) supported_params['bfactor'] = common.parameter( desc='set B factor of atoms', typ=(float,bool) ) def TreatInOutPar(self, set_all_par=False): # just to prevent a useless warning of parameter not used for shelxd process.TreatInOutPar(self,set_all_par) # guess expected number of substr. atoms if needed and possible obj,N=self.inp.GetTargetObjects(self.GetParam('target')) seq_obj=self.inp.Get('sequence', filetype=data.sequence.supported_filetypes) # I think this could fail if multiple models can be inputted for the same crystal # (in which case this might be fixed eg by looping through all such crystals?) if seq_obj and obj and obj['mod'].exp_num_atoms is None and obj['mod'].GetAtomType() is not None: at=obj['mod'].GetAtomType() if not seq_obj.seqstr: seq_obj.GetSequenceString() if at=='SE': obj['mod'].exp_num_atoms = seq_obj.seqstr.count('M') + seq_obj.seqstr.count('U') elif at=='S': obj['mod'].exp_num_atoms = seq_obj.seqstr.count('M') + seq_obj.seqstr.count('C') else: obj['mod'].exp_num_atoms = 2 if seq_obj.monomers_asym is None: matthews=self.AddProcess('matthews', propagate_out=False) matthews.Run() obj['mod'].exp_num_atoms = obj['mod'].exp_num_atoms * seq_obj.monomers_asym if self.GetProg('ecalc'): # this will be used in the emulation mode so that the GUI gets the defaults for substr. cutoff # if too slow then this can be moved to RunPreprocess() fp_obj=self.inp.Get('fsigf',typ='plus') fm_obj=self.inp.Get('fsigf',typ='minus') fa_obj=self.inp.Get('fsigf',typ='fa',filetype='mtz', inp_cont=self.inp.Get('fsigf',typ='average',filetype='mtz',try_convert=False)) if not fa_obj: fa_obj=self.inp.Get('fsigf',typ='delta-anom',filetype='mtz') if not fa_obj: mtzmadmod=self.AddProg('mtzMADmod',propagate_out=False) mtzmadmod.Run(check_bin=True) fa_obj = mtzmadmod.out.Get(filetype='mtz',typ='delta-anom') self.inp.Add(fa_obj) self.inp.Add( self.ExcludeRefs(fa_obj,fp_obj,fm_obj) ) return if not self.GetProg(supported=True): if program.from_name('shelxc',None).CheckBinary(silent=True): self.AddProg('shelxc') else: self.AddProg('afro') # pass shelxd keywords to shelxc # first create copies of crank, substrdet, shelxd and run treatpar to make sure all required keys are passed # this is not ideal but the very idea of shelxc generating a script for shelxd is far from ideal... if self.GetProg('shelxc') and self.GetCrankParent() and self in self.GetCrankParent().processes: cp = self.GetCrankParent() sd_procs = [sd for sd in cp.GetProcesses('substrdet') if sd in cp.processes[cp.processes.index(self):] and sd.GetProg('shelxd')] if sd_procs: crank = copy.copy(cp) crank.processes = cp.processes[:] substrdet = crank.AddProcessCopy( sd_procs[0], propagate_inp=True ) substrdet.TreatInOutPar() shelxd = substrdet.AddProgCopy( sd_procs[0].GetProg('shelxd'), propagate_inp=True ) shelxd.TreatParams() for key in shelxd.key_list: if not self.GetProg('shelxc').GetKey(key): self.GetProg('shelxc').AddToKey(key, shelxd.GetKey(key,allval=True)) # get&set bfactor for afro if self.GetProg('afro') and self.GetParam('bfactor') in (None,True): gcx=self.AddProg('gcx', propagate_out=False) try: gcx.Run() except (gcx.ProgramRunError,): wilsonB = 50. else: wilsonB = gcx.GetStat('wilson_B') self.programs.remove(gcx) self.SetParam('bfactor', wilsonB) # def RunPreprocess(self,*args,**kwargs): # process.RunPreprocess(self,*args,**kwargs) # if self.GetProg('ecalc'): # fp_obj=self.inp.Get('fsigf',typ='plus') # fm_obj=self.inp.Get('fsigf',typ='minus') # fa_obj=self.inp.Get('fsigf',typ='fa',filetype='mtz', inp_cont=self.inp.Get('fsigf',typ='average',filetype='mtz',try_convert=False)) # if not fa_obj: # fa_obj=self.inp.Get('fsigf',typ='delta-anom',filetype='mtz') # if not fa_obj: # mtzmadmod=self.AddProg('mtzMADmod',propagate_out=False) # mtzmadmod.Run(check_bin=True) # fa_obj = mtzmadmod.out.Get(filetype='mtz',typ='delta-anom') # self.inp.Add(fa_obj) # self.ExcludeRefs(fa_obj,fp_obj,fm_obj) def ExcludeRefs(self,fa_obj,fp_obj,fm_obj): # we assume here that plus/minus are included in the same file as fa_obj (true for mtzMADmod but should be generalized) if fp_obj and fm_obj and fa_obj: sft=self.AddProg('sftools') sft_outmtz = 'sft_adjusted.mtz' sft.runname=sft.name+'_sig' sft.SetKey('read', (fa_obj.GetFileName('mtz'),'col','"'+fa_obj.GetLabel('f')+'"','"'+fa_obj.GetLabel('sigf')+'"')) sft.SetKey('read', (fp_obj.GetFileName('mtz'),'col','"'+fp_obj.GetLabel('f')+'"','"'+fp_obj.GetLabel('sigf')+'"')) sft.SetKey('read', (fm_obj.GetFileName('mtz'),'col','"'+fm_obj.GetLabel('f')+'"','"'+fm_obj.GetLabel('sigf')+'"')) sft.SetKey('calc', ('R','col','_minpm_','=','col','"'+fp_obj.GetLabel('f')+'"','col','"'+fm_obj.GetLabel('f')+'"','min')) sft.SetKey('calc', ('R','col','_delpm_','=','col','"'+fp_obj.GetLabel('f')+'"','col','"'+fm_obj.GetLabel('f')+'"','-')) sft.SetKey('calc', ('R','col','_delpm_','=','col','_delpm_','abs')) sft.SetKey('calc', ('R','col','_delratio_','=','col','_delpm_','col','_minpm_','/')) sft.SetKey('select', ('col','_delratio_','<','1')) sft.SetKey('calc', ('R','col','_sigpratio_','=','col','"'+fp_obj.GetLabel('f')+'"','col','"'+fp_obj.GetLabel('sigf')+'"','/')) sft.SetKey('calc', ('R','col','_sigmratio_','=','col','"'+fm_obj.GetLabel('f')+'"','col','"'+fm_obj.GetLabel('sigf')+'"','/')) # sft.SetKey('select', ('col','_sigpratio_','>','2.5')) # sft.SetKey('select', ('col','_sigmratio_','>','2.5')) sft.SetKey('calc', ('R','col','_sigpmratio_','=','col','"'+fp_obj.GetLabel('sigf')+'"','col','"'+fm_obj.GetLabel('sigf')+'"','/')) sft.SetKey('select', ('col','_sigpmratio_','>','0.25')) sft.SetKey('select', ('col','_sigpmratio_','<','4.')) sft.SetKey('write', sft_outmtz) sft.SetKey('quit\nY') sft.Run() fa_obj=sft.out.AddCopy(fa_obj) sft.out.SetFileToChild(fa_obj,sft_outmtz) self.programs.remove(sft) else: common.Warning('Exclusion of reflections could not be done due to F+/F-/FA not supplied.') return fa_obj def PrintActualLogGraph(self): prog=self.GetProg(supported=True) if prog.nick =='shelxc': resol,dano=prog.GetStat("resol"),prog.GetStat("dano") if resol and dano: reso,dsigd = filter(None,resol[0]),filter(None,dano[0]) if self.opened_loggraph: self.GetLogGraphHandle().seek(0,0) self.LGInfo(self.GetCCP4Header()) self.LGInfo(self.GetLogGraphPrograms()) self.LGInfo('\n $TABLE : Anomalous signal to noise vs resolution:') self.LGInfo('$GRAPHS : Dano/Sigdano vs resolution :A:1,2:\n$$') self.LGInfo('Resolution Dano/Sigdano $$ $$') for i in range(min(len(reso),len(dsigd))): self.LGInfo('{0} {1}'.format(reso[i],dsigd[i])) self.LGInfo('$$\n') if self.rv_report is not None: if not hasattr(self,'rv_plot'): self.rv_plot_dano = self.rv_report.Plot( 'Dano/Sigdano vs resolution', "Resolution", "Dano/Sigdano", block="Anomalous signal vs resolution") self.rv_dano1 = self.rv_plot_dano.PlotLine( ["x",reso], ["Dano/Sigdano",dsigd] ) #self.rv_dano1.Data(reso,dsigd,flush=True) # resolution,signal) def RunPostprocess(self,restore=True,*args,**kwargs): self.PrintActualLogGraph() # lets convert to e-values for shelxc (not needed if shelxd follows but can be used otherwise) if self.GetProg('shelxc'): for fa in self.out.GetAll('fsigf',typ='fa'): fa_obj = self.out.Get('fsigf',typ='fa',filetype='mtz', inp_cont=self.inp.Get('fsigf',typ='average',filetype='mtz',xname=fa.xname,try_convert=False)) if fa_obj: ecalc=self.AddProg('ecalc') ecalc.inp.Add(fa_obj) ecalc.Run() self.programs.remove(ecalc) break else: common.Warning('Calculation of E-values failed.') process.RunPostprocess(self,restore,*args,**kwargs)