ProgramsΒΆ
Here we include some programs we developed for analysis of lipid bilayers. Those programs were used in computations for article DMG-alpha - a computational geometry package for multi-molecular systems (submitted).
This script computes volume per lipid and approximate solvent accessible surface of the surface of lipid bilayer.
from pydmga.geometry import OrthogonalGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.io import pdb
from math import sqrt
# this is basic filename of the dataset, withou .pdb extension
filename = "drec-1_590-600ns"
# use heavy atoms only - if True then we skip all hydrogen atoms
hao = True
# old lipid stops at what residue?
end_lipid_at = "UDC"
# new lipid starts at what residue?
start_lipid_at = "KDO"
# we will read from <filename>.pdb
in_file = file("{0}.pdb".format(filename), "r")
# we will write to <filename>.dat
out_file = file("{0}.dat".format(filename), "w")
try:
# this counts current frame, for output
frame = 0
# we collect volume and area per lipid per frame
out_file.write("frame lipid volume area\n")
while True: # iterate over all frames
# here we remember data from pdb for future use
particles = []
while True: # iterate over elements in frame
item = pdb.PDBRecord(pdb.readnext(in_file))
if item.type() == "CRYST1":
# create box geometry periodic in all directions
geometry = OrthogonalGeometry(item["a"], item["b"], item["c"], True, True, True)
# create container for atoms
container = Container(geometry)
if item.type() == "ATOM":
# add new item, if we are using only hao, then skip it if its hydrogen
if (item["atom"][0] != 'H' and item["atom"][1] != 'H') or not hao:
# add perticle to the container
container.add(item.as_particle())
# remember all data for future use
particles.append(item)
if item.type() == "ENDMDL" or item.type() == "TER":
# stop reading this frame
break
# create new Power Diagram
diagram = Diagram(container)
# here we will count which lipid molecule we have
molecule_no = 0
# this is to allow triggering start_lipid_at test for the first iteration
residue = end_lipid_at
#iterate over all cells in this diagram
for i, cell in enumerate(diagram):
# if we found one solvent then skip the rest of atoms - they are solvent
if particles[i]["residue"] == "SOL" or particles[i]["residue"] == "ION":
break
# we change lipid molecule only if we pass from <end_lipid_at> residue to <start_lipid_at> residue
if residue == end_lipid_at and particles[i]["residue"] == start_lipid_at:
# write data to a file if there is a data
if molecule_no: out_file.write("{0} {1} {2} {3}\n".format(frame, molecule_no, volume, area))
molecule_no +=1
area = 0.0
volume = 0.0
# to track when to change residue
residue = particles[i]["residue"]
# find cell volume
volume += cell.volume()
# iterate over all sides to find contact area for this atom
for side in cell.sides:
neighbour = side.neighbour
if particles[neighbour]['residue'] == "SOL" or particles[neighbour]['residue'] == "ION":
# we have solvent <-> particle contact
area += side.area()
# need to write last lipid (we changed frames, not lpipid inside the frame)
out_file.write("{0} {1} {2} {3}\n".format(frame, molecule_no, volume, area))
# count new frame
frame += 1
except StopIteration:
# pdb file ends
pass
This script computes number of contacts lipid to solvent and lipid to ion.
from pydmga.geometry import OrthogonalGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.io import pdb
from math import sqrt
# this is basic filename of the dataset, withou .pdb extension
filename = "drec-1_590-600ns"
# old lipid stops at what residue?
end_lipid_at = "UDC"
# new lipid starts at what residue?
start_lipid_at = "KDO"
# we will read from <filename>.pdb
in_file = file("{0}.pdb".format(filename), "r")
# we will write to <filename>.dat
out_file = file("{0}.dat".format(filename), "w")
try:
# this counts current frame, for output
frame = 0
# we collect number of water in contact and number of ions in contact
out_file.write("frame lipid nb_h2o nb_ion\n")
while True: # iterate over all frames
# here we remember data from pdb for future use
particles = []
while True: # iterate over elements in frame
item = pdb.PDBRecord(pdb.readnext(in_file))
if item.type() == "CRYST1":
# create box geometry periodic in all directions
geometry = OrthogonalGeometry(item["a"], item["b"], item["c"], True, True, True)
# create container for atoms
container = Container(geometry)
if item.type() == "ATOM":
# add new item, for simplicity we use only heavy atoms
if (item["atom"][0] != 'H' and item["atom"][1] != 'H'):
# add perticle to the container
container.add(item.as_particle())
# remember all data for future use
particles.append(item)
if item.type() == "ENDMDL" or item.type() == "TER":
# stop reading this frame
break
# create new Power Diagram
diagram = Diagram(container)
# here we will count which lipid molecule we have
molecule_no = 0
# this is to allow triggering start_lipid_at test for the first iteration
residue = end_lipid_at
#iterate over all cells in this diagram
for i, cell in enumerate(diagram):
# break from iteration when SOL/ION part of file reached
if particles[i]["residue"] == "SOL" or particles[i]["residue"] == "ION": break
# we change lipid molecule only if we pass from <end_lipid_at> residue to <start_lipid_at> residue
if residue == end_lipid_at and particles[i]["residue"] == start_lipid_at:
# write data to a file if there is a data
if molecule_no: out_file.write("{0} {1} {2} {3}\n".format(frame, molecule_no, len(nb_h2o), len(nb_ion)))
molecule_no +=1
# here we hold water particles in contact
nb_h2o = []
# here we hold ions in contact
nb_ion = []
# to track when to change residue
residue = particles[i]["residue"]
# iterate over all sides to find neighbours
for side in cell.sides:
n = side.neighbour
if particles[n]['residue'] == "SOL" and not n in nb_h2o: nb_h2o.append(n)
if particles[n]['residue'] == "ION" and not n in nb_ion: nb_ion.append(n)
# need to write last lipid (we changed frames, not lpipid inside the frame)
out_file.write("{0} {1} {2} {3}\n".format(frame, molecule_no, len(nb_h2o), len(nb_ion)))
# count new frame
frame += 1
except StopIteration:
# pdb file ends
pass
This script computes 2D voronoi diagram on the plane for all atoms in GN2 groups and cumulates them to get area per lipid.
from pydmga.geometry import OrthogonalGeometry
from pydmga.geometry import CastOrthoPlaneGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.io import pdb
filename = "drec-1_590-600ns"
groups = ["GN2"]
end_lipid_at = "UDC"
start_lipid_at = "KDO"
in_file = file("{0}.pdb".format(filename), "r")
out_file = file("{0}.dat".format(filename), "w")
try:
frame = 0
out_file.write("frame lipid sa\n")
while True:
lipid_by_atom = []
residue = ''
molecule_id = 0
while True:
item = pdb.PDBRecord(pdb.readnext(in_file))
if item.type() == "CRYST1":
(a, b, c) = (item["a"], item["b"], item["c"])
box_geometry = OrthogonalGeometry(a, b, c, True, True, True)
geometry_up = CastOrthoPlaneGeometry(box_geometry, (0, 0, 1))
container_up = Container(geometry_up)
if item.type() == "ATOM":
if residue == end_lipid_at and item['residue'] == start_lipid_at:
molecule_id += 1
if item["residue"] in groups and item["atom"][0] != 'H' and item["atom"][1] != 'H':
if item['z'] > c / 2:
container_up.add(item.as_particle())
lipid_by_atom.append(molecule_id)
residue = item["residue"]
if item.type() == "ENDMDL" or item.type() == "TER":
break
diagram_up = Diagram(container_up)
all_area = 0.0
molecule_id = lipid_by_atom[0]
sa = 0.0
for i, cell in enumerate(diagram_up):
if molecule_id != lipid_by_atom[i]:
out_file.write("{0} {1} {2}\n".format(frame, molecule_id, sa))
sa = 0.0
molecule_id = lipid_by_atom[i]
for side in cell.sides:
neighbour = side.neighbour
if geometry_up.on_boundary(neighbour):
sa += side.area()
out_file.write("{0} {1} {2}\n".format(frame, molecule_id, sa))
frame += 1
except StopIteration:
pass
This script computes 2D voronoi diagram on the plane for geometric centers of GN2 groups.
from pydmga.geometry import OrthogonalGeometry
from pydmga.geometry import CastOrthoPlaneGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.io import pdb
filename = "drec-1_590-600ns"
gc = False
groups = ["GN2"]
end_lipid_at = "UDC"
start_lipid_at = "KDO"
in_file = file("{0}.pdb".format(filename), "r")
out_file = file("{0}.dat".format(filename), "w")
try:
frame = 0
out_file.write("frame lipid sa\n")
while True:
lipid_by_atom = []
residue = ''
this_group = []
molecule_id = 0
while True:
item = pdb.PDBRecord(pdb.readnext(in_file))
if item.type() == "CRYST1":
(a, b, c) = (item["a"], item["b"], item["c"])
box_geometry = OrthogonalGeometry(a, b, c, True, True, True)
geometry_up = CastOrthoPlaneGeometry(box_geometry, (0, 0, 1))
container_up = Container(geometry_up)
if item.type() == "ATOM":
if residue == end_lipid_at and item['residue'] == start_lipid_at:
molecule_id += 1
x = sum([p[0] for p in this_group]) / len(this_group)
y = sum([p[1] for p in this_group]) / len(this_group)
z = sum([p[2] for p in this_group]) / len(this_group)
if z > c / 2:
container_up.add(container_up.size(), x, y, z, 1.0)
lipid_by_atom.append(molecule_id)
this_group = []
if item["residue"] in groups and item["atom"][0] != 'H' and item["atom"][1] != 'H':
this_group.append(item.as_coords())
residue = item["residue"]
if item.type() == "ENDMDL" or item.type() == "TER":
x = sum([p[0] for p in this_group]) / len(this_group)
y = sum([p[1] for p in this_group]) / len(this_group)
z = sum([p[2] for p in this_group]) / len(this_group)
if z > c / 2:
container_up.add(container_up.size(), x, y, z, 1.0)
lipid_by_atom.append(molecule_id)
break
diagram_up = Diagram(container_up)
for i, cell in enumerate(diagram_up):
molecule_id = lipid_by_atom[i]
for side in cell.sides:
neighbour = side.neighbour
if geometry_up.on_boundary(neighbour):
out_file.write("{0} {1} {2}\n".format(frame, molecule_id, side.area()))
frame += 1
except StopIteration:
pass
This program was used in the computation of water molecules trapped inside lipid bilayer in poster presentation at Bioinformatics Conference 2013, SocBIN/BIT13.
import sys
sys.path.append('../')
from pydmga.geometry import OrthogonalGeometry
from pydmga.geometry import CastSphereGeometry
from pydmga.geometry import CastCylinderGeometry
from pydmga.container import Container
from pydmga.diagram import Diagram
from pydmga.render import geometry
from pydmga.render import graphics
from pydmga.render import render
from random import random
#from dmga.io.pdb import PDBReader
solvent_from = -1
data = []
import fileinput
i = 0
cryst_line = ""
for line in fileinput.input():
tag = line[0:6].strip()
#print tag
if (tag == 'CRYST1'):
cryst_line = line
box_x = float(line[6:15].strip())
box_y = float(line[15:24].strip())
box_z = float(line[24:33].strip())
tilt_x = float(line[33:40].strip())
tilt_y = float(line[40:47].strip())
tilt_z = float(line[47:54].strip())
group = line[55:66]
z_value = int(line[66:70].strip())
#print "we have box:", (box_x, box_y, box_z), (tilt_x, tilt_y, tilt_z), group, z_value
if (tag == 'ATOM'):
id = int(line[6:11].strip())
atom_name = line[12:16]
alt_loc = line[16:17]
res_name = line[17:20]
chain_id = line[21:22]
seq_no = int(line[22:27])
x = float(line[30:38].strip())
y = float(line[38:46].strip())
z = float(line[46:54].strip())
r = float(line[54:60].strip())
# drop rest...
if (atom_name[1] != 'H'):
data.append( (i, x, y, z, r) )
if (solvent_from < 0 and res_name.strip() == 'SOL'):
solvent_from = i
i += 1
#print "We have atom:", id, atom_name, alt_loc, res_name, chain_id, seq_no, (x, y, z, r)
geometry = OrthogonalGeometry(box_x, box_y, box_z, True, True, True)
container = Container(geometry)
#print "Data lenght:", solvent_from, len(data)
container.add(data)
diagram = Diagram(container)
#display = render.Render()
colors = []
surface_by_part = []
surface_by_part_half = []
mark = [False for i in range(len(data)+1)]
print len(mark)
for i in range(solvent_from, len(data)):
if (mark[i]): continue
queue = [i]
#print "new part at", i
mark[i] = True
surface_by_part.append([])
surface_by_part_half.append([])
while(len(queue) > 0):
i = queue.pop()
cell = diagram.get_cell(i)
for side in cell.sides:
neighbour = side.neighbour
polygon = side.as_coords()
if (neighbour < solvent_from):
if (neighbour < 50 * 160):
surface_by_part_half[-1].append(0)
else:
surface_by_part_half[-1].append(1)
surface_by_part[-1].append(polygon) ##add new polygon to this part
else:
if (not(mark[neighbour])):
mark[neighbour] = True
queue.append(neighbour)
outer = surface_by_part[0]
surface_by_part[0] = []
surface_by_part.insert(0, [])
for i in range(len(outer)):
surface_by_part[surface_by_part_half[0][i]].append(outer[i])
import sys
filename = ".".join(sys.argv[1].split(".")[0:-1])
print "filename is", filename
pa_id = 100000000
#print "There are", len(surface_by_part), "parts of water in this system."
f = None
for i, part in enumerate(surface_by_part):
file_count = 0
for polygon in part:
if (pa_id + len(polygon) >= 9999):
if (f is not None): f.close()
file_count += 1
f = open("/home/robson/{0}_{1:03d}_{2:03d}.pdb".format(filename, i, file_count), "w")
f.write(cryst_line)
pa_id = 0
#print last vertex of polygon
to_connect = []
for v in polygon:
pa_id += 1
to_connect.append(pa_id)
if (i == 0):
if v[2] < 30.0:
v = (v[0], v[1], v[2] + box_z)
if (i == 1):
if v[2] > 30.0:
v = (v[0], v[1], v[2] - box_z)
f.write("HETATM{0: 5d} PS1 PSD P 1 {1: 8.3f}{2: 8.3f}{3: 8.3f} 0.00 0.00 PSDOPS \n".format(pa_id, v[0], v[1], v[2]))
to_connect.insert(0, pa_id)
for j in range(len(polygon)):
f.write("CONECT{0: 5d}{1: 5d}\n".format(to_connect[j], to_connect[j+1]))
pa_id = 100000000
if not f.closed:
f.close()
#set connect_cutoff,0.0
#set connect_mode,1
#display.add_outline_polygon(polygon, color)
#display.run()