Talk:Poisson–Boltzmann profile for an ion channel

From Wikiversity
Jump to navigation Jump to search

Notes for Oxonians[edit source]

General announcements[edit source]

This resource was used as part of the lecture and practical session on electrostatics calculations for the Wellcome Trust and Systems Biology Doctoral Training Centre postgraduate course on computational biochemistry, University of Oxford (Michaelmas 2007). The lecture is on WebLearn at http://tinyurl.com/282dfd (restricted access; a similar lecture without access restriction is in the resource near the top).

Feel free to go to lunch as soon as you have entered your result! I hope you all had fun. Please fill in the feedback form!

There will be no debrief session after lunch. But as the results come in to the table below, think about the possible sources of error and the general shape of the profile. See virtual debrief session below. To get in touch with me by email, look me up at University of Oxford contact search. – Kaihsu 12:40, 5 December 2007 (UTC)[reply]

APBS issues[edit source]

APBS is installed in /usr/apps/apbs-0.5.1-ia32/bin/apbs. You can make an alias for it in the bash shell:

alias apbs=/usr/apps/apbs-0.5.1-ia32/bin/apbs

Then just type apbs, followed by the parameters, to use it. Or, you can add it to your path (the environment variable PATH).

Use this input file: nAChR.in. Do not copy the input file from the PDF file: it is truncated!

PDB2PQR settings[edit source]

Use all the default settings on the PDB2PQR server: do not change anything! This ensures that your colleagues' results are comparable. In case PDB2PQR does not work, get this one from WebLearn: 1OED.pqr.

Results[edit source]

You might want to report your results here, in addition to plotting it with your colleagues:

x y z energy/kJ·mol−1 name
65 65 15 −0.018 Sumeet
65 65 30 −0.065
65 65 35 −1.2747 Phil
65 65 40 −0.357 Ayse
65 65 45 −0.5045 Maria
65 65 50 −1.205 Gareth
65 65 51.25 −1.114 Matt G
65 65 53.75 −2.880 Matt G
65 65 55 −3.33 Ahmed
65 65 60 0.859
65 65 65 23.76 Kit
65 65 66.25 26.82 Seamus
65 65 67.5 18.76 Kit
65 65 68.75 22.22 Seamus
65 65 70 20.41 Anna
65 65 72.5 28.60 Anna
65 65 75 25.45 Alex
65 65 80 20.00 Armin
65 65 82.5 21.94 Sumeet
65 65 85 24.15 Sumeet
65 65 87 6.36 Yaseen
65 65 90 −9.64 Marton
65 65 95 −12.8 Ben
65 65 100 −9.49 Ali
65 65 105 −3.79 Muhan
65 65 110 −1.837
65 65 115 −0.95 Christian
65 65 117.5 −0.84 Christian
65 65 120 −0.47 Christina
65 65 125

File:ACh.jpg

Virtual debrief session[edit source]

If someone is keen, they can post a plot of the results (using gnuplot, xmgrace, or any other plotting software) as they come in.

What are the possible sources of error in this profile? How do we estimate the error? Has anybody tried to (carefully and knowledgeably) change the grid-point numbers?

What are the implications of this profile? (To start, which side is extracellular and which cytoplasmic? 30 Å or 120 Å?) Think the distinction between thermodynamics and kinetics.

We have a set of scripts in our lab to automate the pore-scanning process in parallel on a compute cluster. Get in touch if you would like to take a look (and possibly improve it!).

Kaihsu 12:40, 5 December 2007 (UTC)[reply]

Code for automation[edit source]

I wrote some Python code to automate this process. The job submission requires a queuing system called Grid Engine. Copyright © 2008 Kaihsu Tai. Moral rights asserted (why?). Hereby licensed under either GFDL or GNU General Public License at your option.

Obviously, you will have to change the parameters in the driver-functions (main()) to fit the purposes of this tutorial. That is the exercise for the reader! – Kaihsu 16:12, 23 December 2008 (UTC)[reply]

Also, you will need to generate the PQR file, along with a file containing a list of the coordinates for the sample points. The script submits the job to a queueing system called Grid Engine, but you can submit the job by hand if you do not have this installed. – Kaihsu 11:04, 9 January 2009 (UTC)[reply]

Additionally, you will need to change the ion species. – Kaihsu 17:03, 30 January 2009 (UTC)[reply]

See also Two-dimensional Poisson–Boltzmann profile for an ion channel. – Kaihsu 13:07, 30 April 2009 (UTC)[reply]

placeion.py[edit source]

#!/usr/bin/python
import os

