accessible surface area calculation improvements

src/mol-model-props/computed/accessible-surface-area.ts

@@ -0,0 +1,64 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { ParamDefinition as PD } from '../../mol-util/param-definition'
+import { ShrakeRupleyComputationParams, AccessibleSurfaceArea } from './accessible-surface-area/shrake-rupley';
+import { Structure, CustomPropertyDescriptor } from '../../mol-model/structure';
+import { Task } from '../../mol-task';
+import { idFactory } from '../../mol-util/id-factory';
+
+
+const nextAccessibleSurfaceAreaId = idFactory()
+
+export namespace ComputedAccessibleSurfaceArea {
+    export type Property = {
+        id: number
+        asa: AccessibleSurfaceArea
+    }
+
+    export function get(structure: Structure): Property | undefined {
+        return structure.inheritedPropertyData.__ComputedAccessibleSurfaceArea__;
+    }
+    function set(structure: Structure, prop: Property) {
+        (structure.inheritedPropertyData.__ComputedAccessibleSurfaceArea__ as Property) = prop;
+    }
+
+    export function createAttachTask(params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
+        return (structure: Structure) => Task.create('Compute Accessible Surface Area', async ctx => {
+            if (get(structure)) return true;
+            return await attachFromCifOrCompute(structure, params)
+        });
+    }
+
+    export const Descriptor = CustomPropertyDescriptor({
+        isStatic: true,
+        name: 'molstar_accessible_surface_area',
+        // TODO `cifExport` and `symbol`
+    });
+
+    export async function attachFromCifOrCompute(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) {
+        if (structure.customPropertyDescriptors.has(Descriptor)) return true;
+
+        const compAccessibleSurfaceArea = await computeAccessibleSurfaceArea(structure, params)
+
+        structure.customPropertyDescriptors.add(Descriptor);
+        set(structure, compAccessibleSurfaceArea);
+        return true;
+    }
+}
+
+export const AccessibleSurfaceAreaComputationParams = {
+    ...ShrakeRupleyComputationParams
+}
+export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams
+export type AccessibleSurfaceAreaComputationProps = PD.Values<AccessibleSurfaceAreaComputationParams>
+
+async function computeAccessibleSurfaceArea(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps>): Promise<ComputedAccessibleSurfaceArea.Property> {
+    const p = { ...PD.getDefaultValues(AccessibleSurfaceAreaComputationParams), params };
+
+    const accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure, p);
+    return { id: nextAccessibleSurfaceAreaId(), asa: accessibleSurfaceArea }
+}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts

