Meshing a protein using PyGAMer

This tutorial will walk you through the process of using PyGAMer to generate a tetrahedral mesh of the volume surrounding a protein. The protein of interest, PDBID 2JHO, is that of sperm whale myoglobin at 1.4Å resolution.

# Load in PyGAMer and numpy
import pygamer
print("PyGAMer version:", pygamer.__version__())
import numpy as np
# We will use threevis to visualize the mesh in the notebook
import threevis

PyGAMer version: v2.0.2-dev-16-g36333b7

# Mesh the protein of interest
mesh = pygamer.readPDB_molsurf('../data/2jho.pdb')

When using PyGAMer it’s important to initialize the orientation of the simplices. Notably compute_orientation() will try to apply a self consistent set of orientations. In many cases, you may wish to ensure that the mesh normals are outward facing. This can be achieved by calling correctNormals() which essentially computes the volume of the bounded surface and flips the normals if the volume is negative.

# Compute the normal orientation
components, orientable, manifold = mesh.compute_orientation()
mesh.correctNormals()
print(F"The mesh has {components} components, is"
      F" {'orientable' if orientable else 'non-orientable'}, and is"
      F" {'manifold' if manifold else 'non-manifold'}.")

The mesh has 2 components, is orientable, and is manifold.

Note that this mesh has two components. This is because there is some empty space inside of the protein which however is not solvent acessible. Therefore when the surface is extracted, this void space is encapsulated by a separate mesh. You can split the two surfaces accordingly as follows.

meshes = mesh.splitSurfaces()
for i, m in enumerate(meshes):
    print(F"Mesh {i} is {m.getVolume()} A^3 in volume.")
# Keep only the larger mesh
mesh = meshes[0]

Mesh 0 is 19248.321147885723 A^3 in volume.
Mesh 1 is 51.60539976197604 A^3 in volume.

# You can show the two meshes like follows:
protverts, protedges, protfaces = meshes[0].to_ndarray()
bverts, bedges, bfaces = meshes[1].to_ndarray()
ctx = threevis.Context(width=640, height=480)
colors = np.zeros((len(bfaces), 3))
colors[:,0] = np.ones(len(bfaces))
ctx.draw_faces(bverts, bfaces, colors=threevis.FaceAttribute(colors))
ctx.draw_edges(protverts, protedges)
ctx.display()

Conditioning the protein mesh

# Set selection of all vertices to True so smooth will operate on them.
for v in mesh.vertexIDs:
    v.data().selected = True
# Apply 5 iterations of smoothing
mesh.smooth(max_iter=5, preserve_ridges=True, verbose=True)

Initial Quality: Min Angle = 18.2858, Max Angle = 139.106, # smaller-than-15 = 0, # larger-than-165 = 0
Iteration 1:
Min Angle = 17.0361, Max Angle = 129.799, # smaller-than-15 = 0, # larger-than-165 = 0
Iteration 2:
Min Angle = 17.0849, Max Angle = 125.303, # smaller-than-15 = 0, # larger-than-165 = 0
Iteration 3:
Min Angle = 17.5223, Max Angle = 125.131, # smaller-than-15 = 0, # larger-than-165 = 0
Iteration 4:
Min Angle = 18.0895, Max Angle = 125.635, # smaller-than-15 = 0, # larger-than-165 = 0
Iteration 5:
Min Angle = 18.6719, Max Angle = 125.444, # smaller-than-15 = 0, # larger-than-165 = 0

print(F"The mesh currently has {mesh.nVertices} vertices, \
{mesh.nEdges} edges, and {mesh.nFaces} faces.")

The mesh currently has 42590 vertices, 127770 edges, and 85180 faces.

Iterative Coarsening

Depending on your problem, you may wish to coarsen or decimate the mesh. Basically reduce the number of vertices. GAMer offers several functions to help you with this goal.

for i in range(5):
    # Coarsen dense regions of the mesh
    mesh.coarse_dense(rate = 2, numiter = 3)
    # Coarsen flat regions of the mesh
    mesh.coarse_flat(rate = 0.1, numiter = 3)
    mesh.smooth(max_iter = 3, preserve_ridges = True, verbose = False)
    print(F"Iteration {i}: {mesh.nVertices} vertices, \
{mesh.nEdges} edges, and {mesh.nFaces} faces.")

Iteration 0: 23769 vertices, 71307 edges, and 47538 faces.
Iteration 1: 15651 vertices, 46953 edges, and 31302 faces.
Iteration 2: 11197 vertices, 33591 edges, and 22394 faces.
Iteration 3: 8217 vertices, 24651 edges, and 16434 faces.
Iteration 4: 6218 vertices, 18654 edges, and 12436 faces.

# Smooth the mesh again.
mesh.smooth(max_iter = 10, preserve_ridges = True, verbose = True)

