diff --git a/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts b/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts
index 6ab4c79a204b6956f6be4facf1aad27d70f6ec0b..2d648a90b733750dcecba57ee43fc621846cc119 100644
--- a/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts
+++ b/src/mol-model/structure/model/properties/utils/accessible-surface-area.ts
@@ -1,39 +1,117 @@
-import { ResidueIndex } from 'mol-model/structure';
-import { MoleculeType } from 'mol-model/structure/model/types';
-import { GridLookup3D } from 'mol-math/geometry';
import { Vec3 } from 'mol-math/linear-algebra';
-import { SortedArray } from 'mol-data/int';
-import { ElementIndex } from 'mol-model/structure/model/indexing';
import { AtomicHierarchy, AtomicConformation } from '../atomic';
+import { ElementSymbol, VdwRadii } from '../../types';
+import { skipUntil } from 'rxjs/operators';
+import { atomicHet } from 'mol-model/structure/query/queries/internal';
+import { compile } from 'mol-script/runtime/query/compiler';
+
+const defaultNumberOfPoints = 960;
+const defaultProbeSize = 1.4;
+const trigonalCarbonVdw = 1.76;
+const tetrahedralCarbonVdw = 1.87;
+const trigonalNitrogenVdw = 1.65;
+const tetrahedralNitrogenVdw = 1.50;
+/** deviating radii from the definition in types.ts */
+const oxygenVdw = 1.4;
+const sulfurVdw = 1.85;
export function computeModelASA(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
- const { lookup3d, proteinResidues } = calcAtomicTraceLookup3D(hierarchy, conformation)
- const backboneIndices = calcBackboneAtomIndices(hierarchy, proteinResidues)
-
- const residueCount = proteinResidues.length
- const numberOfSpherePoints = 960
+ const numberOfSpherePoints = defaultNumberOfPoints;
const ctx: ASAContext = {
- probeSize: 1.4,
+ probeSize: defaultProbeSize,
spherePoints: generateSpherePoints(numberOfSpherePoints),
cons: 4.0 * Math.PI / numberOfSpherePoints,
hierarchy,
- proteinResidues
- }
+ conformation
+ }
- calculateASA(ctx)
+ calculateASA(ctx);
}
function calculateASA(ctx: ASAContext) {
+ console.log(ctx.hierarchy)
+ // 1. extract all heavy atoms
+ // 2. assign radius to all heavy atoms - depends on element
+ // 3. for each residue
+ // 3a. calculate the individual ASA of each atom
+ // 3b. sum up
+ // 3c. (optionally) normalize by maximum value expected for a particular amino acid - needs lookup of max values
+
+ // const { type_symbol } = ctx.hierarchy.atoms;
+ // console.log(type_symbol.value(0));
+
+ const { atoms, residues, residueAtomSegments, derived } = ctx.hierarchy;
+ const { type_symbol } = atoms;
+ const radii: number[] = [];
+
+ // TODO can be more elegant by obtaining residue info once and then using offset to navigate to next residue
+ for (let i = 0; i < type_symbol.rowCount; ++i) {
+ // skip hydrogen atoms
+ if (type_symbol.value(i) === 'H' || type_symbol.value(i) === 'D' || type_symbol.value(i) === 'T')
+ continue;
+
+ // determine group this atom belongs to
+ const groupIndex = residueAtomSegments.index[i];
+
+ // skip entities not part of a peptide chain
+ if (derived.residue.moleculeType[groupIndex] !== 4)
+ continue;
+
+ const atomId = atoms.label_atom_id.value(i);
+ const compId = residues.label_comp_id.value(groupIndex);
+ const element = type_symbol.value(i);
+ radii[radii.length] = determineRadius(atomId, element, compId);
+ }
+}
+/**
+ * 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. Thus this method cannot be used in a structure that contains Hydrogens!
+ */
+function determineRadius(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':
+ if (atomId === 'C' || atomId === 'CE1' || atomId === 'CE2' || atomId === 'CE3' ||
+ atomId === 'CH2' || atomId === 'CZ' || atomId === 'CZ2' || atomId === 'CZ3') {
+ return trigonalCarbonVdw;
+ } else if (atomId === 'CA' || atomId === 'CB' || atomId === 'CE' || atomId === 'CG1' ||
+ atomId === 'CG2') {
+ return tetrahedralCarbonVdw;
+ } else if (compId === 'PHE' || compId === 'TRP' || compId === 'TYR' || compId === 'HIS' ||
+ compId === 'ASP' || compId === 'ASN') {
+ return trigonalCarbonVdw;
+ } else if (compId === 'PRO' || compId === 'LYS' || compId === 'ARG' || compId === 'MET' ||
+ compId === 'ILE' || compId === 'LEU') {
+ return tetrahedralCarbonVdw;
+ } else if (compId === 'GLU' || compId === 'GLN') {
+ return atomId === 'CD' ? trigonalCarbonVdw : tetrahedralCarbonVdw;
+ }
+ default:
+ return (VdwRadii as any)[atomId];
+ }
}
+/** 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
-
- return points
+ 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), y, Math.sin(phi) * r] as Vec3;
+ }
+ return points;
}
interface ASAContext {
@@ -42,54 +120,5 @@ interface ASAContext {
cons: number,
hierarchy: AtomicHierarchy,
- proteinResidues: SortedArray<ResidueIndex>
-}
-
-interface BackboneAtomIndices {
- cIndices: ArrayLike<ElementIndex | -1>
- hIndices: ArrayLike<ElementIndex | -1>
- oIndices: ArrayLike<ElementIndex | -1>
- nIndices: ArrayLike<ElementIndex | -1>
-}
-
-function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
- const { x, y, z } = conformation;
- const { moleculeType, traceElementIndex } = hierarchy.derived.residue
- const indices: number[] = []
- const _proteinResidues: number[] = []
- for (let i = 0, il = moleculeType.length; i < il; ++i) {
- if (moleculeType[i] === MoleculeType.protein) {
- indices[indices.length] = traceElementIndex[i]
- _proteinResidues[_proteinResidues.length] = i
- }
- }
- const lookup3d = GridLookup3D({ x, y, z, indices: SortedArray.ofSortedArray(indices) }, 4);
- const proteinResidues = SortedArray.ofSortedArray<ResidueIndex>(_proteinResidues)
- return { lookup3d, proteinResidues }
-}
-
-
-function calcBackboneAtomIndices(hierarchy: AtomicHierarchy, proteinResidues: SortedArray<ResidueIndex>): BackboneAtomIndices {
- const residueCount = proteinResidues.length
- const { index } = hierarchy
-
- const c = new Int32Array(residueCount)
- const h = new Int32Array(residueCount)
- const o = new Int32Array(residueCount)
- const n = new Int32Array(residueCount)
-
- for (let i = 0, il = residueCount; i < il; ++i) {
- const rI = proteinResidues[i]
- c[i] = index.findAtomOnResidue(rI, 'C')
- h[i] = index.findAtomOnResidue(rI, 'H')
- o[i] = index.findAtomOnResidue(rI, 'O')
- n[i] = index.findAtomOnResidue(rI, 'N')
- }
-
- return {
- cIndices: c as unknown as ArrayLike<ElementIndex | -1>,
- hIndices: h as unknown as ArrayLike<ElementIndex | -1>,
- oIndices: o as unknown as ArrayLike<ElementIndex | -1>,
- nIndices: n as unknown as ArrayLike<ElementIndex | -1>,
- }
+ conformation: AtomicConformation
}
\ No newline at end of file