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);

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

@@ -20,7 +20,8 @@ namespace SecondaryStructure {
     export type Element = None | Helix | Sheet
 
     export interface None {
-        kind: 'none'
+        kind: 'none',
+        flags?: SecondaryStructureType // TODO should this level of detail be added to non-defined secondary structure elements?
     }
 
     export interface Helix {

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

@@ -15,12 +15,25 @@ 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 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
+ */
+
+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) {
+export function computeModelDSSP(hierarchy: AtomicHierarchy,
+    conformation: AtomicConformation,
+    oldDefinition = true,
+    oldOrdering = true) {
+    // console.log(`calculating secondary structure elements using ${ oldDefinition ? 'old' : 'revised'} definition and ${ oldOrdering ? 'old' : 'revised'} ordering of secondary structure elements`)
     const { lookup3d, proteinResidues } = calcAtomicTraceLookup3D(hierarchy, conformation)
     const backboneIndices = calcBackboneAtomIndices(hierarchy, proteinResidues)
     const hbonds = calcBackboneHbonds(hierarchy, conformation, proteinResidues, backboneIndices, lookup3d)
@@ -28,42 +41,153 @@ export function computeModelDSSP(hierarchy: AtomicHierarchy, conformation: Atomi
     const residueCount = proteinResidues.length
     const flags = new Uint32Array(residueCount)
 
+    const torsionAngles = calculateDihedralAngles(hierarchy, conformation, proteinResidues, backboneIndices)
+
+    const ladders: Ladder[] = []
+    const bridges: Bridge[] = []
+
+    const getResidueFlag = oldOrdering ? getOriginalResidueFlag : getUpdatedResidueFlag
+    const getFlagName = oldOrdering ? getOriginalFlagName : getUpdatedFlagName
+
     const ctx: DSSPContext = {
+        oldDefinition,
+        oldOrdering,
+        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])
+        // const kind = mapToKind(assign)
+        // TODO is this expected behavior? elements will be strictly split depending on 'winning' flag
+        if (elements.length === 0 || // check ought to fail at very start
+            // elements[elements.length - 1].kind !== kind || // new element if overall kind changed
+            flag !== (elements[elements.length - 1] as SecondaryStructure.Helix | SecondaryStructure.Sheet).flags) { // exact 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
     }
+
+    // console.log(keys)
+    // console.log(elements)
+    // printDSSPString(flags, getFlagName)
+
     return secondaryStructure
 }
 
+// function printDSSPString(flags: Uint32Array, getFlagName: (f: DSSPType) => String) {
+//     let out = ''
+//     for (let i = 0, il = flags.length; i < il; ++i) {
+//         const f = DSSPType.create(flags[i])
+//         out += getFlagName(f)
+//     }
+//     console.log(out)
+// }
+
+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 {
+        return {
+            kind: 'none',
+            flags: getResidueFlag(flag)
+        }
+    }
+}
+
+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 {
+            return 'none'
+        }
+}
+
 interface DSSPContext {
+    /** Whether to use the old DSSP convention for the annotation of turns and helices, causes them to be two residues shorter */
+    oldDefinition: boolean,
+    /** alpha-helices are preferred over 3-10 helix */
+    oldOrdering: boolean,
+    getResidueFlag: (f: DSSPType) => SecondaryStructureType,
+    getFlagName: (f: DSSPType) => String,
+
     hierarchy: AtomicHierarchy
     proteinResidues: SortedArray<ResidueIndex>
     /** flags for each residue */
     flags: Uint32Array
+    hbonds: DsspHbonds,
 
-    hbonds: DsspHbonds
+    torsionAngles: { phi: number[], psi: number[], omega: number[] },
+    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
+}
+
+interface Bridge {
+    partner1: number,
+    partner2: number,
+    type: BridgeType
 }
 
 type DSSPType = BitFlags<DSSPType.Flag>
@@ -82,6 +206,9 @@ namespace DSSPType {
         T3 = 0x80,
         T4 = 0x100,
         T5 = 0x200,
+        T3S = 0x400, // marks 3-turn start
+        T4S = 0x800,
+        T5S = 0x1000
     }
 }
 
