diff --git a/src/mol-model/structure/structure/accessible-surface-area.ts b/src/mol-model/structure/structure/accessible-surface-area.ts
new file mode 100644
index 0000000000000000000000000000000000000000..5735b2df5ec4663b18693b5c09bd80eb09e1dfdf
--- /dev/null
+++ b/src/mol-model/structure/structure/accessible-surface-area.ts
@@ -0,0 +1,321 @@
+/**
+ * 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 atoms vdw radii
+ const trigonalCarbonVdw = 1.76;
+ const tetrahedralCarbonVdw = 1.87;
+ const trigonalNitrogenVdw = 1.65;
+ const tetrahedralNitrogenVdw = 1.50;
+ // deviating radii from definition in types.ts
+ const oxygenVdw = 1.40;
+ const sulfurVdw = 1.85;
+ // Chothia's nucleotide atom vdw radii
+ const nucCarbonVdw = 1.80;
+ const nucNitrogenVdw = 1.60;
+ const nucPhosphorusVdw = 1.40;
+ export const missingValue = -1.0;
+
+ /**
+ * 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 {
+ atomRadius: ctx.atomRadius!,
+ 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?: Float32Array, // TODO there are only 5-10 unique values in this array - rather than storing values, a int pointing to a dictionary will be far more memory efficient
+ 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();
+
+ for (let aI = 0; aI < atomCount; ++aI) {
+ if (aI % 10000 === 0) {
+ console.log(`calculating accessible surface area, current: ${ aI }, max: ${ atomCount }`);
+ if (ctx.rtctx.shouldUpdate) {
+ await ctx.rtctx.update({ message: 'calculating accessible surface area', current: aI, max: atomCount });
+ }
+ }
+
+ const radius1 = atomRadius![aI];
+ if (radius1 === missingValue) 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 = atomRadius![bI];
+ if (aI === bI || radius2 === missingValue) 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];
+ const 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 cutoff3 = (atomRadius![nI] + probeSize!) * (atomRadius![nI] + 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 Float32Array(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] = missingValue;
+ continue;
+ }
+
+ const residueType = moleculeType[rI];
+ // skip non-polymer groups
+ if (!ctx.params.nonPolymer) {
+ if (!isPolymer(residueType)) {
+ ctx.atomRadius[aI] = missingValue;
+ continue;
+ }
+ }
+
+ const atomId = label_atom_id.value(aI);
+ const residueId = label_comp_id.value(rI);
+ if (isNucleic(residueType)) {
+ ctx.atomRadius[aI] = determineRadiusNucl(atomId, element, residueId);
+ } else if (residueType === MoleculeType.protein) {
+ ctx.atomRadius[aI] = determineRadiusAmino(atomId, element, residueId);
+ } else {
+ ctx.atomRadius[aI] = VdwRadius(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 oxygenVdw;
+ case 'S':
+ return sulfurVdw;
+ case 'N':
+ return atomId === 'NZ' ? tetrahedralNitrogenVdw : trigonalNitrogenVdw;
+ case 'C':
+ switch (atomId) {
+ case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
+ return trigonalCarbonVdw;
+ case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
+ return tetrahedralCarbonVdw;
+ default:
+ switch (compId) {
+ case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
+ return trigonalCarbonVdw;
+ case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
+ return tetrahedralCarbonVdw;
+ case 'GLU': case 'GLN':
+ return atomId === 'CD' ? trigonalCarbonVdw : tetrahedralCarbonVdw;
+ }
+ }
+ }
+ return VdwRadius(element);
+ }
+
+ function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
+ switch (element) {
+ case 'C':
+ return nucCarbonVdw;
+ case 'N':
+ return nucNitrogenVdw;
+ case 'P':
+ return nucPhosphorusVdw;
+ case 'O':
+ return oxygenVdw;
+ }
+ return VdwRadius(element);
+ }
+
+ 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 * tetrahedralCarbonVdw // 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 atomRadius?: ArrayLike<number>
+ readonly accessibleSurfaceArea?: ArrayLike<number>
+ readonly relativeAccessibleSurfaceArea?: ArrayLike<number>
+ buried(index: number): boolean
+}
+
+export { AccessibleSurfaceArea }
\ No newline at end of file
diff --git a/src/mol-model/structure/structure/unit.ts b/src/mol-model/structure/structure/unit.ts
index a66e32e6299003abbd7045ac4b07a1ecae6220c8..28869515a8ab573b618a463aaaf92f5b3b84927c 100644
--- a/src/mol-model/structure/structure/unit.ts
+++ b/src/mol-model/structure/structure/unit.ts
@@ -18,8 +18,6 @@ import { IntMap, SortedArray } from 'mol-data/int';
import { hash2, hashFnv32a } from 'mol-data/util';
import { getAtomicPolymerElements, getCoarsePolymerElements, getAtomicGapElements, getCoarseGapElements } from './util/polymer';
import { getNucleotideElements } from './util/nucleotide';
-import { computeAccessibleSurfaceArea } from './unit/accessible-surface-area/compute';
-import { AccessibleSurfaceArea } from './unit/accessible-surface-area/data';
/**
* A building block of a structure that corresponds to an atomic or
@@ -193,12 +191,6 @@ namespace Unit {
return this.props.nucleotideElements.ref;
}
- get accessibleSurfaceArea() {
- if (this.props.accessibleSurfaceArea.ref) return this.props.accessibleSurfaceArea.ref;
- this.props.accessibleSurfaceArea.ref = computeAccessibleSurfaceArea(this);
- return this.props.accessibleSurfaceArea.ref;
- }
-
getResidueIndex(elementIndex: StructureElement.UnitIndex) {
return this.model.atomicHierarchy.residueAtomSegments.index[this.elements[elementIndex]];
}
@@ -223,7 +215,6 @@ namespace Unit {
polymerElements: ValueRef<SortedArray<ElementIndex> | undefined>
gapElements: ValueRef<SortedArray<ElementIndex> | undefined>
nucleotideElements: ValueRef<SortedArray<ElementIndex> | undefined>
- accessibleSurfaceArea: ValueRef<AccessibleSurfaceArea | undefined>
}
function AtomicProperties(): AtomicProperties {
@@ -233,8 +224,7 @@ namespace Unit {
rings: ValueRef.create(void 0),
polymerElements: ValueRef.create(void 0),
gapElements: ValueRef.create(void 0),
- nucleotideElements: ValueRef.create(void 0),
- accessibleSurfaceArea: ValueRef.create(void 0)
+ nucleotideElements: ValueRef.create(void 0)
};
}
diff --git a/src/mol-model/structure/structure/unit/accessible-surface-area/compute.ts b/src/mol-model/structure/structure/unit/accessible-surface-area/compute.ts
deleted file mode 100644
index 873bffead0a8758a0dd7a7aa9032766f4aed9edc..0000000000000000000000000000000000000000
--- a/src/mol-model/structure/structure/unit/accessible-surface-area/compute.ts
+++ /dev/null
@@ -1,265 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import Unit from '../../unit'
-import { Vec3 } from 'mol-math/linear-algebra';
-import { AccessibleSurfaceAreaComputationParams, AccessibleSurfaceArea, SolventAccessibility } from './data';
-import { isHydrogen, getElementIdx } from '../links/common'; // TODO these functions are relevant for many tasks: move them somewhere actually common
-import { ElementSymbol, MaxAsa, DefaultMaxAsa, isPolymer, isNucleic, MoleculeType } from 'mol-model/structure/model/types';
-import { VdwRadius } from 'mol-model/structure/model/properties/atomic/measures';
-import { ParamDefinition as PD } from 'mol-util/param-definition'
-
-// Chothia's amino acid atoms vdw radii
-const trigonalCarbonVdw = 1.76;
-const tetrahedralCarbonVdw = 1.87;
-const trigonalNitrogenVdw = 1.65;
-const tetrahedralNitrogenVdw = 1.50;
-/** deviating radii from definition in types.ts */
-const oxygenVdw = 1.40;
-const sulfurVdw = 1.85;
-// Chothia's nucleotide atom vdw radii
-const nucCarbonVdw = 1.80;
-const nucNitrogenVdw = 1.60;
-const nucPhosphorusVdw = 1.40;
-const missingAccessibleSurfaceAreaValue = -1.0;
-
-interface AccessibleSurfaceAreaContext {
- unit: Unit.Atomic,
- params: PD.Values<AccessibleSurfaceAreaComputationParams>,
- spherePoints: Vec3[],
- cons: number,
- maxLookupRadius: number,
- atomRadius?: Float32Array,
- accessibleSurfaceArea?: Float32Array,
- relativeAccessibleSurfaceArea?: Float32Array,
- buried?: Uint8Array
-}
-
-/**
- * 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).
- */
-function computeAccessibleSurfaceArea(unit: Unit.Atomic, params?: PD.Values<AccessibleSurfaceAreaComputationParams>): AccessibleSurfaceArea {
- if (!params) params = PD.getDefaultValues(AccessibleSurfaceAreaComputationParams);
-
- // TODO non-polymer flag is currently useless as hetatms are located in different units - aim is not to color them, but to compute values correctly - relates to correct ASA computation for inter-chain contacts
- console.log(`computing accessible surface area for unit #${ unit.id + 1 } - ${ params.numberOfSpherePoints } points, ${ params.probeSize } probe size, ${ params.nonPolymer ? 'honoring' : 'ignoring'} non-polymer atoms`);
-
- const ctx = initialize(unit, params);
- assignRadiusForHeavyAtoms(ctx);
- computePerResidue(ctx);
- normalizeAccessibleSurfaceArea(ctx);
-
- return {
- atomRadius: ctx.atomRadius!,
- accessibleSurfaceArea: ctx.accessibleSurfaceArea!,
- relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea!,
- buried: ctx.buried!
- };
-}
-
-function normalizeAccessibleSurfaceArea(ctx: AccessibleSurfaceAreaContext) {
- const { residues, derived } = ctx.unit.model.atomicHierarchy;
- const { accessibleSurfaceArea, relativeAccessibleSurfaceArea } = ctx;
-
- 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;
- ctx.buried![i] |= (rasa < ctx.params.buriedRasaThreshold ? SolventAccessibility.Flag.BURIED : SolventAccessibility.Flag.ACCESSIBLE)
- }
-}
-
-/**
- * notes on performance - scenario: compute for first 10 units of 3j3q @ 960 sphere points
- * lookup3d + refinement: ~5000ms
- * naive approach: ~5600ms - higher variance
- */
-function computePerResidue(ctx: AccessibleSurfaceAreaContext) { // runs at roughly 5000 ms
- const { atomRadius, accessibleSurfaceArea, maxLookupRadius, spherePoints, cons } = ctx;
- const { probeSize } = ctx.params;
- const { elements: atoms, residueIndex } = ctx.unit;
- const { x, y, z } = ctx.unit.model.atomicConformation;
- const atomCount = ctx.unit.elements.length;
- const { lookup3d } = ctx.unit;
-
- const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]);
- const a1Pos = Vec3.zero();
- const a2Pos = Vec3.zero();
-
- for (let _aI = 0; _aI < atomCount; ++_aI) {
- // map the atom index of this unit to the actual 'global' atom index
- const aI = atoms[_aI];
- const radii1 = atomRadius![aI];
- if (radii1 === missingAccessibleSurfaceAreaValue) continue;
-
- // find suitable neighbor candidates by lookup
- const { indices, count } = lookup3d.find(x[aI], y[aI], z[aI], maxLookupRadius);
- position(aI, a1Pos);
-
- // refine list by actual criterion
- const cutoff = probeSize + probeSize + radii1;
- const filteredIndicies = []; // TODO might be better to use IntArray here and reuse it - how to find safe upper limit of possible neighborhood count - BioJava mentions 60 as relatively safe upper bound
- for (let ni = 0; ni < count; ni++) {
- const _bI = indices[ni];
- const bI = atoms[_bI];
- const radii2 = atomRadius![bI];
- if (bI === aI || radii2 === missingAccessibleSurfaceAreaValue) continue;
-
- const cutoff2 = (cutoff + radii2) * (cutoff + radii2);
- // accurately check for neighborhood
- position(bI, a2Pos);
- if (Vec3.squaredDistance(a1Pos, a2Pos) < cutoff2) {
- filteredIndicies[filteredIndicies.length] = bI;
- }
- }
-
- // test sphere points
- const r = probeSize + radii1;
- let accessiblePointCount = 0;
- for (let si = 0; si < spherePoints.length; ++si) {
- const spherePoint = spherePoints[si];
- const testPoint = [spherePoint[0] * r + a1Pos[0], spherePoint[1] * r + a1Pos[1], spherePoint[2] * r + a1Pos[2]] as Vec3;
- let accessible = true;
-
- for (let ni = 0; ni < filteredIndicies.length; ++ni) {
- const naI = filteredIndicies[ni];
- position(naI, a2Pos);
- const cutoff3 = (atomRadius![naI] + probeSize) * (atomRadius![naI] + probeSize);
- if (Vec3.squaredDistance(testPoint, a2Pos) < cutoff3) {
- accessible = false;
- break;
- }
- }
-
- if (accessible) ++accessiblePointCount;
- }
-
- const value = cons * accessiblePointCount * r * r;
- accessibleSurfaceArea![residueIndex[aI]] += value;
- // +30% computation by normalizing partial solutions
- // relativeAccessibleSurfaceArea[residueIndex[aI]] += value * (NormalizationFactors as any)[residueIndex[aI]];
- }
-}
-
-function assignRadiusForHeavyAtoms(ctx: AccessibleSurfaceAreaContext) {
- const atomCount = ctx.unit.elements.length;
- const { elements: atoms, residueIndex } = ctx.unit;
- const { residues } = ctx.unit.model.atomicHierarchy;
- const { moleculeType } = ctx.unit.model.atomicHierarchy.derived.residue;
- const { type_symbol, label_atom_id } = ctx.unit.model.atomicHierarchy.atoms;
- const { label_comp_id } = ctx.unit.model.atomicHierarchy.residues;
-
- const residueCount = residues.label_comp_id.rowCount;
- ctx.atomRadius = new Float32Array(atomCount - 1);
- ctx.accessibleSurfaceArea = new Float32Array(residueCount - 1);
- ctx.relativeAccessibleSurfaceArea = new Float32Array(residueCount - 1);
- ctx.buried = new Uint8Array(residueCount - 1);
-
- for (let _aI = 0; _aI < atomCount; ++_aI) {
- const aI = atoms[_aI];
- const raI = residueIndex[aI];
- const aeI = getElementIdx(type_symbol.value(aI)!);
-
- // skip hydrogen atoms
- if (isHydrogen(aeI)) {
- ctx.atomRadius[aI] = missingAccessibleSurfaceAreaValue;
- continue;
- }
-
- // skip non-polymer groups
- if (!ctx.params.nonPolymer) {
- if (!isPolymer(moleculeType[raI])) {
- ctx.atomRadius[aI] = missingAccessibleSurfaceAreaValue;
- continue;
- }
- }
-
- const atomId = label_atom_id.value(aI);
- const element = type_symbol.value(aI);
- const resn = label_comp_id.value(raI)!;
-
- ctx.atomRadius[aI] = isNucleic(moleculeType[raI]) ? determineRadiusNucl(atomId, element, resn) : moleculeType[raI] === MoleculeType.protein ? determineRadiusAmino(atomId, element, resn) : VdwRadius(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 oxygenVdw;
- case 'S':
- return sulfurVdw;
- case 'N':
- return atomId === 'NZ' ? tetrahedralNitrogenVdw : trigonalNitrogenVdw;
- case 'C':
- switch (atomId) {
- case 'C': case 'CE1': case'CE2': case 'CE3': case 'CH2': case 'CZ': case 'CZ2': case 'CZ3':
- return trigonalCarbonVdw;
- case 'CA': case 'CB': case 'CE': case 'CG1': case 'CG2':
- return tetrahedralCarbonVdw;
- default:
- switch (compId) {
- case 'PHE': case 'TRP': case 'TYR': case 'HIS': case 'ASP': case 'ASN':
- return trigonalCarbonVdw;
- case 'PRO': case 'LYS': case 'ARG': case 'MET': case 'ILE': case 'LEU':
- return tetrahedralCarbonVdw;
- case 'GLU': case 'GLN':
- return atomId === 'CD' ? trigonalCarbonVdw : tetrahedralCarbonVdw;
- }
- }
- }
- return VdwRadius(element);
-}
-
-function determineRadiusNucl(atomId: string, element: ElementSymbol, compId: string): number {
- switch (element) {
- case 'C':
- return nucCarbonVdw;
- case 'N':
- return nucNitrogenVdw;
- case 'P':
- return nucPhosphorusVdw;
- case 'O':
- return oxygenVdw;
- }
- return VdwRadius(element);
-}
-
-function initialize(unit: Unit.Atomic, params: PD.Values<AccessibleSurfaceAreaComputationParams>): AccessibleSurfaceAreaContext {
- return {
- unit: unit,
- params: params,
- spherePoints: generateSpherePoints(params.numberOfSpherePoints),
- cons: 4.0 * Math.PI / params.numberOfSpherePoints,
- maxLookupRadius: 1.4 + 1.4 + 1.87 + 1.87
- }
-}
-
-/** 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 { computeAccessibleSurfaceArea, missingAccessibleSurfaceAreaValue }
\ No newline at end of file
diff --git a/src/mol-model/structure/structure/unit/accessible-surface-area/data.ts b/src/mol-model/structure/structure/unit/accessible-surface-area/data.ts
deleted file mode 100644
index 3a50e10e0a080c542c14430e4bf7c7fa11789e2c..0000000000000000000000000000000000000000
--- a/src/mol-model/structure/structure/unit/accessible-surface-area/data.ts
+++ /dev/null
@@ -1,35 +0,0 @@
-/**
- * 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 { BitFlags } from 'mol-util';
-
-export interface AccessibleSurfaceArea {
- readonly atomRadius: ArrayLike<number>,
- readonly accessibleSurfaceArea: ArrayLike<number>,
- readonly relativeAccessibleSurfaceArea: ArrayLike<number>,
- readonly buried: Uint8Array
-}
-
-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 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
- }
-}
-
-
-export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams
\ No newline at end of file
diff --git a/src/mol-model/structure/structure/util/lookup3d.ts b/src/mol-model/structure/structure/util/lookup3d.ts
index 2427ba51e2d0bb44a4200fe4b0384a63f8669664..d63fb8721dcb5c57df5cb326139d0a8a9b12bc04 100644
--- a/src/mol-model/structure/structure/util/lookup3d.ts
+++ b/src/mol-model/structure/structure/util/lookup3d.ts
@@ -5,12 +5,30 @@
*/
import Structure from '../structure'
-import { Lookup3D, GridLookup3D, Result, Box3D, Sphere3D } from 'mol-math/geometry';
+import { Lookup3D, GridLookup3D, Box3D, Sphere3D, Result } from 'mol-math/geometry';
import { Vec3 } from 'mol-math/linear-algebra';
import { computeStructureBoundary } from './boundary';
import { OrderedSet } from 'mol-data/int';
import { StructureUniqueSubsetBuilder } from './unique-subset-builder';
import StructureElement from '../element';
+import Unit from '../unit';
+
+export interface StructureResult extends Result<number> {
+ units: Unit[]
+}
+
+export namespace StructureResult {
+ export function add(result: StructureResult, unit: Unit, index: number, distSq: number) {
+ result.indices[result.count] = index;
+ result.units[result.count] = unit;
+ result.squaredDistances[result.count] = distSq;
+ result.count++;
+ }
+
+ export function create(): StructureResult {
+ return { count: 0, indices: [], units: [], squaredDistances: [] };
+ }
+}
export class StructureLookup3D {
private unitLookup: Lookup3D;
@@ -20,28 +38,27 @@ export class StructureLookup3D {
return this.unitLookup.find(x, y, z, radius);
}
- // TODO: find another efficient way how to implement this instead of using "tuple".
- // find(x: number, y: number, z: number, radius: number): Result<Element.Packed> {
- // Result.reset(this.result);
- // const { units } = this.structure;
- // const closeUnits = this.unitLookup.find(x, y, z, radius);
- // if (closeUnits.count === 0) return this.result;
-
- // for (let t = 0, _t = closeUnits.count; t < _t; t++) {
- // const unit = units[closeUnits.indices[t]];
- // Vec3.set(this.pivot, x, y, z);
- // if (!unit.conformation.operator.isIdentity) {
- // Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse);
- // }
- // const unitLookup = unit.lookup3d;
- // const groupResult = unitLookup.find(this.pivot[0], this.pivot[1], this.pivot[2], radius);
- // for (let j = 0, _j = groupResult.count; j < _j; j++) {
- // Result.add(this.result, Element.Packed.create(unit.id, groupResult.indices[j]), groupResult.squaredDistances[j]);
- // }
- // }
-
- // return this.result;
- // }
+ private result: StructureResult = StructureResult.create();
+ find(x: number, y: number, z: number, radius: number): StructureResult {
+ Result.reset(this.result);
+ const { units } = this.structure;
+ const closeUnits = this.unitLookup.find(x, y, z, radius);
+ if (closeUnits.count === 0) return this.result;
+
+ for (let t = 0, _t = closeUnits.count; t < _t; t++) {
+ const unit = units[closeUnits.indices[t]];
+ Vec3.set(this.pivot, x, y, z);
+ if (!unit.conformation.operator.isIdentity) {
+ Vec3.transformMat4(this.pivot, this.pivot, unit.conformation.operator.inverse);
+ }
+ const unitLookup = unit.lookup3d;
+ const groupResult = unitLookup.find(this.pivot[0], this.pivot[1], this.pivot[2], radius);
+ for (let j = 0, _j = groupResult.count; j < _j; j++) {
+ StructureResult.add(this.result, unit, groupResult.indices[j], groupResult.squaredDistances[j]);
+ }
+ }
+ return this.result;
+ }
findIntoBuilder(x: number, y: number, z: number, radius: number, builder: StructureUniqueSubsetBuilder) {
const { units } = this.structure;
@@ -98,8 +115,6 @@ export class StructureLookup3D {
}
}
-
-
check(x: number, y: number, z: number, radius: number): boolean {
const { units } = this.structure;
const closeUnits = this.unitLookup.find(x, y, z, radius);
diff --git a/src/mol-theme/color.ts b/src/mol-theme/color.ts
index 48eb75cc2eb2c2b5ea76e999d044cfc283c0f66e..cab7d59fcd15c3839d33722eed8d6d929908f7ec 100644
--- a/src/mol-theme/color.ts
+++ b/src/mol-theme/color.ts
@@ -7,7 +7,6 @@
import { Color } from 'mol-util/color';
import { Location } from 'mol-model/location';
import { ColorType } from 'mol-geo/geometry/color-data';
-import { AccessibleSurfaceAreaColorThemeProvider } from './color/accessible-surface-area'
import { CarbohydrateSymbolColorThemeProvider } from './color/carbohydrate-symbol';
import { UniformColorThemeProvider } from './color/uniform';
import { deepEqual } from 'mol-util';
@@ -63,7 +62,6 @@ namespace ColorTheme {
}
export const BuiltInColorThemes = {
- 'accessible-surface-area': AccessibleSurfaceAreaColorThemeProvider,
'carbohydrate-symbol': CarbohydrateSymbolColorThemeProvider,
'chain-id': ChainIdColorThemeProvider,
'cross-link': CrossLinkColorThemeProvider,
diff --git a/src/mol-theme/color/accessible-surface-area.ts b/src/mol-theme/color/accessible-surface-area.ts
deleted file mode 100644
index caa7521d480bf0187bee3e09090cb4385f35a728..0000000000000000000000000000000000000000
--- a/src/mol-theme/color/accessible-surface-area.ts
+++ /dev/null
@@ -1,69 +0,0 @@
-/**
- * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
- */
-
-import { Color } from 'mol-util/color';
-import { Location } from 'mol-model/location';
-import { ColorTheme, LocationColor } from '../color';
-import { ParamDefinition as PD } from 'mol-util/param-definition'
-import { ThemeDataContext } from '../theme';
-import { ColorListOptions, ColorListName, ColorScale } from 'mol-util/color/scale';
-import { StructureElement, Unit } from 'mol-model/structure';
-import { missingAccessibleSurfaceAreaValue } from 'mol-model/structure/structure/unit/accessible-surface-area/compute';
-
-const DefaultColor = Color(0xFFFFFF)
-const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
-
-export const AccessibleSurfaceAreaColorThemeParams = {
- list: PD.ColorScale<ColorListName>('Rainbow', ColorListOptions),
-}
-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
- let scale: ColorScale | undefined = undefined
-
- if (ctx.structure) {
- scale = ColorScale.create({
- listOrName: props.list,
- minLabel: 'Start',
- maxLabel: 'End'
- })
- const scaleColor = scale.color
-
- color = (location: Location): Color => {
- if (StructureElement.isLocation(location)) {
- if (Unit.isAtomic(location.unit)) {
- const value = location.unit.accessibleSurfaceArea.relativeAccessibleSurfaceArea[location.unit.residueIndex[location.element]];
- return value !== missingAccessibleSurfaceAreaValue ? scaleColor(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
-}
\ No newline at end of file
diff --git a/src/tests/browser/render-structure.ts b/src/tests/browser/render-structure.ts
index 983c3b935aa33061b857437147dc9d4e855203e2..7bc1713df05f0c2b532c15a21d2d1ca4b1dec0cc 100644
--- a/src/tests/browser/render-structure.ts
+++ b/src/tests/browser/render-structure.ts
@@ -7,11 +7,17 @@
import './index.html'
import { Canvas3D } from 'mol-canvas3d/canvas3d';
import CIF, { CifFrame } from 'mol-io/reader/cif'
-import { Model, Structure } from 'mol-model/structure';
-import { ColorTheme } from 'mol-theme/color';
+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 { trajectoryFromMmCIF } from 'mol-model-formats/structure/mmcif';
+import { AccessibleSurfaceArea } from 'mol-model/structure/structure/accessible-surface-area';
+import { Color, ColorScale } from 'mol-util/color';
+import { Location } from 'mol-model/location';
+import { ThemeDataContext } from 'mol-theme/theme';
+import { ParamDefinition as PD } from 'mol-util/param-definition';
+import { ColorListName, ColorListOptions } from 'mol-util/color/scale';
const parent = document.getElementById('app')!
parent.style.width = '100%'
@@ -60,23 +66,78 @@ function getCartoonRepr() {
return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams)
}
-async function init() {
- const cif = await downloadFromPdb(/*'3j3q'*/'1hrc')
+let accessibleSurfaceArea: AccessibleSurfaceArea;
+async function init(props = {}) {
+ const cif = await downloadFromPdb(
+ // '3j3q'
+ '1aon'
+ // '1acj'
+ )
const models = await getModels(cif)
-
const structure = await getStructure(models[0])
+
+ // async compute ASA
+ accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure)
+
const cartoonRepr = getCartoonRepr()
- console.time('ASA');
+ // create color theme
cartoonRepr.setTheme({
- color: reprCtx.colorThemeRegistry.create('accessible-surface-area', { structure }),
+ color: AccessibleSurfaceAreaColorTheme(reprCtx, { ...PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams), ...props }),
size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
})
await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
- console.timeEnd('ASA');
canvas3d.add(cartoonRepr)
canvas3d.resetCamera()
}
-init()
\ No newline at end of file
+init()
+
+const DefaultColor = Color(0xFFFFFF)
+const Description = 'Assigns a color based on the relative accessible surface area of a residue.'
+
+export const AccessibleSurfaceAreaColorThemeParams = {
+ list: PD.ColorScale<ColorListName>('Rainbow', ColorListOptions)
+}
+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 = () => DefaultColor
+ const scale = ColorScale.create({
+ listOrName: props.list,
+ minLabel: '0.0 (buried)',
+ maxLabel: '1.0 (exposed)',
+ domain: [0.0, 1.0]
+ })
+ color = (location: Location): Color => {
+ if (StructureElement.isLocation(location)) {
+ if (Unit.isAtomic(location.unit)) {
+ const value = accessibleSurfaceArea.relativeAccessibleSurfaceArea![location.unit.residueIndex[location.element]];
+ return value !== AccessibleSurfaceArea.missingValue ? scale.color(value) : DefaultColor;
+ }
+ }
+
+ return 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
+}
\ No newline at end of file