Difference between revisions of "Python Example Scripts"

From Phaserwiki
(Normal Mode Analysis: Fix syntax for NMA, but this still crashes because the modes and the pdb files are not in one-to-one correspondence.)
(Reading MTZ Files for Molecular Replacement)
Line 9: Line 9:
  
 
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.
 
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.
<pre>
+
 
 
       #beta_blip_data.py
 
       #beta_blip_data.py
 
       from phaser import *
 
       from phaser import *
 
       i = InputMR_DAT()
 
       i = InputMR_DAT()
 
       HKLIN = "beta_blip.mtz"
 
       HKLIN = "beta_blip.mtz"
      F = "Fobs"
 
      SIGF = "Sigma"
 
 
       i.setHKLI(HKLIN)
 
       i.setHKLI(HKLIN)
      i.setLABI_F_SIGF(F,SIGF)
 
 
       i.setMUTE(True)
 
       i.setMUTE(True)
 
       r = runMR_DAT(i)
 
       r = runMR_DAT(i)
Line 35: Line 32:
 
         print "Job exit status FAILURE"
 
         print "Job exit status FAILURE"
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
</pre>
 
  
 
==Automated Molecular Replacement==
 
==Automated Molecular Replacement==

Revision as of 17:14, 25 November 2016

See Python Interface for introduction to running phaser as a python library.

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"
     i.setHKLI(HKLIN)
     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 "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_F_SIGF("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.setCELL6(r.getUnitCell())
        i.setREFL_F_SIGF(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()
          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(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputANO()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL6(r.getUnitCell())
        i.setREFL_F_SIGF(r.getMiller(),r.getFobs(),r.getSigFobs())
        i.setLABI_F_SIGF(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(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputCCA()
        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 = 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(F,SIGF)
      i.setMUTE(True)
      r = runMR_DAT(i)
      if r.Success():
        i = InputNCS()
        i.setSPAC_HALL(r.getSpaceGroupHall())
        i.setCELL6(r.getUnitCell())
        i.setREFL_F_SIGF(r.getMiller(),r.getF(),r.getSIGF())
        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.addNMA_MODE(4)
      i.setPERT_RMS_DIRE("FORWARD")
      i.setMUTE(True)
      r = runNMAXYZ(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_F_SIGF("Fobs","Sigma")
      i.setMUTE(True)
      o = Output()
      redirect_str = StringIO()
      o.setPackagePhenix(file_object=redirect_str)
      r = runMR_DAT(i,o)