@@ -89,7 +216,7 @@ namespace DSSPType {
 const caMaxDist = 7.0;
 
 /** min distance between two C-alpha atoms to check for hbond */
-const caMinDist = 4.0;
+const caMinDist = 0.5;
 
 function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
     const { x, y, z } = conformation;
@@ -141,6 +268,150 @@ 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()
+    // let bends = 0
+
+    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) {
+            // console.log(`found bend at ${ i } with kappa ${ angle }`)
+            flags[i] |= DSSPType.Flag.S
+            // bends++
+        }
+    }
+    // console.log(bends + ' bends')
+}
+
+function calculateDihedralAngles(hierarchy: AtomicHierarchy, conformation: AtomicConformation, proteinResidues: SortedArray<ResidueIndex>, backboneIndices: BackboneAtomIndices): { phi: number[], psi: number[], omega: number[] } {
+    const { cIndices, nIndices } = backboneIndices
+    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 cPosPrev = Vec3.zero()
+    const caPosPrev = Vec3.zero()
+    const nPosPrev = Vec3.zero()
+    const cPos = Vec3.zero()
+    const caPos = Vec3.zero()
+    const nPos = Vec3.zero()
+    const cPosNext = Vec3.zero()
+    const caPosNext = Vec3.zero()
+    const nPosNext = Vec3.zero()
+
+    const phi: number[] = []
+    const psi: number[] = []
+    const omega: number[] = []
+
+    for (let i = 0; i < residueCount - 1; ++i) {
+        const oPIprev = i - 1
+        const oRIprev = proteinResidues[i - 1]
+        const oPI = i
+        const oRI = proteinResidues[i]
+        const oPInext = i + 1
+        const oRInext = proteinResidues[i + 1]
+
+        const cAtomPrev = cIndices[oPIprev]
+        const caAtomPrev = traceElementIndex[oRIprev]
+        const nAtomPrev = nIndices[oPIprev]
+        const cAtom = cIndices[oPI]
+        const caAtom = traceElementIndex[oRI]
+        const nAtom = nIndices[oPI]
+        const cAtomNext = cIndices[oPInext]
+        const caAtomNext = traceElementIndex[oRInext]
+        const nAtomNext = nIndices[oRInext]
+
+        // ignore C-terminal residue as acceptor
+        // if (index.findAtomOnResidue(oRI, 'OXT') !== -1) continue
+
+        position(cAtomPrev, cPosPrev)
+        position(caAtomPrev, caPosPrev)
+        position(nAtomPrev, nPosPrev)
+        position(cAtom, cPos)
+        position(caAtom, caPos)
+        position(nAtom, nPos)
+        position(cAtomNext, cPosNext)
+        position(caAtomNext, caPosNext)
+        position(nAtomNext, nPosNext)
+
+        const tPhi = calculatePhi(cPosPrev, nPos, caPos, cPos)
+        const tPsi = calculatePsi(nPos, caPos, cPos, nPosNext)
+        const tOmega = calculateOmega(caPos, cPos, nPosNext, caPosNext)
+
+        // console.log(`phi:  ${ tPhi }, psi: ${ tPsi }, omega: ${ tOmega }`)
+
+        phi[phi.length] = tPhi
+        psi[psi.length] = tPsi
+        omega[omega.length] = tOmega
+    }
+
+    return { phi, psi, omega };
+}
+
+// angle calculations return NaN upon missing atoms
+
+function calculatePhi(c: Vec3, nNext: Vec3, caNext: Vec3, cNext: Vec3) {
+    return angle(c, nNext, caNext, cNext);
+}
+
+function calculatePsi(n: Vec3, ca: Vec3, c: Vec3, nNext: Vec3) {
+    return angle(n, ca, c, nNext);
+}
+
+function calculateOmega(ca: Vec3, c: Vec3, nNext: Vec3, caNext: Vec3) {
+    return angle(ca, c, nNext, caNext);
+}
+
 function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConformation, proteinResidues: SortedArray<ResidueIndex>, backboneIndices: BackboneAtomIndices, lookup3d: GridLookup3D): DsspHbonds {
     const { cIndices, hIndices, nIndices, oIndices } = backboneIndices
     const { index } = hierarchy
@@ -163,7 +434,7 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf
     const cPosPrev = Vec3.zero()
     const oPosPrev = Vec3.zero()
 
-    const caMinDistSq = caMinDist * caMinDist
+    // const caMinDistSq = caMinDist * caMinDist
 
     for (let i = 0, il = proteinResidues.length; i < il; ++i) {
         const oPI = i
@@ -183,10 +454,10 @@ 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
+            // if (squaredDistances[j] < caMinDistSq) continue
 
             const nPI = indices[j]
 
@@ -221,6 +492,7 @@ function calcBackboneHbonds(hierarchy: AtomicHierarchy, conformation: AtomicConf
             const e = calcHbondEnergy(oPos, cPos, nPos, hPos)
             if (e > hbondEnergyCutoff) continue
 
+            // console.log(`detected hbond between ${ oPI } and ${ nPI } with energy ${ e }`)
             oAtomResidues[oAtomResidues.length] = oPI
             nAtomResidues[nAtomResidues.length] = nPI
             energies[energies.length] = e
@@ -245,8 +517,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,37 +526,47 @@ 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 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, useOriginal = false) {
-    const getResidueFlag = useOriginal ? getOriginalResidueFlag : getUpdatedResidueFlag
+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>
 }
 
@@ -303,6 +585,8 @@ const Q = -27.888
 /** cutoff for hbonds in kcal/mol, must be lower to be consider as an hbond */
 const hbondEnergyCutoff = -0.5
 
+const hbondEnergyMinimal = -9.9
+
 /**
  * E = Q * (1/r(ON) + l/r(CH) - l/r(OH) - l/r(CN))
  */
@@ -312,9 +596,19 @@ function calcHbondEnergy(oPos: Vec3, cPos: Vec3, nPos: Vec3, hPos: Vec3) {
     const distCN = Vec3.distance(cPos, nPos)
     const distON = Vec3.distance(oPos, nPos)
 
+    if (distOH < caMinDist || distCH < caMinDist || distCN < caMinDist || distON < caMinDist) {
+        return hbondEnergyMinimal
+    }
+
     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
 }
 
 /**
@@ -328,28 +622,43 @@ function assignTurns(ctx: DSSPContext) {
     const { proteinResidues, hbonds, flags, hierarchy } = ctx
     const { chains, residueAtomSegments, chainAtomSegments } = hierarchy
     const { label_asym_id } = chains
+    // let turns = 0, turnOpenings = 0
 
-    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]]
+    const turnFlag = [DSSPType.Flag.T3S, DSSPType.Flag.T4S, DSSPType.Flag.T5S, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5]
 
-        // TODO should take sequence gaps into account
-        for (let k = 3; k <= 5; ++k) {
-            if (i + k >= proteinResidues.length) continue
+    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]]
 
-            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
+            // console.log(`${ i } has ${ idx + 3}-turn opening: ${ hbonds.getDirectedEdgeIndex(i, i + idx + 3) !== -1 }`)
+            if (hbonds.getDirectedEdgeIndex(i, i + idx + 3) !== -1) {
+                // console.log(`${ idx + 3 }-turn opening at ${ i } to ${ i + idx + 3 }`)
+                flags[i] |= turnFlag[idx + 3] | turnFlag[idx]
+                // turnOpenings++
+                if (ctx.oldDefinition) {
+                    for (let k = 1; k < idx + 3; ++k) {
+                        flags[i + k] |= turnFlag[idx + 3] | DSSPType.Flag.T
+                        // turns++
+                    }
+                } else {
+                    for (let k = 0; k <= idx + 3; ++k) {
+                        flags[i + k] |= turnFlag[idx + 3] | DSSPType.Flag.T
+                        // turns++
+                    }
+                }
             }
         }
     }
+
+    // console.log(`${ turns } turns, ${ turnOpenings } turn openings`)
 }
 
 /**
@@ -369,7 +678,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 +694,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] = { partner1: Math.min(i, j), partner2: Math.max(i, j), type: BridgeType.PARALLEL }
             }
 
             // Parallel Bridge(i, j) =: [Hbond(j - 1, i) and Hbond(i, j + 1)]
@@ -393,6 +704,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] = { partner1: Math.min(j, i), partner2: Math.max(j, i), type: BridgeType.PARALLEL }
             }
 
             // Antiparallel Bridge(i, j) =: [Hbond(i, j) and Hbond(j, i)]
@@ -401,6 +713,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] = { partner1: Math.min(j, i), partner2: Math.max(j, i), type: BridgeType.ANTI_PARALLEL }
             }
 
             // Antiparallel Bridge(i, j) =: [Hbond(i - 1, j + 1) and Hbond(j - 1, i + l)]
@@ -409,9 +722,14 @@ 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] = { partner1: Math.min(j, i), partner2: Math.max(j, i), type: BridgeType.ANTI_PARALLEL }
             }
         }
     }
+
+    bridges.sort((a, b) => a.partner1 > b.partner1 ? 1 : a.partner1 < b.partner1 ? -1 : 0)
+    // console.log(`${ bridges.length } bridges`)
+    // bridges.forEach(b => console.log(b))
 }
 
 /**
@@ -428,21 +746,77 @@ 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.oldOrdering ? [4, 3, 5] : [3, 4, 5]
+    // const count = [0, 0, 0, 0, 0, 0]
+    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])
+
+            // TODO rework to elegant solution which will not break instantly
+            if (ctx.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
+                    // console.log('skipping in favor of more preferable helix already present')
+                    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
+                    // console.log('skipping in favor of more preferable helix already present')
+                    continue
+                }
+            }
 
-        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]
+            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
+                // console.log(`valid hit for ${n} from ${i} to ${i+n}`)
+                if (ctx.oldDefinition) {
+                    for (let k = 0; k < n; k++) {
+                        flags[i + k] |= helixFlag[n]
+                        // count[n]++
+                    }
+                } else {
+                    for (let k = -1; k <= n; k++) {
+                        flags[i + k] |= helixFlag[n]
+                        // count[n]++
+                    }
                 }
             }
         }
     }
+    // console.log(`${ count[3] + count[4] + count[5] } helix elements - ${ count[3] } 3-10, ${ count[4] } alpha, ${ count[5] } pi`)
+}
+
+function angle(a: Vec3, b: Vec3, c: Vec3, d: Vec3) {
+    const ab = Vec3.zero();
+    const cb = Vec3.zero();
+    const bc = Vec3.zero();
+    const dc = Vec3.zero();
+
+    Vec3.sub(ab, a, b);
+    Vec3.sub(cb, c, b);
+    Vec3.sub(bc, b, c);
+    Vec3.sub(dc, d, c);
+
+    const abc = Vec3.zero();
+    const bcd = Vec3.zero();
+
+    Vec3.cross(abc, ab, cb);
+    Vec3.cross(bcd, bc, dc);
+
+    const angle = Vec3.angle(abc, bcd) * 360.0 / (2 * Math.PI);
+    const cross = Vec3.zero();
+    Vec3.cross(cross, abc, bcd);
+    const value = Vec3.dot(cb, cross);
+
+    return value > 0 ? angle : -angle;
 }
 
 /**
@@ -451,23 +825,146 @@ function assignHelices(ctx: DSSPContext) {
  * Type: E
  */
 function assignLadders(ctx: DSSPContext) {
-    // TODO
+    const { bridges, ladders } = ctx
+
+    // let ext = 0
+    // 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)) {
+                // ext++
+                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
+            }
+        }
+    }
+    // console.log(`${ ext } extension operations`)
+
+    // 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)) {
+                // console.log(`connecting ladders ${ ladderIndex1 } and ${ ladderIndex2 } with inbetween bulge`)
+                ladder1.nextLadder = ladderIndex2
+                ladder2.previousLadder = ladderIndex1
+            }
+        }
+    }
+
+    // console.log(`${ ladders.length } ladders`)
 }
 
 /**
- * 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
+
+    // console.log(`${ bridge.partner1 }-${ bridge.partner2 } ${ ladder.firstEnd }`)
+    // 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) {
+            // console.log(ladder)
+            return true
+        }
+    } else {
+        if (bridge.partner2 === ladder.secondStart - 1) {
+            // console.log(ladder)
+            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,73 @@
+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

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) {

src/tests/browser/render-structure.ts

@@ -65,9 +65,10 @@ async function init() {
     const cif = await downloadFromPdb('3j3q')
     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 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()
 
@@ -76,7 +77,7 @@ async function init() {
         size: reprCtx.sizeThemeRegistry.create('uniform', { structure })
     })
     await cartoonRepr.createOrUpdate({ ...CartoonRepresentationProvider.defaultValues, quality: 'auto' }, structure).run()
-    
+
     canvas3d.add(cartoonRepr)
     canvas3d.resetCamera()
 }