Initial Quality: Min Angle = 12.2065, Max Angle = 153.945, # smaller-than-15 = 2, # larger-than-165 = 0
Iteration 1:
Min Angle = 11.3888, Max Angle = 155.785, # smaller-than-15 = 2, # larger-than-165 = 0
Iteration 2:
Min Angle = 10.9815, Max Angle = 156.645, # smaller-than-15 = 3, # larger-than-165 = 0
Iteration 3:
Min Angle = 10.9626, Max Angle = 156.641, # smaller-than-15 = 3, # larger-than-165 = 0
Iteration 4:
Min Angle = 11.2803, Max Angle = 155.911, # smaller-than-15 = 4, # larger-than-165 = 0
Iteration 5:
Min Angle = 11.8625, Max Angle = 155.123, # smaller-than-15 = 4, # larger-than-165 = 0
Iteration 6:
Min Angle = 11.6953, Max Angle = 156.412, # smaller-than-15 = 4, # larger-than-165 = 0
Iteration 7:
Min Angle = 10.9964, Max Angle = 157.257, # smaller-than-15 = 4, # larger-than-165 = 0
Iteration 8:
Min Angle = 10.5348, Max Angle = 157.765, # smaller-than-15 = 4, # larger-than-165 = 0
Iteration 9:
Min Angle = 10.2606, Max Angle = 158.013, # smaller-than-15 = 3, # larger-than-165 = 0
Iteration 10:
Min Angle = 10.1356, Max Angle = 158.053, # smaller-than-15 = 3, # larger-than-165 = 0

# Set boundary markers of the mesh to 23
for faceID in mesh.faceIDs:
    faceID.data().marker = 23

# Get the root metadata
gInfo = mesh.getRoot()
gInfo.ishole = True    # Don't mesh the inside of
gInfo.marker = -1

# Center mesh at 0,0,0
center, radius = mesh.getCenterRadius()
mesh.translate(-center)

Now let’s construct a bounding box to represent the cytosol around the protein

# Generate a surrounding box
box = pygamer.surfacemesh.cube(5)

# Set box boundary markers to 50
for faceID in box.faceIDs:
    faceID.data().marker = 50

# Get and set box metadata
gInfo = box.getRoot()
gInfo.ishole = False
gInfo.marker = 5

# Scale the box by 2 times the radius of the protein mesh
box.scale(radius*2)

# Draw the protein and surrounding box together
protverts, _, protfaces = mesh.to_ndarray()
boxverts, boxedges, _ = box.to_ndarray()
ctx = threevis.Context(width=640, height=480)
colors = np.zeros((len(protfaces), 3))
colors[:,0] = np.ones(len(protfaces))
ctx.draw_faces(protverts, protfaces, colors=threevis.FaceAttribute(colors))
ctx.draw_edges(boxverts, boxedges)
ctx.display()

You can write these meshes to file if you wish.

pygamer.writeOFF('2jho.off', mesh)

pygamer.writeOFF('box.off', box)

Tetrahedralizing

# Construct a list of meshes to pass into TetGen
meshes = [mesh, box]

# Call tetgen to tetrahedralize. The string represents command line arguments passed directly to TetGen.
tetmesh = pygamer.makeTetMesh(meshes, "q1.3/10O8/7AC")

Number of vertices: 7756
Number of Faces: 15508
Number of Regions: 1
Number of Holes: 1
Region point: Tensor({0}:-19.242; {1}:10.4768; {2}:5.83092)
Region point: Tensor({0}:-45.3196; {1}:17.2646; {2}:-17.2646)

Tetrahedral meshes can be written out to several formats. Here’s an example to write out in Dolfin XML format: pygamer.writeDolfin('outputTetMesh.xml', tetmesh).

Visualizing the tetrahedralized result

To give a sense of the result, we can cut a slice through the resulting tetrahedral mesh. Here we are also eliminating the internal simplices which would be occluded.

toRemove = []
for vertexID in tetmesh.vertexIDs:
    if vertexID.data()[1] > -1:
        toRemove.append(vertexID)
for v in toRemove:
    tetmesh.removeVertex(v)

# Extract the surface mesh
surfmesh = tetmesh.extractSurface()
# Note that the resulting mesh isn't guaranteed to be pseudo-manifold
# So correctNormals may produce an error to stderr, but this is ok for the purposes of visualization
surfmesh.correctNormals()
v, e, f = surfmesh.to_ndarray()

# Draw the protein and surrounding box together
ctx = threevis.Context(width=640, height=480)

rgb = np.ones((len(f), 3))

for i, face in enumerate(surfmesh.faceIDs):
    if face.data().marker == 23:
        rgb[i,0] = 1
        rgb[i,1] = 0
        rgb[i,2] = 0
    elif face.data().marker == 50:
        rgb[i,0] = 0
        rgb[i,1] = 0
        rgb[i,2] = 1

colors = threevis.FaceAttribute(rgb)

ctx.draw_faces(v, f, colors=colors)
ctx.draw_edges(v, e)
ctx.display()

Since we have designated the protein surface (red) as a hole TetGen has not placed tetrahedral inside of the volume. The gray region shows a slice through the tetrahedralized volume. The surface of the box is shown in blue.