diff --git a/src/mol-model-formats/structure/mmcif/parser.ts b/src/mol-model-formats/structure/mmcif/parser.ts index 934de9467bd082f45449c2f4e56f173d1f495d8d..1230ad95424d314524f358ed72e8b110b8c3559c 100644 --- a/src/mol-model-formats/structure/mmcif/parser.ts +++ b/src/mol-model-formats/structure/mmcif/parser.ts @@ -19,7 +19,7 @@ import { createAssemblies } from './assembly'; import { getAtomicHierarchyAndConformation } from './atomic'; import { ComponentBond } from './bonds'; import { getIHMCoarse, EmptyIHMCoarse, IHMData } from './ihm'; -import { getSecondaryStructureMmCif } from './secondary-structure'; +import { getSecondaryStructure } from './secondary-structure'; import { getSequence } from './sequence'; import { sortAtomSite } from './sort'; import { StructConn } from './bonds/struct_conn'; @@ -199,7 +199,7 @@ function createStandardModel(format: mmCIF_Format, atom_site: AtomSite, entities coarseHierarchy: coarse.hierarchy, coarseConformation: coarse.conformation, properties: { - secondaryStructure: getSecondaryStructureMmCif(format.data, atomic.hierarchy), + secondaryStructure: getSecondaryStructure(format.data, atomic.hierarchy, atomic.conformation), ...formatData }, customProperties: new CustomProperties(), @@ -225,7 +225,7 @@ function createModelIHM(format: mmCIF_Format, data: IHMData, formatData: FormatD coarseHierarchy: coarse.hierarchy, coarseConformation: coarse.conformation, properties: { - secondaryStructure: getSecondaryStructureMmCif(format.data, atomic.hierarchy), + secondaryStructure: getSecondaryStructure(format.data, atomic.hierarchy, atomic.conformation), ...formatData }, customProperties: new CustomProperties(), diff --git a/src/mol-model-formats/structure/mmcif/secondary-structure.ts b/src/mol-model-formats/structure/mmcif/secondary-structure.ts index 3c0def78910279781a359199a4a06a0e841dad20..0ba80546f6aa6d92d472d0a541f4dd7aed16671b 100644 --- a/src/mol-model-formats/structure/mmcif/secondary-structure.ts +++ b/src/mol-model-formats/structure/mmcif/secondary-structure.ts @@ -7,10 +7,20 @@ import { mmCIF_Database as mmCIF, mmCIF_Database } from 'mol-io/reader/cif/schema/mmcif' import { SecondaryStructureType } from 'mol-model/structure/model/types'; -import { AtomicHierarchy } from 'mol-model/structure/model/properties/atomic'; +import { AtomicHierarchy, AtomicConformation } from 'mol-model/structure/model/properties/atomic'; import { SecondaryStructure } from 'mol-model/structure/model/properties/seconday-structure'; import { Column } from 'mol-data/db'; import { ChainIndex, ResidueIndex } from 'mol-model/structure/model/indexing'; +import { computeSecondaryStructure } from 'mol-model/structure/model/properties/utils/secondary-structure'; + +// TODO add parameter to allow forcing computation +export function getSecondaryStructure(data: mmCIF_Database, hierarchy: AtomicHierarchy, conformation: AtomicConformation): SecondaryStructure { + if (!data.struct_conf._rowCount && !data.struct_sheet_range._rowCount) { + return computeSecondaryStructure(hierarchy, conformation) + } else { + return getSecondaryStructureMmCif(data, hierarchy) + } +} export function getSecondaryStructureMmCif(data: mmCIF_Database, hierarchy: AtomicHierarchy): SecondaryStructure { const map: SecondaryStructureMap = new Map(); diff --git a/src/mol-model/structure/model/properties/utils/atomic-derived.ts b/src/mol-model/structure/model/properties/utils/atomic-derived.ts index f2e1833042eff059d102d7d25dd1d860eac8bce7..ffdfc911494992f9602bc818f4da7e3beb563940 100644 --- a/src/mol-model/structure/model/properties/utils/atomic-derived.ts +++ b/src/mol-model/structure/model/properties/utils/atomic-derived.ts @@ -45,8 +45,8 @@ export function getAtomicDerivedData(data: AtomicData, index: AtomicIndex, chemi return { residue: { - traceElementIndex: traceElementIndex as unknown as ArrayLike<ElementIndex>, - directionElementIndex: directionElementIndex as unknown as ArrayLike<ElementIndex>, + traceElementIndex: traceElementIndex as unknown as ArrayLike<ElementIndex>, // TODO maybe -1 + directionElementIndex: directionElementIndex as unknown as ArrayLike<ElementIndex>, // TODO maybe -1 moleculeType: moleculeType as unknown as ArrayLike<MoleculeType>, } } diff --git a/src/mol-model/structure/model/properties/utils/secondary-structure.ts b/src/mol-model/structure/model/properties/utils/secondary-structure.ts new file mode 100644 index 0000000000000000000000000000000000000000..2d01032f84d00aed1d95c70edac7eb58e4c7ed2c --- /dev/null +++ b/src/mol-model/structure/model/properties/utils/secondary-structure.ts @@ -0,0 +1,466 @@ +/** + * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author Alexander Rose <alexander.rose@weirdbyte.de> + */ + +import { SecondaryStructure } from 'mol-model/structure/model/properties/seconday-structure'; +import { ResidueIndex } from 'mol-model/structure'; +import { MoleculeType, SecondaryStructureType } from 'mol-model/structure/model/types'; +import { Vec3 } from 'mol-math/linear-algebra'; +import { GridLookup3D } from 'mol-math/geometry'; +import { SortedArray } from 'mol-data/int'; +import { IntAdjacencyGraph } from 'mol-math/graph'; +import { BitFlags } from 'mol-util'; +import { ElementIndex } from 'mol-model/structure/model/indexing'; +import { AtomicHierarchy, AtomicConformation } from '../atomic'; + +export function computeSecondaryStructure(hierarchy: AtomicHierarchy, conformation: AtomicConformation): SecondaryStructure { + // TODO use Zhang-Skolnik for CA alpha only parts or for coarse parts with per-residue elements + return computeModelDSSP(hierarchy, conformation) +} + +export function computeModelDSSP(hierarchy: AtomicHierarchy, conformation: AtomicConformation) { + const { lookup3d, proteinResidues } = calcAtomicTraceLookup3D(hierarchy, conformation) + const backboneIndices = calcBackboneAtomIndices(hierarchy, proteinResidues) + const hbonds = calcBackboneHbonds(hierarchy, conformation, proteinResidues, backboneIndices, lookup3d) + + const residueCount = proteinResidues.length + const flags = new Uint32Array(residueCount) + + const ctx: DSSPContext = { + hierarchy, + proteinResidues, + flags, + hbonds + } + + assignBends(ctx) + assignTurns(ctx) + assignHelices(ctx) + assignBridges(ctx) + assignLadders(ctx) + assignSheets(ctx) + + const assignment = getDSSPAssignment(flags) + + const type = new Uint32Array(hierarchy.residues._rowCount) as unknown as SecondaryStructureType[] + for (let i = 0, il = proteinResidues.length; i < il; ++i) { + type[proteinResidues[i]] = assignment[i] + } + + const secondaryStructure: SecondaryStructure = { + type, + key: [], // TODO + elements: [] // TODO + } + return secondaryStructure +} + +interface DSSPContext { + hierarchy: AtomicHierarchy + proteinResidues: SortedArray<ResidueIndex> + /** flags for each residue */ + flags: Uint32Array + + hbonds: DsspHbonds +} + +interface DSSPType extends BitFlags<DSSPType.Flag> { } +namespace DSSPType { + export const is: (t: DSSPType, f: Flag) => boolean = BitFlags.has + export const create: (f: Flag) => DSSPType = BitFlags.create + export const enum Flag { + _ = 0x0, + H = 0x1, + B = 0x2, + E = 0x4, + G = 0x8, + I = 0x10, + S = 0x20, + T = 0x40, + T3 = 0x80, + T4 = 0x100, + T5 = 0x200, + } +} + +/** max distance between two C-alpha atoms to check for hbond */ +const caMaxDist = 9.0; + +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) }); + const proteinResidues = SortedArray.ofSortedArray<ResidueIndex>(_proteinResidues) + return { lookup3d, proteinResidues } +} + +interface BackboneAtomIndices { + cIndices: ArrayLike<ElementIndex | -1> + hIndices: ArrayLike<ElementIndex | -1> + oIndices: ArrayLike<ElementIndex | -1> + nIndices: ArrayLike<ElementIndex | -1> +} + +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>, + } +} + +type DsspHbonds = IntAdjacencyGraph<{ readonly energies: ArrayLike<number> }> + +function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConformation, proteinResidues: SortedArray<ResidueIndex>, backboneIndices: BackboneAtomIndices, lookup3d: GridLookup3D): DsspHbonds { + const { cIndices, hIndices, nIndices, oIndices } = backboneIndices + const { index } = hierarchy + const { x, y, z } = conformation + const { traceElementIndex } = hierarchy.derived.residue + + const residueCount = proteinResidues.length + const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i]) + + const oAtomResidues: number[] = []; + const nAtomResidues: number[] = []; + const energies: number[] = []; + + const oPos = Vec3.zero() + const cPos = Vec3.zero() + const caPos = Vec3.zero() + const nPos = Vec3.zero() + const hPos = Vec3.zero() + + const cPosPrev = Vec3.zero() + const oPosPrev = Vec3.zero() + + for (let i = 0, il = proteinResidues.length; i < il; ++i) { + const oPI = i + const oRI = proteinResidues[i] + + const oAtom = oIndices[oPI] + const cAtom = cIndices[oPI] + const caAtom = traceElementIndex[oRI] + + // continue if residue is missing O or C atom + if (oAtom === -1 || cAtom === -1) continue + + // ignore C-terminal residue as acceptor + if (index.findAtomOnResidue(oRI, 'OXT') !== -1) continue + + position(oAtom, oPos) + position(cAtom, cPos) + position(caAtom, caPos) + + const { indices, count } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist) + + for (let j = 0; j < count; ++j) { + const nPI = indices[j] + + // ignore bonds within a residue or to prev or next residue, TODO take chain border into account + if (nPI === oPI || nPI - 1 === oPI || nPI + 1 === oPI) continue + + const nAtom = nIndices[nPI] + if (nAtom === -1) continue + + position(nAtom, nPos) + + const hAtom = hIndices[nPI] + if (hAtom === -1) { + // approximate calculation of H position, TODO factor out + if (nPI === 0) continue + const nPIprev = nPI - 1 + + const oAtomPrev = oIndices[nPIprev] + const cAtomPrev = cIndices[nPIprev] + if (oAtomPrev === -1 || cAtomPrev === -1) continue + + position(oAtomPrev, oPosPrev) + position(cAtomPrev, cPosPrev) + + Vec3.sub(hPos, cPosPrev, oPosPrev) + const dist = Vec3.distance(oPosPrev, cPosPrev) + Vec3.scaleAndAdd(hPos, nPos, hPos, 1 / dist) + } else { + position(hAtom, hPos) + } + + const e = calcHbondEnergy(oPos, cPos, nPos, hPos) + if (e > hbondEnergyCutoff) continue + + oAtomResidues[oAtomResidues.length] = oPI + nAtomResidues[nAtomResidues.length] = nPI + energies[energies.length] = e + } + } + + return buildHbondGraph(residueCount, oAtomResidues, nAtomResidues, energies); +} + +function buildHbondGraph(residueCount: number, oAtomResidues: number[], nAtomResidues: number[], energies: number[]) { + const builder = new IntAdjacencyGraph.DirectedEdgeBuilder(residueCount, oAtomResidues, nAtomResidues); + const _energies = new Float32Array(builder.slotCount); + + for (let i = 0, _i = builder.edgeCount; i < _i; i++) { + builder.addNextEdge(); + builder.assignProperty(_energies, energies[i]); + } + + return builder.createGraph({ energies }); +} + +/** Original priority: H,B,E,G,I,T,S */ +function getOriginalResidueFlag(f: DSSPType) { + if (DSSPType.is(f, DSSPType.Flag.H)) return SecondaryStructureType.SecondaryStructureDssp.H + if (DSSPType.is(f, DSSPType.Flag.B)) return SecondaryStructureType.SecondaryStructureDssp.B + if (DSSPType.is(f, DSSPType.Flag.E)) return SecondaryStructureType.SecondaryStructureDssp.E + if (DSSPType.is(f, DSSPType.Flag.G)) return SecondaryStructureType.SecondaryStructureDssp.G + if (DSSPType.is(f, DSSPType.Flag.I)) return SecondaryStructureType.SecondaryStructureDssp.I + if (DSSPType.is(f, DSSPType.Flag.T)) return SecondaryStructureType.SecondaryStructureDssp.T + if (DSSPType.is(f, DSSPType.Flag.S)) return SecondaryStructureType.SecondaryStructureDssp.S + return SecondaryStructureType.Flag.None +} + +/** Version 2.1.0 priority: I,H,B,E,G,T,S */ +function getUpdatedResidueFlag(f: DSSPType) { + if (DSSPType.is(f, DSSPType.Flag.I)) return SecondaryStructureType.SecondaryStructureDssp.I + if (DSSPType.is(f, DSSPType.Flag.H)) return SecondaryStructureType.SecondaryStructureDssp.H + if (DSSPType.is(f, DSSPType.Flag.B)) return SecondaryStructureType.SecondaryStructureDssp.B + if (DSSPType.is(f, DSSPType.Flag.E)) return SecondaryStructureType.SecondaryStructureDssp.E + if (DSSPType.is(f, DSSPType.Flag.G)) return SecondaryStructureType.SecondaryStructureDssp.G + if (DSSPType.is(f, DSSPType.Flag.T)) return SecondaryStructureType.SecondaryStructureDssp.T + if (DSSPType.is(f, DSSPType.Flag.S)) return SecondaryStructureType.SecondaryStructureDssp.S + return SecondaryStructureType.Flag.None +} + +// function geFlagName(f: DSSPType) { +// if (DSSPType.is(f, DSSPType.Flag.I)) return 'I' +// if (DSSPType.is(f, DSSPType.Flag.H)) return 'H' +// if (DSSPType.is(f, DSSPType.Flag.B)) return 'B' +// if (DSSPType.is(f, DSSPType.Flag.E)) return 'E' +// if (DSSPType.is(f, DSSPType.Flag.G)) return 'G' +// if (DSSPType.is(f, DSSPType.Flag.T)) return 'T' +// if (DSSPType.is(f, DSSPType.Flag.S)) return 'S' +// return '-' +// } + +function getDSSPAssignment(flags: Uint32Array, useOriginal = false) { + const getResidueFlag = useOriginal ? getOriginalResidueFlag : getUpdatedResidueFlag + const type = new Uint32Array(flags.length) + for (let i = 0, il = flags.length; i < il; ++i) { + const f = DSSPType.create(flags[i]) + // console.log(i, geFlagName(f)) + type[i] = getResidueFlag(f) + } + return type as unknown as ArrayLike<SecondaryStructureType> +} + +/** + * Constant for electrostatic energy in kcal/mol + * f * q1 * q2 + * Q = -332 * 0.42 * 0.20 + * + * f is the dimensional factor + * + * q1 and q2 are partial charges which are placed on the C,O + * (+q1,-q1) and N,H (-q2,+q2) + */ +const Q = -27.888 + +/** cutoff for hbonds in kcal/mol, must be lower to be consider as an hbond */ +const hbondEnergyCutoff = -0.5 + +/** + * E = Q * (1/r(ON) + l/r(CH) - l/r(OH) - l/r(CN)) + */ +function calcHbondEnergy(oPos: Vec3, cPos: Vec3, nPos: Vec3, hPos: Vec3) { + const distOH = Vec3.distance(oPos, hPos) + const distCH = Vec3.distance(cPos, hPos) + const distCN = Vec3.distance(cPos, nPos) + const distON = Vec3.distance(oPos, nPos) + + const e1 = Q / distOH - Q / distCH + const e2 = Q / distCN - Q / distON + return e1 + e2 +} + +/** + * The basic turn pattern is a single H bond of type (i, i + n). + * We assign an n-turn at residue i if there is an H bond from CO(i) to NH(i + n), + * i.e., “n-turn(i)=: Hbond(i, i + n), n = 3, 4, 5.” + * + * Type: T + */ +function assignTurns(ctx: DSSPContext) { + const { proteinResidues, hbonds, flags, hierarchy } = ctx + const { chains, residueAtomSegments, chainAtomSegments } = hierarchy + const { label_asym_id } = chains + + const turnFlag = [ 0, 0, 0, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5 ] + + for (let i = 0, il = proteinResidues.length; i < il; ++i) { + const rI = proteinResidues[i] + const cI = chainAtomSegments.index[residueAtomSegments.offsets[rI]] + + // TODO should take sequence gaps into account + for (let k = 3; k <= 5; ++k) { + if (i + k >= proteinResidues.length) continue + + const rN = proteinResidues[i + k] + const cN = chainAtomSegments.index[residueAtomSegments.offsets[rN]] + // check if on same chain + if (!label_asym_id.areValuesEqual(cI, cN)) continue + + // check if hbond exists + if (hbonds.getDirectedEdgeIndex(i, i + k) !== -1) { + flags[i] |= turnFlag[k] | DSSPType.Flag.T + } + } + } +} + +/** + * Two nonoverlapping stretches of three residues each, i - 1, i, i + 1 and j - 1, j, j + 1, + * form either a parallel or antiparallel bridge, depending on which of + * two basic patterns is matched. We assign a bridge between residues i and j + * if there are two H bonds characteristic of P-structure; in particular, + * + * Parallel Bridge(i, j) =: + * [Hbond(i - 1, j) and Hbond(j, i + 1)] or + * [Hbond(j - 1, i) and Hbond(i, j + 1)] + * + * Antiparallel Bridge(i, j) =: + * [Hbond(i, j) and Hbond(j, i)] or + * [Hbond(i - 1, j + 1) and Hbond(j - 1, i + l)] + * + * Type: B + */ +function assignBridges(ctx: DSSPContext) { + const { proteinResidues, hbonds, flags } = ctx + + const { offset, b } = hbonds + let i: number, j: number + + for (let k = 0, kl = proteinResidues.length; k < kl; ++k) { + for (let t = offset[k], _t = offset[k + 1]; t < _t; t++) { + const l = b[t] + if (k > l) continue + + // Parallel Bridge(i, j) =: [Hbond(i - 1, j) and Hbond(j, i + 1)] + i = k + 1 // k is i - 1 + j = l + if (i !== j && hbonds.getDirectedEdgeIndex(j, i + 1) !== -1) { + flags[i] |= DSSPType.Flag.B + flags[j] |= DSSPType.Flag.B + } + + // Parallel Bridge(i, j) =: [Hbond(j - 1, i) and Hbond(i, j + 1)] + i = k + j = l - 1 // l is j + 1 + if (i !== j && hbonds.getDirectedEdgeIndex(j - 1, i) !== -1) { + flags[i] |= DSSPType.Flag.B + flags[j] |= DSSPType.Flag.B + } + + // Antiparallel Bridge(i, j) =: [Hbond(i, j) and Hbond(j, i)] + i = k + j = l + if (i !== j && hbonds.getDirectedEdgeIndex(j, i) !== -1) { + flags[i] |= DSSPType.Flag.B + flags[j] |= DSSPType.Flag.B + } + + // Antiparallel Bridge(i, j) =: [Hbond(i - 1, j + 1) and Hbond(j - 1, i + l)] + i = k + 1 + j = l - 1 + if (i !== j && hbonds.getDirectedEdgeIndex(j - 1, i + 1) !== -1) { + flags[i] |= DSSPType.Flag.B + flags[j] |= DSSPType.Flag.B + } + } + } +} + +/** + * A minimal helix is defined by two consecutive n-turns. + * For example, a 4-helix, of minimal length 4 from residues i to i + 3, + * requires 4-turns at residues i - 1 and i, + * + * 3-helix(i,i + 2)=: [3-turn(i - 1) and 3-turn(i)] + * 4-helix(i,i + 3)=: [4-turn(i - 1) and 4-turn(i)] + * 5-helix(i,i + 4)=: [5-turn(i - 1) and 5-turn(i)] + * + * Type: G (n=3), H (n=4), I (n=5) + */ +function assignHelices(ctx: DSSPContext) { + const { proteinResidues, flags } = ctx + + const turnFlag = [ 0, 0, 0, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5 ] + const helixFlag = [ 0, 0, 0, DSSPType.Flag.G, DSSPType.Flag.H, DSSPType.Flag.I ] + + for (let i = 1, il = proteinResidues.length; i < il; ++i) { + const fI = DSSPType.create(flags[i]) + const fI1 = DSSPType.create(flags[i - 1]) + + for (let k = 3; k <= 5; ++k) { + if (DSSPType.is(fI, turnFlag[k]) && DSSPType.is(fI1, turnFlag[k])) { + for (let l = 0; l < k; ++l) { + flags[i + l] |= helixFlag[k] + } + } + } + } +} + +/** + * ladder=: set of one or more consecutive bridges of identical type + * + * Type: E + */ +function assignLadders(ctx: DSSPContext) { + // TODO +} + +/** + * sheet=: set of one or more ladders connected by shared residues + * + * Type: E + */ +function assignSheets(ctx: DSSPContext) { + // TODO +} + +/** + * Bend(i) =: [angle ((CW - Ca(i - 2)),(C"(i + 2) - C"(i))) > 70"] + * + * Type: S + */ +function assignBends(ctx: DSSPContext) { + // TODO +} \ No newline at end of file diff --git a/src/mol-model/structure/model/types.ts b/src/mol-model/structure/model/types.ts index 8d8be296dbf079b0a1e8aa3333d793afc7797940..a550f8f32516575bf3c926ed5a2728e759829e8c 100644 --- a/src/mol-model/structure/model/types.ts +++ b/src/mol-model/structure/model/types.ts @@ -261,10 +261,6 @@ export const IonNames = new Set([ export interface SecondaryStructureType extends BitFlags<SecondaryStructureType.Flag> { } export namespace SecondaryStructureType { - export const Helix = ['h', 'g', 'i'] - export const Sheet = ['e', 'b'] - export const Turn = ['s', 't', 'l', ''] - export const is: (ss: SecondaryStructureType, f: Flag) => boolean = BitFlags.has export const create: (fs: Flag) => SecondaryStructureType = BitFlags.create @@ -375,7 +371,6 @@ export namespace SecondaryStructureType { I: Flag.Helix | Flag.HelixPi, // PI-helix E: Flag.Beta | Flag.BetaSheet, // Extended conformation B: Flag.Beta | Flag.BetaStrand, // Isolated bridge - b: Flag.Beta | Flag.BetaStrand, // Isolated bridge T: Flag.Turn, // Turn C: Flag.NA, // Coil (none of the above) } diff --git a/src/mol-model/structure/util.ts b/src/mol-model/structure/util.ts index 40b5c13dac895a5a43e62a51aee243f932886054..3f65f5f16755cb2ce569e0940f73e07d07dc4c7a 100644 --- a/src/mol-model/structure/util.ts +++ b/src/mol-model/structure/util.ts @@ -52,6 +52,18 @@ export function residueLabel(model: Model, rI: number) { return `${label_asym_id.value(cI)} ${label_comp_id.value(rI)} ${label_seq_id.value(rI)}` } +export function elementLabel(model: Model, index: ElementIndex) { + const { atoms, residues, chains, residueAtomSegments, chainAtomSegments } = model.atomicHierarchy + const { label_atom_id } = atoms + const { auth_seq_id, auth_comp_id } = residues + const { auth_asym_id } = chains + + const residueIndex = residueAtomSegments.index[index] + const chainIndex = chainAtomSegments.index[residueIndex] + + return `[${auth_comp_id.value(residueIndex)}]${auth_seq_id.value(residueIndex)}:${auth_asym_id.value(chainIndex)}.${label_atom_id.value(index)}` +} + // const centerPos = Vec3.zero() // const centerMin = Vec3.zero() // export function getCenterAndRadius(centroid: Vec3, unit: Unit, indices: ArrayLike<number>) { diff --git a/src/tests/browser/render-structure.ts b/src/tests/browser/render-structure.ts new file mode 100644 index 0000000000000000000000000000000000000000..c601c504ffb595d69ffd5f5eef9daaf37a75312e --- /dev/null +++ b/src/tests/browser/render-structure.ts @@ -0,0 +1,83 @@ +/** + * 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 CIF, { CifFrame } from 'mol-io/reader/cif' +import { Model, Structure } from 'mol-model/structure'; +import { ColorTheme } 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 { computeModelDSSP } from 'mol-model/structure/model/properties/utils/secondary-structure'; + +const parent = document.getElementById('app')! +parent.style.width = '100%' +parent.style.height = '100%' + +const canvas = document.createElement('canvas') +canvas.style.width = '100%' +canvas.style.height = '100%' +parent.appendChild(canvas) + +const canvas3d = Canvas3D.create(canvas, parent) +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); + 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) +} + +async function init() { + const cif = await downloadFromPdb('1d66') + const models = await getModels(cif) + console.time('computeModelDSSP') + const secondaryStructure = computeModelDSSP(models[0].atomicHierarchy, models[0].atomicConformation) + console.timeEnd('computeModelDSSP') + ;(models[0].properties as any).secondaryStructure = secondaryStructure + const structure = await getStructure(models[0]) + const cartoonRepr = getCartoonRepr() + + cartoonRepr.setTheme({ + color: reprCtx.colorThemeRegistry.create('secondary-structure', { structure }), + size: reprCtx.sizeThemeRegistry.create('uniform', { structure }) + }) + await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run() + + canvas3d.add(cartoonRepr) + canvas3d.resetCamera() +} + +init() \ No newline at end of file diff --git a/webpack.config.js b/webpack.config.js index 3a2db252989711386f98d744c3c559fb6858991a..756a94fc275d4ebdb70689f34c154c7ca24bd4a2 100644 --- a/webpack.config.js +++ b/webpack.config.js @@ -88,8 +88,9 @@ module.exports = [ createApp('model-server-query'), createBrowserTest('font-atlas'), - createBrowserTest('render-text'), + createBrowserTest('render-mesh'), createBrowserTest('render-shape'), createBrowserTest('render-spheres'), - createBrowserTest('render-mesh') + createBrowserTest('render-structure'), + createBrowserTest('render-text'), ] \ No newline at end of file