Merge pull request #12 from JonStargaryen/dssp

fixes to DSSP impl, detection of bends and sheets

src/mol-math/linear-algebra/3d/vec3.ts

@@ -414,6 +414,7 @@ namespace Vec3 {
     }
 
     const angleTempA = zero(), angleTempB = zero();
+    /** Computes the angle between 2 vectors, reports in rad. */
     export function angle(a: Vec3, b: Vec3) {
         copy(angleTempA, a);
         copy(angleTempB, b);
@@ -433,6 +434,28 @@ namespace Vec3 {
         }
     }
 
+    const tmp_dh_ab = Vec3.zero();
+    const tmp_dh_cb = Vec3.zero();
+    const tmp_dh_bc = Vec3.zero();
+    const tmp_dh_dc = Vec3.zero();
+    const tmp_dh_abc = Vec3.zero();
+    const tmp_dh_bcd = Vec3.zero();
+    const tmp_dh_cross = Vec3.zero();
+    /** Computes the dihedral angles of 4 related atoms. phi: C, N+1, CA+1, C+1 - psi: N, CA, C, N+1 - omega: CA, C, N+1, CA+1 */
+    export function dihedralAngle(a: Vec3, b: Vec3, c: Vec3, d: Vec3) {
+        Vec3.sub(tmp_dh_ab, a, b);
+        Vec3.sub(tmp_dh_cb, c, b);
+        Vec3.sub(tmp_dh_bc, b, c);
+        Vec3.sub(tmp_dh_dc, d, c);
+
+        Vec3.cross(tmp_dh_abc, tmp_dh_ab, tmp_dh_cb);
+        Vec3.cross(tmp_dh_bcd, tmp_dh_bc, tmp_dh_dc);
+
+        const angle = Vec3.angle(tmp_dh_abc, tmp_dh_bcd) * 360.0 / (2 * Math.PI);
+        Vec3.cross(tmp_dh_cross, tmp_dh_abc, tmp_dh_bcd);
+        return Vec3.dot(tmp_dh_cb, tmp_dh_cross) > 0 ? angle : -angle;
+    }
+
     /**
      * Returns whether or not the vectors have exactly the same elements in the same position (when compared with ===)
      */

src/mol-math/linear-algebra/_spec/vec3.spec.ts

@@ -0,0 +1,14 @@
+import { Vec3 } from '../3d'
+
+describe('vec3', () => {
+    const vec1 = [ 1, 2, 3 ] as Vec3;
+    const vec2 = [ 2, 3, 1 ] as Vec3;
+    const orthVec1 = [ 0, 1, 0 ] as Vec3;
+    const orthVec2 = [ 1, 0, 0 ] as Vec3;
+
+    it('angle calculation', () => {
+        expect(Vec3.angle(vec1, vec1) * 360 / (2 * Math.PI)).toBe(0.0);
+        expect(Vec3.angle(orthVec1, orthVec2) * 360 / (2 * Math.PI)).toBe(90.0);
+        expect(Vec3.angle(vec1, vec2)).toBeCloseTo(0.666946);
+    });
+})

src/mol-model/structure/model/properties/seconday-structure.ts

