BestDatabaseSequenceMapping superposition

CHANGELOG.md

@@ -12,6 +12,7 @@ Note that since we don't clearly distinguish between a public and private interf
 - Fix #178: ``IndexPairBonds`` for non-single residue structures (bug due to atom reordering).
 - Add volumetric color smoothing for MolecularSurface and GaussianSurface representations (#173)
 - Fix nested 3d grid lookup that caused results being overwritten in non-covalent interactions computation.
+- Basic implementation of ``BestDatabaseSequenceMapping`` (parse from CIF, color theme, superposition).
 
 ## [v2.0.5] - 2021-04-26
 

src/mol-model-props/sequence/themes/best-database-mapping.ts

@@ -21,7 +21,7 @@ const Description = 'Assigns a color based on best dababase sequence mapping.';
 const globalAccessionMap = new Map<string, number>();
 
 export const BestDatabaseSequenceMappingColorThemeParams = {
-    ...getPaletteParams({ type: 'colors', colorList: 'set-3' }),
+    ...getPaletteParams({ type: 'colors', colorList: 'set-1' }),
 };
 export type BestDatabaseSequenceMappingColorThemeParams = typeof BestDatabaseSequenceMappingColorThemeParams
 export function getBestDatabaseSequenceMappingColorThemeParams(ctx: ThemeDataContext) {
@@ -35,7 +35,7 @@ export function BestDatabaseSequenceMappingColorTheme(ctx: ThemeDataContext, pro
             const mapping = BestDatabaseSequenceMapping.Provider.get(m).value;
             if (!mapping) continue;
             for (const acc of mapping.accession) {
-                if (!acc) continue;
+                if (!acc || globalAccessionMap.has(acc)) continue;
                 globalAccessionMap.set(acc, globalAccessionMap.size);
             }
         }

src/mol-model/structure/structure/util/superposition-db-mapping.ts

@@ -0,0 +1,159 @@
+/**
+ * Copyright (c) 2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Segmentation } from '../../../../mol-data/int';
+import { MinimizeRmsd } from '../../../../mol-math/linear-algebra/3d/minimize-rmsd';
+import { BestDatabaseSequenceMapping } from '../../../../mol-model-props/sequence/best-database-mapping';
+import { ElementIndex } from '../../model/indexing';
+import { Structure } from '../structure';
+import { Unit } from '../unit';
+
+export interface AlignmentResult {
+    transform: MinimizeRmsd.Result,
+    pivot: number,
+    other: number
+}
+
+export function alignAndSuperposeWithBestDatabaseMapping(structures: Structure[]): AlignmentResult[] {
+    const indexMap = new Map<string, IndexEntry>();
+
+    for (let i = 0; i < structures.length; i++) {
+        buildIndex(structures[i], indexMap, i);
+    }
+
+    const index = Array.from(indexMap.values());
+
+    // TODO: support non-first structure pivots
+    const pairs = findPairs(structures.length, index);
+
+    const ret: AlignmentResult[] = [];
+    for (const p of pairs) {
+        const [a, b] = getPositionTables(index, p.i, p.j, p.count);
+        const transform = MinimizeRmsd.compute({ a, b });
+        ret.push({ transform, pivot: p.i, other: p.j });
+    }
+
+    return ret;
+}
+
+function getPositionTables(index: IndexEntry[], pivot: number, other: number, N: number) {
+    const xs = MinimizeRmsd.Positions.empty(N);
+    const ys = MinimizeRmsd.Positions.empty(N);
+
+    let o = 0;
+    for (const { pivots } of index) {
+        const a = pivots[pivot];
+        const b = pivots[other];
+        if (!a || !b) continue;
+
+        const l = Math.min(a[2] - a[1], b[2] - b[1]);
+
+        for (let i = 0; i < l; i++) {
+            let eI = (a[1] + i) as ElementIndex;
+            xs.x[o] = a[0].conformation.x(eI);
+            xs.y[o] = a[0].conformation.y(eI);
+            xs.z[o] = a[0].conformation.z(eI);
+
+            eI = (b[1] + i) as ElementIndex;
+            ys.x[o] = b[0].conformation.x(eI);
+            ys.y[o] = b[0].conformation.y(eI);
+            ys.z[o] = b[0].conformation.z(eI);
+            o++;
+        }
+    }
+
+    return [xs, ys];
+}
+
+function findPairs(N: number, index: IndexEntry[]) {
+    const pairwiseCounts: number[][] = [];
+    for (let i = 0; i < N; i++) {
+        pairwiseCounts[i] = [];
+        for (let j = 0; j < N; j++) pairwiseCounts[i][j] = 0;
+    }
+
+    for (const { pivots } of index) {
+        for (let i = 0; i < N; i++) {
+            if (!pivots[i]) continue;
+
+            const lI = pivots[i]![2] - pivots[i]![1];
+
+            for (let j = i + 1; j < N; j++) {
+                if (!pivots[j]) continue;
+
+                const lJ = pivots[j]![2] - pivots[j]![1];
+                pairwiseCounts[i][j] = pairwiseCounts[i][j] + Math.min(lI, lJ);
+            }
+        }
+    }
+
+    const ret: { i: number, j: number, count: number }[] = [];
+
+    for (let j = 1; j < N; j++) {
+        ret[j - 1] = { i: 0, j, count: pairwiseCounts[0][j] };
+    }
+
+    // TODO: support non-first structure pivots
+    // for (let i = 0; i < N - 1; i++) {
+    //     let max = 0, maxJ = i;
+    //     for (let j = i + 1; j < N; j++) {
+    //         if (pairwiseCounts[i][j] > max) {
+    //             maxJ = j;
+    //             max = pairwiseCounts[i][j];
+    //         }
+    //     }
+
+    //     ret[i] = { i, j: maxJ, count: max };
+    // }
+
+    return ret;
+}
+
+interface IndexEntry {
+    key: string,
+    pivots: { [i: number]: [unit: Unit.Atomic, start: ElementIndex, end: ElementIndex] | undefined }
+}
+
+function buildIndex(structure: Structure, index: Map<string, IndexEntry>, sI: number) {
+    for (const unit of structure.units) {
+        if (unit.kind !== Unit.Kind.Atomic) continue;
+
+        const { elements, model } = unit;
+        const { offsets: residueOffset } = model.atomicHierarchy.residueAtomSegments;
+
+        const map = BestDatabaseSequenceMapping.Provider.get(model).value;
+        if (!map) return;
+
+        const { dbName, accession, num } = map;
+
+        const chainsIt = Segmentation.transientSegments(unit.model.atomicHierarchy.chainAtomSegments, elements);
+        const residuesIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueAtomSegments, elements);
+
+        while (chainsIt.hasNext) {
+            const chainSegment = chainsIt.move();
+            residuesIt.setSegment(chainSegment);
+            while (residuesIt.hasNext) {
+                const residueSegment = residuesIt.move();
+                const eI = elements[residueSegment.start];
+                const rI = residueOffset[eI];
+
+                if (!dbName[rI]) continue;
+
+                const key = `${dbName[rI]}-${accession[rI]}-${num[rI]}`;
+
+                if (!index.has(key)) {
+                    index.set(key, { key, pivots: { [sI]: [unit, eI, elements[residueSegment.end]] } });
+                } else {
+                    const entry = index.get(key)!;
+
+                    if (!entry.pivots[sI]) {
+                        entry.pivots[sI] = [unit, eI, elements[residueSegment.end]];
+                    }
+                }
+            }
+        }
+    }
+}

src/mol-plugin-ui/structure/superposition.tsx

@@ -20,6 +20,9 @@ import { ParameterControls } from '../controls/parameters';
 import { stripTags } from '../../mol-util/string';
 import { StructureSelectionHistoryEntry } from '../../mol-plugin-state/manager/structure/selection';
 import { ToggleSelectionModeButton } from './selection';
+import { alignAndSuperposeWithBestDatabaseMapping } from '../../mol-model/structure/structure/util/superposition-db-mapping';
+import { PluginCommands } from '../../mol-plugin/commands';
+import { BestDatabaseSequenceMapping } from '../../mol-model-props/sequence/best-database-mapping';
 
 export class StructureSuperpositionControls extends CollapsableControls {
     defaultState() {
@@ -55,6 +58,7 @@ const SuperpositionTag = 'SuperpositionTransform';
 type SuperpositionControlsState = {
     isBusy: boolean,
     action?: 'byChains' | 'byAtoms' | 'options',
+    canUseDb?: boolean,
     options: StructureSuperpositionOptions
 }
 
@@ -68,9 +72,10 @@ interface AtomsLociEntry extends LociEntry {
     atoms: StructureSelectionHistoryEntry[]
 };
 
-export class SuperpositionControls extends PurePluginUIComponent<{}, SuperpositionControlsState> {
+export class SuperpositionControls extends PurePluginUIComponent<{ }, SuperpositionControlsState> {
     state: SuperpositionControlsState = {
         isBusy: false,
+        canUseDb: false,
         action: undefined,
         options: DefaultStructureSuperpositionOptions
     }
@@ -87,6 +92,10 @@ export class SuperpositionControls extends PurePluginUIComponent<{}, Superpositi
         this.subscribe(this.plugin.behaviors.state.isBusy, v => {
             this.setState({ isBusy: v });
         });
+
+        this.subscribe(this.plugin.managers.structure.hierarchy.behaviors.selection, sel => {
+            this.setState({ canUseDb: sel.structures.every(s => !!s.cell.obj?.data && s.cell.obj.data.models.some(m => BestDatabaseSequenceMapping.Provider.isApplicable(m)) ) });
+        });
     }
 
     get selection() {
@@ -164,6 +173,25 @@ export class SuperpositionControls extends PurePluginUIComponent<{}, Superpositi
         }
     }
 
+    superposeDb = async () => {
+        const input = this.plugin.managers.structure.hierarchy.behaviors.selection.value.structures;
+
+        const transforms = alignAndSuperposeWithBestDatabaseMapping(input.map(s => s.cell.obj?.data!));
+
+        let rmsd = 0;
+
+        for (const xform of transforms) {
+            await this.transform(input[xform.other].cell, xform.transform.bTransform);
+            rmsd += xform.transform.rmsd;
+        }
+
+        rmsd /= transforms.length - 1;
+
+        this.plugin.log.info(`Superposed ${input.length} structures with avg. RMSD ${rmsd.toFixed(2)}.`);
+        await new Promise(res => requestAnimationFrame(res));
+        PluginCommands.Camera.Reset(this.plugin);
+    };
+
     toggleByChains = () => this.setState({ action: this.state.action === 'byChains' ? void 0 : 'byChains' });
     toggleByAtoms = () => this.setState({ action: this.state.action === 'byAtoms' ? void 0 : 'byAtoms' });
     toggleOptions = () => this.setState({ action: this.state.action === 'options' ? void 0 : 'options' });
@@ -293,6 +321,14 @@ export class SuperpositionControls extends PurePluginUIComponent<{}, Superpositi
         </>;
     }
 
+    superposeByDbMapping() {
+        return <>
+            <Button icon={SuperposeChainsSvg} title='Superpose structures using database mapping.' className='msp-btn msp-btn-block' onClick={this.superposeDb} style={{ marginTop: '1px' }} disabled={this.state.isBusy}>
+                DB
+            </Button>
+        </>;
+    }
+
     private setOptions = (values: StructureSuperpositionOptions) => {
         this.setState({ options: values });
     }
@@ -300,8 +336,9 @@ export class SuperpositionControls extends PurePluginUIComponent<{}, Superpositi
     render() {
         return <>
             <div className='msp-flex-row'>
-                <ToggleButton icon={SuperposeChainsSvg} label='By Chains' toggle={this.toggleByChains} isSelected={this.state.action === 'byChains'} disabled={this.state.isBusy} />
-                <ToggleButton icon={SuperposeAtomsSvg} label='By Atoms' toggle={this.toggleByAtoms} isSelected={this.state.action === 'byAtoms'} disabled={this.state.isBusy} />
+                <ToggleButton icon={SuperposeChainsSvg} label='Chains' toggle={this.toggleByChains} isSelected={this.state.action === 'byChains'} disabled={this.state.isBusy} />
+                <ToggleButton icon={SuperposeAtomsSvg} label='Atoms' toggle={this.toggleByAtoms} isSelected={this.state.action === 'byAtoms'} disabled={this.state.isBusy} />
+                {this.state.canUseDb && this.superposeByDbMapping()}
                 <ToggleButton icon={TuneSvg} label='' title='Options' toggle={this.toggleOptions} isSelected={this.state.action === 'options'} disabled={this.state.isBusy} style={{ flex: '0 0 40px', padding: 0 }} />
             </div>
             {this.state.action === 'byChains' && this.addByChains()}