|
|
@@ -0,0 +1,351 @@
|
|
|
+/**
|
|
|
+ * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
|
|
|
+ *
|
|
|
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
|
|
|
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
|
|
|
+ */
|
|
|
+
|
|
|
+import { Structure, StructureElement, StructureProperties } from '../../mol-model/structure';
|
|
|
+import { Task, RuntimeContext } from '../../mol-task';
|
|
|
+import { CentroidHelper } from '../../mol-math/geometry/centroid-helper';
|
|
|
+import { AccessibleSurfaceAreaParams } from '../../mol-model-props/computed/accessible-surface-area';
|
|
|
+import { Vec3 } from '../../mol-math/linear-algebra';
|
|
|
+import { getElementMoleculeType } from '../../mol-model/structure/util';
|
|
|
+import { MoleculeType } from '../../mol-model/structure/model/types';
|
|
|
+import { AccessibleSurfaceArea } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley';
|
|
|
+import { ParamDefinition as PD } from '../../mol-util/param-definition';
|
|
|
+import { MembraneOrientation } from './prop';
|
|
|
+import '../../mol-util/polyfill';
|
|
|
+
|
|
|
+const LARGE_CA_THRESHOLD = 5000;
|
|
|
+export const MEMBRANE_STORAGE_KEY = 'MEMBRANE_STORAGE_KEY';
|
|
|
+
|
|
|
+interface TMDETContext {
|
|
|
+ structure: Structure,
|
|
|
+
|
|
|
+ numberOfSpherePoints: number,
|
|
|
+ stepSize: number,
|
|
|
+ minThickness: number,
|
|
|
+ maxThickness: number,
|
|
|
+ asaCutoff: number,
|
|
|
+ adjust: number,
|
|
|
+
|
|
|
+ offsets: ArrayLike<number>,
|
|
|
+ exposed: ArrayLike<number>,
|
|
|
+ hydrophobic: ArrayLike<boolean>,
|
|
|
+ centroid: Vec3,
|
|
|
+ extent: number,
|
|
|
+ large: boolean
|
|
|
+};
|
|
|
+
|
|
|
+export const TMDETParams = {
|
|
|
+ numberOfSpherePoints: PD.Numeric(140, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }),
|
|
|
+ stepSize: PD.Numeric(1, { min: 0.25, max: 4, step: 0.25 }, { description: 'Thickness of membrane slices that will be tested' }),
|
|
|
+ minThickness: PD.Numeric(20, { min: 10, max: 30, step: 1}, { description: 'Minimum membrane thickness used during refinement' }),
|
|
|
+ maxThickness: PD.Numeric(40, { min: 30, max: 50, step: 1}, { description: 'Maximum membrane thickness used during refinement' }),
|
|
|
+ asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Relative ASA cutoff above which residues will be considered' }),
|
|
|
+ adjust: PD.Numeric(14, { min: 0, max: 30, step: 1 }, { description: 'Minimum length of membrane-spanning regions (original values: 14 for alpha-helices and 5 for beta sheets). Set to 0 to not optimize membrane thickness.' }),
|
|
|
+};
|
|
|
+export type TMDETParams = typeof TMDETParams
|
|
|
+export type TMDETProps = PD.Values<TMDETParams>
|
|
|
+
|
|
|
+/**
|
|
|
+ * Implements:
|
|
|
+ * Membrane positioning for high- and low-resolution protein structures through a binary classification approach
|
|
|
+ * Guillaume Postic, Yassine Ghouzam, Vincent Guiraud, and Jean-Christophe Gelly
|
|
|
+ * Protein Engineering, Design & Selection, 2015, 1–5
|
|
|
+ * doi: 10.1093/protein/gzv063
|
|
|
+ */
|
|
|
+export function computeTMDET(structure: Structure, props: TMDETProps) {
|
|
|
+ return Task.create('Compute Membrane Orientation', async runtime => {
|
|
|
+ return await calculate(runtime, structure, props);
|
|
|
+ });
|
|
|
+}
|
|
|
+
|
|
|
+// avoiding namespace lookup improved performance in Chrome (Aug 2020) <<<<<<<<<<<<<<<<< WTH????
|
|
|
+const v3add = Vec3.add;
|
|
|
+const v3clone = Vec3.clone;
|
|
|
+const v3dot = Vec3.dot;
|
|
|
+const v3normalize = Vec3.normalize;
|
|
|
+const v3scale = Vec3.scale;
|
|
|
+const v3set = Vec3.set;
|
|
|
+const v3squaredDistance = Vec3.squaredDistance;
|
|
|
+const v3sub = Vec3.sub;
|
|
|
+const v3zero = Vec3.zero;
|
|
|
+
|
|
|
+const centroidHelper = new CentroidHelper();
|
|
|
+async function initialize(structure: Structure, props: TMDETProps, accessibleSurfaceArea: AccessibleSurfaceArea): Promise<TMDETContext> {
|
|
|
+ const l = StructureElement.Location.create(structure);
|
|
|
+ const { label_atom_id, label_comp_id, x, y, z } = StructureProperties.atom;
|
|
|
+ const asaCutoff = props.asaCutoff / 100;
|
|
|
+ centroidHelper.reset();
|
|
|
+
|
|
|
+ const offsets = new Array<number>();
|
|
|
+ const exposed = new Array<number>();
|
|
|
+ const hydrophobic = new Array<boolean>();
|
|
|
+
|
|
|
+ const vec = v3zero();
|
|
|
+ for (let i = 0, il = structure.units.length; i < il; ++i) {
|
|
|
+ const unit = structure.units[i];
|
|
|
+ const { elements } = unit;
|
|
|
+ l.unit = unit;
|
|
|
+
|
|
|
+ for (let j = 0, jl = elements.length; j < jl; ++j) {
|
|
|
+ const eI = elements[j];
|
|
|
+ l.element = eI;
|
|
|
+
|
|
|
+ // consider only amino acids
|
|
|
+ if (getElementMoleculeType(unit, eI) !== MoleculeType.Protein) {
|
|
|
+ continue;
|
|
|
+ }
|
|
|
+
|
|
|
+ // only CA is considered for downstream operations
|
|
|
+ if (label_atom_id(l) !== 'CA' && label_atom_id(l) !== 'BB') {
|
|
|
+ continue;
|
|
|
+ }
|
|
|
+
|
|
|
+ // original ANVIL only considers canonical amino acids
|
|
|
+ if (!MaxAsa[label_comp_id(l)]) {
|
|
|
+ continue;
|
|
|
+ }
|
|
|
+
|
|
|
+ // while iterating use first pass to compute centroid
|
|
|
+ v3set(vec, x(l), y(l), z(l));
|
|
|
+ centroidHelper.includeStep(vec);
|
|
|
+
|
|
|
+ // keep track of offsets and exposed state to reuse
|
|
|
+ offsets.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
|
|
|
+ if (AccessibleSurfaceArea.getValue(l, accessibleSurfaceArea) / MaxAsa[label_comp_id(l)] > asaCutoff) {
|
|
|
+ exposed.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
|
|
|
+ hydrophobic.push(isHydrophobic(label_comp_id(l)));
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // calculate centroid and extent
|
|
|
+ centroidHelper.finishedIncludeStep();
|
|
|
+ const centroid = v3clone(centroidHelper.center);
|
|
|
+ for (let k = 0, kl = offsets.length; k < kl; k++) {
|
|
|
+ setLocation(l, structure, offsets[k]);
|
|
|
+ v3set(vec, x(l), y(l), z(l));
|
|
|
+ centroidHelper.radiusStep(vec);
|
|
|
+ }
|
|
|
+ const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq);
|
|
|
+
|
|
|
+ return {
|
|
|
+ ...props,
|
|
|
+ structure,
|
|
|
+
|
|
|
+ offsets,
|
|
|
+ exposed,
|
|
|
+ hydrophobic,
|
|
|
+ centroid,
|
|
|
+ extent,
|
|
|
+ large: offsets.length > LARGE_CA_THRESHOLD
|
|
|
+ };
|
|
|
+}
|
|
|
+
|
|
|
+export async function calculate(runtime: RuntimeContext, structure: Structure, params: TMDETProps): Promise<MembraneOrientation> {
|
|
|
+ // can't get away with the default 92 points here
|
|
|
+ const asaProps = { ...PD.getDefaultValues(AccessibleSurfaceAreaParams), probeSize: 4.0, traceOnly: true, numberOfSpherePoints: 184 };
|
|
|
+ const accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure, asaProps).runInContext(runtime);
|
|
|
+
|
|
|
+ const ctx = await initialize(structure, params, accessibleSurfaceArea);
|
|
|
+
|
|
|
+ const normalVector = v3zero();
|
|
|
+ const center = v3zero();
|
|
|
+
|
|
|
+ // localStorage vs sessionStorage
|
|
|
+ const membrane: MembraneOrientation = JSON.parse(
|
|
|
+ window.localStorage.getItem(MEMBRANE_STORAGE_KEY)!
|
|
|
+ );
|
|
|
+ window.console.debug('membrane object from localStorage:', membrane);
|
|
|
+
|
|
|
+ window.console.debug('normal vector:', membrane.normalVector);
|
|
|
+ window.console.debug('plain point 1:', membrane.planePoint1);
|
|
|
+ window.console.debug('plain point 2:', membrane.planePoint2);
|
|
|
+
|
|
|
+ v3sub(normalVector, membrane.planePoint1, membrane.planePoint2);
|
|
|
+ v3normalize(normalVector, normalVector);
|
|
|
+ v3add(center, membrane.planePoint1, membrane.planePoint2);
|
|
|
+ v3scale(center, center, 0.5);
|
|
|
+ window.console.debug('calculated center:', center);
|
|
|
+
|
|
|
+ const candidate: MembraneCandidate = {
|
|
|
+ normalVector: membrane.normalVector,
|
|
|
+ planePoint1: membrane.planePoint1,
|
|
|
+ planePoint2: membrane.planePoint2,
|
|
|
+ stats: HphobHphil.initial(ctx) //null // TODO: WTH?
|
|
|
+ };
|
|
|
+ const extent = adjustExtent(ctx, candidate, center);
|
|
|
+
|
|
|
+ window.console.debug('result of "tmdet / calculate":', {
|
|
|
+ planePoint1: membrane.planePoint1,
|
|
|
+ planePoint2: membrane.planePoint2,
|
|
|
+ normalVector,
|
|
|
+ centroid: center,
|
|
|
+ radius: extent
|
|
|
+ });
|
|
|
+
|
|
|
+ return {
|
|
|
+ planePoint1: membrane.planePoint1,
|
|
|
+ planePoint2: membrane.planePoint2,
|
|
|
+ normalVector,
|
|
|
+ centroid: center,
|
|
|
+ radius: extent
|
|
|
+ };
|
|
|
+}
|
|
|
+
|
|
|
+interface MembraneCandidate {
|
|
|
+ planePoint1: Vec3,
|
|
|
+ planePoint2: Vec3,
|
|
|
+ stats: HphobHphil,
|
|
|
+ normalVector?: Vec3,
|
|
|
+ spherePoint?: Vec3,
|
|
|
+ qmax?: number
|
|
|
+}
|
|
|
+
|
|
|
+namespace MembraneCandidate {
|
|
|
+ export function initial(c1: Vec3, c2: Vec3, stats: HphobHphil): MembraneCandidate {
|
|
|
+ return {
|
|
|
+ planePoint1: c1,
|
|
|
+ planePoint2: c2,
|
|
|
+ stats
|
|
|
+ };
|
|
|
+ }
|
|
|
+
|
|
|
+ export function scored(spherePoint: Vec3, planePoint1: Vec3, planePoint2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate {
|
|
|
+ const normalVector = v3zero();
|
|
|
+ v3sub(normalVector, centroid, spherePoint);
|
|
|
+ return {
|
|
|
+ planePoint1,
|
|
|
+ planePoint2,
|
|
|
+ stats,
|
|
|
+ normalVector,
|
|
|
+ spherePoint,
|
|
|
+ qmax
|
|
|
+ };
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+/** Filter for membrane residues and calculate the final extent of the membrane layer */
|
|
|
+function adjustExtent(ctx: TMDETContext, membrane: MembraneCandidate, centroid: Vec3): number {
|
|
|
+ const { offsets, structure } = ctx;
|
|
|
+ const { normalVector, planePoint1, planePoint2 } = membrane;
|
|
|
+ const l = StructureElement.Location.create(structure);
|
|
|
+ const testPoint = v3zero();
|
|
|
+ const { x, y, z } = StructureProperties.atom;
|
|
|
+
|
|
|
+ const d1 = -v3dot(normalVector!, planePoint1);
|
|
|
+ const d2 = -v3dot(normalVector!, planePoint2);
|
|
|
+ const dMin = Math.min(d1, d2);
|
|
|
+ const dMax = Math.max(d1, d2);
|
|
|
+ let extent = 0;
|
|
|
+
|
|
|
+ for (let k = 0, kl = offsets.length; k < kl; k++) {
|
|
|
+ setLocation(l, structure, offsets[k]);
|
|
|
+ v3set(testPoint, x(l), y(l), z(l));
|
|
|
+ if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
|
|
|
+ const dsq = v3squaredDistance(testPoint, centroid);
|
|
|
+ if (dsq > extent) extent = dsq;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return Math.sqrt(extent);
|
|
|
+}
|
|
|
+
|
|
|
+export function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean {
|
|
|
+ const d1 = -v3dot(normalVector, planePoint1);
|
|
|
+ const d2 = -v3dot(normalVector, planePoint2);
|
|
|
+ return _isInMembranePlane(testPoint, normalVector, Math.min(d1, d2), Math.max(d1, d2));
|
|
|
+}
|
|
|
+
|
|
|
+function _isInMembranePlane(testPoint: Vec3, normalVector: Vec3, min: number, max: number): boolean {
|
|
|
+ const d = -v3dot(normalVector, testPoint);
|
|
|
+ return d > min && d < max;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+interface HphobHphil {
|
|
|
+ hphob: number,
|
|
|
+ hphil: number
|
|
|
+}
|
|
|
+
|
|
|
+namespace HphobHphil {
|
|
|
+ export function initial(ctx: TMDETContext): HphobHphil {
|
|
|
+ const { exposed, hydrophobic } = ctx;
|
|
|
+ let hphob = 0;
|
|
|
+ let hphil = 0;
|
|
|
+ for (let k = 0, kl = exposed.length; k < kl; k++) {
|
|
|
+ if (hydrophobic[k]) {
|
|
|
+ hphob++;
|
|
|
+ } else {
|
|
|
+ hphil++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return { hphob, hphil };
|
|
|
+ }
|
|
|
+
|
|
|
+ const testPoint = v3zero();
|
|
|
+ export function sliced(ctx: TMDETContext, stepSize: number, spherePoint: Vec3, diam: Vec3, diamNorm: number): HphobHphil[] {
|
|
|
+ const { exposed, hydrophobic, structure } = ctx;
|
|
|
+ const { units, serialMapping } = structure;
|
|
|
+ const { unitIndices, elementIndices } = serialMapping;
|
|
|
+ const sliceStats: HphobHphil[] = [];
|
|
|
+ for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
|
|
|
+ sliceStats[sliceStats.length] = { hphob: 0, hphil: 0 };
|
|
|
+ }
|
|
|
+
|
|
|
+ for (let i = 0, il = exposed.length; i < il; i++) {
|
|
|
+ const unit = units[unitIndices[exposed[i]]];
|
|
|
+ const elementIndex = elementIndices[exposed[i]];
|
|
|
+ v3set(testPoint, unit.conformation.x(elementIndex), unit.conformation.y(elementIndex), unit.conformation.z(elementIndex));
|
|
|
+ v3sub(testPoint, testPoint, spherePoint);
|
|
|
+ if (hydrophobic[i]) {
|
|
|
+ sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphob++;
|
|
|
+ } else {
|
|
|
+ sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphil++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sliceStats;
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+/** ANVIL-specific (not general) definition of membrane-favoring amino acids */
|
|
|
+const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'TRP', 'VAL']);
|
|
|
+/** Returns true if ANVIL considers this as amino acid that favors being embedded in a membrane */
|
|
|
+export function isHydrophobic(label_comp_id: string): boolean {
|
|
|
+ return HYDROPHOBIC_AMINO_ACIDS.has(label_comp_id);
|
|
|
+}
|
|
|
+
|
|
|
+/** Accessible surface area used for normalization. ANVIL uses 'Total-Side REL' values from NACCESS, from: Hubbard, S. J., & Thornton, J. M. (1993). naccess. Computer Program, Department of Biochemistry and Molecular Biology, University College London, 2(1). */
|
|
|
+export const MaxAsa: { [k: string]: number } = {
|
|
|
+ 'ALA': 69.41,
|
|
|
+ 'ARG': 201.25,
|
|
|
+ 'ASN': 106.24,
|
|
|
+ 'ASP': 102.69,
|
|
|
+ 'CYS': 96.75,
|
|
|
+ 'GLU': 134.74,
|
|
|
+ 'GLN': 140.99,
|
|
|
+ 'GLY': 32.33,
|
|
|
+ 'HIS': 147.08,
|
|
|
+ 'ILE': 137.96,
|
|
|
+ 'LEU': 141.12,
|
|
|
+ 'LYS': 163.30,
|
|
|
+ 'MET': 156.64,
|
|
|
+ 'PHE': 164.11,
|
|
|
+ 'PRO': 119.90,
|
|
|
+ 'SER': 78.11,
|
|
|
+ 'THR': 101.70,
|
|
|
+ 'TRP': 211.26,
|
|
|
+ 'TYR': 177.38,
|
|
|
+ 'VAL': 114.28
|
|
|
+};
|
|
|
+
|
|
|
+function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
|
|
|
+ l.structure = structure;
|
|
|
+ l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];
|
|
|
+ l.element = structure.serialMapping.elementIndices[serialIndex];
|
|
|
+ return l;
|
|
|
+}
|
