Difference between revisions of "Python Example Scripts"

From Phaserwiki
(Translational NCS Analysis: Change setREFL to setREFL_F_SIGF)
(Anisotropy Correction: Change example script to include the new resharpen Boolean, set by default to True)
 
(6 intermediate revisions by 2 users not shown)
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_P3221.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==
 
Example script for automated structure solution of BETA-BLIP
 
Example script for automated structure solution of BETA-BLIP
<pre>
+
 
 
       #beta_blip_auto.py
 
       #beta_blip_auto.py
 
       from phaser import *
 
       from phaser import *
 
       i = InputMR_DAT()
 
       i = InputMR_DAT()
       i.setHKLI("beta_blip.mtz")
+
       i.setHKLI("beta_blip_P3221.mtz")
      i.setLABI_F_SIGF("Fobs","Sigma")
 
 
       i.setHIRES(6.0)
 
       i.setHIRES(6.0)
 
       i.setMUTE(True)
 
       i.setMUTE(True)
Line 50: Line 45:
 
       if r.Success():
 
       if r.Success():
 
         i = InputMR_AUTO()
 
         i = InputMR_AUTO()
         i.setSPAC_HALL(r.getSpaceGroupHall())
+
         i.setREFL_DATA(r.getREFL_DATA())
        i.setCELL6(r.getUnitCell())
 
        i.setREFL_F_SIGF(r.getMiller(),r.getFobs(),r.getSigFobs())
 
 
         i.setROOT("beta_blip_auto")
 
         i.setROOT("beta_blip_auto")
 
         i.addENSE_PDB_ID("beta","beta.pdb",1.0)
 
         i.addENSE_PDB_ID("beta","beta.pdb",1.0)
Line 76: Line 69:
 
         print "Job exit status FAILURE"
 
         print "Job exit status FAILURE"
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
</pre>
 
  
 
==Reading MTZ Files for Experimental Phasing==
 
==Reading MTZ Files for Experimental Phasing==
Line 195: Line 187:
 
Example script script for anisotropy correction of BETA-BLIP data
 
Example script script for anisotropy correction of BETA-BLIP data
 
        
 
        
<pre>
 
 
       #beta_blip_ano.py
 
       #beta_blip_ano.py
 
       from phaser import *
 
       from phaser import *
 
       i = InputMR_DAT()
 
       i = InputMR_DAT()
       HKLIn = "beta_blip.mtz"
+
       HKLIn = "beta_blip_P3221.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 209: Line 197:
 
         i = InputANO()
 
         i = InputANO()
 
         i.setSPAC_HALL(r.getSpaceGroupHall())
 
         i.setSPAC_HALL(r.getSpaceGroupHall())
         i.setCELL6(r.getUnitCell())
+
         i.setREFL_DATA(r.getREFL_DATA())
        i.setREFL_F_SIGF(r.getMiller(),r.getFobs(),r.getSigFobs())
 
        i.setREFL_ID(F,SIGF)
 
        i.setHKLI(HKLIn)
 
 
         i.setROOT("beta_blip_ano")
 
         i.setROOT("beta_blip_ano")
 
         i.setMUTE(True)
 
         i.setMUTE(True)
Line 231: Line 216:
 
           f = r.getF()
 
           f = r.getF()
 
           sigf = r.getSIGF()
 
           sigf = r.getSIGF()
           f_iso = r.getCorrectedF()
+
           f_iso = r.getCorrectedF(True)
           sigf_iso = r.getCorrectedSIGF()
+
           sigf_iso = r.getCorrectedSIGF(True)
           corr = r.getCorrection()
+
           corr = r.getCorrection(True)
 
           nrefl = min(10,hkl.size())
 
           nrefl = min(10,hkl.size())
 
           print "First ", nrefl , " reflections"
 
           print "First ", nrefl , " reflections"
 
           print "%4s %4s %4s %10s %10s %10s %10s %10s" % \
 
           print "%4s %4s %4s %10s %10s %10s %10s %10s" % \
             ("H","K","L",F,SIGF,r.getLaboutF(),r.getLaboutSIGF(),"Corr\'n")
+
             ("H","K","L","F","SIGF",r.getLaboutF(),r.getLaboutSIGF(),"Corr\'n")
 
           for i in range(0,nrefl):
 
           for i in range(0,nrefl):
 
             print "%4d %4d %4d %10.4f %10.4f %10.4f %10.4f %10.4f" % \
 
             print "%4d %4d %4d %10.4f %10.4f %10.4f %10.4f %10.4f" % \
Line 247: Line 232:
 
         print "Job exit status FAILURE"
 
         print "Job exit status FAILURE"
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
 
         print r.ErrorName(), "ERROR :", r.ErrorMessage()
</pre>
 
  
 
==Cell Content Analysis==
 
==Cell Content Analysis==
Line 330: Line 314:
 
       i.setROOT("beta_nma")
 
       i.setROOT("beta_nma")
 
       i.addENSE_PDB_ID("beta","beta.pdb",1.0)
 
       i.addENSE_PDB_ID("beta","beta.pdb",1.0)
       i.setNMAP_MODES([7,10])
+
       i.addNMA_MODE(4)
       i.setNMAP_FORWARD()
+
       i.setPERT_RMS_DIRE("FORWARD")
 
       i.setMUTE(True)
 
       i.setMUTE(True)
       r = runNMA(i)
+
       r = runNMAXYZ(i)
 
       if r.Success():
 
       if r.Success():
 
         print "Normal Mode Analysis"
 
         print "Normal Mode Analysis"

Latest revision as of 10:14, 12 December 2017

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_P3221.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_P3221.mtz")
     i.setHIRES(6.0)
     i.setMUTE(True)
     r = runMR_DAT(i)
     if r.Success():
       i = InputMR_AUTO()
       i.setREFL_DATA(r.getREFL_DATA())
       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_P3221.mtz"
     i.setHKLI(HKLIn)
     i.setMUTE(True)
     r = runMR_DAT(i)
     if r.Success():
       i = InputANO()
       i.setSPAC_HALL(r.getSpaceGroupHall())
       i.setREFL_DATA(r.getREFL_DATA())
       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(True)
         sigf_iso = r.getCorrectedSIGF(True)
         corr = r.getCorrection(True)
         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)