Python Example Scripts

From Phaserwiki
Revision as of 11:34, 16 May 2012 by Airlie (talk | contribs) (Translational NCS Analysis)

Example scripts for the most popular modes of running Phaser.

Reading MTZ Files for Molecular Replacement

Example script for reading data from MTZ file beta_blip.mtz.

Note that by default reflections are sorted into resolution order upon reading, to achieve a performance gain in the molecular replacement routines. If reflections are not being read from an MTZ file with this script, reflections should be pre-sorted into resolution order to achieve the same performance gain. Sorting is turned off with the setSORT(False) function.

      #beta_blip_data.py
      from phaser import *
      i = InputMR_DAT()
      HKLIN = "beta_blip.mtz"
      F = "Fobs"
      SIGF = "Sigma"
      i.setHKLI(HKLIN)
      i.setLABI(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      print r.logfile()
      if r.Success():
        hkl = r.getMiller()
        fobs = r.getF()
        sigma = r.getSIGF()
        nrefl = min(10,hkl.size())
        print "Data read from: " , HKLIN
        print "Spacegroup Name (Hall symbol) = %s (%s)" % \
          (r.getSpaceGroupName(), r.getSpaceGroupHall())
        print "Unitcell = " , r.getUnitCell()
        print "First ", nrefl , " reflections"
        print "%4s %4s %4s %10s %10s" % ("H","K","L",F,SIGF)
        for i in range(0,nrefl):
          print "%4d %4d %4d %10.4f %10.4f" % \
            (hkl[i][0],hkl[i][1],hkl[i][2],fobs[i],sigma[i])
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()


Automated Molecular Replacement

Example script for automated structure solution of BETA-BLIP

      beta_blip_auto.py
      from phaser import *
      i = InputMR_DAT()
      i.setHKLI("beta_blip.mtz")
      i.setLABI("Fobs","Sigma")
      i.setHIRES(6.0)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputMR_AUTO()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL(r.getUnitCell())
        i.setREFL(r.getMiller(),r.getFobs(),r.getSigFobs())
        i.setROOT("beta_blip_auto")
        i.addENSE_PDB_ID("beta","beta.pdb",1.0)
        i.addENSE_PDB_ID("blip","blip.pdb",1.0)
        i.addCOMP_PROT_MW_NUM(28853,1)
        i.addCOMP_PROT_MW_NUM(17522,1)
        i.addSEAR_ENSE_NUM("beta",1)
        i.addSEAR_ENSE_NUM("blip",1)
        i.setMUTE(True)
        del(r)
        r = runMR_AUTO(i)
        if r.Success():
          if r.foundSolutions() :
            print "Phaser has found MR solutions"
            print "Top LLG = %f" % r.getTopLLG()
            print "Top PDB file = %s" % r.getTopPdbFile()
            print "Spacegroup Name (Hall symbol) = %s (%s)" % \
              (r.getSpaceGroupName(), r.getSpaceGroupHall())
          else:
            print "Phaser has not found any MR solutions"
        else:
          print "Job exit status FAILURE"
          print r.ErrorName(), "ERROR :", r.ErrorMessage()
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Reading MTZ Files for Experimental Phasing

Example script for reading SAD data from MTZ file S-insulin.mtz

      insulin_data.py
      from phaser import *
      i = InputEP_DAT()
      HKLIN = "S-insulin.mtz"
      xtalid = "insulin"
      waveid = "cuka"
      i.setHKLI(HKLIN)
      i.addCRYS_ANOM_LABI(xtalid,waveid,"F(+)","SIGF(+)","F(-)","SIGF(-)")
      i.setMUTE(True)
      r = runEP_DAT(i)
      if r.Success():
        hkl = r.getMiller()
        Fpos = r.getFpos(xtalid,waveid)
        Ppos = r.getPpos(xtalid,waveid)
        Fneg = r.getFneg(xtalid,waveid)
        Pneg = r.getPneg(xtalid,waveid)
        print "Data read from: " , HKLIN
        print "Spacegroup Name (Hall symbol) = %s (%s)" % \
         (r.getSpaceGroupName(), r.getSpaceGroupHall())
        print "Unitcell = " , r.getUnitCell()
        nrefl = min(10,hkl.size())
        print "First ", nrefl , " reflections with anomalous differences"
        print "%4s %4s %4s %10s %10s %10s" % ("H","K","L","F(+)","F(-)","D")
        i = 0
        r = 0
        while r < nrefl:
          if Ppos[i] and Pneg[i] :
            D = abs(Fpos[i]-Fneg[i])
            if D > 0
              print "%4d %4d %4d %10.4f %10.4f %10.4f" % \
                (hkl[i][0],hkl[i][1],hkl[i][2],Fpos[i],Fneg[i],D)
              r=r+1
          i=i+1
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Automated Experimental Phasing

Example script for SAD phasing for insulin

      #insulin_sad.py
      from phaser import *
      from cctbx import xray
      i = InputEP_DAT()
      HKLIN = "S-insulin.mtz"
      xtalid = "insulin"
      waveid = "cuka"
      i.setHKLI(HKLIN)
      i.addCRYS_ANOM_LABI(xtalid,waveid,"F(+)","SIGF(+)","F(-)","SIGF(-)")
      i.setMUTE(True)
      r = runEP_DAT(i)
      if r.Success():
        hkl = r.getMiller()
        Fpos = r.getFpos(xtalid,waveid)
        Spos = r.getSIGFpos(xtalid,waveid)
        Ppos = r.getPpos(xtalid,waveid)
        Fneg = r.getFneg(xtalid,waveid)
        Sneg = r.getSIGFneg(xtalid,waveid)
        Pneg = r.getPneg(xtalid,waveid)
        i = InputEP_AUTO()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL(r.getUnitCell())
        i.setCRYS_MILLER(hkl)
        i.addCRYS_ANOM_DATA(xtalid,waveid,Fpos,Spos,Ppos,Fneg,Sneg,Pneg)
        i.setATOM_PDB(xtalid,"S-insulin_hyss.pdb")
        i.setLLGC_CRYS_COMPLETE(xtalid,True)
        i.addLLGC_CRYS_SCAT_ELEMENT(xtalid,"S")
        i.addCOMP_PROT_FASTA_NUM("S-insulin.seq",1.)
        i.setHKLO(False)
        i.setSCRI(False)
        i.setXYZO(False)
        i.setMUTE(True)
        r = runEP_AUTO(i)
        if r.Success():
          print "SAD phasing"
          print "Data read from: " , HKLIN
          print "Data output to : " , r.getMtzFile()
          print "Spacegroup Name (Hall symbol) = %s (%s)" % \
            (r.getSpaceGroupName(), r.getSpaceGroupHall())
          print "Unitcell = " , r.getUnitCell()
          print "LogLikelihood = " , r.getLogLikelihood()
          atom = r.getAtoms(xtalid)
          print atom.size(), " refined atoms"
          print "%5s %10s %10s %10s %10s %10s" % \
            ("atom","x","y","z","occupancy","u-iso")
          for i in range(0,atom.size()):
            print "%5s %10.4f %10.4f %10.4f %10.4f %10.4f" % \
           
(atom[i].scattering_type,atom[i].site[0],atom[i].site[1],atom[i].site[2],atom[i].occupancy,atom[i].u_iso)
          hkl = r.getMiller();
          fwt = r.getFWT()
          phwt = r.getPHWT()
          fom = r.getFOM()
          nrefl = min(10,hkl.size())
          print "First ", nrefl , " reflections"
          print "%4s %4s %4s %10s %10s %10s" % \
            ("H","K","L","FWT","PHWT","FOM")
          for i in range(0,nrefl):
            print "%4d %4d %4d %10.4f %10.4f %10.4f" % \
              (hkl[i][0],hkl[i][1],hkl[i][2],fwt[i],phwt[i],fom[i])
        else:
          print "Job exit status FAILURE"
          print r.ErrorName(), "ERROR :", r.ErrorMessage()
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Anisotropy Correction

Example script script for anisotropy correction of BETA-BLIP data

      #beta_blip_ano.py
      from phaser import *
      i = InputMR_DAT()
      HKLIn = "beta_blip.mtz"
      F = "Fobs"
      SIGF = "Sigma"
      i.setHKLI(HKLIn)
      i.setLABI(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputANO()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL(r.getUnitCell())
        i.setREFL(r.getMiller(),r.getFobs(),r.getSigFobs())
        i.setREFL_ID(F,SIGF)
        i.setHKLI(HKLIn)
        i.setROOT("beta_blip_ano")
        i.setMUTE(True)
        del(r)
        r = runANO(i)
        if r.Success():
          print "Anisotropy Correction"
          print "Data read from: " , HKLIn
          print "Data output to : " , r.getMtzFile()
          print "Spacegroup Name (Hall symbol) = %s (%s)" % \
            (r.getSpaceGroupName(), r.getSpaceGroupHall())
          print "Unitcell = " , r.getUnitCell()
          print "Principal components = " , r.getEigenBs()
          print "Range of principal components = " , r.getAnisoDeltaB()
          print "Wilson Scale = " , r.getWilsonK()
          print "Wilson B-factor = " , r.getWilsonB()
          hkl = r.getMiller();
          f = r.getF()
          sigf = r.getSIGF()
          f_iso = r.getCorrectedF()
          sigf_iso = r.getCorrectedSIGF()
          corr = r.getCorrection()
          nrefl = min(10,hkl.size())
          print "First ", nrefl , " reflections"
          print "%4s %4s %4s %10s %10s %10s %10s %10s" % \
            ("H","K","L",F,SIGF,r.getLaboutF(),r.getLaboutSIGF(),"Corr\'n")
          for i in range(0,nrefl):
            print "%4d %4d %4d %10.4f %10.4f %10.4f %10.4f %10.4f" % \
              (hkl[i][0],hkl[i][1],hkl[i][2],f[i],sigf[i],f_iso[i],sigf_iso[i],corr[i])
        else:
          print "Job exit status FAILURE"
          print r.ErrorName(), "ERROR :", r.ErrorMessage()
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Cell Content Analysis

Example script for cell content analysis of BETA-BLIP

      #beta_blip_cca.py
      from phaser import *
      i = InputMR_DAT()
      HKLIN = "beta_blip.mtz"
      F = "Fobs"
      SIGF = "Sigma"
      i.setHKLI(HKLIN)
      i.setLABI(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputCCA()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL(r.getUnitCell())
        i.addCOMP_PROT_MW_NUM(28853,1)
        i.addCOMP_PROT_MW_NUM(17522,1)
        i.setMUTE(True)
        del(r)
        r = runCCA(i)
        if r.Success():
          print "Cell Content Analysis"
          print "Molecular weight of assembly = " , r.getAssemblyMW()
          print "Best Z value = " , r.getBestZ()
          print "Best VM value = " , r.getBestVM()
          print "Probability of Best VM = " , r.getBestProb()
        else:
          print "Job exit status FAILURE"
          print r.ErrorName(), "ERROR :", r.ErrorMessage()
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Translational NCS Analysis

Example script for translational NCS analysis of BETA-BLIP

      #beta_blip_ncs.py
      from phaser import *
      i = InputMR_DAT()
      HKLIN = "beta_blip.mtz"
      F = "Fobs"
      SIGF = "Sigma"
      i.setHKLI(HKLIN)
      i.setLABI(F,SIGF)
      i.setMUTE(True)
      i.setREFL(r.getMiller(),r.getF(),r.getSIGF())
      r = runMR_DAT(i)
      if r.Success():
        i = InputNCS()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL6(r.getUnitCell())
        i.addCOMP_PROT_MW_NUM(28853,1)
        i.addCOMP_PROT_MW_NUM(17522,1)
        i.setMUTE(True)
        del(r)
        r = runNCS(i)
        if r.Success():
          print "Translational NCS analysis"
          print "Translational NCS present = ", r.hasTNCS()
          if r.hasTNCS():
            print "Translational NCS vecor = ", r.hasTNCS()
          print "Twinning alpha = ", r.getTwinAlpha()
        else:
          print "Job exit status FAILURE"
          print r.ErrorName(), "ERROR :", r.ErrorMessage()
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Normal Mode Analysis

Example script for normal mode analysis of BETA-BLIP. Note that the space group and unit cell are not required, and so the MTZ file does not need to be read to extract these parameters.

      #beta_nma.py
      from phaser import *
      i = InputNMA()
      i.setROOT("beta_nma")
      i.addENSE_PDB_ID("beta","beta.pdb",1.0)
      i.setNMAP_MODES([7,10])
      i.setNMAP_FORWARD()
      i.setMUTE(True)
      r = runNMA(i)
      if r.Success():
        print "Normal Mode Analysis"
        for i in range(0,r.getNum()):
          print "PDB file = ", r.getPdbFile(i)
          displacement = r.getDisplacements(i)
          mode = r.getModes(i)
          for j in range(0,mode.size()):
            print " Mode = " , mode[j], " Displacement = ", displacement[j]
      else:
        print "Job exit status FAILURE"
        print r.ErrorName(), "ERROR :", r.ErrorMessage()

Logfile Handling

Example of how to redirect phaser output to a python string for real-time viewing of output, but not via standard output. Output to standard out is silenced with setMUTE(True).

      beta_blip_logfile.py
      from phaser import *
      from cStringIO import StringIO
      i = InputMR_DAT()
      i.setHKLI("beta_blip.mtz")
      i.setLABI("Fobs","Sigma")
      i.setMUTE(True)
      o = Output()
      redirect_str = StringIO()
      o.setPackagePhenix(file_object=redirect_str)
      r = runMR_DAT(i,o)