Talk:Poisson–Boltzmann profile for an ion channel

Notes for Oxonians

Latest comment: 16 years ago2 comments1 person in discussion

General announcements

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

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

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

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

x/Å	y/Å	z/Å	energy/kJ·mol⁻¹	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

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

Latest comment: 15 years ago4 comments1 person in discussion

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

placeion.py

#!/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

#!/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()