"""
Module for displaying and manipulating a membrane.
Used by Prime panel, System Builder, and thermal_mmgbsa.py
Copyright Schrodinger, LLC. All rights reserved.
"""
import copy
import math
import sys
import numpy
from sklearn.cluster import SpectralClustering
from schrodinger import structure
from schrodinger.structutils import analyze
from schrodinger.structutils import measure
from schrodinger.structutils import transform
array = numpy.array # Make local
try:
from schrodinger.graphics3d import box # for drawing membrane
except ImportError:
box = None
# Constants used to display the membrane as planes:
MEMBRANE_SIZE = 30
MEMBRANE_COLOR = "red"
MEMBRANE_OPACITY = 0.40
vector_magnitude = transform.get_vector_magnitude
get_normalized = transform.get_normalized_vector
[docs]class Membrane_Model:
[docs] def __init__(self, ct=None):
self.ct = ct
# center of the membrane:
self.center = array([0.0, 0.0, 0.0])
# Normalized vector from one membrane "plane" to the other:
self.orientation = array([0.0, 0.0, 1.0])
self.thickness = 44.5
self.membrane_group = box.Group()
return
[docs] def get_vector_atoms_from_internal_coords(self):
# Membrane axis.
coords1 = self.center - 0.5 * self.thickness * self.orientation
coords2 = self.center + 0.5 * self.thickness * self.orientation
return (coords1, coords2)
[docs] def update_internal_coords_to_vector_atoms(self, coords1, coords2):
"""
Update internal coordinates from 2 numpy arrays containing the x,y,z
coordinates of two atoms defining the membrane.
"""
# In case inputs are python lists:
coords1 = array(coords1)
coords2 = array(coords2)
self.center = (coords1 + coords2) / 2.0
diff = coords2 - coords1
self.thickness = math.sqrt(diff[0]**2 + diff[1]**2 + diff[2]**2)
self.orientation = diff / self.thickness
return
[docs] def getCenterOrientationOfAtoms(self, atom_list):
# Determine geometric center of the protein (workspace):
atom_coordinates = []
for ai in atom_list:
atom = self.ct.atom[ai]
atom_coordinates.append(array(atom.xyz))
if len(atom_coordinates) < 2:
raise RuntimeError("Must have at least 2 atoms to add a membrane.")
# If only one atom is present, vector to center is 0,0,0
# which causes a division error and shows no direction anyway.
center = numpy.average(array(atom_coordinates), 0)
# A simple shape-based approach (assumes the membrane normal is the longest).
# FIXME Please document better why distance cutoffs are needed:
min_atom_distance = 15.0
max_atom_distance = 1000.0
# Determine orientation (normalized vector from one dummy atom to another)
# Algorithm: For all atoms that are within distance threshold from center:
# 1. calculate vector from center to that atom
# 2. Add all such vectors (except those that don't contribute)
# 3. Resulting vector is direction for the membrane
orientation_vectors = []
for coor in atom_coordinates:
vector_from_center = coor - center
# distance between the atom and the center:
distance = vector_magnitude(vector_from_center)
if distance > min_atom_distance and distance < max_atom_distance:
orientation_vectors.append(vector_from_center)
if len(orientation_vectors) == 0:
# Ignore the distance cutoffs Ev:70470
for coor in atom_coordinates:
vector_from_center = coor - center
orientation_vectors.append(vector_from_center)
sum_of_vectors = orientation_vectors[0]
for vector in orientation_vectors[1:]:
# Subract if the vector is in opposite direction:
if numpy.dot(sum_of_vectors, vector) > 0.0:
sum_of_vectors += vector
else:
sum_of_vectors -= vector
# Normalize the orientation vector (use normalize?):
orientation = get_normalized(sum_of_vectors)
return center, orientation
[docs] def autoPlaceByMolecule(self, mol_atom_lists):
"""
Auto places the membrane based on the average vector between all
specified molecules (list of atom iterators)
"""
all_atoms = []
orientation_vectors = []
for molatoms in mol_atom_lists:
try:
center, orientation = self.getCenterOrientationOfAtoms(molatoms)
except RuntimeError:
# This molecule has no atoms in proximity to its center (skip)
continue
orientation_vectors.append(orientation)
all_atoms.extend(molatoms)
if len(orientation_vectors) == 0:
msg = "None of the fragments in transmembrane region (ASL) had"
msg += "\nsignificant direction to determine membrane position."
raise RuntimeError(msg)
sum_of_vectors = orientation_vectors[0]
for vector in orientation_vectors[1:]:
# Subract if the vector is in opposite direction:
if numpy.dot(sum_of_vectors, vector) > 0.0:
sum_of_vectors += vector
else:
sum_of_vectors -= vector
# Normalize the orientation vector (use normalize?):
self.orientation = get_normalized(sum_of_vectors)
# Determine geometric center of the protein (workspace):
coordinates = []
for ai in all_atoms:
atom = self.ct.atom[ai]
coordinates.append(array(atom.xyz))
self.center = numpy.average(array(coordinates), 0)
return
[docs] def placeFromExplicitMembrane(self):
"""
Attempt to detect explicit membrane atoms in the structure, and set
self.center and self.orientation based on its orientation.
"""
MEMBRANE_HEAD_ASL = '(res. DPPC,DMPC,POPE,POPC) AND (atom.ele N,O,P)'
head_atoms = analyze.evaluate_asl(self.ct, MEMBRANE_HEAD_ASL)
if len(head_atoms) < 20:
raise RuntimeError('No explicit membrane present')
# Group membrane head atoms by their coordinates, into 2 clusters,
# one for each membrane "side":
head_coordinates = numpy.asarray(
[self.ct.atom[a].xyz for a in head_atoms])
clust = SpectralClustering(n_clusters=2, affinity='nearest_neighbors')
clust_memberships = clust.fit_predict(head_coordinates)
# Split the membrane head atoms into 2 lists, one for each side:
coords1 = []
coords2 = []
for i, membership in enumerate(clust_memberships):
xyz = head_coordinates[i]
if membership == 0:
coords1.append(xyz)
elif membership == 1:
coords2.append(xyz)
else:
raise RuntimeError('Should never happen')
# Orientation of the membrane is caluclated by averaging the vectors that
# are normal to the 2 planes:
orientation1 = measure.fit_plane_to_points(numpy.asarray(coords1))
orientation2 = measure.fit_plane_to_points(numpy.asarray(coords2))
if numpy.dot(orientation1, orientation2) < 0:
# If plane vectors are in opposite directions, reverse one
orientation2 = -orientation2
self.orientation = numpy.mean([orientation1, orientation2], axis=0)
# This is same as averaging the 2 centroids, of 2 planes:
self.center = numpy.mean(head_coordinates, axis=0)
[docs] def writePrimeProperties(self):
"""
Add implicit membrane properties to the structure, based on the
current membrane orientation.
"""
(coor1, coor2) = self.get_vector_atoms_from_internal_coords()
self.ct.property['r_psp_Memb1_x'] = coor1[0]
self.ct.property['r_psp_Memb1_y'] = coor1[1]
self.ct.property['r_psp_Memb1_z'] = coor1[2]
self.ct.property['r_psp_Memb2_x'] = coor2[0]
self.ct.property['r_psp_Memb2_y'] = coor2[1]
self.ct.property['r_psp_Memb2_z'] = coor2[2]
[docs] def findHydrophobicCenter(self):
"""
Returns coordinates of center of mass of all hydrophobic residues
"""
asl = '( res. PHE,LEU,ILE,TYR,TRP,VAL,MET,PRO,CYS,ALA )'
atom_list = analyze.evaluate_asl(self.ct, asl)
coordinates = []
for ai in atom_list:
atom = self.ct.atom[ai]
coordinates.append(array(atom.xyz))
center = numpy.average(array(coordinates), 0)
return center
[docs] def autoPlace(self):
"""
Automatically orient the membrane according to the protein in self.ct
"""
atom_list = list(range(1, self.ct.atom_total + 1))
self.center, self.orientation = self.getCenterOrientationOfAtoms(
atom_list)
[docs] def rotateProteinToMembrane(self):
"""
Translate the protein so that the membrane will be located at the
origin and rotate the protein so that membrane is along the z-axis.
Assumes that center/orientation/thickness are set.
At the end the protein will not have vector atoms.
"""
# Make the membrane center the new origin:
new_origin = transform.get_coords_array_from_list(self.center)
matrix = transform.get_translation_matrix(-1 * new_origin)
transform.transform_structure(self.ct, matrix)
# We need to rotate the workspace so that the membrane orientation
# changes to going up-down:
z_axis_vector = numpy.array(transform.Z_AXIS)
matrix = transform.get_alignment_matrix(self.orientation, z_axis_vector)
transform.transform_structure(self.ct, matrix)
self.center = array([0.0, 0.0, 0.0])
self.orientation = numpy.array([0.0, 0.0, 1.0])
# PANEL-18318: update 3D representation of the membrane box since
# its center and orientation was changed.
self.calculateMembraneBox()
[docs] def calculateMembraneBox(self):
"""
Creates 3D objects for the representation of the membrane box (2 red
squares) in this instance.
The membrane info is taken from center/orientation/thickness
"""
(coor1, coor2) = self.get_vector_atoms_from_internal_coords()
# Perform the vector operation
# We have coordinates for two positions (coor1, coor2) which
# define the membrane planes. We need to figure out a vector
# perpendiclar to the vector between those two points and a
# series of points in the planes in order to draw a polygon
# Assign the size of the planes
size = MEMBRANE_SIZE
# AB is the vector between the two points
AB = coor2 - coor1
# M is the midpoint between the two points
M = coor1 + 0.5 * AB
# N is an arbitrary point we will project onto M to get
# a perpendicular vector:
N = copy.copy(M)
# Setting N to (10,10,10) ensures that the membrane cubes are drawn
# along the X/Y axis IF the AB vector is along the Z axis.
# NOTE that this will not work in rare case where AB vector is
# in this same direction.
N[0] += 10.0
N[1] += 10.0
N[2] += 10.0
MN = N - M
dot1 = numpy.dot(AB, AB) # magnitude of AB squared
dot2 = numpy.dot(AB, MN)
K = dot2 / dot1
# O is the projection of N onto AB. The vector O-N perpendicular
# to A-B
O = M + numpy.dot(K, AB) # noqa: E741
ON = get_normalized(N - O)
# P is perpendicular to both AB and ON
P = get_normalized(numpy.cross(AB, ON))
ABnorm = get_normalized(AB)
head_thickness = (self.thickness - 30.0) / 2.0
coor1offset = coor1 + head_thickness * ABnorm
coor2offset = coor2 - head_thickness * ABnorm
def create_membrane_box(coor, offset):
points = [
coor - size * P, # "origin" point on main plane
coor - size * ON, # 2nd point on main plane
coor + size * ON, # 3rd point on main plane
offset - size * P, # point from the offset plane adj to origin
]
return box.MaestroParallelepiped(points=points,
color=MEMBRANE_COLOR,
opacity=MEMBRANE_OPACITY,
style=box.FILL)
# This ordering is needed because order of points is relevant
plane1_box = create_membrane_box(coor1offset, coor1)
plane2_box = create_membrane_box(coor2, coor2offset)
self.membrane_group.clear()
self.membrane_group.add(plane1_box)
self.membrane_group.add(plane2_box)
return
[docs] def show(self):
"""
Show the membrane group
"""
self.membrane_group.show()
[docs] def hide(self):
"""
Hide the membrane.
"""
self.membrane_group.hide()
[docs] def clear(self):
"""
Remove the 3D objects from the group, which removes them from Maestro's
fit bounds.
"""
# Ev:82397
self.membrane_group.clear()
[docs] def isDefined(self):
"""
Return True if the membrane dimensions are defined.
"""
return bool(self.membrane_group)
[docs] def write_structure(self, filename):
self.ct.write(filename)
return
if __name__ == '__main__':
mode = sys.argv[1]
in_file = sys.argv[2]
out_file = sys.argv[3]
ct = structure.StructureReader.read(in_file)
membrane = Membrane_Model(ct)
if mode == "-add":
print('\nConstructing membrane...\n')
sys.stdout.flush()
parameter_file = sys.argv[4]
try:
coords1 = array([ct.property['r_psp_Memb1_x'], \
ct.property['r_psp_Memb1_y'], \
ct.property['r_psp_Memb1_z']])
coords2 = array([ct.property['r_psp_Memb2_x'], \
ct.property['r_psp_Memb2_y'], \
ct.property['r_psp_Memb2_z']])
except:
print('ERROR: membrane properties not defined in input file ' +
in_file)
print(
'Please use the Setup Membrane panel from within Prime Refinement to'
)
print('define the membrane coordinates for the selected entries.')
sys.exit()
membrane.update_internal_coords_to_vector_atoms(coords1, coords2)
membrane.write_parameter_file(parameter_file)
elif mode == "-remove":
print('\nRemoving membrane...\n')
sys.stdout.flush()
membrane.write_structure(out_file)
#EOF