@@ -0,0 +1,101 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { Task, RuntimeContext } from '../../../mol-task';
+import { BitFlags } from '../../../mol-util';
+import { ParamDefinition as PD } from '../../../mol-util/param-definition'
+import { Vec3 } from '../../../mol-math/linear-algebra';
+import { Structure } from '../../../mol-model/structure';
+import { assignRadiusForHeavyAtoms } from './shrake-rupley/radii';
+import { ShrakeRupleyContext, VdWLookup } from './shrake-rupley/common';
+import { computePerResidue } from './shrake-rupley/per-residue';
+import { normalizeAccessibleSurfaceArea } from './shrake-rupley/normalize';
+
+export const ShrakeRupleyComputationParams = {
+    numberOfSpherePoints: PD.Numeric(92, {}, { description: 'number of sphere points to sample per atom: 92 (original paper), 960 (BioJava), 3000 (EPPIC) - see Shrake A, Rupley JA: Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 1973.' }),
+    probeSize: PD.Numeric(1.4, {}, { description: 'corresponds to the size of a water molecule: 1.4 (original paper), 1.5 (occassionally used)' }),
+    buriedRasaThreshold: PD.Numeric(0.16, { min: 0.0, max: 1.0 }, { description: 'below this cutoff of relative accessible surface area a residue will be considered buried - see: Rost B, Sander C: Conservation and prediction of solvent accessibility in protein families. Proteins 1994.' }),
+    nonPolymer: PD.Boolean(false, { description: 'Include non-polymer atoms in computation.' })
+}
+export type ShrakeRupleyComputationParams = typeof ShrakeRupleyComputationParams
+
+namespace AccessibleSurfaceArea {
+    /**
+     * Adapts the BioJava implementation by Jose Duarte. That implementation is based on the publication by Shrake, A., and
+     * J. A. Rupley. "Environment and Exposure to Solvent of Protein Atoms. Lysozyme and Insulin." JMB (1973).
+     */
+    export function compute(structure: Structure,
+        params: Partial<PD.Values<ShrakeRupleyComputationParams>> = {}) {
+        params = { ...PD.getDefaultValues(ShrakeRupleyComputationParams), ...params };
+        return Task.create('Compute Accessible Surface Area', async rtctx => {
+            return await _compute(rtctx, structure, params);
+        }).run();
+    }
+
+    async function _compute(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>> = {}): Promise<AccessibleSurfaceArea> {
+        const ctx = initialize(rtctx, structure, params);
+
+        assignRadiusForHeavyAtoms(ctx);
+        computePerResidue(ctx);
+        normalizeAccessibleSurfaceArea(ctx);
+
+        return {
+            accessibleSurfaceArea: ctx.accessibleSurfaceArea,
+            relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea,
+            buried: (index: number) => ctx.relativeAccessibleSurfaceArea[index] < 0.16
+        };
+    }
+
+    function initialize(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>>): ShrakeRupleyContext {
+        const { elementCount, polymerResidueCount } = structure;
+
+        return {
+            rtctx: rtctx,
+            structure: structure,
+            probeSize: params.probeSize!,
+            nonPolymer: params.nonPolymer!,
+            spherePoints: generateSpherePoints(params.numberOfSpherePoints!),
+            cons: 4.0 * Math.PI / params.numberOfSpherePoints!,
+            maxLookupRadius: 2 * params.probeSize! + 2 * VdWLookup[2], // 2x probe size + 2x largest VdW
+            atomRadius: new Int8Array(elementCount),
+            accessibleSurfaceArea: new Float32Array(polymerResidueCount),
+            relativeAccessibleSurfaceArea: new Float32Array(polymerResidueCount),
+            updateChunk: 25000
+        }
+    }
+
+    /** Creates a collection of points on a sphere by the Golden Section Spiral algorithm. */
+    function generateSpherePoints(numberOfSpherePoints: number): Vec3[] {
+        const points: Vec3[] = [];
+        const inc = Math.PI * (3.0 - Math.sqrt(5.0));
+        const offset = 2.0 / numberOfSpherePoints;
+        for (let k = 0; k < numberOfSpherePoints; ++k) {
+            const y = k * offset - 1.0 + (offset / 2.0);
+            const r = Math.sqrt(1.0 - y * y);
+            const phi = k * inc;
+            points[points.length] = [Math.cos(phi) * r, y, Math.sin(phi) * r] as Vec3;
+        }
+        return points;
+    }
+
+    export namespace SolventAccessibility {
+        export const is: (t: number, f: Flag) => boolean = BitFlags.has
+        export const create: (f: Flag) => number = BitFlags.create
+        export const enum Flag {
+            _ = 0x0,
+            BURIED = 0x1,
+            ACCESSIBLE = 0x2
+        }
+    }
+}
+
+interface AccessibleSurfaceArea {
+    readonly accessibleSurfaceArea: ArrayLike<number>
+    readonly relativeAccessibleSurfaceArea: ArrayLike<number>
+    buried(index: number): boolean
+}
+
+export { AccessibleSurfaceArea }

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts

@@ -0,0 +1,63 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { RuntimeContext } from '../../../../mol-task';
+import { Structure } from '../../../../mol-model/structure';
+import { Vec3 } from '../../../../mol-math/linear-algebra';
+
+export interface ShrakeRupleyContext {
+    rtctx: RuntimeContext,
+    structure: Structure,
+    spherePoints: Vec3[],
+    probeSize: number,
+    nonPolymer: boolean,
+    cons: number,
+    maxLookupRadius: number,
+    atomRadius: Int8Array,
+    accessibleSurfaceArea: Float32Array,
+    relativeAccessibleSurfaceArea: Float32Array,
+    updateChunk: number
+}
+
+// Chothia's amino acid and nucleotide atom vdw radii
+export const VdWLookup = [
+    -1.0, // 0: missing
+    1.76, // 1: trigonal C
+    1.87, // 2: tetrahedral C
+    1.65, // 3: trigonal N
+    1.50, // 4: tetrahedral N
+    1.40, // 5: O
+    1.85, // 6: S
+    1.80, // 7: C (nucleic)
+    1.60, // 8: N (nucleic)
+    1.40  // 9: P (nucleic)
+] // can still be appended on-the-fly for rare elements like selenium
+
+
+/** Maximum accessible surface area observed for amino acids. Taken from: http://dx.doi.org/10.1371/journal.pone.0080635 */
+export const MaxAsa = {
+    'ALA': 121.0,
+    'ARG': 265.0,
+    'ASN': 187.0,
+    'ASP': 187.0,
+    'CYS': 148.0,
+    'GLU': 214.0,
+    'GLN': 214.0,
+    'GLY': 97.0,
+    'HIS': 216.0,
+    'ILE': 195.0,
+    'LEU': 191.0,
+    'LYS': 230.0,
+    'MET': 203.0,
+    'PHE': 228.0,
+    'PRO': 154.0,
+    'SER': 143.0,
+    'THR': 163.0,
+    'TRP': 264.0,
+    'TYR': 255.0,
+    'VAL': 165.0
+}
+export const DefaultMaxAsa = 121.0

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/normalize.ts

