|
@@ -23,12 +23,102 @@ import { AtomicHierarchy, AtomicConformation } from '../atomic';
|
|
|
* - 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
|
|
|
+
|
|
|
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)
|
|
|
}
|
|
|
|
|
|
+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[] },
|
|
|
+ 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 function computeModelDSSP(hierarchy: AtomicHierarchy,
|
|
|
conformation: AtomicConformation,
|
|
|
oldDefinition = true,
|
|
@@ -75,9 +165,7 @@ export function computeModelDSSP(hierarchy: AtomicHierarchy,
|
|
|
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[] = []
|
|
|
|
|
@@ -86,9 +174,8 @@ export function computeModelDSSP(hierarchy: AtomicHierarchy,
|
|
|
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 || // check ought to fail at very start
|
|
|
- flag !== (elements[elements.length - 1] as SecondaryStructure.Helix | SecondaryStructure.Sheet).flags) { // exact flag changed
|
|
|
- elements[elements.length] = createElement(mapToKind(assign), flags[i], getResidueFlag)
|
|
|
+ if (elements.length === 0 /* would fail at very start */ || flag !== (elements[elements.length - 1] as SecondaryStructure.Helix | SecondaryStructure.Sheet).flags /* flag changed */) {
|
|
|
+ elements[elements.length] = createElement(mapToKind(assign), flags[i], getResidueFlag)
|
|
|
}
|
|
|
keys[i] = elements.length - 1
|
|
|
}
|
|
@@ -133,89 +220,15 @@ function createElement(kind: string, flag: DSSPType.Flag, getResidueFlag: (f: DS
|
|
|
}
|
|
|
|
|
|
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
|
|
|
-}
|
|
|
-
|
|
|
-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
|
|
|
+ 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'
|
|
|
}
|
|
|
}
|
|
|
|
|
|
-/** max distance between two C-alpha atoms to check for hbond */
|
|
|
-const caMaxDist = 9.0;
|
|
|
-
|
|
|
function calcAtomicTraceLookup3D(hierarchy: AtomicHierarchy, conformation: AtomicConformation) {
|
|
|
const { x, y, z } = conformation;
|
|
|
const { moleculeType, traceElementIndex } = hierarchy.derived.residue
|
|
@@ -284,7 +297,6 @@ function assignBends(ctx: DSSPContext) {
|
|
|
const caPosPrev2 = Vec3.zero()
|
|
|
const caPos = Vec3.zero()
|
|
|
const caPosNext2 = Vec3.zero()
|
|
|
- // let bends = 0
|
|
|
|
|
|
const nIndices = ctx.backboneIndices.nIndices
|
|
|
const cPos = Vec3.zero()
|
|
@@ -321,16 +333,14 @@ function assignBends(ctx: DSSPContext) {
|
|
|
|
|
|
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 { index } = hierarchy
|
|
|
const { x, y, z } = conformation
|
|
|
const { traceElementIndex } = hierarchy.derived.residue
|
|
|
|
|
@@ -349,7 +359,6 @@ function calculateDihedralAngles(hierarchy: AtomicHierarchy, conformation: Atomi
|
|
|
|
|
|
const phi: number[] = []
|
|
|
const psi: number[] = []
|
|
|
- // const omega: number[] = []
|
|
|
|
|
|
for (let i = 0; i < residueCount - 1; ++i) {
|
|
|
const oPIprev = i - 1
|
|
@@ -370,7 +379,7 @@ function calculateDihedralAngles(hierarchy: AtomicHierarchy, conformation: Atomi
|
|
|
const nAtomNext = nIndices[oRInext]
|
|
|
|
|
|
// ignore C-terminal residue as acceptor
|
|
|
- // if (index.findAtomOnResidue(oRI, 'OXT') !== -1) continue
|
|
|
+ if (index.findAtomOnResidue(oRI, 'OXT') !== -1) continue
|
|
|
|
|
|
position(cAtomPrev, cPosPrev)
|
|
|
position(caAtomPrev, caPosPrev)
|
|
@@ -551,23 +560,6 @@ function getDSSPAssignment(flags: Uint32Array, getResidueFlag: (f: DSSPType) =>
|
|
|
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))
|
|
|
*/
|