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

/**
 * 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';

/**
 * 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,
    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)

    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,

        torsionAngles,
        backboneIndices,
        conformation,
        ladders,
        bridges
    }

    assignTurns(ctx)
    assignHelices(ctx)
    assignBends(ctx)
    assignBridges(ctx)
    assignLadders(ctx)
    assignSheets(ctx)

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

    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>
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
    }
}

/** 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 = 0.5;

function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
    const { x, y, z } = conformation;
    const { moleculeType, traceElementIndex } = hierarchy.derived.residue
    const indices: number[] = []
    const _proteinResidues: number[] = []
    for (let i = 0, il = moleculeType.length; i < il; ++i) {
        if (moleculeType[i] === MoleculeType.protein) {
            indices[indices.length] = traceElementIndex[i]
            _proteinResidues[_proteinResidues.length] = i
        }
    }
    const lookup3d = GridLookup3D({ x, y, z, indices: SortedArray.ofSortedArray(indices) }, 4);
    const proteinResidues = SortedArray.ofSortedArray<ResidueIndex>(_proteinResidues)
    return { lookup3d, proteinResidues }
}

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> }>

/**
 * 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
    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()

    // const caMinDistSq = caMinDist * caMinDist

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

            // console.log(`detected hbond between ${ oPI } and ${ nPI } with energy ${ e }`)
            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.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
    if (DSSPType.is(f, DSSPType.Flag.S)) return SecondaryStructureType.SecondaryStructureDssp.S
    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.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 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])
        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

const hbondEnergyMinimal = -9.9

/**
 * 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)

    if (distOH < caMinDist || distCH < caMinDist || distCN < caMinDist || distON < caMinDist) {
        return hbondEnergyMinimal
    }

    const e1 = Q / distOH - Q / distCH
    const e2 = Q / distCN - Q / distON
    const e = e1 + e2

    // cap lowest possible energy
    if (e < hbondEnergyMinimal)
        return hbondEnergyMinimal

    return e
}

/**
 * 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
    // let turns = 0, turnOpenings = 0

    const turnFlag = [DSSPType.Flag.T3S, DSSPType.Flag.T4S, DSSPType.Flag.T5S, DSSPType.Flag.T3, DSSPType.Flag.T4, DSSPType.Flag.T5]

    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
            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
            // 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`)
}

/**
 * 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, bridges } = 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
                // 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)]
            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
                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)]
            i = k
            j = l
            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)]
            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
                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))
}

/**
 * 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 = [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]

    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
                }
            }

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

/**
 * ladder=: set of one or more consecutive bridges of identical type
 *
 * Type: E
 */
function assignLadders(ctx: DSSPContext) {
    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`)
}

/**
 * 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 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)
}

/**
 * sheet=: set of one or more ladders connected by shared residues
 *
 * Type: E
 */
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
            }
        }
    }
}