@@ -0,0 +1,32 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { ShrakeRupleyContext } from './common';
+import { isPolymer, MaxAsa, DefaultMaxAsa } from '../../../../mol-model/structure/model/types';
+
+export async function normalizeAccessibleSurfaceArea(ctx: ShrakeRupleyContext) {
+    const updateChunk = ctx.updateChunk / 10;
+    const residueCount = ctx.relativeAccessibleSurfaceArea.length;
+    for (let i = 0; i < residueCount; i += updateChunk) {
+        computeRange(ctx, i, Math.min(i + updateChunk, residueCount));
+    }
+}
+
+function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
+    const { accessibleSurfaceArea, relativeAccessibleSurfaceArea, structure } = ctx;
+    const { residues, derived } = structure.model.atomicHierarchy;
+
+    for (let i = begin; i < end; ++i) {
+        // skip entities not part of a polymer chain
+        if (!ctx.nonPolymer) {
+            if (!isPolymer(derived.residue.moleculeType[i])) continue;
+        }
+
+        const maxAsa = (MaxAsa as any)[residues.label_comp_id.value(i)];
+        const rasa = accessibleSurfaceArea[i] / (maxAsa === undefined ? DefaultMaxAsa : maxAsa);
+        relativeAccessibleSurfaceArea[i] = rasa;
+    }
+}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/per-residue.ts

@@ -0,0 +1,86 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { ShrakeRupleyContext, VdWLookup } from './common';
+import { Vec3 } from '../../../../mol-math/linear-algebra';
+
+export async function computePerResidue(ctx: ShrakeRupleyContext) {
+    const { rtctx, updateChunk, atomRadius } = ctx;
+    for (let i = 0; i < atomRadius.length; i += updateChunk) {
+        computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length));
+
+        if (rtctx.shouldUpdate) {
+            rtctx.update({ message: 'Computing per residue surface accessibility...', current: i, max: atomRadius.length });
+        }
+    }
+    // async yields 4x speed-up
+}
+
+const aPos = Vec3.zero();
+const bPos = Vec3.zero();
+let testPoint = Vec3.zero();
+function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
+    const { structure, atomRadius, accessibleSurfaceArea, spherePoints, cons, maxLookupRadius, probeSize } = ctx;
+        const { model } = structure;
+        const { x, y, z } = model.atomicConformation;
+        const { residueAtomSegments } = model.atomicHierarchy;
+        const { lookup3d } = structure;
+
+    const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
+
+    // console.log(`computing ASA chunk ${ begin } to ${ end }`);
+    for (let aI = begin; aI < end; aI++) {
+        const radius1 = VdWLookup[atomRadius[aI]];
+        if (radius1 === VdWLookup[0]) continue;
+
+        // pre-filter by lookup3d
+        // lookup provides >10x speed-up compared to naive evaluation
+        const { count, units, indices, squaredDistances } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
+        // we could keep track of all found neighbors of each atom, however slow and memory intensive
+
+        // collect neighbors for each atom
+        const cutoff1 = probeSize + probeSize + radius1;
+        const neighbors = [];
+        for (let iI = 0; iI < count; ++iI) {
+            const bI = units[iI].elements[indices[iI]];
+
+            const radius2 = VdWLookup[atomRadius[bI]];
+            if (aI === bI || radius2 === VdWLookup[0]) continue;
+
+            const cutoff2 = (cutoff1 + radius2) * (cutoff1 + radius2);
+            if (squaredDistances[iI] < cutoff2) {
+                neighbors[neighbors.length] = bI;
+            }
+        }
+
+        // for all neighbors: test all sphere points
+        position(aI, aPos);
+        const scalar = probeSize + radius1;
+        let accessiblePointCount = 0;
+        for (let sI = 0; sI < spherePoints.length; ++sI) {
+            const spherePoint = spherePoints[sI];
+            testPoint = [spherePoint[0] * scalar + aPos[0],
+                spherePoint[1] * scalar + aPos[1],
+                spherePoint[2] * scalar + aPos[2]] as Vec3;
+            let accessible = true;
+
+            for (let _nI = 0; _nI < neighbors.length; ++_nI) {
+                const nI = neighbors[_nI];
+                position(nI, bPos);
+                const radius3 = VdWLookup[atomRadius[nI]];
+                const cutoff3 = (radius3 + probeSize) * (radius3 + probeSize);
+                if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) {
+                    accessible = false;
+                    break;
+                }
+            }
+
+            if (accessible) ++accessiblePointCount;
+        }
+
+        accessibleSurfaceArea[residueAtomSegments.index[aI]] += cons * accessiblePointCount * scalar * scalar;
+    }
+}

