Select Git revision
secondary-structure.ts 32.56 KiB
/**
* 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';
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';
import { ParamDefinition as PD } from 'mol-util/param-definition'
/**
* 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,
params: Partial<PD.Values<SecondaryStructureComputationParams>> = {}): SecondaryStructure {
params = { ...PD.getDefaultValues(SecondaryStructureComputationParams), ...params };
// TODO use Zhang-Skolnik for CA alpha only parts or for coarse parts with per-residue 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)
// 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,
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])
// 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: keys,
elements: elements,
dsspString: composeDSSPString(flags, getFlagName)
}
return secondaryStructure
}
function composeDSSPString(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)
}
return 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 if (kind === 'turn' || kind === 'bend') {
return {
kind: 'turn',
flags: getResidueFlag(flag)
}
} else {
return {
kind: 'none'
}
}
}
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'
}
}
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()
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
const { x, y, z } = conformation
const { traceElementIndex } = hierarchy.derived.residue
const residueCount = proteinResidues.length
const position = (i: number, v: Vec3) => Vec3.set(v, x[i], y[i], z[i])
const oAtomResidues: number[] = [];
const nAtomResidues: number[] = [];
const energies: number[] = [];
const oPos = Vec3.zero()
const cPos = Vec3.zero()
const caPos = Vec3.zero()
const nPos = Vec3.zero()
const hPos = Vec3.zero()
const cPosPrev = Vec3.zero()
const oPosPrev = Vec3.zero()
for (let i = 0, il = proteinResidues.length; i < il; ++i) {
const oPI = i
const oRI = proteinResidues[i]
const oAtom = oIndices[oPI]
const cAtom = cIndices[oPI]
const caAtom = traceElementIndex[oRI]
// continue if residue is missing O or C atom
if (oAtom === -1 || cAtom === -1) continue
// ignore C-terminal residue as acceptor
if (index.findAtomOnResidue(oRI, 'OXT') !== -1) continue
position(oAtom, oPos)
position(cAtom, cPos)
position(caAtom, caPos)
const { indices, count } = lookup3d.find(caPos[0], caPos[1], caPos[2], caMaxDist)
for (let j = 0; j < count; ++j) {
const nPI = indices[j]
// ignore bonds within a residue or to prev or next residue, TODO take chain border into account
if (nPI === oPI || nPI - 1 === oPI || nPI + 1 === oPI) continue
const nAtom = nIndices[nPI]
if (nAtom === -1) continue
position(nAtom, nPos)
const hAtom = hIndices[nPI]
if (hAtom === -1) {
// approximate calculation of H position, TODO factor out
if (nPI === 0) continue
const nPIprev = nPI - 1
const oAtomPrev = oIndices[nPIprev]
const cAtomPrev = cIndices[nPIprev]
if (oAtomPrev === -1 || cAtomPrev === -1) continue
position(oAtomPrev, oPosPrev)
position(cAtomPrev, cPosPrev)
Vec3.sub(hPos, cPosPrev, oPosPrev)
const dist = Vec3.distance(oPosPrev, cPosPrev)
Vec3.scaleAndAdd(hPos, nPos, hPos, 1 / dist)
} else {
position(hAtom, hPos)
}
const e = calcHbondEnergy(oPos, cPos, nPos, hPos)
if (e > hbondEnergyCutoff) continue
oAtomResidues[oAtomResidues.length] = oPI
nAtomResidues[nAtomResidues.length] = nPI
energies[energies.length] = e
}
}
return buildHbondGraph(residueCount, oAtomResidues, nAtomResidues, energies);
}
function buildHbondGraph(residueCount: number, oAtomResidues: number[], nAtomResidues: number[], energies: number[]) {
const builder = new IntAdjacencyGraph.DirectedEdgeBuilder(residueCount, oAtomResidues, nAtomResidues);
const _energies = new Float32Array(builder.slotCount);
for (let i = 0, _i = builder.edgeCount; i < _i; i++) {
builder.addNextEdge();
builder.assignProperty(_energies, energies[i]);
}
return builder.createGraph({ energies });
}
/** Original priority: H,B,E,G,I,T,S */
function getOriginalResidueFlag(f: DSSPType) {
if (DSSPType.is(f, DSSPType.Flag.H)) return SecondaryStructureType.SecondaryStructureDssp.H
if (DSSPType.is(f, DSSPType.Flag.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>
}
/**
* E = Q * (1/r(ON) + l/r(CH) - l/r(OH) - l/r(CN))
*/
function calcHbondEnergy(oPos: Vec3, cPos: Vec3, nPos: Vec3, hPos: Vec3) {
const distOH = Vec3.distance(oPos, hPos)
const distCH = Vec3.distance(cPos, hPos)
const distCN = Vec3.distance(cPos, nPos)
const distON = Vec3.distance(oPos, nPos)
const e1 = Q / distOH - Q / distCH
const e2 = Q / distCN - Q / distON
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
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
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
}
}
}
}
}
}
/**
* 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] = new Bridge(i, j, 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] = new Bridge(j, i, 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] = new Bridge(j, i, 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] = new Bridge(j, i, BridgeType.ANTI_PARALLEL)
}
}
}
bridges.sort((a, b) => a.partner1 > b.partner1 ? 1 : a.partner1 < b.partner1 ? -1 : 0)
}
/**
* 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.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])
// 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]
}
}
}
}
}
}
/**
* ladder=: set of one or more consecutive bridges of identical type
*
* Type: E
*/
function assignLadders(ctx: DSSPContext) {
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
}
}
}
}
/**
* 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
// 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)
}
/**
* 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
}
}
}
}