class placeion:
  "preparing job for APBS energy profiling by placing ions"

  def __init__(self):
    self.cglen = [0, 0, 0]
    self.pqrLines = []

  def readPQR(self):
    pqrFile = open(self.pqrName, "r")
    lines = pqrFile.readlines()
    for line in lines:
      if (line[0:4] == "ATOM"):
        self.pqrLines.append(line)
    pqrFile.close()

    # finding the cglen (coarse grid lengths) automatically
    # This depends on correct spacing in the PQR file.
    minDim = [+100000, +100000, +100000]
    maxDim = [-100000, -100000, -100000]
    for line in lines:
      if (line[0:4] == "ATOM"):
        tokens = line.split()
        for i in range(0, 3):
          if float(tokens[i+5]) < minDim[i]:
            minDim[i] = float(tokens[i+5])
          if float(tokens[i+5]) > maxDim[i]:
            maxDim[i] = float(tokens[i+5])
    for i in [0, 1]:
      self.cglen[i] = (maxDim[i] - minDim[i]) + 40
    self.cglen[2] = (maxDim[2] - minDim[2]) + 80
    
  def writePQRs(self):
    # make directory
    if (not os.path.exists(self.jobName)):
      os.mkdir(self.jobName)

    # write protein
    proFile = open(self.jobName + "/pro" + ".pqr", "w")
    for line in self.pqrLines:
      proFile.write(line)
    proFile.close()

    # write complex and ion
    for point in self.points:
      x = point[0]
      y = point[1]
      z = point[2]
      aID = 50000
      rID = 10000
      ionLines = ""
      for zOffset in self.zOffsets:
        # example: "ATOM  50000  NHX SPM 10000      43.624  57.177  58.408  1.0000 2.1300"
        ionLines += self.getPqrLine(aID, "NHX", rID, "SPM", x, y, z+zOffset, +1.0, 2.130)
        # entry for NH4(+) in Table III in
        # Rashin and Honig (1985) J. Phys. Chem. 89(26):5588--5593
        # http://dx.doi.org/10.1021/j100272a006
        aID += 1
        rID += 1

      # write ion
      ionFile = open(self.jobName + "/ion_" + str(z) + ".pqr", "w")
      ionFile.write(ionLines)
      ionFile.close()

      # write complex
      cpxFile = open(self.jobName + "/cpx_" + str(z) + ".pqr", "w")
      for line in self.pqrLines:
        cpxFile.write(line)
      cpxFile.write(ionLines)
      cpxFile.close()

  def getPqrLine(self, aID, aType, rID, rType, x, y, z, q, r):
    # example: "ATOM  50000  NHX SPM 10000      43.624  57.177  58.408  1.0000 2.1300"
    line =  "ATOM  "
    line += "%5d  " % aID
    line += "%-4s" % aType
    line += "%3s " % rType
    line += "%5d    " % rID
    line += "%8.3f%8.3f%8.3f " % (x, y, z)
    line += "%7.4f%7.4f\n" % (q, r)
    return line

  def readPoints(self):
    pointsFile = open(self.pointsName, "r")
    lines = pointsFile.readlines()
    pointsFile.close()

    points = []
    for line in lines:
      tokens = line.split()
      parsed = [float(tokens[0]), float(tokens[1]), float(tokens[2])]
      points.append(parsed)
    self.points = points

  def writeIn(self):
    for point in self.points:
      z = point[2]
      inFile = open(self.jobName + "/job_" + str(z) + ".in", "w")
      inFile.write("""# APBS input file generated by """ + __file__ + """
# for the complex and the ion at z = """ + str(z) + "\n") 

      inFile.write("""
  read
    mol pqr pro.pqr
    mol pqr ion_""" + str(z) + """.pqr
    mol pqr cpx_""" + str(z) + """.pqr
  end
""")

      count = 1
      for what in ["pro", "ion", "cpx"]:
        inFile.write("""
  elec name """ + what + """
    mg-auto
    mol """ + str(count) + """
    dime 97 97 193
    cglen """ + str(self.cglen[0]) + " " + str(self.cglen[1]) + " " + str(self.cglen[2]) + """
    fglen 20 20 20    
    cgcent mol 3
    fgcent mol 2
    # NaCl ionic strength in mol/L
    ion  1 """ + str(self.ionic) + """ 0.95 # sodium ions
    ion -1 """ + str(self.ionic) + """ 1.81 # chloride ions

    lpbe
    bcfl mdh
    pdie 2.0 # protein and faux-lipid
    sdie 78.5 # Eisenberg and Crothers Phys. Chem. book 1979
    srfm smol
    chgm spl2
    srad 1.4
    swin 0.3
    sdens 10.0
    temp 300
    # gamma 0.105 # Uncomment for old versions of APBS: deprecated for APBS 1.0.0
    calcenergy total
    calcforce no
  end
""")
        count += 1

      inFile.write("""
  print energy cpx - ion - pro end
  quit
""")
      inFile.close()

  def writeJob(self):
    for point in self.points:
      z = point[2]
      jobFile = open(self.jobName + "/job_" + str(z) + ".bash", "w")
      jobFile.write("#!/bin/bash\n")
      jobFile.write("""
#$ -N z""" + str(z) + self.jobName + """
#$ -S /bin/bash
#$ -l glibc=2.3,mem_free=500M,mem_total=500M
#$ -cwd
#$ -j y
#$ -r y
echo job running on $HOSTNAME
ulimit -c 64

/sansom/fedpacks/bin/apbs """)
      jobFile.write("job_" + str(z) + ".in > job_" + str(z) + ".out")
      jobFile.close()

    qsubName = "qsub_" + self.jobName + ".bash"
    jobFile = open(qsubName, "w")
    jobFile.write("""#$ -N """ + self.jobName + """
#$ -S /bin/bash
#$ -l glibc=2.3,mem_free=500M,mem_total=500M
#$ -cwd
#$ -j y
#$ -r y

declare -a job

""")
    count = 0
    for point in self.points:
      z = point[2]
      jobFile.write("job[" + str(count+1) + ']="' + self.jobName + "/job_" + str(z) + ".bash" + '"\n')
      count += 1
    jobFile.write("""
run_d=$(dirname ${job[${SGE_TASK_ID}]})
script=$(basename ${job[${SGE_TASK_ID}]})

cd ${run_d} || { echo "Failed to cd ${run_d}. Abort."; exit 1; }
chmod a+x ${script}

exec ./${script}
""")
    jobFile.close()
    return count

  def run(self):
    self.readPQR()
    self.readPoints()
    self.writePQRs()
    self.writeIn()
    count = self.writeJob()
    # qsub -t 1-80 qsub_all.sh
    os.spawnl(
      os.P_WAIT, "/sansom/fedpacks/opt/SGE6/bin/lx24-x86/qsub", "/sansom/fedpacks/opt/SGE6/bin/lx24-x86/qsub",
      "-t", "1-" + str(count), "qsub_" + self.jobName + ".bash"
    )