src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts

@@ -0,0 +1,119 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
+ */
+
+import { ShrakeRupleyContext, VdWLookup } from './common';
+import { getElementIdx, isHydrogen } from '../../../../mol-model/structure/structure/unit/links/common';
+import { isPolymer, isNucleic, MoleculeType, ElementSymbol } from '../../../../mol-model/structure/model/types';
+import { VdwRadius } from '../../../../mol-model/structure/model/properties/atomic';
+
+export async function assignRadiusForHeavyAtoms(ctx: ShrakeRupleyContext) {
+    const { updateChunk, atomRadius } = ctx;
+    for (let i = 0; i < atomRadius.length; i += updateChunk) {
+        computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length));
+    }
+}
+
+function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) {
+    const { structure } = ctx;
+    const { model } = structure;
+    const { atoms: atomInfo, derived, residues, residueAtomSegments } = model.atomicHierarchy;
+    const { label_comp_id } = residues;
+    const { moleculeType } = derived.residue;
+    const { type_symbol, label_atom_id } = atomInfo;
+
+    for (let aI = begin; aI < end; ++aI) {
+        const rI = residueAtomSegments.index[aI];
+        const element = type_symbol.value(aI);
+        const elementIdx = getElementIdx(element);
+        // skip hydrogen atoms
+        if (isHydrogen(elementIdx)) {
+            ctx.atomRadius[aI] = VdWLookup[0];
+            continue;
+        }
+
+        const residueType = moleculeType[rI];
+        // skip non-polymer groups
+        if (!ctx.nonPolymer) {
+            if (!isPolymer(residueType)) {
+                ctx.atomRadius[aI] = VdWLookup[0];
+                continue;
+            }
+        }
+
+        const atomId = label_atom_id.value(aI);
+        let compId = label_comp_id.value(rI);
+
+        // handle modified residues
+        const parentId = model.properties.modifiedResidues.parentId.get(compId);
+        if (parentId !== void 0) compId = parentId;
+
+        if (isNucleic(residueType)) {
+             ctx.atomRadius[aI] = determineRadiusNucl(atomId, element, compId);
+        } else if (residueType === MoleculeType.Protein) {
+            ctx.atomRadius[aI] = determineRadiusAmino(atomId, element, compId);
+        } else {
+            ctx.atomRadius[aI] = handleNonStandardCase(element);
+        }
+    }
+}
+
+/**
+ * Gets the van der Waals radius of the given atom following the values defined by Chothia (1976)
+ * J.Mol.Biol.105,1-14. NOTE: the vdw values defined by the paper assume no Hydrogens and thus "inflates" slightly
+ * the heavy atoms to account for Hydrogens.
+ */
+function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number {
+    switch (element) {
+        case 'O':
+        return 5;
+        case 'S':
+        return 6;
+        case 'N':
+        return atomId === 'NZ' ? 4 : 3;
+        case 'C':
+        switch (atomId) {
+            case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
+            return 1;
+            case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
+            return 2;
+            default:
+            switch (compId) {
+                case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
+                return 1;
+                case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
+                return 2;
+                case 'GLU': case 'GLN':
+                return atomId === 'CD' ? 1 : 2;
+            }
+        }
+    }
+    return handleNonStandardCase(element);
+}
+
+function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
+    switch (element) {
+        case 'C':
+        return 7;
+        case 'N':
+        return 8;
+        case 'P':
+        return 9;
+        case 'O':
+        return 5;
+    }
+    return handleNonStandardCase(element);
+}
+
+function handleNonStandardCase(element: ElementSymbol): number {
+    const radius = VdwRadius(element);
+    let index = VdWLookup.indexOf(radius);
+    if (index === -1) {
+        // add novel value to lookup array
+        index = VdWLookup.length;
+        VdWLookup[index] = radius;
+    }
+    return index;
+}

src/mol-model/structure/structure/accessible-surface-area.ts

