"""
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 schrodinger import structure
from schrodinger.structutils import analyze, measure, transform
from sklearn.cluster import SpectralClustering
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
l = 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)
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 - l * P, # "origin" point on main plane
coor - l * ON, # 2nd point on main plane
coor + l * ON, # 3rd point on main plane
offset - l * 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