Merge pull request #196 from molstar/sequence-mapping

Basic sequence mapping support

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/best-database-mapping.ts

@@ -0,0 +1,95 @@
+/**
+ * Copyright (c) 2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Column } from '../../mol-data/db';
+import { MmcifFormat } from '../../mol-model-formats/structure/mmcif';
+import { CustomPropertyDescriptor } from '../../mol-model/custom-property';
+import { Model } from '../../mol-model/structure';
+import { StructureElement } from '../../mol-model/structure/structure';
+import { CustomModelProperty } from '../common/custom-model-property';
+
+export { BestDatabaseSequenceMapping };
+
+interface BestDatabaseSequenceMapping {
+    readonly dbName: string[],
+    readonly accession: string[],
+    readonly num: number[],
+    readonly residue: string[]
+}
+
+namespace BestDatabaseSequenceMapping {
+    export const Provider: CustomModelProperty.Provider<{}, BestDatabaseSequenceMapping> = CustomModelProperty.createProvider({
+        label: 'Best Database Sequence Mapping',
+        descriptor: CustomPropertyDescriptor({
+            name: 'molstar_best_database_sequence_mapping'
+        }),
+        type: 'static',
+        defaultParams: {},
+        getParams: () => ({}),
+        isApplicable: (data: Model) => MmcifFormat.is(data.sourceData) && data.sourceData.data.frame.categories?.atom_site.fieldNames.indexOf('db_name') >= 0,
+        obtain: async (ctx, data) => {
+            return { value: fromCif(data) };
+        }
+    });
+
+    export function getKey(loc: StructureElement.Location) {
+        const model = loc.unit.model;
+        const data = Provider.get(model).value;
+        if (!data) return '';
+        const eI = loc.unit.elements[loc.element];
+        const rI = model.atomicHierarchy.residueAtomSegments.offsets[eI];
+        return data.accession[rI];
+    }
+
+    export function getLabel(loc: StructureElement.Location) {
+        const model = loc.unit.model;
+        const data = Provider.get(model).value;
+        if (!data) return;
+        const eI = loc.unit.elements[loc.element];
+        const rI = model.atomicHierarchy.residueAtomSegments.offsets[eI];
+        const dbName = data.dbName[rI];
+        if (!dbName) return;
+        return `${dbName} ${data.accession[rI]} ${data.num[rI]} ${data.residue[rI]}`;
+    }
+
+    function fromCif(model: Model): BestDatabaseSequenceMapping | undefined {
+        if (!MmcifFormat.is(model.sourceData)) return;
+
+        const { atom_site } = model.sourceData.data.frame.categories;
+        const db_name = atom_site.getField('db_name');
+        const db_acc = atom_site.getField('db_acc');
+        const db_num = atom_site.getField('db_num');
+        const db_res = atom_site.getField('db_res');
+
+        if (!db_name || !db_acc || !db_num || !db_res) return;
+
+        const { atomSourceIndex } = model.atomicHierarchy;
+        const { count, offsets: residueOffsets } = model.atomicHierarchy.residueAtomSegments;
+        const dbName = new Array<string>(count);
+        const accession = new Array<string>(count);
+        const num = new Array<number>(count);
+        const residue = new Array<string>(count);
+
+        for (let i = 0; i < count; i++) {
+            const row = atomSourceIndex.value(residueOffsets[i]);
+
+            if (db_name.valueKind(row) !== Column.ValueKind.Present) {
+                dbName[row] = '';
+                accession[row] = '';
+                num[row] = 0;
+                residue[row] = '';
+                continue;
+            }
+
+            dbName[row] = db_name.str(row);
+            accession[row] = db_acc.str(row);
+            num[row] = db_num.int(row);
+            residue[row] = db_res.str(row);
+        }
+
+        return { dbName, accession, num, residue };
+    }
+}

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

@@ -0,0 +1,105 @@
+/**
+ * Copyright (c) 2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Location } from '../../../mol-model/location';
+import { Bond, StructureElement, Unit } from '../../../mol-model/structure';
+import { ColorTheme, LocationColor } from '../../../mol-theme/color';
+import { ThemeDataContext } from '../../../mol-theme/theme';
+import { Color } from '../../../mol-util/color';
+import { getPalette, getPaletteParams } from '../../../mol-util/color/palette';
+import { ParamDefinition as PD } from '../../../mol-util/param-definition';
+import { CustomProperty } from '../../common/custom-property';
+import { BestDatabaseSequenceMapping } from '../best-database-mapping';
+
+const DefaultColor = Color(0xFAFAFA);
+const Description = 'Assigns a color based on best dababase sequence mapping.';
+
+// same colors for same accessions
+const globalAccessionMap = new Map<string, number>();
+
+export const BestDatabaseSequenceMappingColorThemeParams = {
+    ...getPaletteParams({ type: 'colors', colorList: 'set-1' }),
+};
+export type BestDatabaseSequenceMappingColorThemeParams = typeof BestDatabaseSequenceMappingColorThemeParams
+export function getBestDatabaseSequenceMappingColorThemeParams(ctx: ThemeDataContext) {
+    return BestDatabaseSequenceMappingColorThemeParams; // TODO return copy
+}
+export function BestDatabaseSequenceMappingColorTheme(ctx: ThemeDataContext, props: PD.Values<BestDatabaseSequenceMappingColorThemeParams>): ColorTheme<BestDatabaseSequenceMappingColorThemeParams> {
+    let color: LocationColor;
+
+    if (ctx.structure) {
+        for (const m of ctx.structure.models) {
+            const mapping = BestDatabaseSequenceMapping.Provider.get(m).value;
+            if (!mapping) continue;
+            for (const acc of mapping.accession) {
+                if (!acc || globalAccessionMap.has(acc)) continue;
+                globalAccessionMap.set(acc, globalAccessionMap.size);
+            }
+        }
+
+        const l = StructureElement.Location.create(ctx.structure);
+        const palette = getPalette(globalAccessionMap.size + 1, props, { valueLabel: i => `${i}` });
+        const colorMap = new Map<string, Color>();
+
+        const getColor = (location: StructureElement.Location) => {
+            const key = BestDatabaseSequenceMapping.getKey(location);
+            if (!key) return DefaultColor;
+
+            if (colorMap.has(key)) return colorMap.get(key)!;
+
+            const color = palette.color(globalAccessionMap.get(key)!);
+            colorMap.set(key, color);
+            return color;
+        };
+
+        color = (location: Location): Color => {
+            if (StructureElement.Location.is(location) && Unit.isAtomic(location.unit)) {
+                return getColor(location);
+            } else if (Bond.isLocation(location)) {
+                l.unit = location.aUnit;
+                l.element = location.aUnit.elements[location.aIndex];
+                return getColor(l);
+            }
+            return DefaultColor;
+        };
+    } else {
+        color = () => DefaultColor;
+    }
+
+    return {
+        factory: BestDatabaseSequenceMappingColorTheme,
+        granularity: 'group',
+        preferSmoothing: true,
+        color,
+        props,
+        description: Description,
+    };
+}
+
+export const BestDatabaseSequenceMappingColorThemeProvider: ColorTheme.Provider<BestDatabaseSequenceMappingColorThemeParams, 'best-sequence-database-mapping'> = {
+    name: 'best-sequence-database-mapping',
+    label: 'Best Database Sequence Mapping',
+    category: ColorTheme.Category.Residue,
+    factory: BestDatabaseSequenceMappingColorTheme,
+    getParams: getBestDatabaseSequenceMappingColorThemeParams,
+    defaultValues: PD.getDefaultValues(BestDatabaseSequenceMappingColorThemeParams),
+    isApplicable: (ctx: ThemeDataContext) => !!ctx.structure?.models.some(m => BestDatabaseSequenceMapping.Provider.isApplicable(m)),
+    ensureCustomProperties: {
+        attach: async (ctx: CustomProperty.Context, data: ThemeDataContext) => {
+            if (!data.structure) return;
+
+            for (const m of data.structure.models) {
+                await BestDatabaseSequenceMapping.Provider.attach(ctx, m, void 0, true);
+            }
+        },
+        detach: (data) => {
+            if (!data.structure) return;
+            for (const m of data.structure.models) {
+                m.customProperties.reference(BestDatabaseSequenceMapping.Provider.descriptor, false);
+            }
+        }
+    }
+};

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

src/mol-plugin/behavior/dynamic/custom-props.ts

@@ -11,5 +11,6 @@ export { AccessibleSurfaceArea } from './custom-props/computed/accessible-surfac
 export { Interactions } from './custom-props/computed/interactions';
 export { SecondaryStructure } from './custom-props/computed/secondary-structure';
 export { ValenceModel } from './custom-props/computed/valence-model';
+export { BestDatabaseSequenceMapping } from './custom-props/sequence/best-database-mapping';
 
 export { CrossLinkRestraint } from './custom-props/integrative/cross-link-restraint';

src/mol-plugin/behavior/dynamic/custom-props/sequence/best-database-mapping.ts

@@ -0,0 +1,69 @@
+/**
+ * Copyright (c) 2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { OrderedSet } from '../../../../../mol-data/int';
+import { BestDatabaseSequenceMapping as BestDatabaseSequenceMappingProp } from '../../../../../mol-model-props/sequence/best-database-mapping';
+import { BestDatabaseSequenceMappingColorThemeProvider } from '../../../../../mol-model-props/sequence/themes/best-database-mapping';
+import { Loci } from '../../../../../mol-model/loci';
+import { StructureElement } from '../../../../../mol-model/structure';
+import { ParamDefinition as PD } from '../../../../../mol-util/param-definition';
+import { PluginBehavior } from '../../../behavior';
+
+export const BestDatabaseSequenceMapping = PluginBehavior.create<{ autoAttach: boolean, showTooltip: boolean }>({
+    name: 'best-sequence-database-mapping-prop',
+    category: 'custom-props',
+    display: { name: 'Best Database Sequence Mapping' },
+    ctor: class extends PluginBehavior.Handler<{ autoAttach: boolean, showTooltip: boolean }> {
+        private provider = BestDatabaseSequenceMappingProp.Provider
+
+        private labelProvider = {
+            label: (loci: Loci): string | undefined => {
+                if (!this.params.showTooltip) return;
+                return bestDatabaseSequenceMappingLabel(loci);
+            }
+        }
+
+        update(p: { autoAttach: boolean, showTooltip: boolean }) {
+            let updated = (
+                this.params.autoAttach !== p.autoAttach ||
+                this.params.showTooltip !== p.showTooltip
+            );
+            this.params.autoAttach = p.autoAttach;
+            this.params.showTooltip = p.showTooltip;
+            this.ctx.customStructureProperties.setDefaultAutoAttach(this.provider.descriptor.name, this.params.autoAttach);
+            return updated;
+        }
+
+        register(): void {
+            this.ctx.customModelProperties.register(this.provider, this.params.autoAttach);
+            this.ctx.representation.structure.themes.colorThemeRegistry.add(BestDatabaseSequenceMappingColorThemeProvider);
+            this.ctx.managers.lociLabels.addProvider(this.labelProvider);
+        }
+
+        unregister() {
+            this.ctx.customModelProperties.unregister(this.provider.descriptor.name);
+            this.ctx.representation.structure.themes.colorThemeRegistry.remove(BestDatabaseSequenceMappingColorThemeProvider);
+            this.ctx.managers.lociLabels.removeProvider(this.labelProvider);
+        }
+    },
+    params: () => ({
+        autoAttach: PD.Boolean(true),
+        showTooltip: PD.Boolean(true)
+    })
+});
+
+//
+
+function bestDatabaseSequenceMappingLabel(loci: Loci): string | undefined {
+    if(loci.kind === 'element-loci') {
+        if (loci.elements.length === 0) return;
+
+        const e = loci.elements[0];
+        const u = e.unit;
+        const se = StructureElement.Location.create(loci.structure, u, u.elements[OrderedSet.getAt(e.indices, 0)]);
+        return BestDatabaseSequenceMappingProp.getLabel(se);
+    }
+}

src/mol-plugin/spec.ts

@@ -120,6 +120,7 @@ export const DefaultPluginSpec = (): PluginSpec => ({
 
         PluginSpec.Behavior(PluginBehaviors.CustomProps.StructureInfo),
         PluginSpec.Behavior(PluginBehaviors.CustomProps.AccessibleSurfaceArea),
+        PluginSpec.Behavior(PluginBehaviors.CustomProps.BestDatabaseSequenceMapping),
         PluginSpec.Behavior(PluginBehaviors.CustomProps.Interactions),
         PluginSpec.Behavior(PluginBehaviors.CustomProps.SecondaryStructure),
         PluginSpec.Behavior(PluginBehaviors.CustomProps.ValenceModel),