bond-order table for standard residues

src/mol-model-formats/structure/mmcif/bonds/struct_conn.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2018 Mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2019 Mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -11,10 +11,11 @@ import { LinkType } from '../../../../mol-model/structure/model/types'
 import { findEntityIdByAsymId, findAtomIndexByLabelName } from '../util'
 import { Column } from '../../../../mol-data/db'
 import { CustomPropertyDescriptor } from '../../../../mol-model/structure';
-import { mmCIF_Database, mmCIF_Schema } from '../../../../mol-io/reader/cif/schema/mmcif';
+import { mmCIF_Database } from '../../../../mol-io/reader/cif/schema/mmcif';
 import { SortedArray } from '../../../../mol-data/int';
 import { CifWriter } from '../../../../mol-io/writer/cif'
 import { ElementIndex, ResidueIndex } from '../../../../mol-model/structure/model/indexing';
+import { getInterBondOrderFromTable } from '../../../../mol-model/structure/model/properties/atomic/bonds';
 
 export interface StructConn {
     getResidueEntries(residueAIndex: ResidueIndex, residueBIndex: ResidueIndex): ReadonlyArray<StructConn.Entry>,
@@ -211,7 +212,7 @@ export namespace StructConn {
             const partners = _ps(i);
             if (partners.length < 2) continue;
 
-            const type = conn_type_id.value(i) as typeof mmCIF_Schema.struct_conn_type.id.T; // TODO workaround for dictionary inconsistency
+            const type = conn_type_id.value(i)
             const orderType = (pdbx_value_order.value(i) || '').toLowerCase();
             let flags = LinkType.Flag.None;
             let order = 1;
@@ -221,23 +222,24 @@ export namespace StructConn {
                 case 'doub': order = 2; break;
                 case 'trip': order = 3; break;
                 case 'quad': order = 4; break;
+                default:
+                    order = getInterBondOrderFromTable(
+                        struct_conn.ptnr1_label_comp_id.value(i),
+                        struct_conn.ptnr1_label_atom_id.value(i),
+                        struct_conn.ptnr2_label_comp_id.value(i),
+                        struct_conn.ptnr2_label_atom_id.value(i)
+                    )
             }
 
             switch (type) {
                 case 'covale':
-                case 'covale_base':
-                case 'covale_phosphate':
-                case 'covale_sugar':
-                case 'modres':
                     flags = LinkType.Flag.Covalent;
                     break;
                 case 'disulf': flags = LinkType.Flag.Covalent | LinkType.Flag.Sulfide; break;
                 case 'hydrog':
-                case 'mismat':
                     flags = LinkType.Flag.Hydrogen;
                     break;
                 case 'metalc': flags = LinkType.Flag.MetallicCoordination; break;
-                case 'saltbr': flags = LinkType.Flag.Ionic; break;
             }
 
             entries.push({ rowIndex: i, flags, order, distance: pdbx_dist_value.value(i), partners });

src/mol-model/structure/model/properties/atomic/bonds.ts

@@ -0,0 +1,90 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { AminoAcidNames, BaseNames } from '../../types';
+
+/**
+ * Map of intra component bond orders in aminoacids and nucleotides assuming standard IUPAC naming.
+ * The key is constructed as `${compId}|${atomId1}|${atomId2}` with `atomId1 < atomId2`.
+ */
+const IntraBondOrderTable = new Map([
+    ['HIS|CD2|CG', 2],
+    ['HIS|CE1|ND1', 2],
+    ['ARG|CZ|NH2', 2],
+    ['PHE|CE1|CZ', 2],
+    ['PHE|CD2|CE2', 2],
+    ['PHE|CD1|CG', 2],
+    ['TRP|CD1|CG', 2],
+    ['TRP|CD2|CE2', 2],
+    ['TRP|CE3|CZ3', 2],
+    ['TRP|CH2|CZ2', 2],
+    ['ASN|CG|OD1', 2],
+    ['GLN|CD|OE1', 2],
+    ['TYR|CD1|CG', 2],
+    ['TYR|CD2|CE2', 2],
+    ['TYR|CE1|CZ', 2],
+    ['ASP|CG|OD1', 2],
+    ['GLU|CD|OE1', 2],
+
+    ['G|C8|N7', 2],
+    ['G|C4|C5', 2],
+    ['G|C2|N3', 2],
+    ['G|C6|O6', 2],
+    ['C|C4|N3', 2],
+    ['C|C5|C6', 2],
+    ['C|C2|O2', 2],
+    ['A|C2|N3', 2],
+    ['A|C6|N1', 2],
+    ['A|C4|C5', 2],
+    ['A|C8|N7', 2],
+    ['U|C5|C6', 2],
+    ['U|C2|O2', 2],
+    ['U|C4|O4', 2],
+
+    ['DG|C8|N7', 2],
+    ['DG|C4|C5', 2],
+    ['DG|C2|N3', 2],
+    ['DG|C6|O6', 2],
+    ['DC|C4|N3', 2],
+    ['DC|C5|C6', 2],
+    ['DC|C2|O2', 2],
+    ['DA|C2|N3', 2],
+    ['DA|C6|N1', 2],
+    ['DA|C4|C5', 2],
+    ['DA|C8|N7', 2],
+    ['DT|C5|C6', 2],
+    ['DT|C2|O2', 2],
+    ['DT|C4|O4', 2]
+])
+
+/**
+ * Get order for bonds in aminoacids and nucleotides assuming standard IUPAC naming
+ */
+export function getIntraBondOrderFromTable (compId: string, atomId1: string, atomId2: string) {
+    [ atomId1, atomId2 ] = atomId1 < atomId2 ? [ atomId1, atomId2 ] : [ atomId2, atomId1 ]
+    if (AminoAcidNames.has(compId) && atomId1 === 'C' && atomId2 === 'O') return 2
+    if (BaseNames.has(compId) && atomId1 === 'OP1' && atomId2 === 'P') return 2
+    return IntraBondOrderTable.get(`${compId}|${atomId1}|${atomId2}`) || 1
+}
+
+/**
+ * Map of inter component bond orders assuming PDBx/mmCIF naming.
+ * The key is constructed as `${compId1}|${compId2}|${atomId1}|${atomId2}` with `compId1 < compId2`.
+ */
+const InterBondOrderTable = new Map([
+    ['LYS|NZ|RET|C15', 2] // Schiff base in Rhodopsin and Bacteriorhodopsin
+])
+
+/**
+ * Get order for bonds between component assuming PDBx/mmCIF naming.
+ */
+export function getInterBondOrderFromTable (compId1: string, atomId1: string, compId2: string, atomId2: string) {
+    if (compId1 > compId2) {
+        [ compId1, compId2 ] = [ compId2, compId1 ];
+        [ atomId1, atomId2 ] = [ atomId2, atomId1 ];
+    }
+    return InterBondOrderTable.get(`${compId1}|${atomId1}|${compId2}|${atomId2}`) || 1
+}

src/mol-model/structure/structure/unit/links/inter-compute.ts

@@ -16,6 +16,7 @@ import { Vec3, Mat4 } from '../../../../../mol-math/linear-algebra';
 import StructureElement from '../../element';
 import { StructConn } from '../../../../../mol-model-formats/structure/mmcif/bonds';
 import { ElementIndex } from '../../../model/indexing';
+import { getInterBondOrderFromTable } from '../../../../../mol-model/structure/model/properties/atomic/bonds';
 
 const MAX_RADIUS = 4;
 
@@ -52,13 +53,18 @@ function findPairLinks(unitA: Unit.Atomic, unitB: Unit.Atomic, props: LinkComput
     const state: PairState = { mapAB: new Map(), mapBA: new Map(), bondedA: UniqueArray.create(), bondedB: UniqueArray.create() };
     let bondCount = 0;
 
-    const { elements: atomsA } = unitA;
+    const { elements: atomsA, residueIndex: residueIndexA } = unitA;
     const { x: xA, y: yA, z: zA } = unitA.model.atomicConformation;
-    const { elements: atomsB } = unitB;
+    const { elements: atomsB, residueIndex: residueIndexB } = unitB;
     const atomCount = unitA.elements.length;
 
-    const { type_symbol: type_symbolA, label_alt_id: label_alt_idA } = unitA.model.atomicHierarchy.atoms;
-    const { type_symbol: type_symbolB, label_alt_id: label_alt_idB } = unitB.model.atomicHierarchy.atoms;
+    // const { type_symbol: type_symbolA, label_alt_id: label_alt_idA } = unitA.model.atomicHierarchy.atoms;
+    // const { type_symbol: type_symbolB, label_alt_id: label_alt_idB } = unitB.model.atomicHierarchy.atoms;
+    const { type_symbol: type_symbolA, label_alt_id: label_alt_idA, label_atom_id: label_atom_idA } = unitA.model.atomicHierarchy.atoms;
+    const { type_symbol: type_symbolB, label_alt_id: label_alt_idB, label_atom_id: label_atom_idB } = unitB.model.atomicHierarchy.atoms;
+    const { label_comp_id: label_comp_idA } = unitA.model.atomicHierarchy.residues;
+    const { label_comp_id: label_comp_idB } = unitB.model.atomicHierarchy.residues;
+
     const { lookup3d } = unitB;
     const structConn = unitA.model === unitB.model && unitA.model.sourceData.kind === 'mmCIF' ? StructConn.get(unitA.model) : void 0;
 
@@ -104,6 +110,8 @@ function findPairLinks(unitA: Unit.Atomic, unitB: Unit.Atomic, props: LinkComput
         const thresholdA = getElementThreshold(aeI);
         const altA = label_alt_idA.value(aI);
         const metalA = MetalsSet.has(aeI);
+        const atomIdA = label_atom_idA.value(aI);
+        const compIdA = label_comp_idA.value(residueIndexA[aI]);
 
         for (let ni = 0; ni < count; ni++) {
             const _bI = indices[ni];
@@ -123,7 +131,12 @@ function findPairLinks(unitA: Unit.Atomic, unitB: Unit.Atomic, props: LinkComput
 
             if (isHa || isHb) {
                 if (dist < props.maxCovalentHydrogenBondingLength) {
-                    addLink(_aI, _bI, 1, LinkType.Flag.Covalent | LinkType.Flag.Computed, state);
+                    addLink(
+                        _aI, _bI,
+                        1, // covalent bonds involving a hydrogen are always of order 1
+                        LinkType.Flag.Covalent | LinkType.Flag.Computed,
+                        state
+                    );
                     bondCount++;
                 }
                 continue;
@@ -135,7 +148,14 @@ function findPairLinks(unitA: Unit.Atomic, unitB: Unit.Atomic, props: LinkComput
                 : beI < 0 ? thresholdA : Math.max(thresholdA, getElementThreshold(beI));
 
             if (dist <= pairingThreshold) {
-                addLink(_aI, _bI, 1, (isMetal ? LinkType.Flag.MetallicCoordination : LinkType.Flag.Covalent) | LinkType.Flag.Computed, state);
+                const atomIdB = label_atom_idB.value(bI);
+                const compIdB = label_comp_idB.value(residueIndexB[bI]);
+                addLink(
+                    _aI, _bI,
+                    getInterBondOrderFromTable(compIdA, compIdB, atomIdA, atomIdB),
+                    (isMetal ? LinkType.Flag.MetallicCoordination : LinkType.Flag.Covalent) | LinkType.Flag.Computed,
+                    state
+                );
                 bondCount++;
             }
         }

src/mol-model/structure/structure/unit/links/intra-compute.ts

@@ -12,6 +12,7 @@ import { IntAdjacencyGraph } from '../../../../../mol-math/graph';
 import { LinkComputationProps, getElementIdx, MetalsSet, getElementThreshold, isHydrogen, getElementPairThreshold, DefaultLinkComputationProps } from './common';
 import { SortedArray } from '../../../../../mol-data/int';
 import { StructConn, ComponentBond } from '../../../../../mol-model-formats/structure/mmcif/bonds';
+import { getIntraBondOrderFromTable } from '../../../../../mol-model/structure/model/properties/atomic/bonds';
 
 function getGraph(atomA: number[], atomB: number[], _order: number[], _flags: number[], atomCount: number): IntraUnitLinks {
     const builder = new IntAdjacencyGraph.EdgeBuilder(atomCount, atomA, atomB);
@@ -50,20 +51,20 @@ function _computeBonds(unit: Unit.Atomic, props: LinkComputationProps): IntraUni
     for (let _aI = 0; _aI < atomCount; _aI++) {
         const aI =  atoms[_aI];
         const raI = residueIndex[aI];
+        const compId = label_comp_id.value(raI);
 
         if (!props.forceCompute && raI !== lastResidue) {
-            const resn = label_comp_id.value(raI)!;
-            if (!!component && component.entries.has(resn)) {
-                componentMap = component.entries.get(resn)!.map;
+            if (!!component && component.entries.has(compId)) {
+                componentMap = component.entries.get(compId)!.map;
             } else {
                 componentMap = void 0;
             }
         }
         lastResidue = raI;
 
-        const componentPairs = componentMap ? componentMap.get(label_atom_id.value(aI)) : void 0;
-
-        const aeI = getElementIdx(type_symbol.value(aI)!);
+        const aeI = getElementIdx(type_symbol.value(aI));
+        const atomIdA = label_atom_id.value(aI)
+        const componentPairs = componentMap ? componentMap.get(atomIdA) : void 0;
 
         const { indices, count, squaredDistances } = query3d.find(x[aI], y[aI], z[aI], MAX_RADIUS);
         const isHa = isHydrogen(aeI);
@@ -126,7 +127,7 @@ function _computeBonds(unit: Unit.Atomic, props: LinkComputationProps): IntraUni
                 if (dist < props.maxCovalentHydrogenBondingLength) {
                     atomA[atomA.length] = _aI;
                     atomB[atomB.length] = _bI;
-                    order[order.length] = 1;
+                    order[order.length] = 1; // covalent bonds involving hydrogen are always of order 1
                     flags[flags.length] = LinkType.Flag.Covalent | LinkType.Flag.Computed;
                 }
                 continue;
@@ -140,7 +141,7 @@ function _computeBonds(unit: Unit.Atomic, props: LinkComputationProps): IntraUni
             if (dist <= pairingThreshold) {
                 atomA[atomA.length] = _aI;
                 atomB[atomB.length] = _bI;
-                order[order.length] = 1;
+                order[order.length] = getIntraBondOrderFromTable(compId, atomIdA, label_atom_id.value(bI));
                 flags[flags.length] = (isMetal ? LinkType.Flag.MetallicCoordination : LinkType.Flag.Covalent) | LinkType.Flag.Computed;
             }
         }