import time import re, math import subprocess from phaser.pipeline import proxies, homology from phaser import sculptor from phaser import AlterOriginSymmate import iotbx.mtz from iotbx import pdb from iotbx import reflection_file_converter import logging import traceback import gzip import string import shutil import glob import scitbx import phaser.test_phaser_parsers import BLASTMRutils import os, posixpath from mmtbx.scaling import xtriage from mmtbx.command_line import cif_as_mtz from iotbx.bioinformatics import seq_sequence_parse import urllib2 import BLASTphaserthread from phaser import SimpleFileProperties import SCOPdistiller from iotbx import bioinformatics import threading import itertools import copy from phaser.utils2 import * def openpdb(pdbfname): pdbfname = RecaseFilename2(pdbfname) return pdb.input(file_name= pdbfname) LOG_FILENAME = "" # Exceptions class BLASTrunnerException(Exception): pass class MyFormatter(logging.Formatter): def format(self, record): #compute s according to record.levelno #for example, by setting self._fmt #according to the levelno, then calling #the superclass to do the actual formatting indentlevel = len(traceback.format_stack()) indentlevel = max(0,indentlevel - 10) indents = "" for i in range(0,indentlevel): indents = indents + " " if record.levelno >= logging.WARN: self._fmt = "%(levelname)s: %(message)s" else: self._fmt = indents + "%(message)s" return logging.Formatter.format( self, record ) # myobj = BLASTnMRfunc.BLASTrunner("1tfx",r"E:\Programs\Putty\plink.exe",r"E:\Programs\Putty\pscp.exe", "rdo20",r"E:\mtz_files",r"E:\Users\oeffner\Phenix\cctbx_build\exe\phaser.exe",r"E:\Programs\ClustalW2\clustalw2.exe", "wibble","","-ssh") # myobj = BLASTnMRfunc.BLASTrunner("1tfx","ssh","scp","rdo20","/win/mtz_files","phenix.phaser","clustalw2","ln -s") class BLASTrunner(): "Pickup user specific settings before running BLAST and MR calculations" def __init__(self,targetpdbid, sshcmd, scpcmd, userid, mtzfolder, phaserexe, clustalw, superposeexe, calclabel, passwd, sshflag=None, errorcode = None): self.targetpdbid =targetpdbid self.topfolder = os.path.abspath(targetpdbid) self.sshcmd = sshcmd self.sshflag = sshflag self.scpcmd = scpcmd self.userid = userid self.mtzfolder = mtzfolder self.phaserexe = phaserexe self.nmaxsols = 50 self.fmaxhours = 48 self.superpose = superposeexe self.clustalw = clustalw self.osdepbool = False self.srsPC = "@192.168.102.50" #self.srsPC = "@zeus.private.cimr.cam.ac.uk" self.blastPC = "@ulysses.private.cimr.cam.ac.uk" self.pdbnameset = set([]) self.nuniqueseq = 0 self.errorcodesum = [0] # if supplied the external errorcode must be an integer array with one element or None if errorcode != None: # assign to errorcode pointer self.errorcodesum = errorcode self.writeBLASTxml = True self.SCOPtxt = "" self.superposelogfname = "superposition.log" self.sculptprotocol = "1" self.targetscattering = "" self.busecurrenttmpl = False self.resolutions = "" self.buseIntensities = False self.refinevrmsdefault = False self.obsoletes = "" self.pretext = "" self.useuniquesequence = False self.mydb = None self.passwd = passwd self.calclabel = calclabel if os.sys.platform == "win32": self.osdepbool = True # compile re expressions for GatherResults() # we search for all ncs copies of a partial model so there'll be at most number of # ncs copies Z and LLG scores in the annotation string self.RFZrecmp = re.compile("RFZ=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.TFZrecmp = re.compile("TFZ=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.LLGrecmp = re.compile("LLG=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.PAKrecmp = re.compile("PAK=([-+]?[0-9]*)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.TEMPLLLGrecmp = re.compile("LLG=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.Tstarrecmp = re.compile("\*T=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.TFZrecmp = re.compile("TFZ=([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.TFZequivrecmp = re.compile("TFZ==([-+]?[0-9]*\.?[0-9]+)", re.DOTALL|re.MULTILINE|re.IGNORECASE) self.WILSON_Brecmp = re.compile(r"^ \s+ Wilson \s+ B-factor:\s+ ( -? \d+\.\d* )", re.MULTILINE | re.VERBOSE) self.WILSON_Scalerecmp = re.compile(r"^ \s+ Wilson \s+ Scale:\s+ ( -? \d+\.\d* )", re.MULTILINE | re.VERBOSE) self.MATTHEWCOEFrecmp = re.compile(r"^\s+Most probable VM for resolution = (-?\d+\.\d*)", re.MULTILINE) self.CPUtimerecmp = re.compile(r"\r*\nCPU Time: [^\n]* \(\s*(\d*\.\d*) secs\)\r*\n", re.MULTILINE) self.baseinputfname = self.topfolder + os.path.sep + "Phaserinput" \ + os.path.sep + self.calclabel + "input.txt" self.logformat = MyFormatter() self.logger=logging.getLogger("rdo20") self.formatter=MyFormatter() # only add logging handler if none is present if len(self.logger.handlers)==0: handler1= logging.StreamHandler() handler1.setFormatter(self.formatter) self.logger.addHandler(handler1) self.logger.setLevel(logging.DEBUG) if self.calclabel != "": # we're not just creating an empty database layout self.scoppath = os.path.join(self.topfolder,"..", "..", "SCOP_MR") if os.path.exists(self.scoppath) == False: if "SCOPDIR" in os.environ: self.scoppath = os.environ["SCOPDIR"] self.scopobj = SCOPdistiller.SCOPdistiller(os.path.join(self.scoppath,"dir.des.scop.txt"), os.path.join(self.scoppath,"dir.cla.scop.txt"), os.path.join(self.scoppath,"PDBsXray1chainASUto2.5A.txt"), os.path.join(self.scoppath,"All_PDB.txt")) def MakeMyfolder(self,folderpath): try: if os.path.exists(folderpath) == False: self.logger.log(logging.INFO, "creating folder: " + posixpath.abspath(folderpath)) os.makedirs(folderpath) else: self.logger.log(logging.INFO, "Folder " + posixpath.abspath(folderpath) + " already exists") except Exception, m: self.failure = True self.message = m finally: return def FetchPDBfile(self, subfolder, modelpdbid): pdbid = modelpdbid[:4] # only use the first 4 letters, not the chain ID pdbfname = subfolder + pdbid.lower() + '.pdb' if os.path.exists(pdbfname) == True: self.logger.log(logging.INFO, "%s already present" %pdbfname) return pdbfname fname = 'pdb' + pdbid.lower() + '.ent.gz' self.logger.log(logging.INFO, "Fetching pdbfile %s" %fname) lstcommand = [self.scpcmd, #'%s@zeus.private.cimr.cam.ac.uk:/srs/srs/data/pdbcurrent/pdbcurrent/%s'\ '%s%s:/srs/srs/data/pdbcurrent/pdbcurrent/%s'\ %(self.userid, self.srsPC, fname), '%s.'% subfolder] if self.passwd != "": lstcommand.insert(1,self.passwd) lstcommand.insert(1,"-pw") command = tuple(lstcommand) self.logger.log(logging.INFO, "Executing command: " + subprocess.list2cmdline(command)) subprocess.Popen(command, shell=self.osdepbool, stdout=subprocess.PIPE).communicate()[0] if os.path.exists(subfolder + fname): self.logger.log(logging.INFO, "Unpacking " + subfolder + fname) fileObj = gzip.GzipFile(subfolder + fname, 'rb') pdbtext = fileObj.read() # decompress the file fileObj.close() os.remove(fname) # delete the compressed file else: self.logger.log(logging.INFO, "pdbfile %s is not present. Trying RCSB instead..." %fname) u = urllib2.urlopen("http://www.pdb.org/pdb/files/"+ pdbid + ".pdb") pdbtext = u.read() self.logger.log(logging.INFO, "Writing " + pdbfname) fileout = open(pdbfname,'w') fileout.write(pdbtext) fileout.close() return pdbfname def PDB2SeqLocal(self, pdbchainid, uselocalDB = True): self.logger.log(logging.INFO, "PDB2SeqLocal(%s, %s)" %(pdbchainid, str(uselocalDB))) if uselocalDB == True: # get the sequence form our local database lstcommand = [self.sshcmd, "%s%s" %(self.userid, self.srsPC), "(source /etc/profile; getz", "-f", "seq", "-sf", "fasta", "'[pdbcurrent:%s]')" % pdbchainid[:4]]# first 4 characters only are the pdbid #lstcommand = ["getz", "-f", "seq", "-sf", "fasta", "'[pdb:%s]')" % pdbchainid[:4]]# first 4 characters only are the pdbid if self.passwd != "": lstcommand.insert(1,self.passwd) lstcommand.insert(1,"-pw") if self.sshflag != None: lstcommand.insert(1,self.sshflag) command = tuple(lstcommand) self.logger.log(logging.INFO, "Executing command: " + subprocess.list2cmdline(command)) fastastr=subprocess.Popen(command,shell=self.osdepbool, stdout=subprocess.PIPE).communicate()[0] if fastastr.find(">") < 0: # our local pdb server doesn't have this sequence return self.PDB2SeqLocal(pdbchainid, False) t = re.findall('>[^>]*', re.split('Welcome to SRS 8.1.1 -------------',fastastr)[-1]) else: # get the sequence from the RCSB u = urllib2.urlopen("http://www.pdb.org/pdb/files/fasta.txt?structureIdList="+ pdbchainid[:4]) fastastr = u.read() t = re.findall('>[^>]*', fastastr) if len(pdbchainid) != 4 and len(pdbchainid) != 6 and len(pdbchainid) != 11: raise BLASTrunnerException, "Is %s supposed to be a pdbid?!" % pdbchainid sequence = "" if len(pdbchainid) == 4: # we only supplied the pdbid, not the chainid sequence = t # sequence now holds an array of sequences for each chain else: if len(pdbchainid) == 6: # chainid in pdbchainid is specified as 1BL8_D chainid = pdbchainid[5].upper() if len(pdbchainid) == 11: # chainid in pdbchainid is specified as 1BL8_ChainD chainid = pdbchainid[10].upper() for i in range(0, len(t)): if uselocalDB == True: # chain id in returned sequences from our local database is the 12 character in sequence preamble # like: # >1DDL_ChainC # MEQDKILAHQASLNTKPSLLPPPVGNPPPVISYPFQITLASLGTEDAADSVSIASNSVLA if t[i][11].upper() == chainid.upper(): sequence = [t[i],] # we expect an array to be returned, so make one else: # chain id in returned sequences from RCSB is the 6 character in sequence preamble # like: # >1TFX:C|PDBID|CHAIN|SEQUENCE # KPDFCFLEEDPGICRGYITRYFYNNQTKQCERFKYGGCLGNMNNFETLEECKNICEDG if t[i][6].upper() == chainid.upper(): sequence = [t[i],] # we expect an array to be returned, so make one if sequence == "": raise BLASTrunnerException, "Could not find sequence for chain %s in pdb %s" \ %(chainid, pdbchainid[:4]) else: for chainseq in sequence: rawseq = seq_sequence_parse(chainseq)[0][0].sequence m = re.findall(r"[xzbju]", rawseq,re.MULTILINE|re.IGNORECASE) if len(m) > 0: if uselocalDB == True: self.logger.log(logging.WARN,"Found illegal letters %s in the sequence. Trying RCSB instead" %m) newseq = self.PDB2SeqLocal(pdbchainid, False) newsequence = [] for newchainseq in newseq: seq = seq_sequence_parse(newchainseq)[0][0].sequence name = seq_sequence_parse(newchainseq)[0][0].name chainlbl = name[5] newsequence.append(">" + name[:4] + "_Chain" + chainlbl + "\n" + seq + "\n") return newsequence else: self.errorcodesum[0] += 1 raise BLASTrunnerException, "Illegal letter, %s, in the sequence\n%s\nfor %s" \ %(m, chainseq, pdbchainid) return sequence def GetLocalBLASTHomologues(self, maxpdbs, uniquetargetchains): sequence = self.PDB2SeqLocal(self.targetpdbid) # When chains are repeated in the pdb file we use a unique list of sequences. # To distinguish identical chains labbeled with different letters from one another # we strip off the title labels untitledseq = [] for chainseq in sequence: print "chainseq= " + str(chainseq) rawseq = seq_sequence_parse(chainseq)[0][0].sequence # Sequence comment has chain label as the last non-whitespace letter. grab it. # 1DDL_ChainC # MEQDKILAHQASLNTKPSLLPPPVGNPPPVISYPFQITLASLGTEDAADSVSIASNSVLA preamble = seq_sequence_parse(chainseq)[0][0].name if preamble.find("|PDBID|CHAIN|SEQUENCE") > -1: # new convention for preamble chainlabel = seq_sequence_parse(chainseq)[0][0].name[5] else: # still using old convetion as in ">2D5B_ChainA" chainlabel = seq_sequence_parse(chainseq)[0][0].name.rstrip()[-1] untitledseq.append([chainlabel, rawseq]) if self.useuniquesequence: # build the unique list uniqueseq = [] uniqueseq.append(["",""]) # dummy element forces definition of this 2dimensional list for i in range(0,len(untitledseq)): # this will use the first chainid as the key for this dictionary and # will thus overwrite entries with different chainids foundit = False for j in range(0,len(uniqueseq)): if untitledseq[i][1] == uniqueseq[j][1]: foundit = True uniqueseq[j][0] = uniqueseq[j][0] + untitledseq[i][0] if not foundit: uniqueseq.append(untitledseq[i]) uniqueseq.remove(["",""]) # get rid of dummy element self.logger.log(logging.INFO, "Sequence list for %s is:" %self.targetpdbid) self.logger.log(logging.INFO, sequence) self.logger.log(logging.INFO, "Omitting duplicate sequences we have:") else: uniqueseq = untitledseq self.logger.log(logging.INFO, uniqueseq) self.logger.log(logging.INFO,"") # pass the uniqueseq on to the external variable uniquetargetchains.extend(uniqueseq) #uniquetargetchains.remove([])#remove the first empty dummy self.nuniqueseq = len(uniqueseq) # store for use in GetCombinedSearchesforTargetChains() myids = [] for chainid, seq in uniqueseq: # loop over the list of [chainid,sequence] if self.writeBLASTxml: # for each unique sequence in the pdb supplied let's find some homologues lstcommand = [self.sshcmd, \ self.userid + self.srsPC, "/srs/blast/blast-2.2.20/bin/blastall"] #lstcommand = ["/srs/blast/blast-2.2.20/bin/blastall"] if self.passwd != "": lstcommand.insert(1,self.passwd) lstcommand.insert(1,"-pw") if self.sshflag != None: lstcommand.insert(1,self.sshflag) command = tuple(lstcommand) self.logger.log(logging.INFO, "Executing command: " + subprocess.list2cmdline(command)) myproxy = proxies.NCBIBlastStandalone(command, self.osdepbool) blastcmdstr = [""] myproxy.submit(program="blastp",database="/srs/srs/data/blast/pdbaa", \ query=seq, maxresults=maxpdbs, expectation = 1000, pcmdstr = blastcmdstr) self.logger.log(logging.INFO, blastcmdstr[0]) # get the BLAST results as xml data xmldata = myproxy.results(jobid=0) n = xmldata.find(" 0: # strip off garbage such as "Welcome to the SRS" xmldata = xmldata[n:] fname = "%s_%s_blast.xml" %(self.targetpdbid, chainid) fileout = open(fname,"w") fileout.write(xmldata) fileout.close() else: fname = "%s_%s_blast.xml" %(self.targetpdbid, chainid) xmldata = open(fname,"r").read() # parse the xml into python structures myparser = homology.get_ncbi_blast_parser() results = myparser(xmldata) seqlength = results.root.query_len length = len(results.root.iterations[0].hits) self.logger.log(logging.INFO,"Chain %s yields %d BLAST hits\n" %(chainid,length)) for hit in results.root.iterations[0].hits: # label each homologue with sequence identity and the chain id of the unique chain identity = hit.hsps[0].identity align_len = hit.hsps[0].length hit_accession = hit.accession hit_length = hit.length seqid = "%4.3f" %(float(identity*100)/min(hit_length, seqlength)) # seqid = float(identity*100)/seqlength # seqid = float(identity*100)/align_len myids.append([hit_accession, [seqid], chainid]) return myids def GetPDBidsFromPDBLocalWithSequenceIds(self, maxseqid, minseqid, maxpdbs, uniquechains = None): if (maxseqid < minseqid or maxseqid < 0 or minseqid <0 or maxseqid > 100 or minseqid > 100): raise BLASTrunnerException, "Check values of maxseqid and/or minseqid" self.logger.log(logging.EXCEPTION,"Check values of maxseqid and/or minseqid") exit(42) # initialise in case of uniquechains not supplied if uniquechains == None: uniquechains = [] myids = self.GetLocalBLASTHomologues(maxpdbs, uniquechains) selectedpdbs = [] self.logger.log(logging.INFO,u"BLAST sequence identities of models are:") for i in range(0,len(myids)): resol = self.GetResolution(myids[i][0]) if (float(myids[i][1][0]) >= minseqid and float(myids[i][1][0]) <= maxseqid): if float(resol) == 0.0 : self.logger.log(logging.WARN, "The " + myids[i][0][:4] + " entry isn't present in the resolution.txt file. Grab the latest from RCSB!") continue if float(resol) < 0.0 : self.logger.log(logging.INFO, myids[i][0][:4] + " at " + resol \ + u" Å was not determined from crystallography and is ignored.") continue self.logger.log(logging.INFO, myids[i][0][:4] + " at " + resol \ + u" Å, seqid: %s target chain: %s is considered" %(myids[i][1][0], myids[i][2])) selectedpdbs.append(myids[i]) else: self.logger.log(logging.INFO, "%s seqid: %s, target chain: %s is binned" \ %(myids[i][0], myids[i][1][0], myids[i][2])) return selectedpdbs def RemovePDBsWithSimilarSeqIdentFromList(self, selectedpdbs, maxseqid, minseqid, myinc): # we're not interested in many pdbs having almost the same seqid so compose # a list of pdbs evenly distributed over seqids # Do this for each of the chains bins = self.PutPDBsInChainIdBins(selectedpdbs) self.logger.log(logging.INFO,"Purging models with similar sequence identity") evenlyspread = [] #inc = 1 # seqid will increase by at least 1% through each element inc = max(1.0, myinc) for n in range(0,len(bins)): imin = minseqid imax = minseqid + inc while imax <= maxseqid: i=0 found = False while found==False and i< len(bins[n]): # no more than one model between imin and imax if float(bins[n][i][1][0]) >= imin and float(bins[n][i][1][0]) 0: n = self.obsoletes.find("\n", m) # got the line now with the obsolete and the new pdb id match = re.search(r"(" + pdbid[:4].upper() + ")\s+([A-Za-z0-9]+)", self.obsoletes[m:n], re.MULTILINE) if match != None and len(match.groups()) > 1: validpdb = match.group(2) if self.resolutions == "": self.resolutions = open(self.scoppath + os.path.sep + "resolution.txt","r").read() # The resolution file looks like below: # # 1XC4 ; 2.8 # 1XC5 ; -1.00 # 1XC6 ; 2.1 # 1XC7 ; 1.83 n = self.resolutions.find(validpdb, 0) if n<0: self.logger.log(logging.WARN,"No resolution found in resolution.txt for %s." %pdbid) return resolution m = self.resolutions.find("\n",n) # we found the line with the resolution now. So extract the resolution now # if resolution is -1.0 then it's NMR data match = re.search(r'('+ validpdb +'\s;\s)([-+]?[0-9]*\.?[0-9]+)', \ self.resolutions[n:m], re.MULTILINE) resolution = "%3.2f" % float(match.group(2)) except Exception, e: self.errorcodesum[0] += 30 self.logger.log(logging.WARN,"No resolution found in %s.pdb" %pdbid) raise BLASTrunnerException, e return resolution def sculptprot(self): return "sculpt" + self.sculptprotocol + "_" if self.sculptprotocol != "" else "" def DoQuickRefinementOfSuperposition(self, searchmodels): # prepare for a quick refinement of superposed models. The output pdb should then # be a solution that can be used as a template sortedmodels = sorted(searchmodels) # to easily locate an already calculated template sortednames = [] for sortedmodel in sortedmodels: sortednames.append( sortedmodel[0]) sortedmodelname = ",".join(sortednames) modelname = [] for model in searchmodels: modelname.append( model[0]) modelsname = ",".join(modelname) # by comparing with other sorted folder names we avoid many redundant calculations # given that template solutions {pdb_a,pdb_b,pdb_c} are the same as {pdb_b, pdb_a,pdb_c} # and similar for all permutations templfolder = os.path.join(self.topfolder, "TemplateSolutions",sortedmodelname) RNProotname = posixpath.join("TemplateSolutions",sortedmodelname, self.calclabel + "RNP{" + modelsname + "}.superposed") RNPsolfname = RNProotname + ".sol" if os.path.exists(RNPsolfname): self.logger.log(logging.INFO, "Not making template solution %s as it was already created as %s" %(modelsname,sortedmodelname)) return True self.MakeMyfolder(templfolder) self.logger.log(logging.INFO, "Refining superposed search model(s) %s on target structure" %modelsname) rnpfnameprefix = templfolder + os.path.sep + self.calclabel + sortedmodelname + "_quickRNPinput" RNPinputfile = open(rnpfnameprefix + ".txt", "w") RNPinputfile.write("MODE MR_RNP\nMACMR ROT ON TRA ON BFAC OFF VRMS ON\n") solstr = "" try: # get models and use their orientations and positions from the superpose output as solutions for model in sortedmodels: targetpdbchainid = model[2] seqid = model[1][0] ensname = model[3][0][0] pdbfname = self.GetSculptFname(targetpdbchainid, model[0]) pdbfname = posixpath.relpath(pdbfname, self.topfolder) RNPinputfile.write("ENSEMBLE %s PDB %s IDENTITY %s\n" %(ensname, pdbfname, seqid)) for ncscopy in model[3]: euler = ncscopy[1] ortho = ncscopy[2] solstr = solstr + "SOLU 6DIM ENSE " + ensname + " EULER %2.3f %2.3f %2.3f " \ % euler + " ORTH %2.3f %2.3f %2.3f BFAC 0.00\n" %ortho RNPinputfile.write("\nSOLU SET Superposed\n" + solstr) RNPinputfile.write("ROOT " + RNProotname + "\n") # this file is concatenated to the one 3 folders above before submitted to phaser RNPinputfile.close() myphaser = BLASTphaserthread.PhaserRunner(self.baseinputfname, rnpfnameprefix, self.calclabel + "RNP", "", self) # run Phaser. Shouldn't take long if models superposed well on the target self.logger.log(logging.INFO, "Now refining superposed model") stdoutput = myphaser.PhaserRun(templfolder) self.errorcodesum[0] += myphaser.errorcodesum phaserlogfile = open(RNProotname + ".log","w") phaserlogfile.write(stdoutput) phaserlogfile.close() # extract eulers and frac from the solution file. These are then the proposed template solution RNPsolfile = open(RNPsolfname,"r") RNPsol = RNPsolfile.read() RNPsolfile.close() n = RNPsol.find("SOLU 6DIM ENSE",0) solutions = RNPsol[n:] tmplfname = templfolder + posixpath.sep + self.calclabel + "template{" + sortedmodelname + "}.txt" templatefile = open(tmplfname,"w") templatefile.write("SOLU TEMPLATE\n" + solutions) templatefile.close() self.logger.log(logging.INFO, "Eulers and fractionals from refinement of superposition: \n" + solutions) except Exception, m: #raise self.errorcodesum[0] += 1 self.logger.log(logging.ERROR, "Problem refining " + modelsname + "\n" + str(m)) return False olddir = os.getcwd() os.chdir(templfolder) # remove superfluous output mtzs = glob.glob("*.mtz") #delf = mtzs + dbgs delf = mtzs for i in range(0,len(delf)): os.remove(delf[i]) os.chdir(olddir) self.logger.log(logging.INFO,"Done DoQuickRefinementOfSuperposition\n") return True def PrepareTemplateSolutions(self, searchmodels, permutations): model = [] targetpdbchainid = [] seqid = [] for pdb_id in searchmodels: model.append( pdb_id[0]) targetpdbchainid.append( pdb_id[2]) seqid.append( pdb_id[1]) self.logger.log(logging.INFO, "\nPrepareTemplateSolutions(%s)" % ",".join(model)) # successively build up template solutions with more and more models for permutedmodels in permutations: lmodel = [] for lmodels in permutedmodels: lmodel.append(lmodels) if self.DoQuickRefinementOfSuperposition(lmodel) == False: self.logger.log(logging.WARN, "\nPrepareTemplateSolutions(%s) failed.\n" \ "No MR calculations with this combined model will be carried out." % ",".join(model)) self.errorcodesum[0] += 1 return False return True def GetEuler_and_Translation_of_Model_to_Target(self, modelpdbfname, targetpdbfname, targetpdbchainid = ""): self.logger.log(logging.DEBUG, "GetEuler_and_Translation_of_Model_to_Target(%s, %s, %s)" \ %(modelpdbfname, targetpdbfname, targetpdbchainid)) if os.path.exists(modelpdbfname) == False: raise BLASTrunnerException, os.path.abspath(modelpdbfname) + " doesn't exist" if os.path.exists(targetpdbfname) == False: raise BLASTrunnerException, os.path.abspath(targetpdbfname) + " doesn't exist" command = (self.superpose, modelpdbfname , targetpdbfname, "-s", targetpdbchainid, "-o", modelpdbfname[:-4] + '_superposed_%s.pdb' %targetpdbchainid) if targetpdbchainid=="": # target chain id not specified command = (self.superpose, posixpath.relpath(modelpdbfname), posixpath.relpath(targetpdbfname), posixpath.relpath(modelpdbfname[:-4] + '_superposed.pdb')) self.logger.log(logging.INFO, "Executing command:") self.logger.log(logging.INFO, subprocess.list2cmdline(command)) output = subprocess.list2cmdline(command) + "\n" \ + subprocess.Popen(command, shell=self.osdepbool, stdout=subprocess.PIPE).communicate()[0] print "command= " + str(command) self.superposelogfname = self.calclabel + "superposition.log" soutput = open(self.superposelogfname,"w") # create a logfile for debugging soutput.write(output) # dump cmdline and superpose output into logfile soutput.close() # Output from superpose looks like below. We extract euler angles and translation vector # # CCP4 format rotation-translation operator # Polar angles (omega,phi,kappa) : 40.601 -90.131 178.534 # Euler angles (alpha,beta,gamma): -91.096 81.193 -90.834 # Orthogonal translation (/Angst): 31.247 52.164 5.947 # #n = output.find("Euler angles (alpha,beta,gamma):",0) #m = output.find("\n",n) try: #match = re.search(r'(Euler angles \(alpha,beta,gamma\):)(\s*)([\w\W]*)', \ # output[n:m], re.MULTILINE) #print output[1400:2000] match = re.search(r"Rx\s+Ry\s+Rz\s+T\r*\n" + r"\s+(-?\d+\.\d+)"*12, output) strmatrix= [] strmatrix.extend(match.groups()[0:3]) strmatrix.extend(match.groups()[4:7]) strmatrix.extend(match.groups()[8:11]) strtrans = (match.groups()[3], match.groups()[7], match.groups()[11]) rotmatrix = [float(s) for s in strmatrix] orthotrans = [float(s) for s in strtrans] eul = scitbx.math.euler_angles_zyz_angles(rotmatrix) eulerangles = [Roundoff(float(s),3) for s in eul] #stringnumbers = tuple(match.group(3).split()) #eulerangles = [float(stringnumbers[0]), float(stringnumbers[1]), float(stringnumbers[2])] self.logger.log(logging.INFO,"model-target Euler angles are: " + str(eulerangles)) #n = output.find("Orthogonal translation (/Angst):",0) #m = output.find("\n",n) #match = re.search(r'(Orthogonal translation \(/Angst\):)(\s*)([\w\W]*)', \ # output[n:m], re.MULTILINE) #stringnumbers = match.group(3).split() # must have floats and not strings as args for uc.fractionalize() #orthotrans = [float(stringnumbers[0]), float(stringnumbers[1]), float(stringnumbers[2])] self.logger.log(logging.INFO,"model-target orthogonal translation is: " + str(orthotrans)) except Exception, m: self.logger.log(logging.ERROR,"Problem with superpose:\n%s" %output) eulerangles = [0.0, 0.0, 0.0] orthotrans = [0.0, 0.0, 0.0] self.errorcodesum[0] += 1 #raise #BLASTrunnerException, m return (tuple(eulerangles), tuple(orthotrans)) def TrimModel(self, targetpdbchainid, pdbs_seqid, residuenumbers): self.logger.log(logging.DEBUG, "TrimModel(%s, %s)" %(targetpdbchainid, pdbs_seqid[0])) # first get the sequence from the target chain we already got targetseqfile = open(posixpath.join(self.topfolder,"Phaserinput", targetpdbchainid[:4]+"_Chain" + targetpdbchainid[5] + ".seq"),"r") targseq = targetseqfile.read() targetseqfile.close() try: # then get sequence of the model chain if os.path.exists("clustinput.txt") == False: sequences = targseq + "\n" + self.PDB2SeqLocal(pdbs_seqid[0])[0] + "\n" fileout = open("clustinput.txt","w") fileout.write(sequences) fileout.close() else: self.logger.log(logging.INFO, "clustinput.txt already present") if os.path.exists("clust.aln") == False: command = (self.clustalw, '-INFILE=clustinput.txt' ,'-ALIGN',\ '-OUTFILE=clust.aln') self.logger.log(logging.INFO, "Executing command:" + subprocess.list2cmdline(command)) subprocess.Popen(command,shell=self.osdepbool, stdout=subprocess.PIPE).communicate()[0] else: self.logger.log(logging.INFO, "clust.aln already present") # get rid of ncs copies or other chains in the model, i.e. preserve only one chain modelfname = pdbs_seqid[0][:4] + ".pdb" sculptarg= ["output.root = tmp" + self.calclabel, "renumber.use = original", \ "input.model.file_name = %s" %modelfname, "input.alignment.file_name = clust.aln", "input.model.selection = chain '%s'" %pdbs_seqid[0][5], "macromolecule.protocols=%d" %int(self.sculptprotocol)] self.logger.log(logging.INFO, "Sculptor.run: %s\n" % subprocess.list2cmdline(sculptarg)) # create file handler which logs even debug messages fh = open(self.calclabel + "sculptor.log","w") from phaser import tbx_utils factory = tbx_utils.PhilArgumentFactory( master_phil = sculptor.PHIL_MULTIPROTOCOL_MASTER ) sculptor.run( [ factory( arg ) for arg in sculptarg ], fh, setup = sculptor.SETUP_MULTIPROTOCOL,) #sculptor.run(sculptarg, fh) # find first and last residue numbers for tagging on to the file name tmpfname = "tmp" + self.calclabel + "_" + pdbs_seqid[0][:4] + "_" + self.sculptprotocol + ".pdb" mypdb = openpdb(tmpfname) root = mypdb.construct_hierarchy() firstrsd = root.models()[0].chains()[0].residue_groups()[0].resseq_as_int() lastrsd = root.models()[0].chains()[0].residue_groups()[-1].resseq_as_int() residuenumbers[0] = "" # To avoid ENSEMBLE name clashes later we note the number of times this model has been used for # the rare but possible case when BLAST suggest a particular model more than once for our target. # If used more than once then tag the number of times this model is occuring to the modelpdbid # Don't do this if the model occurs only once as it makes file and folder names difficult to read sn = 0 while (pdbs_seqid[0] + str(sn)) in self.pdbnameset: sn += 1 if sn > 0: pdbs_seqid[0] = pdbs_seqid[0] + str(sn) residuenumbers[0] = "{%d-%d}" %(firstrsd,lastrsd) # watch out if a model is suggested more than once self.pdbnameset.add(pdbs_seqid[0]) #sculptfname = "sculpt%s_%s.pdb" %(self.sculptprotocol, pdbs_seqid[0]) sculptfname = os.path.basename(self.GetSculptFname(targetpdbchainid, pdbs_seqid[0])) if os.path.exists(sculptfname): os.remove(sculptfname) os.rename(tmpfname, sculptfname) except Exception, m: #raise self.failure = True self.message = m self.logger.log(logging.ERROR, self.message) self.errorcodesum[0] += 1 return "" else: return sculptfname def PreparePartialModel(self, pdbs_seqid): # create partial modeldir and store sculpted model there self.logger.log(logging.INFO, "\nPreparePartialModel(%s)" % pdbs_seqid) targetpdbchainid = pdbs_seqid[2] modeldir = self.targetpdbid + "_" \ + targetpdbchainid + posixpath.sep + pdbs_seqid[0] + posixpath.sep self.MakeMyfolder(modeldir) olddir = os.getcwd() os.chdir(modeldir) # store the model and sculpt files in the modeldir folder modelfname = self.FetchPDBfile("." + posixpath.sep, pdbs_seqid[0]) residuenumbers = [""] sculptedpdb = self.TrimModel(self.targetpdbid + "_" + targetpdbchainid[0], pdbs_seqid, residuenumbers) if not os.path.exists(sculptedpdb): self.logger.log(logging.WARN, "There's an issue with " + pdbs_seqid[0] + ". No MR calculation will be done with that model.") os.chdir(olddir) return False # insert model completeness into pdbs_seqid list # replace the BLAST sequence identity with the ClustalW2 as it is more "true" myalignment = open(os.path.join("clust.aln"),"r").read() clustparser = bioinformatics.clustal_alignment_parse clustaln = clustparser(myalignment)[0] pdbs_seqid[1][0] = "%4.3f" %(clustaln.identity_fraction()*100.0) sculptedpdb = RecaseFilename2(sculptedpdb) modelscattering, mstr = SimpleFileProperties.GetScatteringOfModelInPDBfile(sculptedpdb) targetpdbfname = os.path.abspath(self.topfolder) + os.path.sep + self.targetpdbid + ".pdb" self.logger.log(logging.INFO, mstr) modelcompleteness = "%4.3f" %(modelscattering/self.targetscattering) pdbs_seqid[1].append(modelcompleteness) # store SSM seqid as well for comparison with clustalw seqid SSMSeqIDs = SimpleFileProperties.GetSSMSequenceIdentity(targetpdbfname, sculptedpdb, targetpdbchainid[0]) pdbs_seqid[1].append(Roundoff(SSMSeqIDs[0][0], 4)) clustalalignlen = clustaln.shortest_seqlen() pdbs_seqid[1].append(clustalalignlen) ssmalignlen = SSMSeqIDs[0][6] pdbs_seqid[1].append(int(ssmalignlen)) euler_ortho = [] # account for NCS in the target structure by making same number of model superpositions samechains = len(targetpdbchainid) superposition = [] ensemblename = "MR_" + pdbs_seqid[0] for i in range(0, samechains): # do a superposition of model onto target for creating a template solution eulerorthos = self.GetEuler_and_Translation_of_Model_to_Target(sculptedpdb, targetpdbfname, targetpdbchainid[i]) # Did superpose fail? if eulerorthos[0][0] == 0.0 and eulerorthos[0][1] == 0.0 and eulerorthos[0][2] == 0.0: self.logger.log(logging.WARN, "Superpose has an issue with " + pdbs_seqid[0] + ". No MR calculation will be done with that model.") os.chdir(olddir) return False euler_ortho.append(eulerorthos) superposition.append((ensemblename, euler_ortho[i][0], euler_ortho[i][1])) pdbs_seqid.append(superposition) pdbs_seqid.append(residuenumbers[0]) os.chdir(olddir) return True def GetSculptFname(self, targetchainid, modelpdbidchainid): # locate the trimmed pdb model prefix = self.calclabel + "sculpt_" sculptedpdb = posixpath.join(self.topfolder, self.targetpdbid \ + "_" + targetchainid, modelpdbidchainid[:6], prefix + modelpdbidchainid + "_singlechain.pdb") if os.path.exists(RecaseFilename2(sculptedpdb)) == False: sculptedpdb = posixpath.join(self.topfolder, self.targetpdbid \ + "_" + targetchainid, modelpdbidchainid[:6], prefix + modelpdbidchainid + ".pdb") # and try without calclabel for backwards compatibility prefix = "sculpt_" if os.path.exists(RecaseFilename2(sculptedpdb)) == False: sculptedpdb = posixpath.join(self.topfolder, self.targetpdbid \ + "_" + targetchainid, modelpdbidchainid[:6], prefix + modelpdbidchainid + "_singlechain.pdb") if os.path.exists(RecaseFilename2(sculptedpdb)) == False: sculptedpdb = posixpath.join(self.topfolder, self.targetpdbid \ + "_" + targetchainid, modelpdbidchainid[:6], prefix + modelpdbidchainid + ".pdb") return RecaseFilename2(sculptedpdb) def PrepareInputFiles(self, searchmodels): #import code; code.interact(local=locals()) workdir = [] for pdb_id in searchmodels: workdir.append( pdb_id[0] ) modelsdir = ",".join( workdir ) self.logger.log(logging.INFO, "Preparing MR AUTO for models %s " %modelsdir) mdlfolder = self.topfolder + posixpath.sep + "MR" + posixpath.sep + modelsdir self.MakeMyfolder(mdlfolder) templatestr = "" partstr = "" ensemblestr = "" searchstr = "" models = [] for search in searchmodels: ensemblename = "MR_" + search[0] # locate the trimmed pdb model sculptedpdb = self.GetSculptFname(search[2], search[0]) if os.path.exists(sculptedpdb) == False: self.logger.log(logging.WARN, "Couldn't find %s " %sculptedpdb) ensemblestr = ensemblestr + "ENSEMBLE " + ensemblename + " PDB " \ + posixpath.relpath(sculptedpdb, self.topfolder) + " IDENTITY " + search[1][0] + "\n" # find calculated template solutions by sorting order of pdbids and then adding combined # solutions to each inputfile as in {pdb_a},{pdb_a,pdb_b},{pdb_a,pdb_b,pdb_c},... templmodels = [] modelinputs = [] for search in searchmodels: # get number of domains to search for with this chain ensemblename = "MR_" + search[0] samechains = len(search[2]) searchstr = "SEARCH ENSEMBLE %s NUM %d\n" %(ensemblename, samechains) templmodels.append(search[0]) templates = [] # templatemodels are sorted so we can look them up sortedtemplmodels = sorted(templmodels) for pdb_id in sortedtemplmodels: templates.append(pdb_id) templatename = ",".join(templates) fulltemplatedir = posixpath.join(self.topfolder, "TemplateSolutions", templatename) templatefname = posixpath.join(fulltemplatedir, self.calclabel + "template{" + templatename + "}.txt") templatestr = "" if os.path.exists(templatefname) == True: templatestr = open(templatefname,"r").read() else: self.logger.log(logging.WARN, "Couldn't find %s " %templatefname) continue self.logger.log(logging.INFO, "Using template solution: %s" %templatefname) # using partial solutions to seed the next search. These are the previous or current # template solutions excluding the last model we've added to the search partstr = "" if len(templmodels[:-1]) > 0: # find the solution file for the template solution excluding the one we've added to this search # We need to sort template model order to retrieve correct file name if self.busecurrenttmpl: #print "templatestr= " + templatestr #import code; code.interact(local=locals()) partstr = "SOLU SET\n" for templmodel in templmodels[:-1]: ensname = "MR_" + templmodel coordstrs = re.findall("SOLU 6DIM ENSE.+" + ensname + "\s+.+", templatestr, re.MULTILINE) soluensstrs = re.findall("SOLU ENSE.+" + ensname + "\s+.+", templatestr, re.MULTILINE) partstr += coordstrs[0] + "\n" + soluensstrs[0] + "\n" else: partsolname = ",".join(sorted(templmodels[:-1])) fullpartsoldir = posixpath.join(self.topfolder, "TemplateSolutions", partsolname) partialfname = posixpath.join(fullpartsoldir, self.calclabel + "template{" + partsolname + "}.txt") if os.path.exists(partialfname) == True: tmpstr = open(partialfname,"r").read() partstr = tmpstr.replace("SOLU TEMPLATE","SOLU SET") else: self.logger.log(logging.WARN, "Couldn't find %s " %partialfname) continue self.logger.log(logging.INFO, "Using partial solution: %s" %partialfname) modelsname = ",".join(templmodels) modestr = """ ELLG USE OFF #MACMR ROT ON TRA ON BFAC OFF VRMS ON RESOLUTION HIGH 2.5 RESOLUTION AUTO HIGH 0.0 MODE MR_AUTO SEARCH ORDER AUTO OFF \n""" modelfname = posixpath.join(mdlfolder, self.calclabel + "modelinput{" + modelsname + "}.txt") modelinputfile = open(modelfname,"w") modelinputfile.write(modestr + ensemblestr + searchstr + partstr \ + templatestr + "\n") modelinputfile.close() self.logger.log(logging.INFO, modelfname + " :\n" + modestr + ensemblestr + searchstr + partstr + templatestr) modelinputs.append(modelsname) return modelinputs def SetupLogfile(self, logfname): if len(self.logger.handlers) <=1: handler= logging.FileHandler(logfname,encoding = "UTF-8") handler.setFormatter(self.formatter) self.logger.addHandler(handler) self.logger.setLevel(logging.DEBUG) self.logger.log(logging.INFO, self.pretext) def erf(self, x): # Error function for use below in GatherResults a1 = 0.254829592 a2 = -0.284496736 a3 = 1.421413741 a4 = -1.453152027 a5 = 1.061405429 p = 0.3275911 # Save the sign of x sign = 1 if x < 0: sign = -1 x = abs(x) # A & S 7.1.26 t = 1.0/(1.0 + p*x) y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) return sign*y # calculate fracvar(rmsd, resolution, modelcompleteness) def FracVar(self, rms, dres, modelcompleteness): var1 = 2.0*math.pi*rms/(math.sqrt(3.0)*dres) fracvar = (math.sqrt(3.0)*dres*self.erf(var1)-4.0*math.exp(-var1*var1) *math.sqrt(math.pi)*rms)*9.0*dres*dres*float(modelcompleteness)/( 32.0*math.pow(math.pi,2.5)*rms*rms*rms) return fracvar def GatherResults(self, search, workdir, partmodels, nadditionalmodel, mrsol, suffix, datares, Unitcell, spacegroup, spacegroupnumber, crystalsystemname, ispolarspacegroup, nsolvedit, PDBTagsets, WilsonB, WilsonScale, Rfactors, EDSentry, MatthewsCoef, TargetScattering, ntargetresidues, Allreflections, ASUvolume, lastMRTargetTbl_id): # note which is the last element in the database so UPDATE matches a corresponding INSERT statement #print "mrsol0=" +str(mrsol[0]) if mrsol[2] == None or mrsol[3] == None: return self.mydb.execute("SELECT count( * ) as total_record FROM MRModelTbl") lastMRModel_id = self.mydb.fetchone()[0]+1 self.mydb.execute("SELECT count( * ) as total_record FROM AnnotationTbl") lastAnnotation_id = self.mydb.fetchone()[0] # partmodels is an element like '1IK6_A0,3EXF_A0' present in the list ['1IK6_A0', '1IK6_A0,3EXF_A0'] that is present in the list: # [[['1IK6_A0', ['42.398', '0.289'], 'B', '{1-326}'], ['3EXF_A0', ['24.607', '0.352'], 'A', '{0-361}']], ['1IK6_A0', '1IK6_A0,3EXF_A0']] # the current search model is the last item in partmodellist, i.e. '3EXF_A0' # corresponding seqid and model completeness is in ['3EXF_A0', ['24.607', '0.352'], 'A', '{0-361}'] # Get that by using counts of commas in '1IK6_A0,3EXF_A0' as index for the list: # [['1IK6_A0', ['42.398', '0.289'], 'B', '{1-326}'], ['3EXF_A0', ['24.607', '0.352'], 'A', '{0-361}']] partmodellist = partmodels.split(",") i = len(partmodellist)-1 searchmodel = partmodellist[i] if suffix == "": self.logger.log(logging.INFO,"Databasing results for search %s in %s. %d targets, %d models, %d solutions." %(partmodels, workdir, lastMRTargetTbl_id[0], lastMRModel_id, lastAnnotation_id)) else: self.logger.log(logging.INFO,"Databasing %s results for the search of %s in %s. %d targets, %d models, %d solutions." %(suffix, partmodels, workdir, lastMRTargetTbl_id[0], lastMRModel_id, lastAnnotation_id)) resolution = self.GetResolution(searchmodel) seqid = search[0][i][1][0] modelcompleteness = search[0][i][1][1] SSMseqid = -1 ClustalWlength = -1 SSMlength = -1 if len(search[0][i][1]) > 2: SSMseqid = search[0][i][1][2] ClustalWlength = search[0][i][1][3] SSMlength = search[0][i][1][4] targetchainids = search[0][i][2] targetchainid = targetchainids[0] clustfname = posixpath.join(self.topfolder, self.targetpdbid + "_" + targetchainids, searchmodel[:6], "clust.aln") myalignment = open(RecaseFilename2(clustfname),"r").read() clustparser = bioinformatics.clustal_alignment_parse clustaln = clustparser(myalignment)[0] ClustalWlength = clustaln.shortest_seqlen() # still need to calculate ssmseqid and ssmalignlen on the fly fullmodelcompleteness = 0.0 for j in range(0,i+1): fullmodelcompleteness += float(search[0][j][1][1]) def GetModelResidueNumbers(smodel, chainids): # get number of residues in the model PDB pdbfname = self.GetSculptFname(chainids, smodel) pdbfname = os.path.relpath(pdbfname, os.path.join(self.topfolder)) nmodelresidues = SimpleFileProperties.GetNumberofResiduesinPDBfile(pdbfname) return nmodelresidues nmodelresidues = GetModelResidueNumbers(searchmodel, targetchainids) def GetSSMrms(superposelog): rms = -1.0 if os.path.exists(superposelog): slogfile = open(superposelog,"r") logtxt = slogfile.read() slogfile.close() m = logtxt.find(" at RMSD = ",0) n = logtxt.find("\n",m) # logtxt[m:n] now holds a string like # " at RMSD = 3.255 and alignment length 3" match1 = re.search(r'(at\sRMSD\s=\s*)([0-9]*\.?[0-9]+)(.*)', logtxt[m:n], re.MULTILINE) # old format m = logtxt.find(" r.m.s.d:",0) n = logtxt.find("\n",m) # logtxt[m:n] now holds a string like # " r.m.s.d: 0.7582 (A)" match2 = re.search(r'\s+r\.m\.s\.d:\s+([0-9]*\.?[0-9]+)', logtxt[m:n], re.MULTILINE) # new format try: if match1 != None: if len(match1.groups()) == 3: rms = float(match1.group(2)) else: if len(match2.groups()) <1: rms = float(match2.group(1)) except Exception, m: self.tlogger.log(logging.WARN, "No SSM RSM in %s for %s " %(superposelog, searchmodel)) rms = -1.0 return rms # get RMS value from SSM superpose.log file slog = posixpath.join(self.topfolder, self.targetpdbid + "_" + targetchainids, searchmodel, self.superposelogfname) SSMrms = GetSSMrms(slog) # Chotia & Lesk seqid-rms formula C_L_rms = max(0.8, 0.4*math.exp(1.87*(1.0-float(seqid)/100.0))) # new rms formula obtained from RMS v seqid scatter plots new_rms = 0.0568808 * math.exp(1.52198 * (1.0 - float(seqid)/100.0)) * math.pow(172.53 + nmodelresidues, 1.0/3.0) # use if we decide to predict rms values with our new fit rather than Chothia and Lesk # C_L_rms = new_rms dres = max(2.5,float(datares)) # old calculations defaults to only 2.5A if self.refinevrmsdefault == True: dres = datares # calculate fracvar(resolution, modelcompleteness, seqid) fracvar = self.FracVar(C_L_rms, dres, modelcompleteness) scopfamnames = self.scopobj.GetFamilyNameFromPDB(searchmodel[:4], searchmodel[5]) scopstr = "\t".join(scopfamnames[0]) # count the number of unique scopids this model has been assigned modelscopidset = set([]) for scopid in scopfamnames: if scopid[0] != '': modelscopidset.add(scopid) nmodelscopid = len(modelscopidset) modelscopids = " | ".join([e[0] for e in modelscopidset]) modelclassnames =" | ".join([e[1] for e in modelscopidset]) modelfoldnames =" | ".join([e[2] for e in modelscopidset]) modelsuperfamilynames =" | ".join([e[3] for e in modelscopidset]) modelfamilynames =" | ".join([e[4] for e in modelscopidset]) targetscopidnames = self.scopobj.GetFamilyNameFromPDB(self.targetpdbid, targetchainid) # count the number of unique scopids the target has been assigned scopidset = set([]) for scopid in targetscopidnames: if scopid[0] != '': scopidset.add(scopid) ntargetscopid = len(scopidset) targetscopids = " | ".join([e[0] for e in scopidset]) targetclassnames =" | ".join([e[1] for e in scopidset]) targetfoldnames =" | ".join([e[2] for e in scopidset]) targetsuperfamilynames =" | ".join([e[3] for e in scopidset]) targetfamilynames =" | ".join([e[4] for e in scopidset]) # get the Wilson B-factor from logfile if we haven't got it already from previous logfile path = os.path.join(self.topfolder, "MR", workdir) path = os.path.relpath(path, self.topfolder) logfname = os.path.join(path, self.calclabel + "{" + workdir + "}" + partmodels + ".log") #print "logfname= " + logfname # might have lost the {} in the log file name on unix if not os.path.exists(logfname): logfname = os.path.join(path, self.calclabel + workdir + partmodels + ".log") if os.path.exists(logfname): logfile = file(logfname).read() else: logfile = "" nreflct = int(mrsol[2]) if WilsonB[0] < 0.0: # i.e. it's the first time we pick up results for the target mtz file wba = self.WILSON_Brecmp.findall(logfile) if len(wba) > 0: WilsonB[0] = float(wba[0]) if WilsonScale[0] < 0.0: # i.e. it's the first time we pick up results for the target mtz file wsc = self.WILSON_Scalerecmp.findall(logfile) if len(wsc) > 0: WilsonScale[0] = float(wsc[0]) # get the MatthewsCoefficient from logfile if we haven't got it already from previous logfile if MatthewsCoef[0] < 0.0: # i.e. it's the first time we pick up results for the target mtz file mca = self.MATTHEWCOEFrecmp.findall(logfile) if len(mca) > 0: MatthewsCoef[0] = float(mca[0]) if TargetScattering[0] < 0.0: # i.e. it's the first time we pick up results for the target mtz file os.chdir("Phaserinput") files = glob.glob("*.seq") TargetScattering[0] = 0 for seqfile in files: seqtxt = open(seqfile,"r").read() scat, msg = SimpleFileProperties.GetScatteringFromSequence(seqtxt) TargetScattering[0] += scat if os.path.exists("xtriagecache.txt") == False: mtzfname = glob.glob("*.mtz")[0] xtriagecache = self.GetXtriageResults(mtzfname) open("xtriagecache.txt", "w").write(str(xtriagecache)) else: xtriagecache = eval(open("xtriagecache.txt", "r").read()) (IsigIgrt3res, hires, lowres, noutliers, icerings, twinLtest) = xtriagecache os.chdir(self.topfolder) # note which is the last element in the database so UPDATE matches a corresponding INSERT statement self.mydb.execute("SELECT COUNT(*) AS total_record FROM MRTargetTbl") lastMRTargetTbl_id[0] = self.mydb.fetchone()[0] +1 if suffix == "": self.mydb.execute("INSERT INTO MRTargetTbl (" "TargetPDBid, SCOPid, NumberofUniqueSCOPids, TargetClass, TargetFold, " "TargetSuperfamily, TargetFamily, Unitcell, Spacegroup, Spacegroupnumber, " "Crystalsystem, IsPolarSpacegroup, TargetScattering, NumberofResiduesinTarget, " "WilsonB, WilsonScale, MatthewsCoefficient, Resolution, UsedReflections, " "Allreflections, ASUvolume, Robs, Rall, Rwork, Rfree, EDSentry, " "IsigIgrt3res, hires, lowres, outliers, icerings, twinLtest) " "VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)", (self.targetpdbid, targetscopids, ntargetscopid, targetclassnames, targetfoldnames, targetsuperfamilynames, targetfamilynames, Unitcell, spacegroup, spacegroupnumber, crystalsystemname, ispolarspacegroup, TargetScattering[0], ntargetresidues, WilsonB[0], WilsonScale[0], MatthewsCoef[0], datares, nreflct, Allreflections, Roundoff(ASUvolume,2), Rfactors[0], Rfactors[1], Rfactors[2], Rfactors[3], EDSentry, Roundoff(IsigIgrt3res,3), Roundoff(hires,3), Roundoff(lowres,3), noutliers, icerings, Roundoff(twinLtest,3))) fbval = "" # fallback value when llgs or other values aren't present cputimes = self.CPUtimerecmp.findall(logfile) cputime = -1 if len(cputimes) > 0: cputime = float(cputimes[-1]) if suffix == "": self.mydb.execute("INSERT INTO MRModelTbl (MRTarget_id, Fullmodel, ModelOrder, " "Foundit, Comment, C_L_rms, SSMrms, NewRMSD, Fracvar, PartialModel, " "CPUtime, NumberofResiduesinModel) VALUES " "(?,?,?,?,?,?,?,?,?,?,?,?)", (lastMRTargetTbl_id[0], workdir, nadditionalmodel+1, nsolvedit, "", Roundoff(C_L_rms,4), Roundoff(SSMrms,4), Roundoff(new_rms,4), Roundoff(fracvar,5), searchmodel, cputime, nmodelresidues)) anno = "" annolist = [] if nsolvedit >= 0: for sol in mrsol[0]: anno += sol.ANNOTATION + ", " annolist.append(sol.ANNOTATION) def GetBfacFracvarVrms(sol): vrms = [] fracvarvrms = [] Bfac = [] modelids = [] if len(sol.NEWVRMS.items()) > 0: # assuming there is only one VRMS per model for m,v in enumerate(sol.NEWVRMS.items()): vrms.append(v[1][0]) vrms_ensname = v[0] Bfac.append(sol.KNOWN[m].getBfac() ) modelids.append(sol.KNOWN[m].getModlid()) # extract partial model name '1IRU_K0' from ensemble name 'MR_1IRU_K0' partmodelname = vrms_ensname[3:] # find target chainid corresponding to this partial model name in search[0] ch = "" partmodelcompleteness = "" for s in search[0]: if s[0] == partmodelname: ch=s[2] partmodelcompleteness = s[1][1] fracvarvrms.append(self.FracVar(v[1][0], dres, partmodelcompleteness)) return Bfac, vrms, fracvarvrms, modelids ncomponents = 0 itemplllg = [None] templCorrCoef = mrsol[3][6] #if mrsol[3][6] != None else [(fbval, fbval)] templllg = mrsol[3][3][0] if mrsol[3][3] != None else fbval llgtemplreflct = float(templllg)/nreflct if templllg != fbval else fbval #import code; code.interact(local=locals()) #print "templCorrCoef= " + str(templCorrCoef) if mrsol[1] != None: templanno = mrsol[1][0].ANNOTATION # get the annotation for template solution itemplllg = self.TEMPLLLGrecmp.findall(templanno) # store VRMS for the template solution if present # assuming there is only one VRMS per model if len(mrsol[1][0].NEWVRMS.values()) > 0: if suffix == "": self.mydb.execute("INSERT INTO TemplateRMSTbl (MRTarget_id, MRModel_id, " "RefinedModel, LLGvrms, iLLGtemplate, LLGreflvrms) VALUES (?,?,?,?,?,?)", (lastMRTargetTbl_id[0], lastMRModel_id, partmodels, Roundoff(templllg,3), itemplllg[0], Roundoff(llgtemplreflct,6))) tmplBfac, tmplvrms, tmplfracvarvrms, tmplmodelids = GetBfacFracvarVrms(mrsol[1][0]) for i in range(len(tmplvrms)): for (k,p) in enumerate(partmodellist): #import code; code.interact(local=locals()) if p == tmplmodelids[i][3:]: sqlexp = "UPDATE TemplateRMSTbl SET vrms%d=?, Bfac%d=?, fracvarvrms%d=? \ WHERE MRModel_id=?" %(k,k,k) self.mydb.execute(sqlexp, (Roundoff(tmplvrms[i],4), Roundoff(tmplBfac[i],2), Roundoff(tmplfracvarvrms[i],5), lastMRModel_id)) if templCorrCoef != None: tmplCCglobal = templCorrCoef[0][0] tmplCCensmap = [] tmplCCenspdb = [] tmplmladens = [] tmplllgens = [] knowns = mrsol[1][0].KNOWN for i in range(len(knowns)): if isinstance( templCorrCoef[0][2], tuple ): # using old format # [[0.3576, 0.0, (0.3105, 0.203, 0.8105), (0.2222, 0.136, 1.0067)], # [0.283, 0.0, (0.3056, 0.203, 0.8174), (0.0381, 0.011, 3.4518)]] tmplCCensmap.append(templCorrCoef[0][2+i][0]) tmplCCenspdb.append(templCorrCoef[0][2+i][1]) tmplmladens.append(templCorrCoef[0][2+i][2]) tmplllgens.append(None) else: tmplCCensmap.append(templCorrCoef[0][3+i][0]) tmplCCenspdb.append(templCorrCoef[0][3+i][1]) tmplmladens.append(templCorrCoef[0][3+i][2]) tmplllgens.append(templCorrCoef[0][3+i][3]) #import code; code.interact(local=locals()) for i in range(len(tmplCCensmap)): sqlexp = "UPDATE TemplateRMSTbl SET CCglobal=?, CCmap%d=?, \ CCpdb%d=?, MLAD%d=?, LLGensmbl%d=? WHERE MRModel_id=?" %(i,i,i,i) self.mydb.execute(sqlexp, (tmplCCglobal, tmplCCensmap[i], tmplCCenspdb[i], tmplmladens[i], tmplllgens[i], lastMRModel_id)) # get LLG values from separate text files with higher precision than in the annotation strings LLGreals = mrsol[3][0] LLGVRMSreals = mrsol[3][1] LLGVRMSfullres = mrsol[3][2] CorrCoefs = mrsol[3][4] #if mrsol[3][4] != None else [(fbval,fbval)] CorrCoeffullres = mrsol[3][5] if mrsol[3][5] != None else [(fbval,fbval)] if suffix == "": self.mydb.execute("UPDATE MRModelTbl SET ModelResolution=?, " "SequenceIdentity=?, SSMseqid=?, ModelCompleteness=?, FullModelCompleteness=?, " "ClustalWlength=?, SSMlength=?, Annotation=?, " "ModelSCOPid=?, NumberofUniqueModelSCOPids =?, ModelClass=?, ModelFold=?, " "ModelSuperfamily=?, ModelFamily=? WHERE ROWID=?", (resolution, seqid, SSMseqid, modelcompleteness, fullmodelcompleteness, ClustalWlength, SSMlength, anno, modelscopids, nmodelscopid, modelclassnames, modelfoldnames, modelsuperfamilynames, modelfamilynames, lastMRModel_id)) # label model according to which set of externally supplied PDB codes it belongs to (say NMR, EM) for tagset in PDBTagsets: belongstoset = 0 if searchmodel[:4] in tagset[0]: belongstoset = 1 sqlxpr = "UPDATE MRModelTbl SET %s=%d WHERE ROWID=%d"\ %(tagset[1], belongstoset, lastMRModel_id) if suffix == "": self.mydb.execute(sqlxpr) # populate ELLGtables if it is one of those 500 calculations if mrsol[9] != []: for k, line in enumerate(mrsol[9]): if k==0 or len(line) < 9: continue foundit = eval(line[6])[0] if len(eval(line[6])) > 0 else 0 if suffix == "": self.mydb.execute("INSERT INTO ELLGtbl (MRModel_id, eLLGtarget, " "achievedELLG, LLGachievedres, LLGfullres, TFZfullres, TFZachievedres, "\ + "foundit, vrmsachievedres, achievedresol, foundit2) VALUES (?,?,?,?,?,?,?,?,?,?,?)", (lastMRModel_id, float(line[0]), float(line[1]), float(line[2]), float(line[3]), float(line[4]), float(line[5]), foundit, float(line[7]), float(line[8]), line[9])) llgs = [] maxllgperreflct = -1000 # unlikely to be retained when comparing against llgperreflct maxllgreal = -1000 # if nsolvedit >= 0: ncomponents = len(search[0][i][2]) # if only partially found then count how many components phaser found foundcomponents = 0 if len(mrsol[0]) > 0: for component in mrsol[0][0].KNOWN: if component.getModlid().find(searchmodel) > 0: foundcomponents +=1 if ncomponents != foundcomponents: self.logger.log(logging.WARN,"Solution file has %d components of %s but %d is expected." \ %(foundcomponents, searchmodel, ncomponents)) rfzs = [] tfzs = [] llgs = [] llgfullress = [] paks = [] tfzequiv = [] llgmax = fbval tfzmax = fbval for (j,anns) in enumerate(annolist): rfz = self.RFZrecmp.findall(anns) tfz = self.TFZrecmp.findall(anns) tfzequiv = self.TFZequivrecmp.findall(anns) llg = self.LLGrecmp.findall(anns) pak = self.PAKrecmp.findall(anns) star = self.Tstarrecmp.findall(anns) llgvrmsreal = fbval #llgreal = LLGreals[j] if LLGreals != None and LLGreals != [None] and len(LLGreals) > j else fbval #llgvrmsreal = LLGVRMSreals[j] if LLGVRMSreals != None and len(LLGVRMSreals) > j and LLGVRMSreals != [None] else fbval llgvrmsfullres = LLGVRMSfullres[j] if LLGVRMSfullres != None and LLGVRMSfullres != [None] and len(LLGVRMSfullres) > j else fbval #CC = CorrCoefs[j] if len(CorrCoefs) > j else (fbval, fbval) CCglobal = fbval CCensmap = [] CCenspdb = [] mladens = [] llgsens = [] #knowns = mrsol[0][j].KNOWN if type(CorrCoefs) == list and len(CorrCoefs) > j and len(CorrCoefs[j]) > 0: if isinstance( CorrCoefs[j][2], tuple ): # using old format # [[0.3576, 0.0, (0.3105, 0.203, 0.8105), (0.2222, 0.136, 1.0067)], # [0.283, 0.0, (0.3056, 0.203, 0.8174), (0.0381, 0.011, 3.4518)]] CCglobal = CorrCoefs[j][0] llgvrmsreal = LLGVRMSreals[j] if LLGVRMSreals != None and len(LLGVRMSreals) > j and LLGVRMSreals != [None] else fbval for i in range(len(CorrCoefs[j])-2): CCensmap.append(CorrCoefs[j][2+i][0]) CCenspdb.append(CorrCoefs[j][2+i][1]) mladens.append(CorrCoefs[j][2+i][2]) llgsens.append(None) else: # using new format which includes LLG values # [[0.3576, 0.0, 89.639, (0.3105, 0.203, 0.8105, 45.367), (0.2222, 0.136, 1.0067, 23.764)], # [0.283, 0.0, 37.207, (0.3056, 0.203, 0.8174, 12.976), (0.0381, 0.011, 3.4518, 62.891)]] CCglobal = CorrCoefs[j][0] llgvrmsreal = CorrCoefs[j][2] for i in range(len(CorrCoefs[j])-3): CCensmap.append(CorrCoefs[j][3+i][0]) CCenspdb.append(CorrCoefs[j][3+i][1]) mladens.append(CorrCoefs[j][3+i][2]) llgsens.append(CorrCoefs[j][3+i][3]) #CCfullres = CorrCoeffullres[j] if len(CorrCoeffullres) > j else (fbval, fbval) if tfzequiv != [] and tfzequiv[0]: #llgmax = float(llg[0]) # LLG corresponding to TFZequiv for MR without vrms refinement #tfzmax = float(tfz[0]) # TFZ ----- if suffix == "": self.mydb.execute("UPDATE MRModelTbl SET TFZequiv=? WHERE ROWID=?", (tfzequiv[0], lastMRModel_id)) #llgmax = float(llg[0]) # LLG corresponding to TFZequiv for MR without vrms refinement else: #print "workdir= %s, searchmodel= %s, tfzequiv= %s " %(workdir, searchmodel, tfzequiv) # now write TFZequiv from other solution results to database if provided #sqlexp = "UPDATE MRModelTbl SET TFZequiv%s = 42 WHERE MRTarget_id=? " %suffix #self.mydb.execute(sqlexp, lastMRTargetTbl_id[0]) sqlexp = "UPDATE MRModelTbl SET TFZequiv%s=? WHERE MRTarget_id=? " \ "AND Fullmodel=? AND PartialModel=?" %suffix #print "sqlexp= " + sqlexp self.mydb.execute(sqlexp, (tfzequiv[0], lastMRTargetTbl_id[0], workdir, searchmodel)) #llgmax = float(llg[-1]) # LLG corresponding to TFZequiv for additional RNP if star != []: Tstar = star[0] if Tstar: if suffix == "": self.mydb.execute("UPDATE MRModelTbl SET SolutionRank=? WHERE ROWID=?", ((j+1), lastMRModel_id)) else: sqlexp = "UPDATE MRModelTbl SET SolutionRank%s=? WHERE MRTarget_id=? " \ "AND Fullmodel=? AND PartialModel=?" %suffix self.mydb.execute(sqlexp, ((j+1), lastMRTargetTbl_id[0], workdir, searchmodel)) else: Tstar = 0 rfz1 = float(rfz[0]) if len(rfz) > 0 else fbval rfzs.append(rfz1) tfz1 = float(tfz[0]) if len(tfz) > 0 else fbval tfzs.append(tfz1) pak1 = float(pak[0]) if len(pak) > 0 else fbval paks.append(pak1) if len(llg) > 1: llgs.append(float(llg[-1])) # for now assume last LLG value is with restricted resolution llgfullress.append(float(llg[1])) else: llgfullress.append(fbval) llgs = [fbval] # fallback if value is missing llg = [fbval] # fallback if value is missing if len(llg) > 2: illgvrms = llg[2] else: illgvrms = fbval #llgperreflct = float(llgreal)/float(nreflct) if llgreal != fbval else fbval llgperreflvrms = float(llgvrmsreal)/float(nreflct) if llgvrmsreal != fbval else fbval llgperreflvrmsfullres = float(llgvrmsfullres)/float(Allreflections) if llgvrmsfullres != fbval else fbval llgvrmscombined = fbval llgvrmscombrefl = fbval # assuming the highest llg corresponds to the llg for the solution emitted by phaser # get it now so we can put it into MRModelTbl maxllgperreflct = max(llgperreflvrms, maxllgperreflct) maxllgreal = max(llgvrmsreal, maxllgreal) #import code; code.interact(local=locals()) Bfac, vrms, fracvarvrms, modelids = GetBfacFracvarVrms(mrsol[0][j]) # and compute fracvarvrms for the full resolution vrmsfullres = fbval fracvarvrmsfullres = fbval if len(mrsol[7]) > j and len(mrsol[7][j].NEWVRMS.values()) > 0: Bfacfullres, vrmsfullres, fracvarvrmsfullres, modelids = GetBfacFracvarVrms(mrsol[7][j]) # insert annotations for individual solutions into the second table of the database if suffix == "": self.mydb.execute("INSERT INTO AnnotationTbl (MRTarget_id, MRModel_id, " "Tstar, RFZ, TFZ, PAK, " "CCglobal, LLGfullres_vrms, LLGfullres_refl_vrms, iLLG, iLLGvrms, " "LLGvrms, LLGrefl_vrms, VRMSfullres, " "FracvarVrmsfullres) " "VALUES (?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)", (lastMRTargetTbl_id[0], lastMRModel_id, Tstar, rfz1, tfz1, pak1, Roundoff(CCglobal,5), llgvrmsfullres, Roundoff(llgperreflvrmsfullres,6), llg[0], illgvrms, Roundoff(llgvrmsreal,3), Roundoff(llgperreflvrms,6), Roundoff(vrmsfullres,4), Roundoff(fracvarvrmsfullres,5) ) ) self.mydb.execute("SELECT COUNT(*) as total_record FROM AnnotationTbl") lastAnnotation_id = self.mydb.fetchone()[0] for i in range(len(CCensmap)): sqlexp = "UPDATE AnnotationTbl SET CCmap%d=?, CCpdb%d=?, MLAD%d=?, \ LLGensmbl%d=? WHERE ROWID=?" %(i,i,i,i) self.mydb.execute(sqlexp, (CCensmap[i], CCenspdb[i], mladens[i], llgsens[i], lastAnnotation_id)) for i in range(len(vrms)): sqlexp = "UPDATE AnnotationTbl SET vrms%d=?, Bfac%d=?, fracvarvrms%d=? \ WHERE ROWID=?" %(i,i,i) self.mydb.execute(sqlexp, (Roundoff(vrms[i],4), Roundoff(Bfac[i],2), Roundoff(fracvarvrms[i],5), lastAnnotation_id)) # in case calculations didn't finish make empty entries as needed in the # arrays to preserve formatting in the spreadsheet while len(rfzs) < ncomponents: rfzs.append(None) while len(tfzs) < ncomponents: tfzs.append(None) while len(paks) < ncomponents: paks.append(None) maxllgperreflct = fbval if maxllgperreflct == -1000 else maxllgperreflct maxllgreal = fbval if maxllgreal == -1000 else maxllgreal for k in range(0,ncomponents): if suffix == "": # list LLG/reflections for the intermediate LLGs, not the final LLG that uses all reflections sqlexp = "UPDATE MRModelTbl SET RFZ%d=?, TFZ%d=?, PAK%d=?, LLG%d=?, \ iLLG%d=?, LLGreflection%d=? WHERE ROWID=?" %(k,k,k,k,k,k) self.mydb.execute(sqlexp, (rfzs[k], tfzs[k], paks[k], Roundoff(maxllgreal,3), llgmax, Roundoff(maxllgperreflct,6), lastMRModel_id)) # last LLG is the final LLG. Print only if available, i.e. solution file is present if len(llgs) > 0: self.mydb.execute("UPDATE MRModelTbl SET iLLGfullres%d=? WHERE ROWID=?" %k, (llgfullress[k], lastMRModel_id)) # now write llg from other solution results to database if provided if suffix != "": sqlexp = "UPDATE MRModelTbl SET iLLG%s=? WHERE MRTarget_id=? " \ "AND Fullmodel=? AND PartialModel=?" %suffix self.mydb.execute(sqlexp, (llg1, lastMRTargetTbl_id[0], workdir, searchmodel)) # now write first values from other solution results to database if provided sqlexp = "UPDATE MRModelTbl SET TFZ0%s=? WHERE MRTarget_id=? " \ "AND Fullmodel=? AND PartialModel=?" %suffix self.mydb.execute(sqlexp, (tfz1, lastMRTargetTbl_id[0], workdir, searchmodel)) sqlexp = "UPDATE MRModelTbl SET Foundit%s=? WHERE MRTarget_id=? " \ "AND Fullmodel=? AND PartialModel=?" %suffix self.mydb.execute(sqlexp, (nsolvedit, lastMRTargetTbl_id[0], workdir, searchmodel)) def multiply(self, lists): # cartesian product of a list of bins if not lists: return [] product=zip(lists[0]) for l in lists[1:]: prods=[] for p in product: prods.extend([p+(tmp,) for tmp in l]) product = prods return product def PutPDBsInChainIdBins(self, pdbs_seqids): # Get list of unique chain IDs from the pdbs_seqids list. # pdbs_seqids should already be sorted according to chain IDs chainIDs = [] for i in range(0,len(pdbs_seqids)): if len(chainIDs) == 0: chainIDs.append(pdbs_seqids[i][2]) else: if pdbs_seqids[i][2] != chainIDs[len(chainIDs)-1]: chainIDs.append(pdbs_seqids[i][2]) if self.useuniquesequence: self.logger.log(logging.INFO,"pdb_chainid models refers to %d unique chains" % len(chainIDs)) # make bins of chain labels where each bin holds an array of models for that chain models_chains = [] for i in range(0,len(chainIDs)):#first label each sub array with the chain models_chains.append((chainIDs[i],[])) # sort pdbs_seqids into those bins for i in range(0,len(pdbs_seqids)): for j in range(0,len(models_chains)): if models_chains[j][0] == pdbs_seqids[i][2]: models_chains[j][1].append(pdbs_seqids[i]) #invert indexing and ignore the very first element which is just the chainID bins = [models_chains[i][1] for i in range(0,len(models_chains))] for i in range(0,len(bins)): self.logger.log(logging.INFO,"There are %d models for chain %s" %(len(bins[i]), models_chains[i][0])) return bins def GetCombinedSearchesforTargetChains(self, pdbs_seqids): bins = self.PutPDBsInChainIdBins(pdbs_seqids) myseq = self.PDB2SeqLocal(self.targetpdbid) if len(bins) < self.nuniqueseq: self.logger.log(logging.WARN, "Models were not found for all chains in the target structure!") self.logger.log(logging.INFO,"pdb_chainid,\tClustalW2 seq identity,\tCompleteness of model\t target pdb chainid(s)") for i in range(0, len(pdbs_seqids)): self.logger.log(logging.INFO,"%s \t\t %s \t\t %s \t\t %s" %(pdbs_seqids[i][0], pdbs_seqids[i][1][0], pdbs_seqids[i][1][1], pdbs_seqids[i][2])) # do cartesian product of the bins with one another to get all combinations of models searches = self.multiply(bins) #mfile = open("testlist.txt","w") #mfile.write(str(searches)) #mfile.close() uniquesearchdict = {} #for e in searches: # print e for search in searches: residuenumbers = "" for pdbseqid in search: residuenumbers += pdbseqid[4] uniquesearchdict[residuenumbers] = search uniquesearches = [] for e in uniquesearchdict: uniquesearches.append( uniquesearchdict[e] ) #for e in uniquesearches: # print e # A sorted list of pdbs in the searches is helpful for identifying other searches # that only differs by permuations from this one sortedsearches = [] for search in uniquesearches: sortedsearches.append( sorted(search) ) self.logger.log(logging.INFO,"Combining those models yields the following %s MR calculations:" %len(sortedsearches)) # pretty print as to only display the model ids for i in range(0,len(sortedsearches)): mrcalcstr = sortedsearches[i][0][0] for j in range(1,len(sortedsearches[i])): mrcalcstr = mrcalcstr + ", " + sortedsearches[i][j][0] self.logger.log(logging.INFO,mrcalcstr) return sortedsearches def PermuteSearches(self, search): # make all possible permutations of search orders with these models self.logger.log(logging.INFO,"\nPermuting search orders of models:") # pretty print search models as concatenated PDB codes separated by commas modelsdir = [] for pdb_id in search: modelsdir.append( pdb_id[0]) self.logger.log(logging.INFO, ",".join( modelsdir ) ) # now permute the list lst = list( itertools.permutations( search ) ) permutations = [] for m in lst: permutations.append(copy.deepcopy(list(m))) self.logger.log(logging.INFO,"yields the following search orders:") for line in permutations: # pretty print permuted search order as concatenated PDB codes separated by commas modelsdir = [] for pdb_id in line: modelsdir.append( pdb_id[0]) self.logger.log(logging.INFO, ",".join( modelsdir ) ) """ return a list of permutations like search: 1I1A_D,1NFD_E,1W72_H yields the following permutations: 1I1A_D,1NFD_E,1W72_H 1I1A_D,1W72_H,1NFD_E 1NFD_E,1I1A_D,1W72_H 1NFD_E,1W72_H,1I1A_D 1W72_H,1I1A_D,1NFD_E 1W72_H,1NFD_E,1I1A_D """ return permutations def GetXtriageResults(self, mtzfname): mtzobj = iotbx.mtz.object(file_name = mtzfname) Fcol = [""]; dummy1 = [""]; dummy2 = [""]; dummy3 = [""] IsigIgrt3res = ""; hires = ""; lowres = ""; noutlier_outside_rescut = "" nicerings = ""; twinLtest = "" self.GetMtzColumns(mtzfname,Fcol, dummy1, dummy2, dummy3) lowres = mtzobj.max_min_resolution()[0] hires = mtzobj.max_min_resolution()[1] try: xtrobj = xtriage.run([mtzfname,"--quiet", "scaling.input.xray_data.obs_labels=%s" %Fcol[0]]) IsigIgrt3res = xtrobj.basic_results.data_strength.resolution_cut noutlier_outside_rescut = 0 # get outliers and sort them according to d_spacing ac_outliers = xtrobj.basic_results.basic_data_stats.outlier.acentric_outliers_table.data if ac_outliers != None: zip_outliers = zip(ac_outliers[0], ac_outliers[1], ac_outliers[2], ac_outliers[3], ac_outliers[4]) sortoutliers = sorted(zip_outliers, key=lambda tup: tup[0]) for e in sortoutliers: noutlier_outside_rescut +=1 if e[0] > IsigIgrt3res: break nicerings = 0 for e in xtrobj.basic_results.basic_data_stats.ijsco.ice_ring_bin_location: if e != None: nicerings += 1 twinLtest = xtrobj.twin_results.twin_summary.twin_results.maha_l except Exception, m: self.failure = True self.message = m self.logger.log(logging.ERROR, self.message) self.errorcodesum[0] += 50 return (IsigIgrt3res, hires, lowres, noutlier_outside_rescut, nicerings, twinLtest) def GetMtzColumns(self, fname, Fcol, sigFcol, Icol, sigIcol): try: obj = iotbx.mtz.object(file_name = fname) columns = obj.column_labels() if self.buseIntensities: for i in range(0,len(columns)): columntype = obj.get_column(columns[i]).type() if columntype=='J': Icol[0] = columns[i] for i in range(0,len(columns)): columntype = obj.get_column(columns[i]).type() if columntype=='Q' and columns[i].find(Icol[0]) > 1: sigIcol[0] = columns[i] else: for i in range(0,len(columns)): columntype = obj.get_column(columns[i]).type() # if Fcol[0] != '' and Sigcol[0] != '': # break # found the first columns of types amplitude and sigma so use these # else: if columntype=='F': Fcol[0] = columns[i] for i in range(0,len(columns)): columntype = obj.get_column(columns[i]).type() if columntype=='Q' and columns[i].find(Fcol[0]) > 1: sigFcol[0] = columns[i] except Exception, m: self.failure = True self.message = m self.logger.log(logging.ERROR, self.message) self.errorcodesum[0] += 50 #raise BLASTrunnerException, m def RunMRauto(self, psearches, njobs, minmodelorder): if len(psearches) > 0: ndistinctchains = len(psearches[0][1:][0]) self.logger.log(logging.INFO, "Will now do at most " + str(ndistinctchains* len(psearches)) + " phaser MR calculations in %d threads\n" % njobs) else: return os.chdir(self.topfolder) # Spread evenly searches with low seq ids (the lengthy jobs) onto each thread. # As the searches list begins with models having lowest seq id do this by dishing out # searches to threads like cards to game players. shortlists=[] for i in range(0,njobs): shortlists.append([]) # distribute searches as evenly as possible into the shortlists for i in range(0, len(psearches)): shortlists[i%njobs].append(psearches[i]) mylock = threading.Lock() phaserjobs = [] for ( index, shortlist ) in enumerate( shortlists ): self.logger.log(logging.INFO,"%d. thread runs the following MR jobs:" %(index+1)) for line in shortlist: # create modelsdir as concatenated PDB codes separated by commas modelsdir = [ [e[0] for e in e1] for e1 in line[:-1]][0] self.logger.log(logging.INFO, ",".join( modelsdir ) ) phaserjobs.append(BLASTphaserthread.PhaserRunner( baseinputfname = self.baseinputfname, modelinputfnameprefix = self.calclabel + "modelinput", calclabel = self.calclabel, rootstr = "", parent = self, searches = shortlist, fmaxhours = self.fmaxhours, nmaxsols = self.nmaxsols, mylock = mylock, threadnumber = index+1, maxthreads = njobs, minmodelorder = minmodelorder ) ) self.logger.log(logging.INFO, "\nTo abort the jobs create an empty file named StopThreadsGracefully in %s\n" %self.topfolder) for i in range(0,njobs): phaserjobs[i].start() for i in range(0,njobs): phaserjobs[i].join() self.errorcodesum[0] += phaserjobs[i].errorcodesum self.logger.log(logging.INFO, "Done with RunMRauto") def BLASTnSetupMRruns(self, maxseqid, minseqid, maxpdbs, myinc = 1.0, noRNP=False): mstr = "%s, %s, in BLASTnSetupMRruns, maxseqid: %s, \ minseqid: %s, maxpdbs: %s, similarseqidentmodels: %s, Seed with refined template in multicomponent calculations: %s\n" \ %(self.targetpdbid, self.calclabel,maxseqid, minseqid, maxpdbs, myinc, str(self.busecurrenttmpl)) reposfile = open("BLASTMRcalculations.txt","a+") reposfile.write(mstr) reposfile.close() self.MakeMyfolder(self.topfolder) global LOG_FILENAME if LOG_FILENAME == "": LOG_FILENAME = os.path.join(self.topfolder, self.calclabel + 'BLASTnSetupMRruns.log') self.SetupLogfile(LOG_FILENAME) os.chdir(self.topfolder) uniquechains = [] self.logger.log(logging.INFO, mstr) pdbsseqids = self.GetPDBidsFromPDBLocalWithSequenceIds(maxseqid, minseqid, maxpdbs, uniquechains) self.RemoveQuarantinedPDBs(pdbsseqids) if myinc > 0.0: pdbs_seqids = self.RemovePDBsWithSimilarSeqIdentFromList(pdbsseqids, maxseqid, minseqid, myinc) else: pdbs_seqids = pdbsseqids self.logger.log(logging.INFO,"\nModels to be used are\npdb_id+chain_id, BLAST sequence_id, target chain_id:") mdlist = "\n".join( ["%s, %s, %s"%(a,b,c) for a,b,c in pdbs_seqids]) self.logger.log(logging.INFO,"\n%s\n" %mdlist) self.logger.log(logging.INFO,"There is(are) %d unique chain(s) in the target." %len(uniquechains)) self.FetchPDBfile("." + posixpath.sep, self.targetpdbid) # self.targetmass = self.GetMassOfModelInPDBfile(self.targetpdbid + ".pdb") # self.targetscattering = self.GetScatteringFromSequence(uniquechains[0][1], True) self.targetscattering = 0.0 for seq in uniquechains: scat, mstr= SimpleFileProperties.GetScatteringFromSequence(seq[1], True) self.targetscattering += scat self.MakeMyfolder("Phaserinput") os.chdir("Phaserinput") # Hop down in Phaserinput and store common input files there. myseq = self.PDB2SeqLocal(self.targetpdbid) seqfname =[] # store fasta sequences in folder, say 1tfx/fasta1.seq for i in range(0,len(myseq)): if myseq[i].find("|PDBID|CHAIN|SEQUENCE") > -1: # new convention for preambling sequences seqfname.append(myseq[i][1:5] + "_Chain" + myseq[i][6] + ".seq") else: # old preample convention seqfname.append(myseq[i][1:12] + '.seq')# use preamble in sequence for naming if os.path.exists(seqfname[i]) == True: self.logger.log(logging.INFO, "Using existing sequence file, %s" %seqfname[i]) else: self.logger.log(logging.INFO,"Writing sequence file %s", seqfname[i]) fileout = open(seqfname[i] ,"w") fileout.write(myseq[i]) fileout.close() #localmtzname = self.targetpdbid.lower() + '.mtz' localmtzname = self.targetpdbid + '.mtz' mtzfname = self.mtzfolder + posixpath.sep + localmtzname # Get the mtz file self.logger.log(logging.INFO, "Copying %s to %s" %(mtzfname, os.getcwd())) if os.path.exists(mtzfname) != True: self.logger.log(logging.WARN,"We don't have an mtz file for " + mtzfname + ". Will try downloading from the RCSB...") # mtz file not avialable so download cif file from RCSB u = urllib2.urlopen("http://www.pdb.org/pdb/files/r" + self.targetpdbid + "sf.ent.gz") cifgz = u.read() gzfile = open("tmp.gz",'wb') gzfile.write(cifgz) gzfile.close() self.logger.log(logging.INFO, "Writing tmp.gz") fileObj = gzip.GzipFile("tmp.gz", 'rb') cifdata = fileObj.read() # decompress the file fileObj.close() self.logger.log(logging.INFO, "Writing tmp.cif") fileout = open("tmp.cif",'w') fileout.write(cifdata) fileout.close() os.remove("tmp.gz") # delete the compressed file # convert cif file to mtz # using phenix.cif_as_mtz command = ["tmp.cif", "--symmetry=.." + posixpath.sep + self.targetpdbid + ".pdb", "--output-file-name=" + self.targetpdbid + ".mtz"] self.logger.log(logging.INFO, "Executing command: " + subprocess.list2cmdline(command)) cif_as_mtz.run(command) os.remove("tmp.cif") # delete cif file if os.path.exists(localmtzname) != True and os.path.exists(mtzfname) == True: shutil.copyfile(mtzfname, os.path.join(os.getcwd(), localmtzname)) psearches = [] if os.path.exists(localmtzname) != True: self.logger.log(logging.FATAL, "No mtz file available!") self.errorcodesum[0] = 42 return psearches Fcol = [''] Sigcol = [''] Icol = [''] SigIcol = [''] self.GetMtzColumns(localmtzname,Fcol,Sigcol,Icol, SigIcol) # only intensities avaialable in the mtz file so convert it to amplitudes if Fcol[0] == '' and Icol[0] != '' and Sigcol[0] != '': self.logger.log(logging.INFO, "Only intensities available in %s. Converting to amplitudes..." %localmtzname) mtzwithamplname = self.targetpdbid.lower() + '_I2ampl.mtz' if os.path.exists(mtzwithamplname): os.remove(mtzwithamplname) # to avoid windows python crash in reflection_file_converter.run intensity2amplarg = ["--write-mtz-amplitudes","--mtz-root-label=FOBS","--label=%s"% Icol[0], "--mtz=%s" %mtzwithamplname, "%s" %localmtzname] self.logger.log(logging.INFO, "reflection_file_converter " + subprocess.list2cmdline(intensity2amplarg)) reflection_file_converter.run(args = intensity2amplarg) FcolI2amp = [''] SigcolI2amp = [''] IcolI2amp = [''] # check that none of the intensities are negative due to measurement errors if self.GetMtzColumns(mtzwithamplname,FcolI2amp,SigcolI2amp,IcolI2amp) == False: self.logger.log(logging.INFO, "Some reflections from the conversion %s " "have amplitude values of zero. Must massage intensities..." %mtzwithamplname) mtzwithmssgamplname = self.targetpdbid.lower() + '_mssg_I2ampl.mtz' if os.path.exists(mtzwithmssgamplname): os.remove(mtzwithmssgamplname) # to avoid windows python crash in reflection_file_converter.run # if so then convert again and massage intensities intensity2amplarg = ["--write-mtz-amplitudes","--mtz-root-label=FOBS","--label=%s"% Icol[0], "--mtz=%s" %mtzwithmssgamplname, "--massage-intensities", "%s" %localmtzname] self.logger.log(logging.INFO, "reflection_file_converter " + subprocess.list2cmdline(intensity2amplarg)) reflection_file_converter.run(args = intensity2amplarg) localmtzname = mtzwithmssgamplname else: localmtzname = mtzwithamplname Fcol = [''] Sigcol = [''] self.GetMtzColumns(localmtzname,Fcol,Sigcol,Icol, SigIcol) if (Fcol[0] == '' or Sigcol[0] == '') and (Icol == '' or SigIcol == ''): self.logger.log(logging.FATAL, "Neither intensity nor amplitude columns seem present in %s" %localmtzname) self.errorcodesum[0] += 41 return psearches self.logger.log(logging.INFO, "Writing base input.txt file in " + os.getcwd()) infname =self.calclabel + "input.txt" if os.path.exists(infname): self.logger.log(logging.WARN, infname + " exists and will be overwritten.") inputfile = open(infname, "w") mtzstr = "Phaserinput" + posixpath.sep + localmtzname inputfile.write("HKLIN " + mtzstr + "\n") if self.buseIntensities: inputfile.write("LABIN I = " + Icol[0] + " SIGI = " + SigIcol[0] + "\n") else: inputfile.write("LABIN F = " + Fcol[0] + " SIGF = " + Sigcol[0] + "\n") for i in range(0, len(seqfname)): inputfile.write("COMPOSITION PROTEIN SEQUENCE Phaserinput" \ + posixpath.sep + seqfname[i] + ' NUM 1\n') #inputfile.write("XYZOUT ON ENSEMBLE ON\n") inputfile.close() os.chdir("..") preparedmodels = [] for pdbseqid in pdbs_seqids: if (self.PreparePartialModel(pdbseqid)): preparedmodels.append(pdbseqid) if noRNP: return psearches searches = self.GetCombinedSearchesforTargetChains(preparedmodels) for search in searches: psearch = self.PermuteSearches(search) if self.PrepareTemplateSolutions(search, psearch): # create new list of those models that made it through the sculpting and RNP psearches.extend(psearch) for psearch in psearches: for model in psearch: # don't need the model superposition eulers and translations anymore as they are very verbose model.remove(model[3]) modelinputs = self.PrepareInputFiles(psearch) psearch.append(modelinputs) #psearches format is like # [[[['1KCW_A', ['32.615', '0.358', 75.279, 742, 538], 'A', '{1-814}'], ['2ORX_A', ['34.395', '0.159', 40.94, 314, 149], 'B', '{273-586}']], ['1KCW_A', '1KCW_A,2ORX_A']], # [[['2ORX_A', ['34.395', '0.159', 40.94, 314, 149], 'B', '{273-586}'], ['1KCW_A', ['32.615', '0.358', 75.279, 742, 538], 'A', '{1-814}']], ['2ORX_A', '2ORX_A,1KCW_A']]] # changing it to # [[['1KCW_A', ['32.615', '0.358', 75.279, 742, 538], 'A', '{1-814}'], ['2ORX_A', ['34.395', '0.159', 40.94, 314, 149], 'B', '{273-586}'], ['1KCW_A', '1KCW_A,2ORX_A']], # [['2ORX_A', ['34.395', '0.159', 40.94, 314, 149], 'B', '{273-586}'], ['1KCW_A', ['32.615', '0.358', 75.279, 742, 538], 'A', '{1-814}'], ['2ORX_A', '2ORX_A,1KCW_A']]] # # last element is a list of succesive search model composition. Want this element to # be put in a separate list, i.e. convert list from [[m]] into [[m-1],[1]] #import code; code.interact(local=locals()) psearches = [[[e for e in e1[:-1]], e1[-1]] for e1 in psearches ] os.chdir(self.topfolder) blastout = open(self.calclabel + self.targetpdbid + "_BLASTchain_ids.txt" ,"w") blastout.write(str(psearches)) blastout.close() self.logger.log(logging.INFO, "\nSubfolders now set up for MR calculations") self.logger.log(logging.INFO,self.calclabel + self.targetpdbid + "_BLASTchain_ids.txt lists model pdb_chain ids") self.logger.log(logging.INFO,"with sequence identity and the permuted MR searches.") self.logger.log(logging.INFO,self.calclabel + "input.txt is toplevel MR input to be joined with") self.logger.log(logging.INFO,self.calclabel + "modelinput.txt before each MR calulation.") self.logger.log(logging.INFO,"\nDone with BLASTnSetupMRruns.") return psearches def CollateMRautoResultsFromlistFile(self, DBfname, TaggedPDBfileliststr, OnlyCreateDBcolumns, OtherSolutions): import sqlite3 # compile tab separated tex file to load into a spreadsheet reposfile = open("BLASTMRcalculations.txt","a+") reposfile.write("%s, %s, in CollateMRautoResultsFromlistFile\n" %(self.targetpdbid, self.calclabel)) reposfile.close() LOG_FILENAME = os.path.join(self.topfolder, self.calclabel + "MRResults.log") self.SetupLogfile(LOG_FILENAME) self.logger.log(logging.INFO, "CollateMRautoResultsFromlistFile, pdb: %s, " \ "calclabel: %s\n" %(self.targetpdbid, self.calclabel)) # a list of filenames for externally supplied PDB codes and corresponding labels (say NMR, EM) searchlistfname = os.path.join(self.topfolder, self.calclabel + self.targetpdbid + "_BLASTchain_ids.txt") searchlistfname = RecaseFilename2(searchlistfname) psearches = eval(open(searchlistfname,"r").read()) PDBTagsets = [] if OnlyCreateDBcolumns == False: TaggedPDBfilelst = TaggedPDBfileliststr.split() nlen = int(len(TaggedPDBfilelst)/2) for i in range(0, nlen): fname = TaggedPDBfilelst[i*2] tag = TaggedPDBfilelst[i*2+1] self.logger.log(logging.INFO, "MR models matching PDB ids in %s will be tagged as %s" %(fname,tag)) # create sets of PDB ids from each file mfile = open(fname, "r") mtxt = mfile.read() mfile.close() tpl = (set([]), tag) for pdbid in mtxt.split(): tpl[0].add(pdbid.upper()) PDBTagsets.append(tpl) # being used by GatherResults os.chdir(self.topfolder) olddir = os.getcwd() mtzname = posixpath.join("Phaserinput", self.targetpdbid.lower() + '.mtz') mtzname = RecaseFilename2(mtzname) mtzobj = iotbx.mtz.object(file_name = mtzname) xtalsysname = mtzobj.space_group().crystal_system() spgname = mtzobj.space_group_info().type().lookup_symbol() spgnumber = mtzobj.space_group_info().type().number() ispolarspacegroup = int(AlterOriginSymmate.list_origins(spgnumber)[1]) ucell = mtzobj.crystals()[0].unit_cell_parameters() unitcell = "%3.3f %3.3f %3.3f %3.3f %3.3f %3.3f" %(ucell[0], ucell[1],ucell[2],ucell[3],ucell[4],ucell[5]) # get size of ASU as the unit cell volume divided by the number of symmetry operations ASUvolume = mtzobj.crystals()[0].unit_cell().volume()/len(mtzobj.space_group().all_ops()) Allreflections = len(mtzobj.get_column(mtzobj.column_labels()[0]).extract_values()) # get the Rfactors for the target structure Rfactors = [-1.0, -1.0, -1.0, -1.0] Rfacstr = open(os.path.join(self.scoppath,"Rfactors.txt"),"r").read() m = re.search("("+ self.targetpdbid + ")\t(.*)\t(.*)\t(.*)\t(.*)\t(.*)\n", Rfacstr,re.IGNORECASE) if m != None: r = m.groups() Rfactors = [Roundoff(r[1],3), Roundoff(r[2],3), Roundoff(r[3],3), Roundoff(r[4],3)] # check if EDS can reproduce the Rfactor EDSstr = open(os.path.join(self.scoppath, "eds_holdings.txt"), "r").read() m = re.search("("+ self.targetpdbid + ")\n", EDSstr,re.IGNORECASE) if m != None: EDSreproduce = 1 else: EDSreproduce = 0 # open the database conn = sqlite3.connect(DBfname) self.mydb = conn.cursor() # some settings for increasing speed when inserting values self.mydb.execute('PRAGMA synchronous = OFF') self.mydb.execute('PRAGMA journal_mode = MEMORY') self.mydb.execute("PRAGMA cache_size=200000") # first check for the existence of our table self.mydb.execute('select tbl_name from sqlite_master') list_tables = self.mydb.fetchall() createtable = True for t in list_tables: if t[0] == "MRTargetTbl": createtable = False break if createtable: # Otherwise create the table with columns # First create the main table sqlxprs = "create table MRTargetTbl (p_id integer primary key autoincrement, " \ + "TargetPDBid text, SCOPid text, icerings int, twinLtest real, "\ + "NumberofUniqueSCOPids int, TargetClass text, TargetFold text, "\ + "TargetSuperfamily text, TargetFamily text, Spacegroup text, "\ + "Spacegroupnumber int, IsPolarSpacegroup int, Crystalsystem text, "\ + "IsigIgrt3res real, hires real, lowres real, outliers int, " \ + "Resolution real, Unitcell text, MatthewsCoefficient real, "\ + "NumberofResiduesinTarget int, WilsonB real, WilsonScale real, "\ + "Robs real, Rall real, Rwork real, Rfree real, EDSentry int, "\ + "TargetScattering int, UsedReflections int, Allreflections int, ASUVolume real)" self.mydb.execute(sqlxprs) sqlxprs = "create table MRModelTbl (p_id integer primary key autoincrement," \ + "MRTarget_id int, Fullmodel text, Comment text,"\ + "ModelOrder int, Foundit int, PartialModel text, "\ + "ModelResolution real, SequenceIdentity real, SSMseqid real, "\ + "NumberofResiduesinModel int, SolutionRank int, ClustalWlength int, SSMlength int, "\ + "FullModelCompleteness real, ModelCompleteness real, Annotation text, "\ + "C_L_rms real, SSMrms real, NewRMSD real, TFZequiv real, "\ + "Fracvar real, ModelSCOPid text, "\ + "NumberofUniqueModelSCOPids int, ModelClass text, "\ + "ModelFold text, ModelSuperfamily text, ModelFamily text, "\ + "CPUtime real)" self.mydb.execute(sqlxprs) # add columns labeling each model according to which PDB set it belongs to for tagset in PDBTagsets: sqlxprs = "alter table MRModelTbl add column %s int" %tagset[1] self.mydb.execute(sqlxprs) # then create a separate table of annotations to accommodate all the solutions returned from an # MR calculation. We set MRModel_id to the corresponding p_id of the MRModelTbl table sqlxprs = "create table AnnotationTbl (MRTarget_id int, MRModel_id int, Tstar int, " \ + "RFZ real, TFZ real, PAK real, LLGfullres_vrms real, iLLGvrms real, " \ + "LLGreflection real, LLGrefl_vrms real, iLLG real, CCglobal real, "\ + "LLGvrms real, LLGfullres_refl_vrms real, "\ + "VRMSfullres real, FracvarVrmsfullres real)" self.mydb.execute(sqlxprs) # Likewise create table for template RMS to accommodate for as many ensembles as present in the calculation sqlxprs = "create table TemplateRMSTbl (MRTarget_id int, MRModel_id int, CCglobal real, " \ + "RefinedModel text, LLGvrms real, iLLGtemplate real, LLGreflvrms real)" self.mydb.execute(sqlxprs) # Make an ELLG table for the 500 calculations with different ELLG target values sqlxprs = "create table ELLGtbl (MRModel_id int, eLLGtarget real, achievedELLG real, "\ + "LLGachievedres real, LLGfullres real, TFZfullres real, TFZachievedres real, "\ + "foundit int, foundit2 int, vrmsachievedres real, achievedresol real)" self.mydb.execute(sqlxprs) ncomponents = len(psearches[0][1]) #for j in range(1): for j in range(ncomponents): if createtable: sqlxprs = "alter table MRModelTbl add column RFZ%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column TFZ%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column PAK%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column LLG%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column LLGreflection%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column iLLG%d real" %j self.mydb.execute(sqlxprs) self.mydb.execute("alter table MRModelTbl add column iLLGfullres%d real" %j) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN CCmap%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN CCpdb%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN MLAD%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN LLGensmbl%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN VRMS%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN Bfac%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE AnnotationTbl ADD COLUMN FracvarVrms%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN CCmap%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN CCpdb%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN MLAD%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN LLGensmbl%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN VRMS%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN Bfac%d real" %j self.mydb.execute(sqlxprs) sqlxprs = "ALTER TABLE TemplateRMSTbl ADD COLUMN FracvarVrms%d real" %j self.mydb.execute(sqlxprs) if createtable and len(OtherSolutions) > 1: for s in OtherSolutions[1:]: # first element is empty as that's for the original Strider calculations suffix = s[0] sqlxprs = "alter table MRModelTbl add column TFZequiv%s real" %suffix self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column TFZ0%s real" %suffix self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column iLLG%s real" %suffix self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column Foundit%s real" %suffix self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column SolutionRank%s int" %suffix self.mydb.execute(sqlxprs) sqlxprs = "alter table MRModelTbl add column CCglobal%s real" %suffix self.mydb.execute(sqlxprs) completedrunfname = posixpath.join(self.topfolder, self.calclabel + self.targetpdbid + "_DoneMRJobs.txt") completedruns = open(completedrunfname,"r") doneruns = completedruns.read() completedruns.close() testsol = BLASTMRutils.BLASTMRutils(self.topfolder, self) datares = self.GetResolution(self.targetpdbid) WilsonB = [-1.0] WilsonScale = [-1.0] MatthewsCoef = [-1.0] TargetScattering = [-1.0] lastMRTargetTbl_id = [1] if OnlyCreateDBcolumns == False: os.chdir(posixpath.join(self.topfolder, "Phaserinput")) seqfiles = glob.glob("*.seq") ntargetresidues = 0 for seqfile in seqfiles: seqtxt = open(seqfile,"r").read() ntargetresidues += SimpleFileProperties.GetNumberofResiduesFromSequence(seqtxt) os.chdir(olddir) for s in OtherSolutions: # first element is empty as that's for the original Strider calculations OtherSolfname = s[1] suffix = s[0] if OtherSolfname =="": SetOfResults = set([]) for search in psearches: # create modelsdir as concatenated PDB codes separated by commas previousmodelscomponentstabs = "" modelsdir = [ [e[0] for e in e1] for e1 in search[:-1]][0] workdir = ",".join( modelsdir ) for (i,partmodels) in enumerate(search[1]): if not partmodels in SetOfResults: # don't bother writing results for the same partmodels more than once SetOfResults.add(partmodels) #print "partmodels= " + str(partmodels) bsolved = [False] founddir = [""] if testsol.PhaserRanPartModelsBefore(founddir,partmodels,bsolved, doneruns): # array that'll hold MR solution, template solution, FTF_scrores data, VRMS and LLG float values try: if OtherSolfname != "": nsolvedit = testsol.IsthereaPartialSolution(founddir[0], partmodels, OtherSolfname) else: nsolvedit = testsol.IsthereaPartialSolution(founddir[0], partmodels) testsol.mrsol #print "testsol.mrsol[3][6] = " +str(testsol.mrsol[3][6]) self.GatherResults(search, workdir, partmodels, i, testsol.mrsol, suffix, datares, unitcell, spgname, spgnumber, xtalsysname, ispolarspacegroup, nsolvedit, PDBTagsets, WilsonB, WilsonScale, Rfactors, EDSreproduce, MatthewsCoef, TargetScattering, ntargetresidues, Allreflections, ASUvolume, lastMRTargetTbl_id,) except Exception, m: #self.logger.log(logging.WARN, m) self.logger.log(logging.ERROR, str(m) + "\n" \ + traceback.format_exc()) ncomponents = len(search[0][i][2]) # write database to file conn.commit() self.mydb.close() self.logger.log(logging.INFO, "\nResults written to " + DBfname) def BLASTnMRrun(self, maxseqid, minseqid, maxpdbs, njobs, removesimilarseqidentmodels): reposfile = open("BLASTMRcalculations.txt","a+") reposfile.write("%s, %s, in BLASTnMRrun\n" %(self.targetpdbid, self.calclabel)) reposfile.close() self.MakeMyfolder(self.topfolder) LOG_FILENAME = os.path.join(self.topfolder, self.calclabel + 'BLASTnMRrun.log') self.SetupLogfile(LOG_FILENAME) pdbs_seqids = self.BLASTnSetupMRruns(maxseqid, minseqid, maxpdbs, removesimilarseqidentmodels) minmodelorder = len(pdbs_seqids[0][1]) # only do MR on the last component using the template solution of the previous one self.RunMRauto(pdbs_seqids, njobs, minmodelorder) self.logger.log(logging.INFO,"Done with BlastnMRrun.") def RunMRautoFromlistFile(self, njobs): searchlistfname = self.calclabel + self.targetpdbid + "_BLASTchain_ids.txt" reposfile = open("BLASTMRcalculations.txt","a+") reposfile.write("%s, %s, in RunMRautoFromlistFile, pdb searchlist: %s\n" \ %(self.targetpdbid, self.calclabel, searchlistfname)) reposfile.close() self.MakeMyfolder(self.topfolder) LOG_FILENAME = os.path.join(self.topfolder, self.calclabel + 'RunMRautoFromlistFile.log') self.SetupLogfile(LOG_FILENAME) self.logger.log(logging.INFO, "RunMRautoFromlistFile, pdb: %s, pdb searchlist: %s, \ calclabel: %s\n" %(self.targetpdbid, searchlistfname, self.calclabel)) searchlistfname = os.path.join(self.topfolder, searchlistfname) searchlistfname = PosixSlashFilePath(searchlistfname) searchlistfname = RecaseFilename2(searchlistfname) pdbsseqidfile = open(searchlistfname,"r") txt = pdbsseqidfile.read() pdbsseqidfile.close() #cast text into the expected list format pdbs_seqids = eval(txt) # only do MR on the last component using the template solution of the previous one minmodelorder = len(pdbs_seqids[0][1]) self.RunMRauto(pdbs_seqids, njobs, minmodelorder) self.logger.log(logging.INFO,"Finished RunMRautoFromlistFile on %s" %self.targetpdbid)