Issue #2: transforming membrane instead of structure

src/apps/viewer/index.ts

@@ -22,7 +22,7 @@ require('mol-plugin-ui/skin/light.scss');
 
 export { PLUGIN_VERSION as version } from '../../mol-plugin/version';
 export { setDebugMode, setProductionMode } from '../../mol-util/debug';
-export { loadWithUNITMPMembraneRepresentation } from '../../extensions/tmdet/behavior';
+export { loadWithUNITMPMembraneRepresentation, pdbtmDescriptor } from '../../extensions/tmdet/behavior';
 
 const Extensions = {
     'tmdet-membrane-orientation': PluginSpec.Behavior(TMDETMembraneOrientation)

src/extensions/tmdet/algorithm.ts

@@ -5,39 +5,13 @@
  * @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 { Structure, StructureElement } from '../../mol-model/structure';
 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' }),
@@ -49,194 +23,7 @@ export const TMDETParams = {
 export type TMDETParams = typeof TMDETParams
 export type TMDETProps = PD.Values<TMDETParams>
 
-export function computeTMDET(structure: Structure, props: TMDETProps) {
-    return Task.create('Compute Membrane Orientation', async runtime => {
-        return 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)!
-    );
-    // TODO:
-    window.console.debug('membrane object from localStorage:', membrane);
-
-    v3sub(normalVector, membrane.planePoint1, membrane.planePoint2);
-    v3normalize(normalVector, normalVector);
-    v3add(center, membrane.planePoint1, membrane.planePoint2);
-    v3scale(center, center, 0.5);
-
-    const candidate: MembraneCandidate = {
-        normalVector: membrane.normalVector,
-        planePoint1: membrane.planePoint1,
-        planePoint2: membrane.planePoint2,
-        stats: HphobHphil.initial(ctx)
-    };
-    window.console.debug('candidate:', candidate);
-
-    const extent = adjustExtent(ctx, candidate, center);
-
-    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);
@@ -248,87 +35,3 @@ function _isInMembranePlane(testPoint: Vec3, normalVector: Vec3, min: number, ma
     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;
-}

src/extensions/tmdet/behavior.ts

@@ -25,10 +25,9 @@ import { Color } from '../../mol-util/color';
 import { createStructureRepresentationParams } from '../../mol-plugin-state/helpers/structure-representation-params';
 import { PluginUIContext } from '../../mol-plugin-ui/context';
 import { StateObjectSelector } from "../../mol-state/object";
-import { MEMBRANE_STORAGE_KEY } from '../../extensions/tmdet/algorithm';
 import { Vec3 } from '../../mol-math/linear-algebra';
 import { applyTransformations, getAtomGroupExpression } from './transformation';
-import { PDBTMDescriptor, PDBTMRegion, PMS } from './types';
+import { PDBTMDescriptor, PDBTMRegion, PDBTMTransformationMatrix, PMS } from './types';
 
 type StructureComponentType = StateObjectSelector<
         PMS,
@@ -112,15 +111,14 @@ export const isTransmembrane = StructureSelectionQuery('Residues Embedded in Mem
 
 // //////////////////////////// TMDET VIEWER FUNCTIONS
 
-let membrane: MembraneOrientation;
-
+export let membraneOrientation: MembraneOrientation;
+export let pdbtmDescriptor: PDBTMDescriptor;
 
 export async function loadWithUNITMPMembraneRepresentation(plugin: PluginUIContext, params: any) {
-    const pdbtmDescriptor: any = await downloadRegionDescriptor(plugin, params);
-
-    membrane = createMembraneOrientation(pdbtmDescriptor);
+    // const pdbtmDescriptor: any = await downloadRegionDescriptor(plugin, params);
+    pdbtmDescriptor = await downloadRegionDescriptor(plugin, params);
 
-    localStorage.setItem(MEMBRANE_STORAGE_KEY, JSON.stringify(membrane));
+    membraneOrientation = createMembraneOrientation(pdbtmDescriptor);
 
     // load structure
     await loadStructure(plugin, params, pdbtmDescriptor);
@@ -133,24 +131,35 @@ export async function loadWithUNITMPMembraneRepresentation(plugin: PluginUIConte
     requestAnimationFrame(() => plugin.canvas3d?.requestCameraReset());
 }
 
+function vadd(a: Vec3, b: Vec3): Vec3 {
+    return Vec3.add(Vec3.zero(), a, b);
+}
+
+function vneg(u: Vec3): Vec3 {
+    return Vec3.negate(Vec3.zero(), u);
+}
+
 function createMembraneOrientation(pdbtmDescriptor: PDBTMDescriptor): MembraneOrientation {
-    const normal = pdbtmDescriptor.additional_entry_annotations.membrane.normal;
-    const membraneNormal: Vec3 = Vec3.fromObj(normal);
-
-    const membraneOrientation: MembraneOrientation = {
-        planePoint1: Vec3.fromArray(Vec3.zero(), membraneNormal, 0),
-        planePoint2: Vec3.fromArray(Vec3.zero(), membraneNormal, 0),
-        centroid: Vec3.fromArray(
-            Vec3.zero(), [0, 0, 0], 0
-        ),
+    const membrane = pdbtmDescriptor.additional_entry_annotations.membrane;
+    const mx: PDBTMTransformationMatrix = membrane.transformation_matrix;
+    const translation: Vec3 = Vec3.create(mx.rowx.t, mx.rowy.t, mx.rowz.t);
+    Vec3.negate(translation, translation);
+
+    const membraneNormal: Vec3 = Vec3.fromObj(membrane.normal);
+
+    const result: MembraneOrientation = {
+        planePoint1: vadd(membraneNormal, translation),
+        planePoint2: vadd(vneg(membraneNormal), translation),
+        centroid: translation,
         normalVector: membraneNormal,
 
         // (NOTE: the TMDET extension calculates and sets it during applying preset)
-        radius: pdbtmDescriptor.additional_entry_annotations.membrane.radius
+        radius: membrane.radius
     };
-    membraneOrientation.planePoint2[2] *= -1;
-
-    return membraneOrientation;
+    console.log("createMemOri\n", result);
+    Vec3.add(result.normalVector, result.normalVector, translation);
+    console.log("createMemOri2\n", result);
+    return result;
 }
 
 async function createStructureRepresentation(plugin: PluginUIContext, pdbtmDescriptor: any) {

src/extensions/tmdet/prop.ts

@@ -8,13 +8,14 @@
 import { ParamDefinition as PD } from '../../mol-util/param-definition';
 import { Structure, StructureProperties, Unit } from '../../mol-model/structure';
 import { CustomPropertyDescriptor } from '../../mol-model/custom-property';
-import { TMDETParams, TMDETProps, computeTMDET, isInMembranePlane } from './algorithm';
+import { isInMembranePlane, TMDETParams } from './algorithm';
 import { CustomStructureProperty } from '../../mol-model-props/common/custom-structure-property';
 import { CustomProperty } from '../../mol-model-props/common/custom-property';
 import { Vec3 } from '../../mol-math/linear-algebra';
 import { QuerySymbolRuntime } from '../../mol-script/runtime/query/base';
 import { CustomPropSymbol } from '../../mol-script/language/symbol';
 import { Type } from '../../mol-script/language/type';
+import { membraneOrientation } from './behavior';
 
 export const MembraneOrientationParams = {
     ...TMDETParams
@@ -75,11 +76,8 @@ export const MembraneOrientationProvider: CustomStructureProperty.Provider<Membr
     obtain: async (ctx: CustomProperty.Context, data: Structure, props: Partial<MembraneOrientationProps>) => {
         const p = { ...PD.getDefaultValues(MembraneOrientationParams), ...props };
         window.console.log('obtain:: DEBUG', data.customPropertyDescriptors);
-        return { value: await compute(ctx, data, p) };
+        let result = membraneOrientation;
+        console.log("RESULT of 'obtain:'", result);
+        return { value: result };
     }
 });
-
-async function compute(ctx: CustomProperty.Context, data: Structure, props: Partial<TMDETProps>): Promise<MembraneOrientation> {
-    const p = { ...PD.getDefaultValues(TMDETParams), ...props };
-    return computeTMDET(data, p).runInContext(ctx.runtime);
-}

src/extensions/tmdet/transformation.ts

@@ -9,9 +9,23 @@ import { Color } from '../../mol-util/color';
 import { Expression } from '../../mol-script/language/expression';
 
 export async function applyTransformations(ctx: PluginUIContext, pdbtmDescriptor: PDBTMDescriptor) {
+    const annotations = pdbtmDescriptor.additional_entry_annotations;
+    const transformation = transformationForStateTransform(annotations.membrane.transformation_matrix);
 
-    console.log('BIOMX', pdbtmDescriptor.additional_entry_annotations.biomatrix);
-    pdbtmDescriptor.additional_entry_annotations.biomatrix.matrix_list.forEach(function(mx) {
+    console.log("MEMBRANE TRNASFORMATION:\n\n", Mat4.makeTable(transformation));
+    const structure: StateObjectRef<PMS> = ctx.managers.structure.hierarchy.current.models[0].structures[0].cell;
+    const update: StateBuilder.To<any, any> = ctx.build().to(structure);
+    update.apply(StateTransforms.Model.TransformStructureConformation, {
+        "transform": { name: "matrix", params: { data: transformation, transpose: false } }
+    });
+    update.commit();
+
+    console.log('BIOMX', annotations.biomatrix);
+    if (annotations.biomatrix == undefined) {
+        return;
+    }
+
+    annotations.biomatrix.matrix_list.forEach(function(mx) {
         mx.apply_to_chain_list.forEach(function(chainPair) {
             let id  = chainPair.chain_id;
             let newId = chainPair.new_chain_id;