@@ -1,336 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import Structure from './structure';
-import { Task, RuntimeContext } from '../../../mol-task';
-import { BitFlags } from '../../../mol-util';
-import { ParamDefinition as PD } from '../../../mol-util/param-definition'
-import { Vec3 } from '../../../mol-math/linear-algebra';
-import { isPolymer, ElementSymbol, isNucleic, MoleculeType } from '../model/types';
-import { VdwRadius } from '../model/properties/atomic';
-import { isHydrogen, getElementIdx } from './unit/links/common';
-
-namespace AccessibleSurfaceArea {
-    // Chothia's amino acid and nucleotide atom vdw radii
-    export const VdWLookup = [
-        -1.0, // 0: missing
-        1.76, // 1: trigonal C
-        1.87, // 2: tetrahedral C
-        1.65, // 3: trigonal N
-        1.50, // 4: tetrahedral N
-        1.40, // 5: O
-        1.85, // 6: S
-        1.80, // 7: C (nucleic)
-        1.60, // 8: N (nucleic)
-        1.40  // 9: P (nucleic)
-    ] // can still be appended on-the-fly for rare elements like selenium
-
-    /**
-     * Adapts the BioJava implementation by Jose Duarte. That implementation is based on the publication by Shrake, A., and
-     * J. A. Rupley. "Environment and Exposure to Solvent of Protein Atoms. Lysozyme and Insulin." JMB (1973).
-     */
-    export function compute(structure: Structure,
-        params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>> = {}) {
-        params = { ...PD.getDefaultValues(AccessibleSurfaceAreaComputationParams), ...params };
-        return Task.create('Compute Accessible Surface Area', async rtctx => {
-            return await _compute(rtctx, structure, params);
-        }).run();
-    }
-
-    async function _compute(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>> = {}): Promise<AccessibleSurfaceArea> {
-        const ctx = initialize(rtctx, structure, params);
-
-        assignRadiusForHeavyAtoms(ctx);
-        computePerResidue(ctx);
-        normalizeAccessibleSurfaceArea(ctx);
-
-        return {
-            accessibleSurfaceArea: ctx.accessibleSurfaceArea!,
-            relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea!,
-            buried: (index: number) => ctx.relativeAccessibleSurfaceArea![index] < 0.16 // TODO this doesnt seem super elegant
-        };
-    }
-
-    interface AccessibleSurfaceAreaContext {
-        rtctx: RuntimeContext,
-        structure: Structure,
-        params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>>,
-        spherePoints: Vec3[],
-        cons: number,
-        maxLookupRadius: number,
-        atomRadius?: Int8Array,
-        accessibleSurfaceArea?: Float32Array,
-        relativeAccessibleSurfaceArea?: Float32Array
-    }
-
-    function normalizeAccessibleSurfaceArea(ctx: AccessibleSurfaceAreaContext) {
-        const { accessibleSurfaceArea, relativeAccessibleSurfaceArea, structure } = ctx;
-        const { residues, derived } = structure.model.atomicHierarchy;
-
-        for (let i = 0; i < residues.label_comp_id.rowCount; ++i) {
-            // skip entities not part of a polymer chain
-            if (!ctx.params.nonPolymer) {
-                if (!isPolymer(derived.residue.moleculeType[i])) continue;
-            }
-
-            const maxAsa = (MaxAsa as any)[residues.label_comp_id.value(i)];
-            const rasa = accessibleSurfaceArea![i] / (maxAsa === undefined ? DefaultMaxAsa : maxAsa);
-            relativeAccessibleSurfaceArea![i] = rasa;
-        }
-    }
-
-    async function computePerResidue(ctx: AccessibleSurfaceAreaContext) {
-        const { structure, atomRadius, accessibleSurfaceArea, spherePoints, cons, params, maxLookupRadius } = ctx;
-        const { probeSize } = params;
-        const { model, elementCount: atomCount } = structure;
-        const { x, y, z } = model.atomicConformation;
-        const { residueAtomSegments } = model.atomicHierarchy;
-        const { lookup3d } = structure;
-
-        const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
-        const aPos = Vec3.zero();
-        const bPos = Vec3.zero();
-        let testPoint = Vec3.zero();
-
-        for (let aI = 0; aI < atomCount; ++aI) {
-            if (aI % 10000 === 0) {
-                if (ctx.rtctx.shouldUpdate) {
-                    await ctx.rtctx.update({ message: 'calculating accessible surface area', current: aI, max: atomCount });
-                }
-            }
-
-            const radius1 = VdWLookup[atomRadius![aI]];
-            if (radius1 === VdWLookup[0]) continue;
-
-            // pre-filter by lookup3d
-            // 36275 ms - lookup ~3000 ms
-            const { count, units, indices, squaredDistances } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
-
-            // collect neighbors for each atom
-            const cutoff1 = probeSize! + probeSize! + radius1;
-            const neighbors = [];
-            for (let iI = 0; iI < count; ++iI) {
-                const bI = units[iI].elements[indices[iI]];
-                const radius2 = VdWLookup[atomRadius![bI]];
-                if (aI === bI || radius2 === VdWLookup[0]) continue;
-
-                const cutoff2 = (cutoff1 + radius2) * (cutoff1 + radius2);
-                if (squaredDistances[iI] < cutoff2) {
-                    neighbors[neighbors.length] = bI;
-                }
-            }
-
-            // for all neighbors: test all sphere points
-            position(aI, aPos);
-            const scalar = probeSize! + radius1;
-            let accessiblePointCount = 0;
-            for (let sI = 0; sI < spherePoints.length; ++sI) {
-                const spherePoint = spherePoints[sI];
-                testPoint = [spherePoint[0] * scalar + aPos[0], spherePoint[1] * scalar + aPos[1], spherePoint[2] * scalar + aPos[2]] as Vec3;
-                let accessible = true;
-
-                for (let _nI = 0; _nI < neighbors.length; ++_nI) {
-                    const nI = neighbors[_nI];
-                    position(nI, bPos);
-                    const radius3 = VdWLookup[atomRadius![nI]];
-                    const cutoff3 = (radius3 + probeSize!) * (radius3 + probeSize!);
-                    if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) {
-                        accessible = false;
-                        break;
-                    }
-                }
-
-                if (accessible) ++accessiblePointCount;
-            }
-
-            accessibleSurfaceArea![residueAtomSegments.index[aI]] += cons * accessiblePointCount * scalar * scalar;
-        }
-    }
-
-    function assignRadiusForHeavyAtoms(ctx: AccessibleSurfaceAreaContext) {
-        const { structure } = ctx;
-        const { model, elementCount: atomCount } = structure;
-        const { atoms: atomInfo, derived, residues, residueAtomSegments } = model.atomicHierarchy;
-        const { label_comp_id } = residues;
-        const { moleculeType } = derived.residue;
-        const { type_symbol, label_atom_id } = atomInfo;
-        const residueCount = moleculeType.length;
-
-        // with atom and residue count at hand: initialize arrays
-        ctx.atomRadius = new Int8Array(atomCount - 1);
-        ctx.accessibleSurfaceArea = new Float32Array(residueCount - 1);
-        ctx.relativeAccessibleSurfaceArea = new Float32Array(residueCount - 1);
-
-        for (let aI = 0; aI < atomCount; ++aI) {
-            const rI = residueAtomSegments.index[aI];
-            const element = type_symbol.value(aI);
-            const elementIdx = getElementIdx(element);
-            // skip hydrogen atoms
-            if (isHydrogen(elementIdx)) {
-                ctx.atomRadius[aI] = VdWLookup[0];
-                continue;
-            }
-
-            const residueType = moleculeType[rI];
-            // skip non-polymer groups
-            if (!ctx.params.nonPolymer) {
-                if (!isPolymer(residueType)) {
-                    ctx.atomRadius[aI] = VdWLookup[0];
-                    continue;
-                }
-            }
-
-            const atomId = label_atom_id.value(aI);
-            let compId = label_comp_id.value(rI);
-
-            // handle modified residues
-            const parentId = model.properties.modifiedResidues.parentId.get(compId);
-            if (parentId !== void 0) compId = parentId;
-
-            if (isNucleic(residueType)) {
-                 ctx.atomRadius[aI] = determineRadiusNucl(atomId, element, compId);
-            } else if (residueType === MoleculeType.Protein) {
-                ctx.atomRadius[aI] = determineRadiusAmino(atomId, element, compId);
-            } else {
-                ctx.atomRadius[aI] = handleNonStandardCase(element);
-            }
-        }
-    }
-
-    /**
-     * Gets the van der Waals radius of the given atom following the values defined by Chothia (1976)
-     * J.Mol.Biol.105,1-14. NOTE: the vdw values defined by the paper assume no Hydrogens and thus "inflates" slightly
-     * the heavy atoms to account for Hydrogens.
-     */
-    function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number {
-        switch (element) {
-            case 'O':
-            return 5;
-            case 'S':
-            return 6;
-            case 'N':
-            return atomId === 'NZ' ? 4 : 3;
-            case 'C':
-            switch (atomId) {
-                case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
-                return 1;
-                case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
-                return 2;
-                default:
-                switch (compId) {
-                    case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
-                    return 1;
-                    case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
-                    return 2;
-                    case 'GLU': case 'GLN':
-                    return atomId === 'CD' ? 1 : 2;
-                }
-            }
-        }
-        return handleNonStandardCase(element);
-    }
-
-    function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
-        switch (element) {
-            case 'C':
-            return 7;
-            case 'N':
-            return 8;
-            case 'P':
-            return 9;
-            case 'O':
-            return 5;
-        }
-        return handleNonStandardCase(element);
-    }
-
-    function handleNonStandardCase(element: ElementSymbol): number {
-        const radius = VdwRadius(element);
-        let index = VdWLookup.indexOf(radius);
-        if (index === -1) {
-            // add novel value to lookup array
-            index = VdWLookup.length;
-            VdWLookup[index] = radius;
-        }
-        return index;
-    }
-
-    function initialize(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<AccessibleSurfaceAreaComputationParams>>): AccessibleSurfaceAreaContext {
-        return {
-            rtctx: rtctx,
-            structure: structure,
-            params: params,
-            spherePoints: generateSpherePoints(params.numberOfSpherePoints!),
-            cons: 4.0 * Math.PI / params.numberOfSpherePoints!,
-            maxLookupRadius: 2 * params.probeSize! + 2 * VdWLookup[2] // 2x probe size + 2x largest VdW
-        }
-    }
-
-    /** Creates a collection of points on a sphere by the Golden Section Spiral algorithm. */
-    function generateSpherePoints(numberOfSpherePoints: number): Vec3[] {
-        const points: Vec3[] = [];
-        const inc = Math.PI * (3.0 - Math.sqrt(5.0));
-        const offset = 2.0 / numberOfSpherePoints;
-        for (let k = 0; k < numberOfSpherePoints; ++k) {
-            const y = k * offset - 1.0 + (offset / 2.0);
-            const r = Math.sqrt(1.0 - y * y);
-            const phi = k * inc;
-            points[points.length] = [Math.cos(phi) * r, y, Math.sin(phi) * r] as Vec3;
-        }
-        return points;
-    }
-
-    /** Maximum accessible surface area observed for amino acids. Taken from: http://dx.doi.org/10.1371/journal.pone.0080635 */
-    export const MaxAsa = {
-        'ALA': 121.0,
-        'ARG': 265.0,
-        'ASN': 187.0,
-        'ASP': 187.0,
-        'CYS': 148.0,
-        'GLU': 214.0,
-        'GLN': 214.0,
-        'GLY': 97.0,
-        'HIS': 216.0,
-        'ILE': 195.0,
-        'LEU': 191.0,
-        'LYS': 230.0,
-        'MET': 203.0,
-        'PHE': 228.0,
-        'PRO': 154.0,
-        'SER': 143.0,
-        'THR': 163.0,
-        'TRP': 264.0,
-        'TYR': 255.0,
-        'VAL': 165.0
-    }
-    export const DefaultMaxAsa = 121.0
-
-    export const AccessibleSurfaceAreaComputationParams = {
-        numberOfSpherePoints: PD.Numeric(92, {}, { description: 'number of sphere points to sample per atom: 92 (original paper), 960 (BioJava), 3000 (EPPIC) - see Shrake A, Rupley JA: Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 1973.' }),
-        probeSize: PD.Numeric(1.4, {}, { description: 'corresponds to the size of a water molecule: 1.4 (original paper), 1.5 (occassionally used)' }),
-        buriedRasaThreshold: PD.Numeric(0.16, { min: 0.0, max: 1.0 }, { description: 'below this cutoff of relative accessible surface area a residue will be considered buried - see: Rost B, Sander C: Conservation and prediction of solvent accessibility in protein families. Proteins 1994.' }),
-        nonPolymer: PD.Boolean(false, { description: 'Include non-polymer atoms in computation.' })
-    }
-    export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams
-
-    export namespace SolventAccessibility {
-        export const is: (t: number, f: Flag) => boolean = BitFlags.has
-        export const create: (f: Flag) => number = BitFlags.create
-        export const enum Flag {
-            _ = 0x0,
-            BURIED = 0x1,
-            ACCESSIBLE = 0x2
-        }
-    }
-}
-
-interface AccessibleSurfaceArea {
-    readonly accessibleSurfaceArea?: ArrayLike<number>
-    readonly relativeAccessibleSurfaceArea?: ArrayLike<number>
-    buried(index: number): boolean
-}
-
-export { AccessibleSurfaceArea }