@@ -17,12 +17,17 @@ interface SecondaryStructure {
 }
 
 namespace SecondaryStructure {
-    export type Element = None | Helix | Sheet
+    export type Element = None | Turn | Helix | Sheet
 
     export interface None {
         kind: 'none'
     }
 
+    export interface Turn {
+        kind: 'turn',
+        flags: SecondaryStructureType
+    }
+
     export interface Helix {
         kind: 'helix',
         flags: SecondaryStructureType,

src/mol-model/structure/model/properties/utils/secondary-structure.ts

@@ -2,6 +2,7 @@
  * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
  */
 
 import { SecondaryStructure } from 'mol-model/structure/model/properties/seconday-structure';
@@ -14,13 +15,120 @@ import { IntAdjacencyGraph } from 'mol-math/graph';
 import { BitFlags } from 'mol-util';
 import { ElementIndex } from 'mol-model/structure/model/indexing';
 import { AtomicHierarchy, AtomicConformation } from '../atomic';
+import { ParamDefinition as PD } from 'mol-util/param-definition'
 
-export function computeSecondaryStructure(hierarchy: AtomicHierarchy, conformation: AtomicConformation): SecondaryStructure {
+/**
+ * TODO bugs to fix:
+ * - some turns are not detected correctly: see e.g. pdb:1acj - maybe more than 2 hbonds require some residue to donate electrons
+ * - some sheets are not extended correctly: see e.g. pdb:1acj
+ * - validate new helix definition
+ * - validate new ordering of secondary structure elements
+ */
+
+ /** max distance between two C-alpha atoms to check for hbond */
+const caMaxDist = 9.0;
+
+/**
+ * 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
+/** prevent extremely low hbond energies */
+const hbondEnergyMinimal = -9.9
+
+interface DSSPContext {
+    params: Partial<PD.Values<SecondaryStructureComputationParams>>,
+    getResidueFlag: (f: DSSPType) => SecondaryStructureType,
+    getFlagName: (f: DSSPType) => String,
+
+    hierarchy: AtomicHierarchy
+    proteinResidues: SortedArray<ResidueIndex>
+    /** flags for each residue */
+    flags: Uint32Array
+    hbonds: DsspHbonds,
+
+    torsionAngles: { phi: Float32Array, psi: Float32Array },
+    backboneIndices: BackboneAtomIndices,
+    conformation: AtomicConformation,
+    ladders: Ladder[],
+    bridges: Bridge[]
+}
+
+interface Ladder {
+    previousLadder: number,
+    nextLadder: number,
+    firstStart: number,
+    secondStart: number,
+    secondEnd: number,
+    firstEnd: number,
+    type: BridgeType
+}
+
+const enum BridgeType {
+    PARALLEL = 0x0,
+    ANTI_PARALLEL = 0x1
+}
+
+class Bridge {
+    partner1: number;
+    partner2: number;
+    type: BridgeType;
+
+    constructor(p1: number, p2: number, type: BridgeType) {
+        this.partner1 = Math.min(p1, p2)
+        this.partner2 = Math.max(p1, p2)
+        this.type = type
+    }
+}
+
+type DSSPType = 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,
+        T3S = 0x400, // marks 3-turn start
+        T4S = 0x800,
+        T5S = 0x1000
+    }
+}
+
+export const SecondaryStructureComputationParams = {
+    oldDefinition: PD.Boolean(true, { description: 'Whether to use the old DSSP convention for the annotation of turns and helices, causes them to be two residues shorter' }),
+    oldOrdering: PD.Boolean(true, { description: 'Alpha-helices are preferred over 3-10 helices' })
+}
+export type SecondaryStructureComputationParams = typeof SecondaryStructureComputationParams
+
+export function computeSecondaryStructure(hierarchy: AtomicHierarchy,
+    conformation: AtomicConformation) {
     // 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) {
+export function computeModelDSSP(hierarchy: AtomicHierarchy,
+    conformation: AtomicConformation,
+    params: Partial<PD.Values<SecondaryStructureComputationParams>> = {}): SecondaryStructure {
+    params = { ...PD.getDefaultValues(SecondaryStructureComputationParams), ...params };
+
     const { lookup3d, proteinResidues } = calcAtomicTraceLookup3D(hierarchy, conformation)
     const backboneIndices = calcBackboneAtomIndices(hierarchy, proteinResidues)
     const hbonds = calcBackboneHbonds(hierarchy, conformation, proteinResidues, backboneIndices, lookup3d)
@@ -28,69 +136,103 @@ export function computeModelDSSP(hierarchy: AtomicHierarchy, conformation: Atomi
     const residueCount = proteinResidues.length
     const flags = new Uint32Array(residueCount)
 
+    // console.log(`calculating secondary structure elements using ${ params.oldDefinition ? 'old' : 'revised'} definition and ${ params.oldOrdering ? 'old' : 'revised'} ordering of secondary structure elements`)
+
+    const torsionAngles = calculateDihedralAngles(hierarchy, conformation, proteinResidues, backboneIndices)
+
+    const ladders: Ladder[] = []
+    const bridges: Bridge[] = []
+
+    const getResidueFlag = params.oldDefinition ? getOriginalResidueFlag : getUpdatedResidueFlag
+    const getFlagName = params.oldOrdering ? getOriginalFlagName : getUpdatedFlagName
+
     const ctx: DSSPContext = {
+        params,
+        getResidueFlag,
+        getFlagName,
+
         hierarchy,
         proteinResidues,
         flags,
-        hbonds
+        hbonds,
+
+        torsionAngles,
+        backboneIndices,
+        conformation,
+        ladders,
+        bridges
     }
 
-    assignBends(ctx)
     assignTurns(ctx)
     assignHelices(ctx)
+    assignBends(ctx)
     assignBridges(ctx)
     assignLadders(ctx)
     assignSheets(ctx)
 
-    const assignment = getDSSPAssignment(flags)
-
+    const assignment = getDSSPAssignment(flags, getResidueFlag)
     const type = new Uint32Array(hierarchy.residues._rowCount) as unknown as SecondaryStructureType[]
+    const keys: number[] = []
+    const elements: SecondaryStructure.Element[] = []
+
     for (let i = 0, il = proteinResidues.length; i < il; ++i) {
-        type[proteinResidues[i]] = assignment[i]
+        const assign = assignment[i]
+        type[proteinResidues[i]] = assign
+        const flag = getResidueFlag(flags[i])
+        // TODO is this expected behavior? elements will be strictly split depending on 'winning' flag
+        if (elements.length === 0 /* would fail at very start */ || flag !== (elements[elements.length - 1] as SecondaryStructure.Helix | SecondaryStructure.Sheet | SecondaryStructure.Turn).flags /* flag changed */) {
+            elements[elements.length] = createElement(mapToKind(assign), flags[i], getResidueFlag)
+        }
+        keys[i] = elements.length - 1
     }
 
     const secondaryStructure: SecondaryStructure = {
         type,
-        key: [], // TODO
-        elements: [] // TODO
+        key: keys,
+        elements: elements
     }
+
     return secondaryStructure
 }
 
-interface DSSPContext {
-    hierarchy: AtomicHierarchy
-    proteinResidues: SortedArray<ResidueIndex>
-    /** flags for each residue */
-    flags: Uint32Array
-
-    hbonds: DsspHbonds
+function createElement(kind: string, flag: DSSPType.Flag, getResidueFlag: (f: DSSPType) => SecondaryStructureType): SecondaryStructure.Element {
+    // TODO would be nice to add more detailed information
+    if (kind === 'helix') {
+        return {
+            kind: 'helix',
+            flags: getResidueFlag(flag)
+        } as SecondaryStructure.Helix
+    } else if (kind === 'sheet') {
+        return {
+            kind: 'sheet',
+            flags: getResidueFlag(flag)
+        } as SecondaryStructure.Sheet
+    } else if (kind === 'turn' || kind === 'bend') {
+        return {
+            kind: 'turn',
+            flags: getResidueFlag(flag)
+        }
+    } else {
+        return {
+            kind: 'none'
+        }
+    }
 }
 
-type DSSPType = 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,
+function mapToKind(assignment: SecondaryStructureType.Flag) {
+    if (assignment === SecondaryStructureType.SecondaryStructureDssp.H || assignment === SecondaryStructureType.SecondaryStructureDssp.G || assignment === SecondaryStructureType.SecondaryStructureDssp.I) {
+        return 'helix'
+    } else if (assignment === SecondaryStructureType.SecondaryStructureDssp.B || assignment === SecondaryStructureType.SecondaryStructureDssp.E) {
+        return 'sheet'
+    } else if (assignment === SecondaryStructureType.SecondaryStructureDssp.T) {
+        return 'turn'
+    } else if (assignment === SecondaryStructureType.SecondaryStructureDssp.S) {
+        return 'bend'
+    } else {
+        return 'none'
     }
 }
 
-/** max distance between two C-alpha atoms to check for hbond */
-const caMaxDist = 7.0;
-
-/** min distance between two C-alpha atoms to check for hbond */
-const caMinDist = 4.0;
-
 function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
     const { x, y, z } = conformation;
     const { moleculeType, traceElementIndex } = hierarchy.derived.residue
@@ -141,6 +283,104 @@ function calcBackboneAtomIndices(hierarchy: AtomicHierarchy, proteinResidues: So
 
 type DsspHbonds = IntAdjacencyGraph<{ readonly energies: ArrayLike<number> }>
 
+/**
+ * Bend(i) =: [angle ((CW - Ca(i - 2)),(C"(i + 2) - C"(i))) > 70"]
+ *
+ * Type: S
+ */
+function assignBends(ctx: DSSPContext) {
+    const flags = ctx.flags
+    const { x, y, z } = ctx.conformation
+    const { traceElementIndex } = ctx.hierarchy.derived.residue
+
+    const proteinResidues = ctx.proteinResidues
+    const residueCount = proteinResidues.length
+
+    const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i])
+
+    const caPosPrev2 = Vec3.zero()
+    const caPos = Vec3.zero()
+    const caPosNext2 = Vec3.zero()
+
+    const nIndices = ctx.backboneIndices.nIndices
+    const cPos = Vec3.zero()
+    const nPosNext = Vec3.zero()
+
+    f1: for (let i = 2; i < residueCount - 2; i++) {
+        // check for peptide bond
+        for (let k = 0; k < 4; k++) {
+            let index = i + k - 2
+            position(traceElementIndex[index], cPos)
+            position(nIndices[index + 1], nPosNext)
+            if (Vec3.squaredDistance(cPos, nPosNext) > 6.25 /* max squared peptide bond distance allowed */) {
+                continue f1
+            }
+        }
+
+        const oRIprev2 = proteinResidues[i - 2]
+        const oRI = proteinResidues[i]
+        const oRInext2 = proteinResidues[i + 2]
+
+        const caAtomPrev2 = traceElementIndex[oRIprev2]
+        const caAtom = traceElementIndex[oRI]
+        const caAtomNext2 = traceElementIndex[oRInext2]
+
+        position(caAtomPrev2, caPosPrev2)
+        position(caAtom, caPos)
+        position(caAtomNext2, caPosNext2)
+
+        const caMinus2 = Vec3.zero()
+        const caPlus2 = Vec3.zero()
+
+        Vec3.sub(caMinus2, caPosPrev2, caPos)
+        Vec3.sub(caPlus2, caPos, caPosNext2)
+
+        const angle = Vec3.angle(caMinus2, caPlus2) * 360 / (2 * Math.PI)
+        if (angle && angle > 70.00) {
+            flags[i] |= DSSPType.Flag.S
+        }
+    }
+}
+
+function calculateDihedralAngles(hierarchy: AtomicHierarchy, conformation: AtomicConformation, proteinResidues: SortedArray<ResidueIndex>, backboneIndices: BackboneAtomIndices): { phi: Float32Array, psi: Float32Array } {
+    const { cIndices, nIndices } = 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])
+
+    let cPosPrev = Vec3.zero(), caPosPrev = Vec3.zero(), nPosPrev = Vec3.zero()
+    let cPos = Vec3.zero(), caPos = Vec3.zero(), nPos = Vec3.zero()
+    let cPosNext = Vec3.zero(), caPosNext = Vec3.zero(), nPosNext = Vec3.zero()
+
+    const phi: Float32Array = new Float32Array(residueCount - 1)
+    const psi: Float32Array = new Float32Array(residueCount - 1)
+
+    const cAtomPrev = cIndices[-1], caAtomPrev = traceElementIndex[proteinResidues[-1]], nAtomPrev = nIndices[-1]
+    position(cAtomPrev, cPosPrev), position(caAtomPrev, caPosPrev), position(nAtomPrev, nPosPrev)
+    const cAtom = cIndices[0], caAtom = traceElementIndex[proteinResidues[0]], nAtom = nIndices[0]
+    position(cAtom, cPos), position(caAtom, caPos), position(nAtom, nPos)
+    const cAtomNext = cIndices[1], caAtomNext = traceElementIndex[proteinResidues[1]], nAtomNext = nIndices[1]
+    position(cAtomNext, cPosNext), position(caAtomNext, caPosNext), position(nAtomNext, nPosNext)
+    for (let i = 0; i < residueCount - 1; ++i) {
+        // ignore C-terminal residue as acceptor
+        if (index.findAtomOnResidue(proteinResidues[i], 'OXT') !== -1) continue
+
+        // returns NaN for missing atoms
+        phi[i] = Vec3.dihedralAngle(cPosPrev, nPos, caPos, cPos)
+        psi[i] = Vec3.dihedralAngle(nPos, caPos, cPos, nPosNext)
+
+        cPosPrev = cPos, caPosPrev = caPos, nPosPrev = nPos
+        cPos = cPosNext, caPos = caPosNext, nPos = nPosNext
+
+        position(cIndices[i + 1], cPosNext), position(traceElementIndex[proteinResidues[i + 1]], caPosNext), position(nIndices[i + 1], nPosNext)
+    }
+
+    return { phi, psi };
+}
+
 function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConformation, proteinResidues: SortedArray<ResidueIndex>, backboneIndices: BackboneAtomIndices, lookup3d: GridLookup3D): DsspHbonds {
     const { cIndices, hIndices, nIndices, oIndices } = backboneIndices
     const { index } = hierarchy
@@ -163,8 +403,6 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf
     const cPosPrev = Vec3.zero()
     const oPosPrev = Vec3.zero()
 
-    const caMinDistSq = caMinDist * caMinDist
-
     for (let i = 0, il = proteinResidues.length; i < il; ++i) {
         const oPI = i
         const oRI = proteinResidues[i]
@@ -183,11 +421,9 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf
         position(cAtom, cPos)
         position(caAtom, caPos)
 
-        const { indices, count, squaredDistances } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist)
+        const { indices, count } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist)
 
         for (let j = 0; j < count; ++j) {
-            if (squaredDistances[j] < caMinDistSq) continue
-
             const nPI = indices[j]
 
             // ignore bonds within a residue or to prev or next residue, TODO take chain border into account
@@ -245,8 +481,8 @@ function buildHbondGraph(residueCount: number, oAtomResidues: number[], nAtomRes
 /** 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.B)) return SecondaryStructureType.SecondaryStructureDssp.B
     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
@@ -254,55 +490,50 @@ function getOriginalResidueFlag(f: DSSPType) {
     return SecondaryStructureType.Flag.None
 }
 
+function getOriginalFlagName(f: DSSPType) {
+    if (DSSPType.is(f, DSSPType.Flag.H)) return 'H'
+    if (DSSPType.is(f, DSSPType.Flag.E)) return 'E'
+    if (DSSPType.is(f, DSSPType.Flag.B)) return 'B'
+    if (DSSPType.is(f, DSSPType.Flag.G)) return 'G'
+    if (DSSPType.is(f, DSSPType.Flag.I)) return 'I'
+    if (DSSPType.is(f, DSSPType.Flag.T)) return 'T'
+    if (DSSPType.is(f, DSSPType.Flag.S)) return 'S'
+    return '-'
+}
+
 /** 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.B)) return SecondaryStructureType.SecondaryStructureDssp.B
     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
+function getUpdatedFlagName(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.E)) return 'E'
+    if (DSSPType.is(f, DSSPType.Flag.B)) return 'B'
+    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, getResidueFlag: (f: DSSPType) => SecondaryStructureType) {
     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))
  */
@@ -314,7 +545,13 @@ function calcHbondEnergy(oPos: Vec3, cPos: Vec3, nPos: Vec3, hPos: Vec3) {
 
     const e1 = Q / distOH - Q / distCH
     const e2 = Q / distCN - Q / distON
-    return e1 + e2
+    const e = e1 + e2
+
+    // cap lowest possible energy
+    if (e < hbondEnergyMinimal)
+        return hbondEnergyMinimal
+
+    return e
 }
 
 /**
@@ -329,24 +566,31 @@ function assignTurns(ctx: DSSPContext) {
     const { chains, residueAtomSegments, chainAtomSegments } = hierarchy
     const { label_asym_id } = chains
 
-    const turnFlag = [0, 0, 0, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5]
+    const turnFlag = [DSSPType.Flag.T3S, DSSPType.Flag.T4S, DSSPType.Flag.T5S, 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]]
+    for (let idx = 0; idx < 3; idx++) {
+        for (let i = 0, il = proteinResidues.length - 1; 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]
+            // TODO should take sequence gaps into account
+            const rN = proteinResidues[i + idx + 3]
             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
+            if (hbonds.getDirectedEdgeIndex(i, i + idx + 3) !== -1) {
+                flags[i] |= turnFlag[idx + 3] | turnFlag[idx]
+                if (ctx.params.oldDefinition) {
+                    for (let k = 1; k < idx + 3; ++k) {
+                        flags[i + k] |= turnFlag[idx + 3] | DSSPType.Flag.T
+                    }
+                } else {
+                    for (let k = 0; k <= idx + 3; ++k) {
+                        flags[i + k] |= turnFlag[idx + 3] | DSSPType.Flag.T
+                    }
+                }
             }
         }
     }
@@ -369,7 +613,7 @@ function assignTurns(ctx: DSSPContext) {
  * Type: B
  */
 function assignBridges(ctx: DSSPContext) {
-    const { proteinResidues, hbonds, flags } = ctx
+    const { proteinResidues, hbonds, flags, bridges } = ctx
 
     const { offset, b } = hbonds
     let i: number, j: number
@@ -385,6 +629,8 @@ function assignBridges(ctx: DSSPContext) {
             if (i !== j && hbonds.getDirectedEdgeIndex(j, i + 1) !== -1) {
                 flags[i] |= DSSPType.Flag.B
                 flags[j] |= DSSPType.Flag.B
+                // TODO move to constructor, actually omit object all together
+                bridges[bridges.length] = new Bridge(i, j, BridgeType.PARALLEL)
             }
 
             // Parallel Bridge(i, j) =: [Hbond(j - 1, i) and Hbond(i, j + 1)]
@@ -393,6 +639,7 @@ function assignBridges(ctx: DSSPContext) {
             if (i !== j && hbonds.getDirectedEdgeIndex(j - 1, i) !== -1) {
                 flags[i] |= DSSPType.Flag.B
                 flags[j] |= DSSPType.Flag.B
+                bridges[bridges.length] = new Bridge(j, i, BridgeType.PARALLEL)
             }
 
             // Antiparallel Bridge(i, j) =: [Hbond(i, j) and Hbond(j, i)]
@@ -401,6 +648,7 @@ function assignBridges(ctx: DSSPContext) {
             if (i !== j && hbonds.getDirectedEdgeIndex(j, i) !== -1) {
                 flags[i] |= DSSPType.Flag.B
                 flags[j] |= DSSPType.Flag.B
+                bridges[bridges.length] = new Bridge(j, i, BridgeType.ANTI_PARALLEL)
             }
 
             // Antiparallel Bridge(i, j) =: [Hbond(i - 1, j + 1) and Hbond(j - 1, i + l)]
@@ -409,9 +657,12 @@ function assignBridges(ctx: DSSPContext) {
             if (i !== j && hbonds.getDirectedEdgeIndex(j - 1, i + 1) !== -1) {
                 flags[i] |= DSSPType.Flag.B
                 flags[j] |= DSSPType.Flag.B
+                bridges[bridges.length] = new Bridge(j, i, BridgeType.ANTI_PARALLEL)
             }
         }
     }
+
+    bridges.sort((a, b) => a.partner1 > b.partner1 ? 1 : a.partner1 < b.partner1 ? -1 : 0)
 }
 
 /**
@@ -428,17 +679,41 @@ function assignBridges(ctx: DSSPContext) {
 function assignHelices(ctx: DSSPContext) {
     const { proteinResidues, flags } = ctx
 
-    const turnFlag = [0, 0, 0, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5]
+    const turnFlag = [DSSPType.Flag.T3S, DSSPType.Flag.T4S, DSSPType.Flag.T5S, 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])
+    const helixCheckOrder = ctx.params.oldOrdering ? [4, 3, 5] : [3, 4, 5]
+    for (let ni = 0; ni < helixCheckOrder.length; ni++) {
+        const n = helixCheckOrder[ni]
+
+        for (let i = 1, il = proteinResidues.length - n; i < il; i++) {
+            const fI = DSSPType.create(flags[i])
+            const fI1 = DSSPType.create(flags[i - 1])
+            const fI2 = 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]
+            // TODO rework to elegant solution which will not break instantly
+            if (ctx.params.oldOrdering) {
+                if ((n === 3 && (DSSPType.is(fI, DSSPType.Flag.H) || DSSPType.is(fI2, DSSPType.Flag.H)) || // for 3-10 yield to alpha helix
+                    (n === 5 && ((DSSPType.is(fI, DSSPType.Flag.H) || DSSPType.is(fI, DSSPType.Flag.G)) || (DSSPType.is(fI2, DSSPType.Flag.H) || DSSPType.is(fI2, DSSPType.Flag.G)))))) { // for pi yield to all other helices
+                    continue
+                }
+            } else {
+                if ((n === 4 && (DSSPType.is(fI, DSSPType.Flag.G) || DSSPType.is(fI2, DSSPType.Flag.G)) || // for alpha helix yield to 3-10
+                    (n === 5 && ((DSSPType.is(fI, DSSPType.Flag.H) || DSSPType.is(fI, DSSPType.Flag.G)) || (DSSPType.is(fI2, DSSPType.Flag.H) || DSSPType.is(fI2, DSSPType.Flag.G)))))) { // for pi yield to all other helices
+                    continue
+                }
+            }
+
+            if (DSSPType.is(fI, turnFlag[n]) && DSSPType.is(fI, turnFlag[n - 3]) && // check fI for turn start of proper type
+                DSSPType.is(fI1, turnFlag[n]) && DSSPType.is(fI1, turnFlag[n - 3])) { // check fI1 accordingly
+                if (ctx.params.oldDefinition) {
+                    for (let k = 0; k < n; k++) {
+                        flags[i + k] |= helixFlag[n]
+                    }
+                } else {
+                    for (let k = -1; k <= n; k++) {
+                        flags[i + k] |= helixFlag[n]
+                    }
                 }
             }
         }
@@ -451,23 +726,137 @@ function assignHelices(ctx: DSSPContext) {
  * Type: E
  */
 function assignLadders(ctx: DSSPContext) {
-    // TODO
+    const { bridges, ladders } = ctx
+
+    // create ladders
+    for (let bridgeIndex = 0; bridgeIndex < bridges.length; bridgeIndex++) {
+        const bridge = bridges[bridgeIndex]
+        let found = false
+        for (let ladderIndex = 0; ladderIndex < ladders.length; ladderIndex++) {
+            const ladder = ladders[ladderIndex]
+            if (shouldExtendLadder(ladder, bridge)) {
+                found = true
+                ladder.firstEnd++
+                if (bridge.type === BridgeType.PARALLEL) {
+                    ladder.secondEnd++
+                } else {
+                    ladder.secondStart--
+                }
+            }
+        }
+
+        // no suitable assignment: create new ladder with single bridge as content
+        if (!found) {
+            ladders[ladders.length] = {
+                previousLadder: 0,
+                nextLadder: 0,
+                firstStart: bridge.partner1,
+                firstEnd: bridge.partner1,
+                secondStart: bridge.partner2,
+                secondEnd: bridge.partner2,
+                type: bridge.type
+            }
+        }
+    }
+
+    // connect ladders
+    for (let ladderIndex1 = 0; ladderIndex1 < ladders.length; ladderIndex1++) {
+        const ladder1 = ladders[ladderIndex1]
+        for (let ladderIndex2 = ladderIndex1; ladderIndex2 < ladders.length; ladderIndex2++) {
+            const ladder2 = ladders[ladderIndex2]
+            if (resemblesBulge(ladder1, ladder2)) {
+                ladder1.nextLadder = ladderIndex2
+                ladder2.previousLadder = ladderIndex1
+            }
+        }
+    }
 }
 
 /**
- * sheet=: set of one or more ladders connected by shared residues
- *
- * Type: E
+ * For beta structures, we define: a bulge-linked ladder consists of two ladders or bridges of the same type
+ * connected by at most one extra residue of one strand and at most four extra residues  on the other strand,
+ * all residues in bulge-linked ladders are marked E, including any extra residues.
  */
-function assignSheets(ctx: DSSPContext) {
-    // TODO
+function resemblesBulge(ladder1: Ladder, ladder2: Ladder) {
+    if (!(ladder1.type === ladder2.type && ladder2.firstStart - ladder1.firstEnd < 6 &&
+        ladder1.firstStart < ladder2.firstStart && ladder2.nextLadder === 0)) return false
+
+    if (ladder1.type === BridgeType.PARALLEL) {
+        return bulgeCriterion2(ladder1, ladder2)
+    } else {
+        return bulgeCriterion2(ladder2, ladder1)
+    }
+}
+
+function bulgeCriterion2(ladder1: Ladder, ladder2: Ladder) {
+    return ladder2.secondStart - ladder1.secondEnd > 0 && ((ladder2.secondStart - ladder1.secondEnd < 6 &&
+        ladder2.firstStart - ladder1.firstEnd < 3) || ladder2.secondStart - ladder1.secondEnd < 3)
+}
+
+function shouldExtendLadder(ladder: Ladder, bridge: Bridge): boolean {
+    // in order to extend ladders, same type must be present
+    if (bridge.type !== ladder.type) return false
+
+    // only extend if residue 1 is sequence neighbor with regard to ladder
+    if (bridge.partner1 !== ladder.firstEnd + 1) return false
+
+    if (bridge.type === BridgeType.PARALLEL) {
+        if (bridge.partner2 === ladder.secondEnd + 1) {
+            return true
+        }
+    } else {
+        if (bridge.partner2 === ladder.secondStart - 1) {
+            return true
+        }
+    }
+
+    return false
+}
+
+function isHelixType(f: DSSPType) {
+    return DSSPType.is(f, DSSPType.Flag.G) || DSSPType.is(f, DSSPType.Flag.H) || DSSPType.is(f, DSSPType.Flag.I)
 }
 
 /**
- * Bend(i) =: [angle ((CW - Ca(i - 2)),(C"(i + 2) - C"(i))) > 70"]
+ * sheet=: set of one or more ladders connected by shared residues
  *
- * Type: S
+ * Type: E
  */
-function assignBends(ctx: DSSPContext) {
-    // TODO
+function assignSheets(ctx: DSSPContext) {
+    const { ladders, flags } = ctx
+    for (let ladderIndex = 0; ladderIndex < ladders.length; ladderIndex++) {
+        const ladder = ladders[ladderIndex]
+        for (let lcount = ladder.firstStart; lcount <= ladder.firstEnd; lcount++) {
+            const diff = ladder.firstStart - lcount
+            const l2count = ladder.secondStart - diff
+
+            if (ladder.firstStart !== ladder.firstEnd) {
+                flags[lcount] |= DSSPType.Flag.E
+                flags[l2count] |= DSSPType.Flag.E
+            } else {
+                if (!isHelixType(flags[lcount]) && DSSPType.is(flags[lcount], DSSPType.Flag.E)) {
+                    flags[lcount] |= DSSPType.Flag.B
+                }
+                if (!isHelixType(flags[l2count]) && DSSPType.is(flags[l2count], DSSPType.Flag.E)) {
+                    flags[l2count] |= DSSPType.Flag.B
+                }
+            }
+        }
+
+        if (ladder.nextLadder === 0) continue
+
+        const conladder = ladders[ladder.nextLadder]
+        for (let lcount = ladder.firstStart; lcount <= conladder.firstEnd; lcount++) {
+            flags[lcount] |= DSSPType.Flag.E
+        }
+        if (ladder.type === BridgeType.PARALLEL) {
+            for (let lcount = ladder.secondStart; lcount <= conladder.secondEnd; lcount++) {
+                flags[lcount] |= DSSPType.Flag.E
+            }
+        } else {
+            for (let lcount = conladder.secondEnd; lcount <= ladder.secondStart; lcount++) {
+                flags[lcount] |= DSSPType.Flag.E
+            }
+        }
+    }
 }

src/mol-model/structure/model/properties/utils/secondary-structure.validation

@@ -0,0 +1,75 @@
+compares Mol* port of DSSP (with default parameters) to the BioJava implementation
+
+### pdb:1pga ###
+# turns #
+Mol*: ----------------------TTTTTTTTTTTTTTTT--------TTTT------
+53 turns, 18 openings
+DSSP: ----------------------TTTTTTTTTTTTTTTT--------TTTT------
+53 turns, 18 openings
+
+# bends #
+Mol*: ---------SS---------SSSSSSSSSSSSSSSSSS--------SSS-------
+23 bends
+DSSP: ---------SS---------SSSSSSSSSSSSSSSSSS--------SSS-------
+23 bends
+
+# helices #
+Mol*: ----------------------HHHHHHHHHHHHHHTT--------TTTT------
+44 helix elements - 0 3-10, 44 alpha, 0 pi
+DSSP: ----------------------HHHHHHHHHHHHHHTT--------TTTT------
+44 helix elements - 0 3-10, 44 alpha, 0 pi
+
+# all #
+Mol*: -EEEEEEE-SS-EEEEEEE-SSHHHHHHHHHHHHHHTT---EEEEETTTTEEEEE-
+DSSP: -EEEEEEE-SS-EEEEEEE-SSHHHHHHHHHHHHHHTT---EEEEETTTTEEEEE-
+
+
+### pdb:1bta ###
+# turns #
+Mol*: ------TTT---TTTTTTTTTTTTT--TT----TTTTTTTTTTT-----------TTTTTTTTT--TTTTTTTTTTTTTTT--------
+127 turns, 44 openings
+DSSP: ------TTT---TTTTTTTTTTTTT--TT----TTTTTTTTTTT-----------TTTTTTTTT--TTTTTTTTTTTTTTT--------
+127 turns, 44 openings
+
+# bends #
+Mol*: ------SSS--SSSSSSSSSSSSS---SS---SSSSSSSSSSSSS-SS------SSSSSSSSSSSSSSSSSSSSSSSSSSS--------
+60 bends
+DSSP: ------SSS--SSSSSSSSSSSSS---SS---SSSSSSSSSSSSS-SS------SSSSSSSSSSSSSSSSSSSSSSSSSSS--------
+60 bends
+
+# helices #
+Mol*: ------TTT---HHHHHHHHHHHHT--TT----HHHHHHHHTTT-----------TTHHHHTTT--HHHHHHHHHHHHHTT--------
+100 helix elements - 0 3-10, 100 alpha, 0 pi
+DSSP: ------TTT---HHHHHHHHHHHHT--TT----HHHHHHHHTTT-----------TTHHHHTTT--HHHHHHHHHHHHHTT--------
+100 helix elements - 0 3-10, 100 alpha, 0 pi
+
+# all #
+Mol*: -EEEEETTT--SHHHHHHHHHHHHT--TT---SHHHHHHHHTTTS-SSEEEEEESTTHHHHTTTSSHHHHHHHHHHHHHTT--EEEEE-
+DSSP: -EEEEETTT--SHHHHHHHHHHHHT--TT---SHHHHHHHHTTTS-SSEEEEEESTTHHHHTTTSSHHHHHHHHHHHHHTT--EEEEE-
+
+
+### pdb:1acj ###
+# turns #
+Mol*: -------TT----------TT----------------TTTTT----------------------------TTTT-TTTTTT----------------------------------TTT------TTT-TTTTTTTTT-----------TTTT---TT-------TTTTTTTTTTTTTTTTTTTTT--TT-------TTTTTTTTTTTT-TTTTTT----------TT-TT----TTTTTTTTTTTTTTTT-----TTTTTTTTTT--TTTTTTTTTTT-----------------------TTTTTTTT----------------TTTTTTT-TT--TT------TTTTTTTTTTTTTT--TTTTTTTTTTT--TTTTT-TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT-------------TT----TTT---TTTTTTTTTTTTT-TTT---TTTTTTTTTTTTTTTTTTTT--------------------------------TTTTTTTTTTTTTTTTTTT-
+614 turns, 223 openings
+DSSP: -------TT----------TT----------------TTTTT------------------------------TT-TTTTTT----------------------------------TTT------TTT-TTTTTTTTT-----------TTTT---TT-------TTTTTTTTTTTTTTTTTTTTT--TT-------TTTTTTTTTTTT-TTTTTT----------TT-TT----TTTTTTTTTTTTTTTT-----TTTTTTTTTT--TTTTTTTTTTT-----------------------TTTTTTTT----------------TTTTTTT-TT--TT------TTTTTTTTTTT-TT--TTTTTTTTTTT--TTTTT-TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT-------------TT----TTT---TTTTTTTTTTTTT-TTT---TTTTTTTTTTTTTTTTTTTT--------------------------------TTTTTTTTTTTTTTTTTTT-
+606 turns, 220 openings
+
+# bends #
+Mol*: --S----SS----------SSS----------S----SS-SSS--------SS-----S-----------SSSS-SSSSSSS--S---S----------SS--SS---------SSSS---S-SSSS--SSSSSSS-----------SSSS----SS-SSS-S-SSSSSSSSSSSSSSSSSSSSS--SSS------SSSSSSSSSSSS-SSSSSS-S----SS--SSSSSS---SSSSSSSSSSSSSSS----S-SSSSSSSSSSS-SSSSSSSSSS--SS--SS--S------SSSSSS-SSSSSSS--S--S-------S--SSSSSSSSSSS--SSS-----SSSSSSSSSSS-SS--SSSSSSSSSSS--SSSSS-SSSSSSSSSSSSSSSSS----SSSSSSSSSSSS-----------SS--S-SSS-S-SS-SSSSSS--SS-SSS---SSSSSSSSSSSSSSSSSSSSSSSS---------SSS------SSSS----SSSS-S----SSSSSSSSSS--
+305 bends
+DSSP: --S----SS----------SSS----------S----SS-SSS--------SS-----S-----------SSSS-SSSSSSS--S---S----------SS--SS---------SSSS---S-SSSS--SSSSSSS-----------SSSS----SS-SSS-S-SSSSSSSSSSSSSSSSSSSSS--SSS------SSSSSSSSSSSS-SSSSSS-S----SS--SSSSSS---SSSSSSSSSSSSSSS----S-SSSSSSSSSSS-SSSSSSSSSS--SS--SS--S------SSSSSS-SSSSSSS--S--S-------S--SSSSSSSSSSS--SSS-----SSSSSSSSSSS-SS--SSSSSSSSSSS--SSSSS-SSSSSSSSSSSSSSSSS----SSSSSSSSSSSS-----------SS--S-SSS-S-SS-SSSSSS--SS-SSS---SSSSSSSSSSSSSSSSSSSSSSSS---------SSS------SSSS----SSSS-S----SSSSSSSSSS--
+305 bends
+
+# helices #
+Mol*: -------TT----------TT----------------GGGTT----------------------------TTTT-HHHHTT----------------------------------TTT------GGG-THHHHHHHT-----------HHHH---TT-------HHHHHHHHHHHHHHHHGGGGT--TT-------THHHHHHHHHHH-HHHHTT----------TT-TT----HHHHHHHHHHHHHHTT-----HHHHHHHHHH--HHHHHHHHGGG-----------------------HHHHHHHT----------------HHHHHHH-TT--TT------HHHHHHHHHHHTTT--HHHHHHHHHHH--GGGTT-HHHHHHHHHHHHHHHHTHHHHHHHHHHHHTT-------------TT----GGG---TTTTHHHHTTGGG-GGG---HHHHHHHHHHHHHHHHHHHT--------------------------------TTHHHHHHHHTHHHHHHHH-
+523 helix elements - 27 3-10, 496 alpha, 0 pi
+DSSP: -------TT----------TT----------------GGGTT------------------------------TT-HHHHTT----------------------------------TTT------GGG-THHHHHHHT-----------HHHH---TT-------HHHHHHHHHHHHHHHHGGGGT--TT-------THHHHHHHHHHH-HHHHTT----------TT-TT----HHHHHHHHHHHHHHTT-----HHHHHHHHHH--HHHHHHHHGGG-----------------------HHHHHHHT----------------HHHHHHH-TT--TT------HHHHHHHHHHH-TT--HHHHHHHHHHH--GGGTT-HHHHHHHHHHHHHHHHTHHHHHHHHHHHHTT-------------TT----GGG---TTTTHHHHTTGGG-GGG---HHHHHHHHHHHHHHHHHHHT--------------------------------TTHHHHHHHHTHHHHHHHH-
+523 helix elements - 27 3-10, 496 alpha, 0 pi
+
+# all #
+Mol*: --SEEEETTEEEE-EEEEETTEEEEEEEEEE-EE---GGGTTS--EE----SSEEE--S---B-------TTTT-HHHHTTS--S-B-S---EEEEEE-SS--SSEEEEEEE--STTT---S-SGGG-THHHHHHHT-EEEE-----SHHHH---TT-SSS-S-HHHHHHHHHHHHHHHHGGGGTEEEEEEEEEEETHHHHHHHHHHH-HHHHTT-SEEEEES--TTSTTSEEEHHHHHHHHHHHHHHTT---S-HHHHHHHHHHS-HHHHHHHHGGG-SS--SS--S--EEE-SSSSSS-HHHHHHHT-S--S-EEEEEESB-SHHHHHHHSTT--TTS-----HHHHHHHHHHHTTT--HHHHHHHHHHH--GGGTT-HHHHHHHHHHHHHHHHTHHHHHHHHHHHHTTSS-EEEEEE----TT--S-GGG-SBTTTTHHHHTTGGG-GGG---HHHHHHHHHHHHHHHHHHHTSSSS---------SSS-EEEEESSSS--EEESTTHHHHHHHHTHHHHHHHH-
+DSSP: --SEEEETTEEEE-EEEEETTEEEEEEEEEE-EE---GGGTTS--EE----SSEEE--S---B-------SSTT-HHHHTTS--S-B-S---EEEEEE-SS--SSEEEEEEE--STTT---S-SGGG-THHHHHHHT-EEEE-----SHHHH---TT-SSS-S-HHHHHHHHHHHHHHHHGGGGTEEEEEEEEEEETHHHHHHHHHHH-HHHHTT-SEEEEES--TTSTTS-EEHHHHHHHHHHHHHHTT---S-HHHHHHHHHHS-HHHHHHHHGGG-SS--SS--S---EE-SSSSSS-HHHHHHHT-S--S-EEEEEESB-SHHHHHHHSTT--TTS-----HHHHHHHHHHH-TT--HHHHHHHHHHH--GGGTT-HHHHHHHHHHHHHHHHTHHHHHHHHHHHHTTSS-EEEEEE----TT--S-GGG-SBTTTTHHHHTTGGG-GGG---HHHHHHHHHHHHHHHHHHHTSSSS---------SSS-EEEEESSSS--EEESTTHHHHHHHHTHHHHHHHH-
+
+TODO fix mismatches                                                   e.g. here
+TODO move to spec.ts once tests are running

src/mol-model/structure/model/types.ts

@@ -275,41 +275,42 @@ export namespace SecondaryStructureType {
         DoubleHelix = 0x1,
         Helix = 0x2,
         Beta = 0x4,
-        Turn = 0x8,
+        Bend = 0x8,
+        Turn = 0x10,
 
         // category variant
-        LeftHanded = 0x10,  // helix
-        RightHanded = 0x20,
+        LeftHanded = 0x20,  // helix
+        RightHanded = 0x40,
 
-        ClassicTurn = 0x40,  // turn
-        InverseTurn = 0x80,
+        ClassicTurn = 0x80,  // turn
+        InverseTurn = 0x100,
 
         // sub-category
-        HelixOther = 0x100,  // protein
-        Helix27 = 0x200,
-        Helix3Ten = 0x400,
-        HelixAlpha = 0x800,
-        HelixGamma = 0x1000,
-        HelixOmega = 0x2000,
-        HelixPi = 0x4000,
-        HelixPolyproline = 0x8000,
-
-        DoubleHelixOther = 0x10000,  // nucleic
-        DoubleHelixZ = 0x20000,
-        DoubleHelixA = 0x40000,
-        DoubleHelixB = 0x80000,
-
-        BetaOther = 0x100000,  // protein
-        BetaStrand = 0x200000,  // single strand
-        BetaSheet = 0x400000,  // multiple hydrogen bonded strands
-        BetaBarell = 0x800000,  // closed series of sheets
-
-        TurnOther = 0x1000000,  // protein
-        Turn1 = 0x2000000,
-        Turn2 = 0x4000000,
-        Turn3 = 0x8000000,
-
-        NA = 0x10000000,  // not applicable/available
+        HelixOther = 0x200,  // protein
+        Helix27 = 0x400,
+        Helix3Ten = 0x800,
+        HelixAlpha = 0x1000,
+        HelixGamma = 0x2000,
+        HelixOmega = 0x4000,
+        HelixPi = 0x8000,
+        HelixPolyproline = 0x10000,
+
+        DoubleHelixOther = 0x20000,  // nucleic
+        DoubleHelixZ = 0x40000,
+        DoubleHelixA = 0x80000,
+        DoubleHelixB = 0x100000,
+
+        BetaOther = 0x200000,  // protein
+        BetaStrand = 0x400000,  // single strand
+        BetaSheet = 0x800000,  // multiple hydrogen bonded strands
+        BetaBarell = 0x1000000,  // closed series of sheets
+
+        TurnOther = 0x2000000,  // protein
+        Turn1 = 0x4000000,
+        Turn2 = 0x8000000,
+        Turn3 = 0x10000000,
+
+        NA = 0x20000000,  // not applicable/available
     }
 
     export const SecondaryStructureMmcif: { [value: string]: number } = {
@@ -386,7 +387,7 @@ export namespace SecondaryStructureType {
         G: Flag.Helix | Flag.Helix3Ten,  // 3-helix (310 helix)
         I: Flag.Helix | Flag.HelixPi,  // 5 helix (pi-helix)
         T: Flag.Turn,  // hydrogen bonded turn
-        S: Flag.Turn,  // bend
+        S: Flag.Bend,  // bend
     }
 }
 

src/mol-theme/color/secondary-structure.ts

@@ -22,6 +22,8 @@ const SecondaryStructureColors = ColorMap({
     'betaTurn': 0x6080FF,
     'betaStrand': 0xFFC800,
     'coil': 0xFFFFFF,
+    'bend': 0x66D8C9 /* biting original color used 0x00FF00 */,
+    'turn': 0x00B266,
 
     'dna': 0xAE00FE,
     'rna': 0xFD0162,
@@ -53,8 +55,10 @@ export function secondaryStructureColor(unit: Unit, element: ElementIndex): Colo
         return SecondaryStructureColors.alphaHelix
     } else if (SecondaryStructureType.is(secStrucType, SecondaryStructureType.Flag.Beta)) {
         return SecondaryStructureColors.betaStrand
+    } else if (SecondaryStructureType.is(secStrucType, SecondaryStructureType.Flag.Bend)) {
+        return SecondaryStructureColors.bend
     } else if (SecondaryStructureType.is(secStrucType, SecondaryStructureType.Flag.Turn)) {
-        return SecondaryStructureColors.coil
+        return SecondaryStructureColors.turn
     } else {
         const moleculeType = getElementMoleculeType(unit, element)
         if (moleculeType === MoleculeType.DNA) {