def main():
  myP = placeion()
  myP.pqrName = "/sansom/s66/kaihsu/works/oxford/physiome/magnesium/H011.pqr"
  myP.pointsName = "/sansom/s66/kaihsu/works/oxford/physiome/magnesium/samplepoints011.dat"
  myP.ionic = 0.15

  # in angstroms
  # set myP.zOffsets to be [-6.3, 0.0, +5.0] for 'butane'  end towards cytoplasm
  # set myP.zOffsets to be [-5.0, 0.0, +6.3] for 'propane' end towards cytoplasm
  myP.zOffsets = [-6.3, 0.0, +5.0] 
  myP.jobName = "H011cytoBut"

  myP.run()

if __name__ == "__main__":
  main()

analyze.py[edit source]

#!/usr/bin/python
import Gnuplot, os

class analyze:
  "analyze APBS energy profiling results"

  def readPoints(self):
    pointsFile = open(self.pointsName, "r")
    lines = pointsFile.readlines()
    pointsFile.close()

    points = []
    for line in lines:
      tokens = line.split()
      parsed = [float(tokens[0]), float(tokens[1]), float(tokens[2])]
      points.append(parsed)
    self.points = points

  def accumulate(self):
    self.readPoints()
    self.zE = []
    for point in self.points:
      z = point[2]
      outName = self.jobName + "/job_" + str(z) + ".out"
      lines = ""
      if (os.path.exists(outName)):
        outFile = open(outName, "r")
        lines = outFile.readlines()
        outFile.close()
      for line in lines:
        if (line[0:18] == "  Local net energy"):
          self.zE.append([z, float(line.split()[6])])
    
  def write(self):
    outName = self.jobName + ".dat"
    outFile = open(outName, "w")
    outFile.write("# z/angstrom E/(kJ/mol)\n")
    for each in self.zE:
      outFile.write("%8.3f %8.3e\n" % (each[0], each[1]))
    outFile.close()

  def plot(self):
    outName = self.jobName + ".dat"
    plotName = self.jobName + ".eps"

    # initialize
    g = Gnuplot.GnuplotProcess()
    cmd  = "set terminal postscript eps colour\n"
    cmd += 'set output "' + plotName + '"\n'
    cmd += """set style data lines
set xlabel "z / nm"
set ylabel "energy / (kJ/mol)"
"""
    cmd += 'plot "' + outName + '" using ($1/10):($2) title "' + self.jobName + '"\n'

    # do it
    g(cmd)

  def run(self):
    self.accumulate()
    self.write()
    self.plot()

def main():
  myP = analyze()
  myP.pointsName = "/sansom/s66/kaihsu/works/oxford/physiome/magnesium/samplepoints011.dat"
  myP.jobName = "H011cytoBut"
  myP.run()

if __name__ == "__main__":
  main()