src/mol-theme/color/accessible-surface-area.ts

@@ -0,0 +1,71 @@
+/**
+ * Copyright (c) 2019 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 { ThemeDataContext } from '../theme';
+import { ColorTheme, LocationColor } from '../color';
+import { ParamDefinition as PD } from '../../mol-util/param-definition'
+import { ColorScale, Color } from '../../mol-util/color';
+import { Unit, StructureElement } from '../../mol-model/structure';
+import { Location } from '../../mol-model/location';
+import { ComputedAccessibleSurfaceArea } from '../../mol-model-props/computed/accessible-surface-area';
+import { ColorListName, ColorListOptionsScale } from '../../mol-util/color/lists';
+import { VdWLookup } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley/common';
+
+const DefaultColor = Color(0xFFFFFF)
+const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
+
+export const AccessibleSurfaceAreaColorThemeParams = {
+    list: PD.ColorList<ColorListName>('rainbow', ColorListOptionsScale)
+}
+export type AccessibleSurfaceAreaColorThemeParams = typeof AccessibleSurfaceAreaColorThemeParams
+export function getAccessibleSurfaceAreaColorThemeParams(ctx: ThemeDataContext) {
+    return AccessibleSurfaceAreaColorThemeParams // TODO return copy
+}
+export function AccessibleSurfaceAreaColorTheme(ctx: ThemeDataContext, props: PD.Values<AccessibleSurfaceAreaColorThemeParams>): ColorTheme<AccessibleSurfaceAreaColorThemeParams> {
+    let color: LocationColor
+
+    const scale = ColorScale.create({
+        listOrName: props.list,
+        minLabel: 'buried',
+        maxLabel: 'exposed',
+        domain: [0.0, 1.0]
+    })
+
+    const accessibleSurfaceArea = ctx.structure ? ComputedAccessibleSurfaceArea.get(ctx.structure)!.asa : undefined
+
+    if (accessibleSurfaceArea) {
+        color = (location: Location): Color => {
+            if (StructureElement.Location.is(location)) {
+                if (Unit.isAtomic(location.unit)) {
+                    const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
+                    return value !== VdWLookup[0] /* signals missing value */ ? scale.color(value) : DefaultColor;
+                }
+            }
+
+            return DefaultColor
+        }
+    } else {
+        color = () => DefaultColor
+    }
+
+    return {
+        factory: AccessibleSurfaceAreaColorTheme,
+        granularity: 'group',
+        color,
+        props,
+        description: Description,
+        legend: scale ? scale.legend : undefined
+    }
+}
+
+export const AccessibleSurfaceAreaColorThemeProvider: ColorTheme.Provider<AccessibleSurfaceAreaColorThemeParams> = {
+    label: 'Accessible Surface Area',
+    factory: AccessibleSurfaceAreaColorTheme,
+    getParams: getAccessibleSurfaceAreaColorThemeParams,
+    defaultValues: PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams),
+    isApplicable: (ctx: ThemeDataContext) => !!ctx.structure
+}

