diff --git a/src/mol-model-props/computed/accessible-surface-area.ts b/src/mol-model-props/computed/accessible-surface-area.ts index 9a184af4d9392c0e764c43aac4967993d4e5b50e..b5d228eca175c9c2ef30b1389591594be8c12393 100644 --- a/src/mol-model-props/computed/accessible-surface-area.ts +++ b/src/mol-model-props/computed/accessible-surface-area.ts @@ -7,7 +7,7 @@ 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 { Task, RuntimeContext } from '../../mol-task'; import { idFactory } from '../../mol-util/id-factory'; @@ -27,9 +27,13 @@ export namespace ComputedAccessibleSurfaceArea { } 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) + return (structure: Structure) => attachTask(structure, params) + } + + export function attachTask(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) { + return Task.create('Compute Accessible Surface Area', async ctx => { + if (get(structure)) return; + return await attachFromCifOrCompute(ctx, structure, params) }); } @@ -39,14 +43,13 @@ export namespace ComputedAccessibleSurfaceArea { // TODO `cifExport` and `symbol` }); - export async function attachFromCifOrCompute(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) { - if (structure.customPropertyDescriptors.has(Descriptor)) return true; + export async function attachFromCifOrCompute(ctx: RuntimeContext, structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps> = {}) { + if (structure.customPropertyDescriptors.has(Descriptor)) return; - const compAccessibleSurfaceArea = await computeAccessibleSurfaceArea(structure, params) + const compAccessibleSurfaceArea = await computeAccessibleSurfaceArea(ctx, structure, params) structure.customPropertyDescriptors.add(Descriptor); set(structure, compAccessibleSurfaceArea); - return true; } } @@ -56,9 +59,9 @@ export const AccessibleSurfaceAreaComputationParams = { export type AccessibleSurfaceAreaComputationParams = typeof AccessibleSurfaceAreaComputationParams export type AccessibleSurfaceAreaComputationProps = PD.Values<AccessibleSurfaceAreaComputationParams> -async function computeAccessibleSurfaceArea(structure: Structure, params: Partial<AccessibleSurfaceAreaComputationProps>): Promise<ComputedAccessibleSurfaceArea.Property> { +async function computeAccessibleSurfaceArea(ctx: RuntimeContext, 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 } + const asa = await AccessibleSurfaceArea.compute(structure, p).runInContext(ctx); + return { id: nextAccessibleSurfaceAreaId(), asa } } \ No newline at end of file diff --git a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts index b9c407791d761b35a5e12fb074b12a98da4cf249..c6742a3480b6f5f923c3b6567387e4ce99c7ac7a 100644 --- a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts +++ b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley.ts @@ -2,68 +2,69 @@ * 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 { Task, RuntimeContext } from '../../../mol-task'; -import { BitFlags } from '../../../mol-util'; +// 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'; +import { ShrakeRupleyContext, VdWLookup, MaxAsa, DefaultMaxAsa } from './shrake-rupley/common'; +import { computeArea } from './shrake-rupley/area'; 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.' }) + numberOfSpherePoints: PD.Numeric(92, { min: 12, max: 120, step: 1 }, { 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, { min: 0.1, max: 4, step: 0.01 }, { 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 as occluders.' }) } export type ShrakeRupleyComputationParams = typeof ShrakeRupleyComputationParams +export type ShrakeRupleyComputationProps = PD.Values<ShrakeRupleyComputationParams> + +// TODO +// - add back buried and relative asa 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(); + export function compute(structure: Structure, props: Partial<ShrakeRupleyComputationProps> = {}) { + const p = { ...PD.getDefaultValues(ShrakeRupleyComputationParams), ...props }; + return Task.create('Compute Accessible Surface Area', async runtime => { + return await calculate(runtime, structure, p); + }); } - async function _compute(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>> = {}): Promise<AccessibleSurfaceArea> { - const ctx = initialize(rtctx, structure, params); + async function calculate(runtime: RuntimeContext, structure: Structure, props: ShrakeRupleyComputationProps): Promise<AccessibleSurfaceArea> { + const ctx = initialize(structure, props); assignRadiusForHeavyAtoms(ctx); - computePerResidue(ctx); - normalizeAccessibleSurfaceArea(ctx); + await computeArea(runtime, ctx); + const { accessibleSurfaceArea, serialResidueIndex } = ctx return { - accessibleSurfaceArea: ctx.accessibleSurfaceArea, - relativeAccessibleSurfaceArea: ctx.relativeAccessibleSurfaceArea, - buried: (index: number) => ctx.relativeAccessibleSurfaceArea[index] < 0.16 + serialResidueIndex, + accessibleSurfaceArea }; } - function initialize(rtctx: RuntimeContext, structure: Structure, params: Partial<PD.Values<ShrakeRupleyComputationParams>>): ShrakeRupleyContext { - const { elementCount, polymerResidueCount } = structure; + function initialize(structure: Structure, props: ShrakeRupleyComputationProps): ShrakeRupleyContext { + const { elementCount, atomicResidueCount } = structure; + const { probeSize, nonPolymer, numberOfSpherePoints } = props 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 + structure, + probeSize, + nonPolymer, + spherePoints: generateSpherePoints(numberOfSpherePoints!), + scalingConstant: 4.0 * Math.PI / numberOfSpherePoints!, + maxLookupRadius: 2 * props.probeSize + 2 * VdWLookup[2], // 2x probe size + 2x largest VdW + atomRadiusType: new Int8Array(elementCount), + serialResidueIndex: new Int32Array(elementCount), + accessibleSurfaceArea: new Float32Array(atomicResidueCount) } } @@ -76,26 +77,31 @@ namespace AccessibleSurfaceArea { 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; + points[points.length] = Vec3.create(Math.cos(phi) * r, y, Math.sin(phi) * r); } 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 - } + // 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 + // } + // } + + /** Get relative area for a given component id */ + export function normalize(compId: string, asa: number) { + const maxAsa = MaxAsa[compId] || DefaultMaxAsa; + return asa / maxAsa } } interface AccessibleSurfaceArea { + readonly serialResidueIndex: 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-props/computed/accessible-surface-area/shrake-rupley/area.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/area.ts new file mode 100644 index 0000000000000000000000000000000000000000..51c7159b326e87750b286d0021fb1936fbb0584d --- /dev/null +++ b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/area.ts @@ -0,0 +1,103 @@ +/** + * 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 { ShrakeRupleyContext, VdWLookup } from './common'; +import { Vec3 } from '../../../../mol-math/linear-algebra'; +import { RuntimeContext } from '../../../../mol-task'; +import { StructureProperties, StructureElement, Structure } from '../../../../mol-model/structure/structure'; + +// TODO +// - iterate over units and elements +// - avoid using serial-element index whenever possible +// - calculate atomRadiusType only for invariant units +// - factor serialResidueIndex out + +const updateChunk = 5000 +export async function computeArea(runtime: RuntimeContext, ctx: ShrakeRupleyContext) { + const { atomRadiusType: atomRadius } = ctx; + for (let i = 0; i < atomRadius.length; i += updateChunk) { + if (runtime.shouldUpdate) { + await runtime.update({ message: 'Computing per residue surface accessibility...', current: i, max: atomRadius.length }); + } + + computeRange(ctx, i, Math.min(i + updateChunk, atomRadius.length)); + } +} + +const aPos = Vec3(); +const bPos = Vec3(); +const testPoint = Vec3(); +const aLoc = StructureElement.Location.create() +const bLoc = StructureElement.Location.create() + +function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) { + l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]] + l.element = structure.serialMapping.elementIndices[serialIndex] + return l +} + +function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) { + const { structure, atomRadiusType, serialResidueIndex, accessibleSurfaceArea, spherePoints, scalingConstant, maxLookupRadius, probeSize } = ctx; + const { x, y, z } = StructureProperties.atom + const { lookup3d, serialMapping, unitIndexMap } = structure; + const { cumulativeUnitElementCount } = serialMapping + + for (let aI = begin; aI < end; ++aI) { + const radius1 = VdWLookup[atomRadiusType[aI]]; + if (radius1 === VdWLookup[0]) continue; + + setLocation(aLoc, structure, aI) + const aX = x(aLoc) + const aY = y(aLoc) + const aZ = z(aLoc) + + // pre-filter by lookup3d (provides >10x speed-up compared to naive evaluation) + const { count, units, indices, squaredDistances } = lookup3d.find(aX, aY, aZ, maxLookupRadius); + + // collect neighbors for each atom + const cutoff1 = probeSize + probeSize + radius1; + const neighbors = []; // TODO reuse + for (let iI = 0; iI < count; ++iI) { + const bUnit = units[iI] + StructureElement.Location.set(bLoc, bUnit, bUnit.elements[indices[iI]]) + const bI = cumulativeUnitElementCount[unitIndexMap.get(bUnit.id)] + indices[iI] + + const radius2 = VdWLookup[atomRadiusType[bI]]; + if (StructureElement.Location.areEqual(aLoc, bLoc) || 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 + Vec3.set(aPos, aX, aY, aZ) + const scale = probeSize + radius1; + let accessiblePointCount = 0; + for (let sI = 0; sI < spherePoints.length; ++sI) { + Vec3.scaleAndAdd(testPoint, aPos, spherePoints[sI], scale) + let accessible = true; + + for (let _nI = 0; _nI < neighbors.length; ++_nI) { + const nI = neighbors[_nI]; + setLocation(bLoc, structure, nI) + Vec3.set(bPos, x(bLoc), y(bLoc), z(bLoc)) + const radius3 = VdWLookup[atomRadiusType[nI]]; + const cutoff3 = (radius3 + probeSize) * (radius3 + probeSize); + if (Vec3.squaredDistance(testPoint, bPos) < cutoff3) { + accessible = false; + break; + } + } + + if (accessible) ++accessiblePointCount; + } + + accessibleSurfaceArea[serialResidueIndex[aI]] += scalingConstant * accessiblePointCount * scale * scale; + } +} \ No newline at end of file diff --git a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts index e2360a917a9112eb88a12b96be0f31c1eb7cba38..c1f78e42b6d0cf0dcb3ee282567cdde1187aab3c 100644 --- a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts +++ b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/common.ts @@ -4,25 +4,22 @@ * @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, + scalingConstant: number, maxLookupRadius: number, - atomRadius: Int8Array, - accessibleSurfaceArea: Float32Array, - relativeAccessibleSurfaceArea: Float32Array, - updateChunk: number + atomRadiusType: Int8Array, + serialResidueIndex: Int32Array, + accessibleSurfaceArea: Float32Array } -// Chothia's amino acid and nucleotide atom vdw radii +/** Chothia's amino acid and nucleotide atom vdw radii */ export const VdWLookup = [ -1.0, // 0: missing 1.76, // 1: trigonal C @@ -36,9 +33,8 @@ export const VdWLookup = [ 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 = { +export const MaxAsa: { [k: string]: number } = { 'ALA': 121.0, 'ARG': 265.0, 'ASN': 187.0, diff --git a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/normalize.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/normalize.ts deleted file mode 100644 index 9331942c6e9b1da0415d63d901e82d7dd8c18c93..0000000000000000000000000000000000000000 --- a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/normalize.ts +++ /dev/null @@ -1,32 +0,0 @@ -/** - * 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; - } -} \ No newline at end of file diff --git a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/per-residue.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/per-residue.ts deleted file mode 100644 index 5cf13e6aa60841f6cee19fe4670e1db74f90ee98..0000000000000000000000000000000000000000 --- a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/per-residue.ts +++ /dev/null @@ -1,86 +0,0 @@ -/** - * 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; - } -} \ No newline at end of file diff --git a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts index 421c34f0fc138eda4fedbe11278fcb9a3173ac39..6ad39b5f9dbddbc1289eaf9ba2edba66a2326677 100644 --- a/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts +++ b/src/mol-model-props/computed/accessible-surface-area/shrake-rupley/radii.ts @@ -2,61 +2,75 @@ * 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 { 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'; +import { StructureElement, StructureProperties } from '../../../../mol-model/structure/structure'; +import { getElementMoleculeType } from '../../../../mol-model/structure/util'; -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)); - } -} +const l = StructureElement.Location.create() -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; - } +export function assignRadiusForHeavyAtoms(ctx: ShrakeRupleyContext) { + const { label_comp_id, key } = StructureProperties.residue + const { type_symbol, label_atom_id } = StructureProperties.atom + const { structure, atomRadiusType, serialResidueIndex } = ctx; + + let prevResidueIdx = 0 + let residueIdx = 0 + let serialResidueIdx = -1 + + for (let i = 0, m = 0, il = structure.units.length; i < il; ++i) { + const unit = structure.units[i] + const { elements } = unit + l.unit = unit + + prevResidueIdx = -1 + + for (let j = 0, jl = elements.length; j < jl; ++j) { + const eI = elements[j] + const mj = m + j + + l.element = eI + residueIdx = key(l) + + if (prevResidueIdx !== residueIdx) ++serialResidueIdx + prevResidueIdx = residueIdx - const residueType = moleculeType[rI]; - // skip non-polymer groups - if (!ctx.nonPolymer) { - if (!isPolymer(residueType)) { - ctx.atomRadius[aI] = VdWLookup[0]; + const element = type_symbol(l); + const elementIdx = getElementIdx(element); + + // skip hydrogen atoms + if (isHydrogen(elementIdx)) { + atomRadiusType[mj] = VdWLookup[0]; + serialResidueIndex[mj] = -1 continue; } - } - const atomId = label_atom_id.value(aI); - let compId = label_comp_id.value(rI); + const residueType = getElementMoleculeType(unit, eI) + // skip non-polymer groups + if (!ctx.nonPolymer && !isPolymer(residueType)) { + atomRadiusType[mj] = VdWLookup[0]; + serialResidueIndex[mj] = -1 + continue; + } - // handle modified residues - const parentId = model.properties.modifiedResidues.parentId.get(compId); - if (parentId !== void 0) compId = parentId; + const atomId = label_atom_id(l); + const compId = label_comp_id(l); - 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); + if (isNucleic(residueType)) { + atomRadiusType[mj] = determineRadiusNucl(atomId, element, compId); + } else if (residueType === MoleculeType.Protein) { + atomRadiusType[mj] = determineRadiusAmino(atomId, element, compId); + } else { + atomRadiusType[mj] = handleNonStandardCase(element); + } + serialResidueIndex[mj] = serialResidueIdx } + m += elements.length } } @@ -68,41 +82,37 @@ function computeRange(ctx: ShrakeRupleyContext, begin: number, end: number) { function determineRadiusAmino(atomId: string, element: ElementSymbol, compId: string): number { switch (element) { case 'O': - return 5; + return 5; case 'S': - return 6; + return 6; case 'N': - return atomId === 'NZ' ? 4 : 3; + 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; + 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; + case 'C': return 7; + case 'N': return 8; + case 'P': return 9; + case 'O': return 5; } return handleNonStandardCase(element); } diff --git a/src/mol-theme/color.ts b/src/mol-theme/color.ts index 9288607d6d919b675f76a7ddd2b686f63c5bbb55..11f037f15b064c83482991d61a9738c549b9fca2 100644 --- a/src/mol-theme/color.ts +++ b/src/mol-theme/color.ts @@ -33,6 +33,7 @@ import { ModelIndexColorThemeProvider } from './color/model-index'; import { OccupancyColorThemeProvider } from './color/occupancy'; import { OperatorNameColorThemeProvider } from './color/operator-name'; import { OperatorHklColorThemeProvider } from './color/operator-hkl'; +import { AccessibleSurfaceAreaColorThemeProvider } from './color/accessible-surface-area'; export type LocationColor = (location: Location, isSecondary: boolean) => Color @@ -74,6 +75,7 @@ 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 index 2cf09cc62030703a5c9dfe7c8f1ad03232234667..dea5769d384f1d7bcfd558e7a36823defc5b61ed 100644 --- a/src/mol-theme/color/accessible-surface-area.ts +++ b/src/mol-theme/color/accessible-surface-area.ts @@ -15,7 +15,7 @@ import { ComputedAccessibleSurfaceArea } from '../../mol-model-props/computed/ac 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 DefaultColor = Color(0xFAFAFA) const Description = 'Assigns a color based on the relative accessible surface area of a residue.' export const AccessibleSurfaceAreaColorThemeParams = { @@ -35,17 +35,20 @@ export function AccessibleSurfaceAreaColorTheme(ctx: ThemeDataContext, props: PD domain: [0.0, 1.0] }) - const accessibleSurfaceArea = ctx.structure ? ComputedAccessibleSurfaceArea.get(ctx.structure)!.asa : undefined + const accessibleSurfaceArea = ctx.structure ? ComputedAccessibleSurfaceArea.get(ctx.structure.root)?.asa : undefined + + if (accessibleSurfaceArea && ctx.structure) { + const { structure } = ctx + const { getSerialIndex } = structure.root.serialMapping + const { relativeAccessibleSurfaceArea, serialResidueIndex } = accessibleSurfaceArea - 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; + const rSI = serialResidueIndex[getSerialIndex(location.unit, location.element)] + return rSI === -1 ? DefaultColor : scale.color(relativeAccessibleSurfaceArea[rSI]) } } - return DefaultColor } } else { @@ -67,5 +70,5 @@ export const AccessibleSurfaceAreaColorThemeProvider: ColorTheme.Provider<Access factory: AccessibleSurfaceAreaColorTheme, getParams: getAccessibleSurfaceAreaColorThemeParams, defaultValues: PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams), - isApplicable: (ctx: ThemeDataContext) => !!ctx.structure + isApplicable: (ctx: ThemeDataContext) => !!ctx.structure, } \ No newline at end of file diff --git a/src/tests/browser/render-asa.ts b/src/tests/browser/render-asa.ts deleted file mode 100644 index 7f5c8d1348ee4956e51e3a089cf89b8b4284ecfc..0000000000000000000000000000000000000000 --- a/src/tests/browser/render-asa.ts +++ /dev/null @@ -1,150 +0,0 @@ -/** - * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info. - * - * @author Alexander Rose <alexander.rose@weirdbyte.de> - */ - -import './index.html' -import { Canvas3D } from '../../mol-canvas3d/canvas3d'; -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 { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif'; -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 { 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%' -parent.style.height = '100%' - -const canvas = document.createElement('canvas') -parent.appendChild(canvas) -resizeCanvas(canvas, parent) - -const canvas3d = Canvas3D.fromCanvas(canvas) -canvas3d.animate() - - -async function parseCif(data: string|Uint8Array) { - const comp = CIF.parse(data); - const parsed = await comp.run(); - if (parsed.isError) throw parsed; - return parsed.result; -} - -async function downloadCif(url: string, isBinary: boolean) { - const data = await fetch(url); - return parseCif(isBinary ? new Uint8Array(await data.arrayBuffer()) : await data.text()); -} - -async function downloadFromPdb(pdb: string) { - // const parsed = await downloadCif(`https://files.rcsb.org/download/${pdb}.cif`, false); - const parsed = await downloadCif(`https://webchem.ncbr.muni.cz/ModelServer/static/bcif/${pdb}`, true); - return parsed.blocks[0]; -} - -async function getModels(frame: CifFrame) { - return await trajectoryFromMmCIF(frame).run(); -} - -async function getStructure(model: Model) { - return Structure.ofModel(model); -} - -const reprCtx = { - colorThemeRegistry: ColorTheme.createRegistry(), - sizeThemeRegistry: SizeTheme.createRegistry() -} -function getCartoonRepr() { - return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams) -} - -let accessibleSurfaceArea: AccessibleSurfaceArea; -async function init(props = {}) { - const cif = await downloadFromPdb( - // '3j3q' - // '1aon' - // '1acj' - // '1pga' - // '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() - - // create color theme - cartoonRepr.setTheme({ - color: AccessibleSurfaceAreaColorTheme(reprCtx, { ...PD.getDefaultValues(AccessibleSurfaceAreaColorThemeParams), ...props }), - size: reprCtx.sizeThemeRegistry.create('uniform', { structure }) - }) - await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run() - - canvas3d.add(cartoonRepr) - canvas3d.resetCamera() -} - -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.ColorList<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.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 - } - - 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 e85346bd90ffa49f271d60b86870d8dca20d3708..98632439c9262d54f7cde91e57869dc921163194 100644 --- a/src/tests/browser/render-structure.ts +++ b/src/tests/browser/render-structure.ts @@ -15,8 +15,15 @@ import { trajectoryFromMmCIF } from '../../mol-model-formats/structure/mmcif'; import { MolecularSurfaceRepresentationProvider } from '../../mol-repr/structure/representation/molecular-surface'; import { BallAndStickRepresentationProvider } from '../../mol-repr/structure/representation/ball-and-stick'; import { GaussianSurfaceRepresentationProvider } from '../../mol-repr/structure/representation/gaussian-surface'; -import { ComputedSecondaryStructure } from '../../mol-model-props/computed/secondary-structure'; +// import { ComputedSecondaryStructure } from '../../mol-model-props/computed/secondary-structure'; import { resizeCanvas } from '../../mol-canvas3d/util'; +import { Representation } from '../../mol-repr/representation'; +import { throttleTime } from 'rxjs/operators'; +import { MarkerAction } from '../../mol-util/marker-action'; +import { EveryLoci } from '../../mol-model/loci'; +import { lociLabel } from '../../mol-theme/label'; +import { InteractionsRepresentationProvider } from '../../mol-repr/structure/representation/interactions'; +import { ComputedInteractions } from '../../mol-model-props/computed/interactions'; const parent = document.getElementById('app')! parent.style.width = '100%' @@ -29,6 +36,33 @@ resizeCanvas(canvas, parent) const canvas3d = Canvas3D.fromCanvas(canvas) canvas3d.animate() +const info = document.createElement('div') +info.style.position = 'absolute' +info.style.fontFamily = 'sans-serif' +info.style.fontSize = '16pt' +info.style.bottom = '20px' +info.style.right = '20px' +info.style.color = 'white' +parent.appendChild(info) + +let prevReprLoci = Representation.Loci.Empty +canvas3d.input.move.pipe(throttleTime(100)).subscribe(({x, y}) => { + const pickingId = canvas3d.identify(x, y) + let label = '' + if (pickingId) { + const reprLoci = canvas3d.getLoci(pickingId) + label = lociLabel(reprLoci.loci) + if (!Representation.Loci.areEqual(prevReprLoci, reprLoci)) { + canvas3d.mark(prevReprLoci, MarkerAction.RemoveHighlight) + canvas3d.mark(reprLoci, MarkerAction.Highlight) + prevReprLoci = reprLoci + } + } else { + canvas3d.mark({ loci: EveryLoci }, MarkerAction.RemoveHighlight) + prevReprLoci = Representation.Loci.Empty + } + info.innerHTML = label +}) async function parseCif(data: string|Uint8Array) { const comp = CIF.parse(data); @@ -43,8 +77,7 @@ async function downloadCif(url: string, isBinary: boolean) { } async function downloadFromPdb(pdb: string) { - // const parsed = await downloadCif(`https://files.rcsb.org/download/${pdb}.cif`, false); - const parsed = await downloadCif(`https://webchem.ncbr.muni.cz/ModelServer/static/bcif/${pdb}`, true); + const parsed = await downloadCif(`https://models.rcsb.org/${pdb}.bcif`, true); return parsed.blocks[0]; } @@ -57,6 +90,7 @@ async function getStructure(model: Model) { } const reprCtx = { + webgl: canvas3d.webgl, colorThemeRegistry: ColorTheme.createRegistry(), sizeThemeRegistry: SizeTheme.createRegistry() } @@ -64,6 +98,10 @@ function getCartoonRepr() { return CartoonRepresentationProvider.factory(reprCtx, CartoonRepresentationProvider.getParams) } +function getInteractionRepr() { + return InteractionsRepresentationProvider.factory(reprCtx, InteractionsRepresentationProvider.getParams) +} + function getBallAndStickRepr() { return BallAndStickRepresentationProvider.factory(reprCtx, BallAndStickRepresentationProvider.getParams) } @@ -80,34 +118,54 @@ async function init() { const cif = await downloadFromPdb('1crn') const models = await getModels(cif) const structure = await getStructure(models[0]) - console.time('computeDSSP') - await ComputedSecondaryStructure.attachFromCifOrCompute(structure) - console.timeEnd('computeDSSP'); + // console.time('computeDSSP') + // await ComputedSecondaryStructure.attachFromCifOrCompute(structure) + // console.timeEnd('computeDSSP'); + + // console.time('computeValenceModel') + // await ComputedValenceModel.attachFromCifOrCompute(structure) + // console.timeEnd('computeValenceModel'); + // console.log(ComputedValenceModel.get(structure)) + + console.time('computeInteractions') + await ComputedInteractions.attachFromCifOrCompute(structure) + console.timeEnd('computeInteractions'); + console.log(ComputedInteractions.get(structure)) const show = { - cartoon: false, + cartoon: true, + interaction: true, ballAndStick: true, - molecularSurface: true, + molecularSurface: false, gaussianSurface: false, } const cartoonRepr = getCartoonRepr() + const interactionRepr = getInteractionRepr() const ballAndStickRepr = getBallAndStickRepr() const molecularSurfaceRepr = getMolecularSurfaceRepr() const gaussianSurfaceRepr = getGaussianSurfaceRepr() if (show.cartoon) { cartoonRepr.setTheme({ - color: reprCtx.colorThemeRegistry.create('secondary-structure', { structure }), + color: reprCtx.colorThemeRegistry.create('element-symbol', { structure }), size: reprCtx.sizeThemeRegistry.create('uniform', { structure }) }) await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run() } + if (show.interaction) { + interactionRepr.setTheme({ + color: reprCtx.colorThemeRegistry.create('interaction-type', { structure }), + size: reprCtx.sizeThemeRegistry.create('uniform', { structure }) + }) + await interactionRepr.createOrUpdate({ ...InteractionsRepresentationProvider.defaultValues, quality: 'auto' }, structure).run() + } + if (show.ballAndStick) { ballAndStickRepr.setTheme({ - color: reprCtx.colorThemeRegistry.create('secondary-structure', { structure }), - size: reprCtx.sizeThemeRegistry.create('uniform', { structure }) + color: reprCtx.colorThemeRegistry.create('element-symbol', { structure }), + size: reprCtx.sizeThemeRegistry.create('uniform', { structure }, { value: 1 }) }) await ballAndStickRepr.createOrUpdate({ ...BallAndStickRepresentationProvider.defaultValues, quality: 'auto' }, structure).run() } @@ -133,6 +191,7 @@ async function init() { } if (show.cartoon) canvas3d.add(cartoonRepr) + if (show.interaction) canvas3d.add(interactionRepr) if (show.ballAndStick) canvas3d.add(ballAndStickRepr) if (show.molecularSurface) canvas3d.add(molecularSurfaceRepr) if (show.gaussianSurface) canvas3d.add(gaussianSurfaceRepr) diff --git a/webpack.config.js b/webpack.config.js index 4f8dbf1867aac1050fbf1c5fb09f3add1c103870..f33a385636c5aff1423777bb51873a7d60824ddb 100644 --- a/webpack.config.js +++ b/webpack.config.js @@ -98,7 +98,6 @@ module.exports = [ createApp('model-server-query'), createBrowserTest('font-atlas'), - createBrowserTest('render-asa'), createBrowserTest('marching-cubes'), createBrowserTest('render-lines'), createBrowserTest('render-mesh'),