src/tests/browser/render-asa.ts

@@ -5,20 +5,21 @@
  */
 
 import './index.html'
-import { resizeCanvas } from '../../mol-canvas3d/util';
 import { Canvas3D } from '../../mol-canvas3d/canvas3d';
-import { CIF, CifFrame } from '../../mol-io/reader/cif';
-import { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif';
+import { CifFrame, CIF } from '../../mol-io/reader/cif'
 import { Model, Structure, StructureElement, Unit } from '../../mol-model/structure';
 import { ColorTheme, LocationColor } from '../../mol-theme/color';
 import { SizeTheme } from '../../mol-theme/size';
 import { CartoonRepresentationProvider } from '../../mol-repr/structure/representation/cartoon';
-import { AccessibleSurfaceArea } from '../../mol-model/structure/structure/accessible-surface-area';
+import { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif';
 import { Color, ColorScale } from '../../mol-util/color';
-import { ColorListName, ColorListOptions } from '../../mol-util/color/lists';
+import { Location } from '../../mol-model/location';
 import { ThemeDataContext } from '../../mol-theme/theme';
 import { ParamDefinition as PD } from '../../mol-util/param-definition';
-import { Location } from '../../mol-model/location';
+import { AccessibleSurfaceArea } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley';
+import { VdWLookup } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley/common';
+import { resizeCanvas } from '../../mol-canvas3d/util';
+import { ColorListName, ColorListOptions } from '../../mol-util/color/lists';
 
 const parent = document.getElementById('app')!
 parent.style.width = '100%'
@@ -73,14 +74,16 @@ async function init(props = {}) {
         // '1aon'
         // '1acj'
         // '1pga'
-        '1brr'
-        // '1hrc'
+        // '1brr'
+        '1hrc'
     )
     const models = await getModels(cif)
     const structure = await getStructure(models[0])
 
     // async compute ASA
+    console.time('computeASA')
     accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure)
+    console.timeEnd('computeASA')
 
     const cartoonRepr = getCartoonRepr()
 
@@ -116,11 +119,12 @@ export function AccessibleSurfaceAreaColorTheme(ctx: ThemeDataContext, props: PD
         maxLabel: '1.0 (exposed)',
         domain: [0.0, 1.0]
     })
+
     color = (location: Location): Color => {
         if (StructureElement.Location.is(location)) {
             if (Unit.isAtomic(location.unit)) {
-                const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea![location.unit.residueIndex[location.element]];
-                return value !== AccessibleSurfaceArea.VdWLookup[0] /* signals missing value */ ? scale.color(value) : DefaultColor;
+                const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
+                return value !== VdWLookup[0] /* signals missing value */ ? scale.color(value) : DefaultColor;
             }
         }