From 8176052b4b63ee2fc6ddc53c3d14cb62f1cf81f1 Mon Sep 17 00:00:00 2001 From: Alan Cleary Date: Fri, 16 Aug 2019 15:14:59 -0600 Subject: [PATCH] Added selectors for aligning search result micro tracks to the micro tracks reducer, thus deprecating the alignment service. This work included refactoring the Smith-Waterman and Repeat alignment algorithms. Also, a basic workflow that utilizes the new selectors to dynamically generate micro synteny components in GoldenLayout was put in place to help with test driven development of the reducers. --- client/package-lock.json | 13 + client/package.json | 1 + .../src/app/components/gene/gene.component.ts | 73 ++-- .../app/components/gene/micro.component.ts | 36 +- .../app/components/multi/multi.component.ts | 93 +---- .../search/search-params.component.ts | 13 +- .../app/components/search/search.component.ts | 110 +----- .../app/directives/golden-layout.directive.ts | 6 + .../src/app/models/mixins/alignment.model.ts | 1 + client/src/app/services/alignment.service.ts | 126 ------- client/src/app/services/index.ts | 3 - .../src/app/services/micro-tracks.service.ts | 23 +- .../app/store/effects/micro-tracks.effects.ts | 2 +- .../store/effects/pairwise-blocks.effects.ts | 2 +- .../store/reducers/micro-tracks.reducer.ts | 83 ++++- .../src/assets/js/gcv/alignment/alignment.ts | 38 -- client/src/assets/js/gcv/alignment/index.ts | 6 +- client/src/assets/js/gcv/alignment/merge.ts | 126 ------- .../js/gcv/alignment/models/alignment.ts | 4 + .../assets/js/gcv/alignment/models/index.ts | 5 + .../alignment/models/internal-alignment.ts | 4 + .../js/gcv/alignment/models/interval.ts | 1 + .../assets/js/gcv/alignment/models/scores.ts | 6 + .../js/gcv/alignment/models/traceback.ts | 6 + client/src/assets/js/gcv/alignment/repeat.ts | 349 ++++++++++-------- .../assets/js/gcv/alignment/smith-waterman.ts | 237 +++++++----- .../src/assets/js/gcv/alignment/trackify.ts | 92 ----- .../gcv/alignment/utils/alignment-interval.ts | 22 ++ .../js/gcv/alignment/utils/compute-score.ts | 13 + .../assets/js/gcv/alignment/utils/index.ts | 4 + .../gcv/alignment/utils/intervals-to-sets.ts | 62 ++++ .../gcv/alignment/utils/merge-alignments.ts | 167 +++++++++ .../assets/js/gcv/common/filter-by-index.ts | 11 + client/src/assets/js/gcv/common/index.ts | 4 + client/src/assets/js/gcv/common/matrix.ts | 21 ++ client/src/assets/js/gcv/common/set-option.ts | 12 + client/src/assets/js/gcv/common/sum.ts | 8 + client/src/assets/js/gcv/graph/msa-hmm.ts | 5 + 38 files changed, 898 insertions(+), 890 deletions(-) delete mode 100644 client/src/app/services/alignment.service.ts delete mode 100644 client/src/assets/js/gcv/alignment/alignment.ts delete mode 100644 client/src/assets/js/gcv/alignment/merge.ts create mode 100644 client/src/assets/js/gcv/alignment/models/alignment.ts create mode 100644 client/src/assets/js/gcv/alignment/models/index.ts create mode 100644 client/src/assets/js/gcv/alignment/models/internal-alignment.ts create mode 100644 client/src/assets/js/gcv/alignment/models/interval.ts create mode 100644 client/src/assets/js/gcv/alignment/models/scores.ts create mode 100644 client/src/assets/js/gcv/alignment/models/traceback.ts delete mode 100644 client/src/assets/js/gcv/alignment/trackify.ts create mode 100644 client/src/assets/js/gcv/alignment/utils/alignment-interval.ts create mode 100644 client/src/assets/js/gcv/alignment/utils/compute-score.ts create mode 100644 client/src/assets/js/gcv/alignment/utils/index.ts create mode 100644 client/src/assets/js/gcv/alignment/utils/intervals-to-sets.ts create mode 100644 client/src/assets/js/gcv/alignment/utils/merge-alignments.ts create mode 100644 client/src/assets/js/gcv/common/filter-by-index.ts create mode 100644 client/src/assets/js/gcv/common/matrix.ts create mode 100644 client/src/assets/js/gcv/common/set-option.ts create mode 100644 client/src/assets/js/gcv/common/sum.ts diff --git a/client/package-lock.json b/client/package-lock.json index 7f82e0c3..a47b1304 100644 --- a/client/package-lock.json +++ b/client/package-lock.json @@ -7931,6 +7931,14 @@ "minimist": "0.0.8" } }, + "mnemonist": { + "version": "0.30.0", + "resolved": "https://registry.npmjs.org/mnemonist/-/mnemonist-0.30.0.tgz", + "integrity": "sha512-g9rbX4ug7am8oW3jnM7adaFaj5vv5PeOaMPNUwjKQQaTyY+qPQHTmUF59/ZOfQdtTknkjA7+7RhtmL2C0mtwPA==", + "requires": { + "obliterator": "^1.5.0" + } + }, "move-concurrently": { "version": "1.0.1", "resolved": "https://registry.npmjs.org/move-concurrently/-/move-concurrently-1.0.1.tgz", @@ -8499,6 +8507,11 @@ "isobject": "^3.0.1" } }, + "obliterator": { + "version": "1.5.0", + "resolved": "https://registry.npmjs.org/obliterator/-/obliterator-1.5.0.tgz", + "integrity": "sha512-dENe0UviDf8/auXn0bIBKwCcUr49khvSBWDLlszv/ZB2qz1VxWDmkNKFqO2nfmve7hQb/QIDY7+rc7K3LdJimQ==" + }, "obuf": { "version": "1.1.2", "resolved": "https://registry.npmjs.org/obuf/-/obuf-1.1.2.tgz", diff --git a/client/package.json b/client/package.json index 46eb5b1f..4a663b36 100644 --- a/client/package.json +++ b/client/package.json @@ -36,6 +36,7 @@ "d3": "^5.8.2", "golden-layout": "^1.5.9", "jquery": "^3.3.1", + "mnemonist": "^0.30.0", "rxjs": "^6.4.0", "split.js": "^1.5.10", "tippy.js": "^4.0.1", diff --git a/client/src/app/components/gene/gene.component.ts b/client/src/app/components/gene/gene.component.ts index 901726e1..fd6a55ec 100644 --- a/client/src/app/components/gene/gene.component.ts +++ b/client/src/app/components/gene/gene.component.ts @@ -1,7 +1,10 @@ // Angular -import { Component, ViewChild } from "@angular/core"; +import { AfterViewInit, Component, OnDestroy, ViewChild } from "@angular/core"; +import { Subject } from "rxjs"; +import { takeUntil } from "rxjs/operators"; // app import { GoldenLayoutDirective } from "../../directives"; +import { MicroTracksService } from "../../services"; import { MacroComponent } from "./macro.component"; import { MicroComponent } from "./micro.component"; import { PlotComponent } from "./plot.component"; @@ -11,10 +14,12 @@ import { PlotComponent } from "./plot.component"; styleUrls: ["./gene.component.scss"], templateUrl: "./gene.component.html", }) -export class GeneComponent { +export class GeneComponent implements OnDestroy, AfterViewInit { @ViewChild(GoldenLayoutDirective) goldenLayoutDirective; + private _destroy: Subject = new Subject(); + layoutComponents = [ {component: MacroComponent, name: "macro"}, {component: MicroComponent, name: "micro"}, @@ -23,26 +28,50 @@ export class GeneComponent { layoutConfig = { content: [{ type: "column", - content: [ - { - type: "component", - componentName: "macro", - isClosable: false - }, - { - type: "component", - componentName: "micro", - componentState: { - inputs: {key: "value"}, - outputs: { - plot: (() => { - this.goldenLayoutDirective.addItem({type: "component", componentName: "plot"}); - }) - }, - }, - isClosable: false - } - ] + content: [] }] }; + + constructor(private _microTracksService: MicroTracksService) { } + + // Angular hooks + + ngOnDestroy(): void { + this._destroy.next(true); + this._destroy.complete(); + } + + ngAfterViewInit(): void { + this._microTracksService.clusterIDs + .pipe(takeUntil(this._destroy)) + .subscribe((IDs) => { + this._setLayoutConfig(IDs); + }); + } + + // private + + private _setLayoutConfig(clusterIDs): void { + const microConfig = (id) => { + return { + type: "component", + componentName: "micro", + componentState: { + inputs: {cluster: id}, + outputs: { + plot: (() => { + this.goldenLayoutDirective.addItem({type: "component", componentName: "plot"}); + }) + }, + }, + isClosable: false + }; + }; + // TODO: remove existing components before adding new ones + clusterIDs.forEach((id) => { + this.goldenLayoutDirective.addItem(microConfig(id)); + }); + } + + // public } diff --git a/client/src/app/components/gene/micro.component.ts b/client/src/app/components/gene/micro.component.ts index 59bda9cf..49e4e6c1 100644 --- a/client/src/app/components/gene/micro.component.ts +++ b/client/src/app/components/gene/micro.component.ts @@ -1,7 +1,11 @@ // Angular + dependencies -import { AfterViewInit, Component, ElementRef, EventEmitter, OnDestroy, Output, - ViewChild } from "@angular/core"; +import { AfterViewInit, Component, ElementRef, EventEmitter, Input, OnDestroy, + Output, ViewChild } from "@angular/core"; +import { Subject } from "rxjs"; +import { takeUntil } from "rxjs/operators"; +// app import { GCV } from "../../../assets/js/gcv"; +import { MicroTracksService } from "../../services"; @Component({ selector: "micro", @@ -12,7 +16,7 @@ import { GCV } from "../../../assets/js/gcv"; parameters
- micro viewer + micro viewer: {{cluster}}
macro legend @@ -22,20 +26,40 @@ import { GCV } from "../../../assets/js/gcv"; }) export class MicroComponent implements AfterViewInit, OnDestroy { + @Input() cluster: number; + @Input() colors: any; // D3 color function @Output() plot = new EventEmitter(); @ViewChild("container") container: ElementRef; + private _destroy: Subject = new Subject(); + + constructor(private _microTracksService: MicroTracksService) { } + + // Angular hooks + ngAfterViewInit() { - //const viewer = new GCV.visualization.Micro(this.container.nativeElement); - //this.container.nativeElement.innerHTML = "micro-synteny viewer"; + // fetch own data because injected components don't have change detection + this._microTracksService.getCluster(this.cluster) + .pipe(takeUntil(this._destroy)) + .subscribe(this._draw); } ngOnDestroy() { - //console.log('destroyed'); + this._destroy.next(true); + this._destroy.complete(); } + // public + spawnPlot() { this.plot.emit(); } + + // private + + _draw(): void { + //const viewer = new GCV.visualization.Micro(this.container.nativeElement); + //this.container.nativeElement.innerHTML = "micro-synteny viewer"; + } } diff --git a/client/src/app/components/multi/multi.component.ts b/client/src/app/components/multi/multi.component.ts index 00519d92..f40d0d8b 100644 --- a/client/src/app/components/multi/multi.component.ts +++ b/client/src/app/components/multi/multi.component.ts @@ -12,7 +12,7 @@ import { AppConfig } from "../../app.config"; import { Alert, Family, Gene, Group, MacroTracks, MicroTracks } from "../../models"; import { ClusterMixin, DrawableMixin, PointMixin } from "../../models/mixins"; import { microTracksOperator, multiMacroTracksOperator } from "../../operators"; -import { AlignmentService, FilterService, InterAppCommunicationService, +import { FilterService, InterAppCommunicationService, MacroTracksService, MicroTracksService } from "../../services"; import { AlertComponent } from "../shared/alert.component"; @@ -79,8 +79,7 @@ export class MultiComponent implements AfterViewInit, OnDestroy, OnInit { // track which families have been checked private familyGlyphsSubject = new BehaviorSubject<{[key: string]: string}>({}); - constructor(private alignmentService: AlignmentService, - private resolver: ComponentFactoryResolver, + constructor(private resolver: ComponentFactoryResolver, private filterService: FilterService, private communicationService: InterAppCommunicationService, private macroTracksService: MacroTracksService, @@ -158,94 +157,6 @@ export class MultiComponent implements AfterViewInit, OnDestroy, OnInit { return glyphs; }), takeUntil(this.destroy)); - - const glyphedAlignedTracks: Observable> = combineLatest( - this.alignmentService.alignedMicroTracks, - this.familyGlyphs) - .pipe( - map(([tracks, glyphs]) => { - return { - // ensure the orphan family is present - families: [...tracks.families, {id: "", name: ""}].map((f) => { - if (glyphs.hasOwnProperty(f.id)) { - const family = Object.create(f); - family.glyph = glyphs[f.id]; - return family; - } - return f; - }), - groups: tracks.groups.map((t) => { - const track = Object.create(t); - track.genes = t.genes.map((g) => { - if (glyphs.hasOwnProperty(g.family)) { - const gene = Object.create(g); - gene.glyph = glyphs[g.family]; - return gene; - } - return g; - }); - return track; - }), - }; - })); - - this.alignmentService.alignedMicroTracks - .pipe(takeUntil(this.destroy)) - .subscribe((tracks) => { - this.hideLeftSlider(); - this._setFamilyGlyph(null, null); - }); - - glyphedAlignedTracks - .pipe( - withLatestFrom(this.microTracksService.routeParams), - takeUntil(this.destroy)) - .subscribe(([tracks, route]) => { - this._onAlignedMicroTracks(tracks, route); - }); - - const filteredMicroTracks = - combineLatest( - glyphedAlignedTracks, - this.filterService.regexpAlgorithm, - this.filterService.orderAlgorithm) - .pipe(microTracksOperator({prefix: (t) => "group " + t.cluster + " - "})); - - filteredMicroTracks - .pipe(takeUntil(this.destroy)) - .subscribe((tracks) => { - // add gene tooltips - tracks.groups.forEach((group) => { - group.genes = group.genes.map((gene) => { - const g = Object.create(gene); - g.htmlAttributes = { - "data-tippy-content": function (g) { - return ` -
-
-
${g.name} (${g.family})
- ${g.fmin} - ${g.fmax} -
-
- `; - } - }; - return g; - }); - }); - this.microTracks = tracks as MicroTracks; - }); - - // subscribe to macro track data - combineLatest( - this.macroTracksService.multiMacroTracks - .pipe(filter((tracks) => tracks !== undefined)), - filteredMicroTracks) - .pipe(multiMacroTracksOperator()) - .pipe(takeUntil(this.destroy)) - .subscribe((tracks) => { - this._onMacroTracks(tracks); - }); } // public diff --git a/client/src/app/components/search/search-params.component.ts b/client/src/app/components/search/search-params.component.ts index ed5aba66..3165e8d0 100644 --- a/client/src/app/components/search/search-params.component.ts +++ b/client/src/app/components/search/search-params.component.ts @@ -7,7 +7,7 @@ import { takeUntil } from "rxjs/operators"; import { AppConfig } from "../../app.config"; import { ALIGNMENT_ALGORITHMS } from "../../algorithms"; import { AlignmentParams, BlockParams, QueryParams } from "../../models"; -import { AlignmentService, MacroTracksService, MicroTracksService } from "../../services"; +import { MacroTracksService, MicroTracksService } from "../../services"; @Component({ selector: "search-params", @@ -38,8 +38,7 @@ export class SearchParamsComponent implements OnDestroy, OnInit { private destroy: Subject; // constructor - constructor(private alignmentService: AlignmentService, - private fb: FormBuilder, + constructor(private fb: FormBuilder, private macroTracksService: MacroTracksService, private microTracksService: MicroTracksService) { this.destroy = new Subject(); @@ -65,9 +64,9 @@ export class SearchParamsComponent implements OnDestroy, OnInit { // initialize alignment group and subscribe to store updates const defaultAlignment = new AlignmentParams(); this.alignmentGroup = this.fb.group(defaultAlignment.formControls()); - this.alignmentService.alignmentParams - .pipe(takeUntil(this.destroy)) - .subscribe((params) => this._updateGroup(this.alignmentGroup, params)); + //this.alignmentService.alignmentParams + // .pipe(takeUntil(this.destroy)) + // .subscribe((params) => this._updateGroup(this.alignmentGroup, params)); // submit the updated form this.blockGroup.markAsDirty(); @@ -95,7 +94,7 @@ export class SearchParamsComponent implements OnDestroy, OnInit { }); // submit alignment params this._submitGroup(this.alignmentGroup, (params) => { - this.alignmentService.updateParams(params); + //this.alignmentService.updateParams(params); }); } else { this.invalid.emit(); diff --git a/client/src/app/components/search/search.component.ts b/client/src/app/components/search/search.component.ts index 3fe95fba..24e5ef83 100644 --- a/client/src/app/components/search/search.component.ts +++ b/client/src/app/components/search/search.component.ts @@ -12,7 +12,7 @@ import { AppConfig } from "../../app.config"; import { Alert, Family, Gene, Group, MacroTracks, MicroTracks } from "../../models"; import { DrawableMixin } from "../../models/mixins"; import { macroTracksOperator, microTracksOperator, plotsOperator } from "../../operators"; -import { AlignmentService, FilterService, InterAppCommunicationService, +import { FilterService, InterAppCommunicationService, MacroTracksService, MicroTracksService, PlotsService } from "../../services"; import { AlertComponent } from "../shared/alert.component"; import { PlotViewerComponent } from "../viewers/plot.component"; @@ -106,8 +106,7 @@ export class SearchComponent implements AfterViewInit, OnDestroy, OnInit { private macroTrackObservable: Observable<[MacroTracks, any]>; // Observable<[MacroTracks, MicroTracks]>; private macroSub: any; - constructor(private alignmentService: AlignmentService, - private resolver: ComponentFactoryResolver, + constructor(private resolver: ComponentFactoryResolver, private filterService: FilterService, private communicationService: InterAppCommunicationService, private macroTracksService: MacroTracksService, @@ -195,117 +194,12 @@ export class SearchComponent implements AfterViewInit, OnDestroy, OnInit { }), takeUntil(this.destroy)); - const glyphedAlignedTracks = combineLatest( - this.alignmentService.alignedMicroTracks, - this.familyGlyphs) - .pipe( - map(([tracks, glyphs]) => { - return { - // ensure the orphan family is present - families: [...tracks.families, {id: "", name: ""}].map((f) => { - if (glyphs.hasOwnProperty(f.id)) { - const family = Object.create(f); - family.glyph = glyphs[f.id]; - return family; - } - return f; - }), - groups: tracks.groups.map((t) => { - const track = Object.create(t); - track.genes = t.genes.map((g) => { - if (glyphs.hasOwnProperty(g.family)) { - const gene = Object.create(g); - gene.glyph = glyphs[g.family]; - return gene; - } - return g; - }); - return track; - }), - }; - })); - this.microTracksService.microTracks .pipe(takeUntil(this.destroy)) .subscribe((tracks) => { this.hideLeftSlider(); this._setFamilyGlyph(null, null); }); - - glyphedAlignedTracks - .pipe( - withLatestFrom( - this.microTracksService.routeParams, - this.microTracksService.microTracks - .pipe(map((tracks) => tracks.groups.length))), - takeUntil(this.destroy)) - .subscribe(([tracks, route, numReturned]) => { - this._onAlignedMicroTracks(tracks as MicroTracks, route, numReturned); - }); - - const filteredMicroTracks = - combineLatest( - glyphedAlignedTracks, - this.filterService.regexpAlgorithm, - this.filterService.orderAlgorithm) - .pipe(microTracksOperator({skipFirst: true})); - - filteredMicroTracks - .pipe(takeUntil(this.destroy)) - .subscribe((tracks) => { - // add gene tooltips - tracks.groups.forEach((group) => { - group.genes = group.genes.map((gene) => { - const g = Object.create(gene); - g.htmlAttributes = { - "data-tippy-content": function (g) { - return ` -
-
-
${g.name} (${g.family})
- ${g.fmin} - ${g.fmax} -
-
- `; - } - }; - return g; - }); - }); - this.microTracks = tracks as MicroTracks; - }); - - // subscribe to macro track data - combineLatest( - this.macroTracksService.macroTracks, - filteredMicroTracks, - this.macroConfigObservable) - .pipe(macroTracksOperator()) - .pipe( - withLatestFrom(this.microTracksService.routeParams), - filter(([tracks, route]) => route.gene !== undefined), - takeUntil(this.destroy)) - .subscribe(([tracks, route]) => this._onMacroTracks(tracks)); - - // subscribe to micro-plots changes - combineLatest(this.plotsService.localPlots, filteredMicroTracks) - .pipe(takeUntil(this.destroy)) - .pipe(plotsOperator()) - .subscribe((plots) => this.microPlots = plots); - this.plotsService.selectedLocalPlot - .pipe( - filter((plot) => plot !== null), - takeUntil(this.destroy)) - .subscribe((plot) => { - this.selectedLocalPlot = plot; - }); - this.plotsService.selectedGlobalPlot - .pipe( - filter((plot) => plot !== null), - takeUntil(this.destroy)) - .subscribe((plot) => { - this.selectedGlobalPlot = plot; - }); } // public diff --git a/client/src/app/directives/golden-layout.directive.ts b/client/src/app/directives/golden-layout.directive.ts index a040e0f5..dcd1d880 100644 --- a/client/src/app/directives/golden-layout.directive.ts +++ b/client/src/app/directives/golden-layout.directive.ts @@ -25,6 +25,8 @@ export class GoldenLayoutDirective implements AfterContentInit { private _el: ElementRef, private _injector: Injector) { } + // Angular hooks + ngAfterContentInit() { // set the initial layout configuration if (this.config === undefined) { @@ -67,6 +69,8 @@ export class GoldenLayoutDirective implements AfterContentInit { } + // private + private _createComponent(component, $container, inputs, outputs): ComponentRef { const factory = this._componentFactoryResolver .resolveComponentFactory(component); @@ -93,6 +97,8 @@ export class GoldenLayoutDirective implements AfterContentInit { componentRef.destroy(); } + // public + addItem(itemConfig) { this._layout.root.contentItems[0].addChild(itemConfig); } diff --git a/client/src/app/models/mixins/alignment.model.ts b/client/src/app/models/mixins/alignment.model.ts index 67de2c35..b83ee082 100644 --- a/client/src/app/models/mixins/alignment.model.ts +++ b/client/src/app/models/mixins/alignment.model.ts @@ -1,3 +1,4 @@ export interface AlignmentMixin { + alignment: number[]; score: number; } diff --git a/client/src/app/services/alignment.service.ts b/client/src/app/services/alignment.service.ts deleted file mode 100644 index 009cae0e..00000000 --- a/client/src/app/services/alignment.service.ts +++ /dev/null @@ -1,126 +0,0 @@ -// Angular -import { Injectable } from "@angular/core"; -import { Observable } from "rxjs"; -// store -import { Store } from "@ngrx/store"; -import * as routerActions from "../store/actions/router.actions"; -import * as fromRoot from "../store/reducers"; -//import * as fromAlignedMicroTracks from "../store/reducers/aligned-micro-tracks.store"; -import * as fromRouter from "../store/reducers/router.reducer"; -// app -import { GCV } from "../../assets/js/gcv"; -import { ALIGNMENT_ALGORITHMS } from "../algorithms"; -import { AlignmentParams, Group, MicroTracks } from "../models"; -import { AlignmentMixin, ClusterMixin } from "../models/mixins"; - -@Injectable() -export class AlignmentService { - alignedMicroTracks: Observable; - alignmentParams: Observable; - - private algorithms = ALIGNMENT_ALGORITHMS; - private algorithmIDs = this.algorithms.map((a) => a.id); - - constructor(private store: Store) { - // initialize observables - //this.alignedMicroTracks = this.store.select(fromAlignedMicroTracks.getAlignedMicroTracks); - this.alignmentParams = this.store.select(fromRouter.getMicroAlignmentParams); - } - - updateParams(params: AlignmentParams): void { - const path = []; - const query = Object.assign({}, params); - this.store.dispatch(new routerActions.Go({path, query})); - } - - getPairwiseReference(track: Group): Observable { - return Observable.create((observer) => { - // create modifiable extension of the query track - const reference = Object.assign({}, track); - reference.genes = []; - for (let i = 0; i < track.genes.length; i++) { - const g = Object.create(track.genes[i]); - g.x = i; - g.y = 0; - reference.genes.push(g); - } - // push the "aligned" query track to the store; - observer.next(reference); - observer.complete(); - }); - } - - getPairwiseAlignment( - query: Group, - tracks: MicroTracks, - params: AlignmentParams - ): Observable { - return Observable.create((observer) => { - // create modifiable extension of tracks - let alignedTracks = new MicroTracks() as MicroTracks<{}, AlignmentMixin, AlignmentMixin>; - alignedTracks.families = tracks.families; - for (const group of tracks.groups) { - alignedTracks.groups.push({ - ...group, - genes: group.genes.map(g => Object.create(g)), - } as Group & AlignmentMixin); - } - // get the alignment algorithm - const algorithmID = this.algorithmIDs.indexOf(params.algorithm); - const algorithm = this.algorithms[algorithmID].algorithm; - // get the algorithm options - const options = Object.assign({}, { - accessor: (g) => { - if (g.reversed) { - return (g.strand === -1 ? "+" : "-") + g.family; - } - return (g.strand === -1 ? "-" : "+") + g.family; - }, - reverse: (s) => { - const r = []; - for (const g of s) { - const g2 = Object.create(Object.getPrototypeOf(g)); - g2.reversed = true; - r.push(g2); - } - r.reverse(); - return r; - }, - scores: Object.assign({}, params), - suffixScores: true, - ignore: ["+", "-"], - }); - // perform the alginments - const alignments = []; - for (let i = 0; i < alignedTracks.groups.length; ++i) { - const result = alignedTracks.groups[i]; - const al = algorithm(query.genes, result.genes, options); - // save tracks that meet the threshold - for (const a of al) { - if (a.score >= options.scores.threshold) { - a.track = Object.assign({}, result); - alignments.push(a); - } - } - } - // convert the alignments into tracks - alignedTracks = GCV.alignment.trackify(query, alignedTracks, alignments); - // merge tracks from same alignment set and remove residual suffix scores - alignedTracks = GCV.alignment.merge(alignedTracks); - for (let i = 1; i < alignedTracks.groups.length; ++i) { - const genes = alignedTracks.groups[i].genes; - for (const g of genes) { - delete g.score; - } - } - // TODO: move to standalone filter - alignedTracks.groups = alignedTracks.groups.filter((g, i) => { - return i === 0 || g.score >= params.score; - }); - // push the aligned tracks to the store; - observer.next(alignedTracks); - // complete the observable - observer.complete(); - }); - } -} diff --git a/client/src/app/services/index.ts b/client/src/app/services/index.ts index faf9a1dc..7e15c08f 100644 --- a/client/src/app/services/index.ts +++ b/client/src/app/services/index.ts @@ -1,4 +1,3 @@ -import { AlignmentService } from "./alignment.service"; import { ChromosomeService } from "./chromosome.service"; import { ClusteringService } from "./clustering.service"; import { DetailsService } from "./details.service"; @@ -12,7 +11,6 @@ import { PlotsService } from "./plots.service"; import { TourService } from "./tour.service"; export const services: any[] = [ - AlignmentService, ChromosomeService, ClusteringService, DetailsService, @@ -26,7 +24,6 @@ export const services: any[] = [ TourService, ]; -export * from "./alignment.service"; export * from "./chromosome.service"; export * from "./clustering.service"; export * from "./details.service"; diff --git a/client/src/app/services/micro-tracks.service.ts b/client/src/app/services/micro-tracks.service.ts index 29d4d379..65cfa320 100644 --- a/client/src/app/services/micro-tracks.service.ts +++ b/client/src/app/services/micro-tracks.service.ts @@ -5,16 +5,17 @@ import { Observable, combineLatest, empty, merge, onErrorResumeNext, throwError } from "rxjs"; import { catchError, filter, map, take } from "rxjs/operators"; // store -import { Store } from "@ngrx/store"; +import { Store, select } from "@ngrx/store"; import * as routerActions from "../store/actions/router.actions"; import * as fromRoot from "../store/reducers"; +import * as fromMicroTracks from "../store/reducers/micro-tracks.reducer"; //import * as fromMacroChromosome from "../store/reducers/macro-chromosome.store"; import * as fromRouter from "../store/reducers/router.reducer"; //import * as fromSearchQueryTrack from "../store/reducers/search-query-track.store"; // app import { AppConfig } from "../app.config"; import { Group, MicroTracks, QueryParams, Track } from "../models"; -import { PointMixin } from "../models/mixins"; +import { AlignmentMixin, ClusterMixin, PointMixin } from "../models/mixins"; import { HttpService } from "./http.service"; @Injectable() @@ -24,8 +25,11 @@ export class MicroTracksService extends HttpService { routeParams: Observable; searchQueryTrack: Observable; + clusterIDs: Observable; + constructor(private _http: HttpClient, private store: Store) { super(_http); + this.clusterIDs = store.select(fromMicroTracks.getClusterIDs); // initialize observables //this.microTracks = store.select(fromMicroTracks.getMicroTracks); this.queryParams = store.select(fromRouter.getMicroQueryParams); @@ -58,9 +62,9 @@ export class MicroTracksService extends HttpService { }; return this._makeRequest>(serverID, "microSearch", body).pipe( map((tracks) => { - this._removeQuery(query, tracks); - this._mergeOverlappingTracks(tracks); - this._parseTracks(serverID, tracks); + //this._removeQuery(query, tracks); + //this._mergeOverlappingTracks(tracks); + //this._parseTracks(serverID, tracks); return tracks; }), catchError((error) => throwError(error)), @@ -108,6 +112,15 @@ export class MicroTracksService extends HttpService { this.store.dispatch(new routerActions.Go({path: [url, { routeParam: 1 }]})); } + // returns all the aligned micro-tracks (selected and search result) belonging + // to a specific cluster + getCluster(id: number): Observable<(Track | ClusterMixin | AlignmentMixin)[]> + { + return this.store.pipe( + select(fromMicroTracks.getAlignedMicroTrackCluster(id)) + ); + } + // adds the server id the track came from to the track and its genes private _parseTrack(source: string, track: Group, i: number = -1): void { track.source = source; diff --git a/client/src/app/store/effects/micro-tracks.effects.ts b/client/src/app/store/effects/micro-tracks.effects.ts index 8c80daab..71240e9a 100644 --- a/client/src/app/store/effects/micro-tracks.effects.ts +++ b/client/src/app/store/effects/micro-tracks.effects.ts @@ -47,7 +47,7 @@ export class MicroTracksEffects { // TODO: update so it only gets tracks that haven't been fetched already, e.g. // if a source is (de)selected @Effect() - consensusSearch = combineLatest( + consensusSearch$ = combineLatest( this.store.select( fromMicroTracks.getClusteredAndAlignedSelectedMicroTracks), this.store.select(fromRouter.getMicroQueryParamSources) diff --git a/client/src/app/store/effects/pairwise-blocks.effects.ts b/client/src/app/store/effects/pairwise-blocks.effects.ts index f9240556..4d7d4a34 100644 --- a/client/src/app/store/effects/pairwise-blocks.effects.ts +++ b/client/src/app/store/effects/pairwise-blocks.effects.ts @@ -43,7 +43,7 @@ export class PairwiseBlocksEffects { // gets blocks for selected chromosomes from selected sources @Effect() - getSelected = combineLatest( + getSelected$ = combineLatest( this.store.select(fromChromosome.getSelectedChromosomes), this.store.select(fromRouter.getMicroQueryParamSources), this.store.select( diff --git a/client/src/app/store/reducers/micro-tracks.reducer.ts b/client/src/app/store/reducers/micro-tracks.reducer.ts index 9be8b6fb..c15e2e99 100644 --- a/client/src/app/store/reducers/micro-tracks.reducer.ts +++ b/client/src/app/store/reducers/micro-tracks.reducer.ts @@ -31,8 +31,9 @@ import * as microTrackActions from "../actions/micro-tracks.actions"; // app import * as clusterfck from "../../../assets/js/clusterfck"; import { GCV } from "../../../assets/js/gcv"; -import { ClusteringParams, Gene, Track } from "../../models"; -import { ClusterMixin } from "../../models/mixins"; +import { ALIGNMENT_ALGORITHMS } from "../../algorithms"; +import { AlignmentParams, ClusteringParams, Gene, Track } from "../../models"; +import { AlignmentMixin, ClusterMixin } from "../../models/mixins"; declare var Object: any; // because TypeScript doesn't support Object.values @@ -242,11 +243,24 @@ export const getClusteredSelectedMicroTracks = createSelector( } ); +export const getClusterIDs = createSelector( + getClusteredSelectedMicroTracks, + (tracks: (Track | ClusterMixin)[]): number[] => { + const IDs = tracks.map((t: ClusterMixin) => t.cluster); + const uniqueIDs = new Set(IDs); + return Array.from(uniqueIDs); + }, +); + // performs a multiple sequence alignment of the tracks in each cluster and // outputs the aligned tracks and their consensus sequence export const getClusteredAndAlignedSelectedMicroTracks = createSelector( getClusteredSelectedMicroTracks, - (tracks: (Track | ClusterMixin)[]) => { + (tracks: (Track | ClusterMixin)[]): + { + consensuses: string[][], + tracks: (Track | ClusterMixin | AlignmentMixin)[] + } => { // group tracks by cluster const clusterer = (accumulator, track) => { if (!(track.cluster in accumulator)) { @@ -289,3 +303,66 @@ export const getClusteredAndAlignedSelectedMicroTracks = createSelector( return clusteredAlignments; }, ); + +export const getClusteredAndAlignedSearchMicroTracks = createSelector( + getMicroTracksState, + getClusteredAndAlignedSelectedMicroTracks, + fromRouter.getMicroAlignmentParams, + (state: State, {consensuses, tracks}, params: AlignmentParams): + (Track | ClusterMixin | AlignmentMixin)[] => { + // get selected alignment algorithm + const algorithmIDs = ALIGNMENT_ALGORITHMS.map((a) => a.id); + const algorithmID = algorithmIDs.indexOf(params.algorithm); + const algorithm = ALIGNMENT_ALGORITHMS[algorithmID].algorithm; + // creates an alignment mixin for the given track + const mixin = (track, alignment) => { + const t = Object.create(track); + t.alignment = alignment.alignment; + t.score = alignment.score; + return t; + }; + // returns one or more alignments (depending on the alignment algorithm) for + // the given sequence relative to the given reference + const aligner = (reference, sequence) => { + const options = { + omit: new Set(""), + scores: Object.assign({}, params), + }; + const alignments = algorithm(reference, sequence); + // TODO: combine repeat tracks + return alignments; + }; + // aligns each track to its cluster's consensus sequence, creates an + // alignment mixin for each track's alignments, and return a flattened array + const reducer = (accumulator, track) => { + const cluster = (track as ClusterMixin).cluster; + const consensus = consensuses[cluster]; + const families = (track as Track).families; + const alignments = aligner(consensus, families); + const trackAlignments = alignments.map((a) => mixin(track, a)); + accumulator.push(...trackAlignments); + return accumulator; + }; + // align the tracks + const searchTracks = Object.values(state.entities); + const alignedTracks = searchTracks.reduce(reducer, []); + return alignedTracks; + }, +); + +export const getAllAlignedAndClusteredMicroTrackCluster = createSelector( + getClusteredAndAlignedSelectedMicroTracks, + getClusteredAndAlignedSearchMicroTracks, + ({consensuses, tracks}, searchTracks): + (Track | ClusterMixin | AlignmentMixin)[] => { + return tracks.concat(searchTracks); + }, +); + +export const getAlignedMicroTrackCluster = (id: number) => createSelector( + getClusteredAndAlignedSelectedMicroTracks, + ({consensuses, tracks}) => { + const cluster = tracks.filter((t: ClusterMixin) => t.cluster === id); + return cluster; + }, +); diff --git a/client/src/assets/js/gcv/alignment/alignment.ts b/client/src/assets/js/gcv/alignment/alignment.ts deleted file mode 100644 index f46650e6..00000000 --- a/client/src/assets/js/gcv/alignment/alignment.ts +++ /dev/null @@ -1,38 +0,0 @@ -/** - * Uses an accessor to compute a score by comparing the given elements. - * @param {generic} a - The first element. - * @param {generic} b - The second element. - * @param {function} accessor - The accessor used to compare elements. - * @param {object} scores - The contains the match and mismatch scores. - * @return {int} - The computed score. - */ -export function computeScore(a, b, accessor, scores, ignore) { - a = accessor(a); - b = accessor(b); - if (a === b && ignore.indexOf(a) === -1) { - return scores.match; - } - return scores.mismatch; -} - -/** - * Creates a matrix with the given initial value. - * @param {number} m - Number of columns. - * @param {number} n - Number of rows. - * @param {number} initial - The initial value to put in each cell. - * @return {Array} - An m x n array with initial value in each cell. - */ -export function matrix(m, n, initial) { - let a; - let i; - let j; - const mat = []; - for (i = 0; i < m; i += 1) { - a = []; - for (j = 0; j < n; j += 1) { - a[j] = initial; - } - mat[i] = a; - } - return mat; -} diff --git a/client/src/assets/js/gcv/alignment/index.ts b/client/src/assets/js/gcv/alignment/index.ts index bf3fc63e..8d7ed9b1 100644 --- a/client/src/assets/js/gcv/alignment/index.ts +++ b/client/src/assets/js/gcv/alignment/index.ts @@ -1,4 +1,2 @@ -export { merge } from "./merge"; -export { default as repeat } from "./repeat"; -export { default as smithWaterman } from "./smith-waterman"; -export { default as trackify } from "./trackify"; +export { repeat } from "./repeat"; +export { smithWaterman } from "./smith-waterman"; diff --git a/client/src/assets/js/gcv/alignment/merge.ts b/client/src/assets/js/gcv/alignment/merge.ts deleted file mode 100644 index f4b510f1..00000000 --- a/client/src/assets/js/gcv/alignment/merge.ts +++ /dev/null @@ -1,126 +0,0 @@ -/** - * Merges overlapping tracks from the same group to maximize alignment score. - * @param {object} tracks - The micro-synteny viewer tracks to be merged. - */ -export function merge(tracks) { -// make a copy of the data (tracks) - const mergedTracks = {families: tracks.families, groups: []}; - // groups tracks by group id - const groups = {}; - for (let i = 0; i < tracks.groups.length; i++) { // skip first (query) track - const track = tracks.groups[i]; - groups[track.id] = groups[track.id] || []; - groups[track.id].push(track); - } - // try to merge each partition - for (const id in groups) { - if (!groups.hasOwnProperty(id)) { - continue; - } - const groupTracks = groups[id]; - const merged = []; // which tracks have been merged into another - // iterate pairs of tracks to see if one is a sub-inversion of the other - for (let j = 0; j < groupTracks.length; j++) { - if (merged.indexOf(j) !== -1) { - continue; - } - const jTrack = groupTracks[j]; - let jLevels = 0; - const jIds = jTrack.genes.map((g) => g.id); - for (let k = j + 1; k < groupTracks.length; k++) { - if (merged.indexOf(j) !== -1 || merged.indexOf(k) !== -1) { - continue; - } - const kTrack = groupTracks[k]; - let kLevels = 0; - const kIds = kTrack.genes.map((g) => g.id); - // compute the intersection - const overlap = jIds.filter((jId) => { - return kIds.indexOf(jId) !== -1; - }); - if (overlap.length > 0) { - // j is the inversion - if (kIds.length > jIds.length) { - // get index list - const indices = overlap.map((jId) => { - return kIds.indexOf(jId); - }); - // compute the score of the inverted sequence before inverting - const min = Math.min.apply(null, indices); - let max = Math.max.apply(null, indices); - const startGene = kTrack.genes[min]; - const endGene = kTrack.genes[max]; - const score = endGene.score - startGene.score; - // perform the inversion if it will improve the super-track's score - if (jTrack.score > score) { - merged.push(j); - // increment the level counter - kLevels++; - // perform the inversion - jTrack.genes.reverse(); - const args = [min, max - min + 1]; - const geneArgs = args.concat(jTrack.genes); - Array.prototype.splice.apply(kTrack.genes, geneArgs); - // adjust inversion scores and y coordinates - max = min + jTrack.genes.length; - const pred = (min > 0) ? kTrack.genes[min - 1].score : 0; - for (let l = min; l < max; l++) { - kTrack.genes[l].score += pred; - kTrack.genes[l].y = kLevels; - } - // adjust post-inversion scores - const adjustment = jTrack.score - score; - for (let l = max; l < kTrack.genes.length; l++) { - kTrack.genes[l].score += adjustment; - } - kTrack.score += adjustment; - } - // k is the inversion - } else if (jIds.length >= kIds.length) { - // get index list - const indices = overlap.map((jId) => { - return jIds.indexOf(jId); - }); - // compute the score of the inverted sequence before inverting - const min = Math.min.apply(null, indices); - let max = Math.max.apply(null, indices); - const startGene = jTrack.genes[min]; - const endGene = jTrack.genes[max]; - const score = endGene.score - startGene.score; - // perform the inversion if it will improve the super-track's score - if (kTrack.score > score) { - merged.push(k); - // increment the level counter - jLevels++; - // perform the inversion - kTrack.genes.reverse(); - const args = [min, max - min + 1]; - const geneArgs = args.concat(kTrack.genes); - const idArgs = args.concat(kIds); - Array.prototype.splice.apply(jTrack.genes, geneArgs); - Array.prototype.splice.apply(jIds, idArgs); - // adjust inversion scores and y coordinates - max = min + kTrack.genes.length; - const pred = (min > 0) ? jTrack.genes[min - 1].score : 0; - for (let l = min; l < max; l++) { - jTrack.genes[l].score += pred; - jTrack.genes[l].y = jLevels; - } - // adjust post-inversion scores - const adjustment = kTrack.score - score; - for (let l = max; l < jTrack.genes.length; l++) { - jTrack.genes[l].score += adjustment; - } - jTrack.score += adjustment; - } - } - } - } - // add the track if it wasn't merged during its iteration - if (merged.indexOf(j) === -1) { - mergedTracks.groups.push(jTrack); - } - } - } - return mergedTracks; -} diff --git a/client/src/assets/js/gcv/alignment/models/alignment.ts b/client/src/assets/js/gcv/alignment/models/alignment.ts new file mode 100644 index 00000000..bc118dbf --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/alignment.ts @@ -0,0 +1,4 @@ +export class Alignment { + alignment: (number|null)[]; + score: number; +} diff --git a/client/src/assets/js/gcv/alignment/models/index.ts b/client/src/assets/js/gcv/alignment/models/index.ts new file mode 100644 index 00000000..bd321799 --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/index.ts @@ -0,0 +1,5 @@ +export { Alignment } from "./alignment"; +export { InternalAlignment } from "./internal-alignment"; +export { Interval } from "./interval"; +export { Scores } from "./scores"; +export { Traceback } from "./traceback"; diff --git a/client/src/assets/js/gcv/alignment/models/internal-alignment.ts b/client/src/assets/js/gcv/alignment/models/internal-alignment.ts new file mode 100644 index 00000000..cedcaa4a --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/internal-alignment.ts @@ -0,0 +1,4 @@ +export class InternalAlignment { + coordinates: (number|null)[]; + scores: (number|null)[]; +} diff --git a/client/src/assets/js/gcv/alignment/models/interval.ts b/client/src/assets/js/gcv/alignment/models/interval.ts new file mode 100644 index 00000000..4eaa1edd --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/interval.ts @@ -0,0 +1 @@ +export type Interval = [number, number]; diff --git a/client/src/assets/js/gcv/alignment/models/scores.ts b/client/src/assets/js/gcv/alignment/models/scores.ts new file mode 100644 index 00000000..90e98cf7 --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/scores.ts @@ -0,0 +1,6 @@ +export class Scores { + match: number; + mismatch: number; + gap: number; + threshold?: number; +} diff --git a/client/src/assets/js/gcv/alignment/models/traceback.ts b/client/src/assets/js/gcv/alignment/models/traceback.ts new file mode 100644 index 00000000..41b2e9fa --- /dev/null +++ b/client/src/assets/js/gcv/alignment/models/traceback.ts @@ -0,0 +1,6 @@ +export enum Traceback { + FIRST, + DIAGONAL, + LEFT, + UP +}; diff --git a/client/src/assets/js/gcv/alignment/repeat.ts b/client/src/assets/js/gcv/alignment/repeat.ts index 10c6eefd..4bdf7ea6 100644 --- a/client/src/assets/js/gcv/alignment/repeat.ts +++ b/client/src/assets/js/gcv/alignment/repeat.ts @@ -1,176 +1,203 @@ -import { computeScore, matrix } from "./alignment"; +// library +import { matrix, setOption, sum } from "../common"; +import { Alignment, InternalAlignment, Scores, Traceback } from "./models"; +import { computeScore, mergeAlignments } from "./utils"; -/** The Repeat algorithm by Durbin et al. */ -export default function repeat(sequence, reference, options) { - /** - * The alignment algorithm. - * @param {Array} seq - The sequence being aligned to. - * @param {Array} ref - The sequence being aligned. - * @return {object} - The hidden iframe. - */ - function align(seq, ref) { - // populate the matrix - const rows = seq.length + 1; // first item is at index 1 - const cols = ref.length + 1; // first item is at index 1 - const a = matrix(cols, rows, 0); - const choice = [0, 0, 0, 0]; - for (let i = 1; i < cols; i++) { // first column is all 0"s - // handle unmatched regions and ends of matches - const scores = a[i - 1].map((score, j) => { - return (j > 0) ? score - options.scores.threshold : score; - }); - a[i][0] = Math.max.apply(null, scores); - // handle starts of matches and extensions - for (let j = 1; j < rows; j++) { - choice[0] = a[i][0]; - choice[1] = a[i - 1][j - 1] + computeScore( - ref[i - 1], seq[j - 1], options.accessor, options.scores, options.ignore, - ); - // adding because passed value should be negative - choice[2] = a[i - 1][j] + options.scores.gap; - choice[3] = a[i][j - 1] + options.scores.gap; - a[i][j] = Math.max.apply(null, choice); - } - } - let str = ""; - for (let j = 0; j < rows; j++) { - for (let i = 0; i < cols; i++) { - str += "\t" + a[i][j].toString(); +/** + * Given a matrix, a column index, and a threshold, the function returns the + * value for the first cell for the column and a pointer to the cell in the + * previous column that the value was derived from. + * @param {number[][]} m - The score matrix. + * @param {number} i - The index of the column to compute the first cell for. + * @param {number} threshold - The threshold to use when computing the first + * cell. + * @return {Array[2]>[2]} - An array containing the + * computed value for the first cell of the specified column and an array + * containing a pointer to the cell in the preceding column that the value was + * derived from. + */ +function computeFirstCellOfColumn(m: number[][], i: number, threshold: number): +[number, [number, number]] { + const values = m[i-1].map((s) => s-threshold); + values[0] = m[i-1][0]; + const v = Math.max(...values); + const k = values.indexOf(v); + return [v, [i-1, k]]; +} + + +/** + * Aligns the given sequence to the given reference using the repeat algorithm. + * Repeat blocks are detected within the reference and alignments are reported + * as arrays that assign sequence coordinates to elements in the reference. + * @param {Array} seq - The sequence to be aligned to the reference. + * @param {Array} ref - The reference to align the sequence to. + * @param {Scores} scores - An object defining the (mis)match, gap, and + * threshold values to be used during scoring. + * @param {Set} omit - A set of elements to be considered unmatchable. + * @return {InternalAlignment[]} - An array containing an alignment object for + * each aligned block. Each object has a "coordinates" attribute that is an + * array describing each reference element's aligned position in the sequence + * (null if it wasn't aligned) and a "scores" attribute that is an array + * describing what each element in the reference contributes to the alignment's + * score (null if the element wasn't aligned). + * Ex: input: seq: ["A", "B", "C", "D", "E"] + * ref: ["A", "B", "A", "B", "F", "G", "H", "C", "D"] + * scores: {match: 5, mismath: 0, gap: -1, threshold: 10} + * output: [ + * { + * coordinates: [0, 1, null, null, null, null, null, null, null], + * scores: [5, 5, null, null, null, null, null, null, null] + * }, + * { + * coordinates: [null, null, 0, 1, 1.25, 1.5, 1.75, 2, 3], + * scores: [null, null, 5, 5, -1, -1, -1, 5, 5] + * } + * ] + */ +function align( + seq: T[], + ref: T[], + scores: Scores, + omit: Set=new Set): InternalAlignment[] +{ + + // construct score and traceback matrices + const cols = ref.length + 1; // first item is at index 1 + const rows = seq.length + 1; // ditto + const m = matrix(cols, rows, 0); // scores + const t = matrix(cols, rows, [0, 0]); // traceback + for (let i = 1; i < cols; i++) { + // handle unmatched regions and ends of matches + [m[i][0], t[i][0]] = computeFirstCellOfColumn(m, i, scores.threshold); + // handle starts of matches and extensions + for (let j = 1; j < rows; j++) { + const choices = [ + m[i][0], + m[i-1][j-1] + computeScore(ref[i-1], seq[j-1], scores, omit), + m[i-1][j] + scores.gap, + m[i][j-1] + scores.gap + ]; + m[i][j] = Math.max(...choices); + switch (choices.indexOf(m[i][j])) { + case Traceback.FIRST: + t[i][j] = [i, 0]; + break; + case Traceback.DIAGONAL: + t[i][j] = [i-1, j-1]; + break; + case Traceback.LEFT: + t[i][j] = [i-1, j]; + break; + case Traceback.UP: + t[i][j] = [i, j-1]; + break; } - str += "\n"; } - // traceback - make a track for each qualified path in the matrix - let i = cols - 1; // start in the extra cell - let j = 0; // rows - const alignments = []; - let index = -1; - let saving = false; - let length = 0; - while (!(i === 0 && j === 0)) { - if (j === 0) { - if (saving && length < 2) { - alignments.pop(); - index--; - length = 0; - } - saving = false; - const max = Math.max.apply(null, a[i]); - const jMax = a[i].lastIndexOf(max); - // start a new alignment only if j is a match and the alignment's - // score meets the threshold - // TODO: why is the matrix not being used? - const ref_f = options.accessor(ref[i - 1]); - const seq_f = options.accessor(seq[jMax - 1]); - if (jMax > 0 && i > 0 && ref_f === seq_f && - options.ignore.indexOf(ref_f) === -1 && max >= options.scores.threshold) { - length = 1; - saving = true; - j = jMax; - alignments.push({sequence: [], reference: [], score: max}); - index++; - // pad seq with the genes not traversed by the alignment - for (let k = seq.length - 1; k >= j; k--) { - alignments[index].sequence.unshift(seq[k]); - alignments[index].reference.unshift(null); - } - // add the starting match - alignments[index].sequence.unshift(seq[j - 1]); - alignments[index].reference.unshift(ref[i - 1]); - if (options.suffixScores) { - alignments[index].reference[0].score = a[i][j]; - } - } else { - // try starting an alignment in the next column - i--; - } - } else if (i === 0) { - j = 0; - } else { - // diag, up, left - const scores = [a[i - 1][j - 1], a[i][j - 1], a[i - 1][j]]; - const max = Math.max.apply(null, scores); - // stop alignment if a 0 cell was reached - if (max === 0) { - // add any missing genes to seq - if (saving) { - for (let k = j - 1; k > 0; k--) { - alignments[index].sequence.unshift(seq[k - 1]); - alignments[index].reference.unshift(null); - } - } - // get ready for a new alignment - i--; - j = 0; - } else { - switch (scores.lastIndexOf(max)) { - // diag - case 0: - j--; - i--; - // no alignments happen in the first row or column - if (saving && j > 0 && i > 0) { - alignments[index].sequence.unshift(seq[j - 1]); - alignments[index].reference.unshift((ref[i - 1])); - if (options.suffixScores) { - alignments[index].reference[0].score = a[i][j]; - } - length++; - } - break; - // up - case 1: - j--; - if (saving && j > 0) { - alignments[index].sequence.unshift(seq[j - 1]); - alignments[index].reference.unshift(null); - } - break; - // left - case 2: - i--; - if (saving && i > 0) { - alignments[index].sequence.unshift(null); - alignments[index].reference.unshift(ref[i - 1]); - if (options.suffixScores) { - alignments[index].reference[0].score = a[i][j]; - } - } - break; - } + } + + // construct alignments via traceback + const alignments: InternalAlignment[] = []; + let a = {coordinates: [], scores: []}; + let insertion = 0; + let [, [i, j]] = computeFirstCellOfColumn(m, cols, scores.threshold); + while (!(i === 0 && j === 0)) { + let [i2, j2] = t[i][j]; + // start new alignment + if (j2 > j) { + a = { + coordinates: Array(ref.length-i2).fill(null), + scores: Array(ref.length-i2).fill(null) + }; + // insertion + } else if (j2 === j && j !== 0) { + insertion += 1; + // (mis)match + } else if (j2 === j-1 && i2 === i-1) { + // backfill insertion + if (insertion > 0) { + let step = 1/(insertion+1); + for (let k = insertion-1; k >= 0; k--) { + const x = j + (k+1)*step; + a.coordinates.unshift(x); + a.scores.unshift(m[i][j+k+1]-m[i][j+k]); } + insertion = 0; } + a.coordinates.unshift(j) + a.scores.unshift(m[i][j]-m[i2][j2]); } - if (saving && length < 2) { - alignments.pop(); - index--; - length = 0; + // end alignment + if (j > 0 && j2 === 0) { + if (a.coordinates.filter((x) => x !== null).length > 2) { + const fill = Array(ref.length-a.coordinates.length).fill(null); + a.coordinates.unshift(...fill); + a.scores.unshift(...fill); + alignments.push(a); + } + insertion = 0; } - return alignments; + [i, j] = [i2, j2]; } + return alignments; +} + +/** The Repeat local alignment algorithm by Durbin et al. + * @param {Array} reference - The sequence to be aligned to. + * @param {Array} sequence - The sequence to align to the reference. + * @param {object} options - Optional parameters. + * @return {Alignment[]} - Local alignments in the form of coordinates + * describing which element in reference each element in sequence aligned to. In + * the case of an insertion, the value will be a decimal between the flanking + * match states. Non-aligned genes get a null coordinate value. + * Ex: [0, 1, 2, 3, 4, 5] // a perfect alignment + * [0, 0.5, 1, 2, 3, 4] // an insertion + * [0, 0.3, 0.6, 1, 2 ] // sequential insertions + * [0, 1, 2, 5, 6, 7] // deletions + * [null, null, 3, 4, 5, null] // ends weren't aligned + */ +export function repeat( + reference: T[], + sequence: T[], + options: any={}): Alignment[] +{ + // parse optional parameters - options = options || {}; - options.accessor = options.accessor || ((e) => e); - options.scores = options.scores || {}; - options.scores.match = options.scores.match || 5; - options.scores.mismatch = options.scores.mismatch || 0; - options.scores.gap = options.scores.gap || -1; - options.scores.threshold = options.scores.threshold || 10; - options.suffixScores = options.suffixScores || false; - options.ignore = options.ignore || []; - options.reverse = options.reverse || ((s) => { - const r = s.slice(); - r.reverse(); - return r; - }); + options = Object.assign({}, options); + options.scores = Object.assign({}, options.scores); + setOption(options.scores, "match", 5); + setOption(options.scores, "mismatch", 0); + setOption(options.scores, "gap", -1); + setOption(options.scores, "threshold", 10); + setOption(options, "omit", new Set()); + setOption(options, "reversals", true); + setOption(options, "inversions", true); // perform forward and reverse alignments - const forwards = align(sequence, reference); - const reverseReference = options.reverse(reference); - const reverses = align(sequence, reverseReference); + const forward = sequence; + const forwardAlignments = + align(reference, forward, options.scores, options.omit); + const reverse = (options.reversals || options.inversions) ? + [...forward].reverse() : []; + const reverseAlignments = + align(reference, reverse, options.scores, options.omit); + // reverse the reverse alignment orderings and the alignments themselves + reverseAlignments + .reverse() + .forEach(({coordinates, scores}) => { + coordinates.reverse(); + scores.reverse(); + }); + + // merge alignments + const alignments = mergeAlignments( + forwardAlignments, + reverseAlignments, + options.reversals, + options.inversions) + .map((a) => ({alignment: a.coordinates, score: sum(a.scores)})); - // clone each object in the arrays and flip the strand for each selected gene - const output = forwards.concat(reverses); - return output; + return alignments; } diff --git a/client/src/assets/js/gcv/alignment/smith-waterman.ts b/client/src/assets/js/gcv/alignment/smith-waterman.ts index bab333cd..c3df9ff1 100644 --- a/client/src/assets/js/gcv/alignment/smith-waterman.ts +++ b/client/src/assets/js/gcv/alignment/smith-waterman.ts @@ -1,112 +1,157 @@ -import { computeScore, matrix } from "./alignment"; +// library +import { matrix, setOption, sum } from "../common"; +import { Alignment, InternalAlignment, Scores, Traceback } from "./models"; +import { computeScore, mergeAlignments } from "./utils"; + /** - * The Smith-Waterman algorithm. - * @param {Array} sequence - The sequence to be aligned. - * @param {Array} reference - The sequence to be aligned to. - * @param {object} options - Optional parameters. - * @return {Array} - The aligned sequences and score. + * Aligns the given sequence to the given reference using the Smith-Waterman + * alignment algorithm. + * @param {Array} seq - The sequence to be aligned to the refernece. + * @param {Array} ref - The reference to align the sequence to. + * @param {Scores} scores - An object defining the (mis)match and gap values to + * be used during scoring. + * @param {Set} omit - A set of elements to be considered unmatchable. + * @return {InternalAlignment} - An alignment object with a "coordinates" + * attribute that is an array describing each reference element's aligned + * position in the sequence (null if it wasn't aligned) and a "scores" attribute + * that is an array describing what each element in the reference contributes to + * the alignment's score (null if the element wasn't aligned). + * Ex: TODO */ -export default function smithWaterman(sequence, reference, options) { +function align( + seq: T[], + ref: T[], + scores: Scores, + omit: Set=new Set): InternalAlignment +{ - /** - * The actual alignment algorithm. - * @param {Array} seq - The sequence being aligned to. - * @param {Array} ref - The sequence being aligned. - * @return {object} - The hidden iframe. - */ - function align(seq, ref) { - // initialize letiables - const rows = ref.length + 1; - const cols = seq.length + 1; - const a = matrix(rows, cols, 0); - let i = 0; - let j = 0; - const choice = [0, 0, 0, 0]; - // populate the matrix - let max = 0; - let iMax = 0; - let jMax = 0; - for (i = 1; i < rows; i++) { - for (j = 1; j < cols; j++) { - choice[0] = 0; - choice[1] = a[i - 1][j - 1] + computeScore( - ref[i - 1], seq[j - 1], options.accessor, options.scores, options.ignore, - ); - choice[2] = a[i - 1][j] + options.scores.gap; - choice[3] = a[i][j - 1] + options.scores.gap; - a[i][j] = Math.max.apply(null, choice); - if (a[i][j] >= max) { - max = a[i][j]; - iMax = i; - jMax = j; - } + // construct score and traceback matrices + const cols = ref.length + 1; // first item is at index 1 + const rows = seq.length + 1; // ditto + const m = matrix(cols, rows, 0); // scores + const t = matrix(cols, rows, [0, 0]); // traceback + let max = 0; + let maxCell = [0, 0]; + for (let i = 1; i < cols; i++) { + for (let j = 1; j < rows; j++) { + const choices = [ + 0, + m[i-1][j-1] + computeScore(ref[i-1], seq[j-1], scores, omit), + m[i-1][j] + scores.gap, + m[i][j-1] + scores.gap + ]; + m[i][j] = Math.max(...choices); + switch (choices.indexOf(m[i][j])) { + case Traceback.FIRST: + // points to default [0, 0] + break; + case Traceback.DIAGONAL: + t[i][j] = [i-1, j-1]; + break; + case Traceback.LEFT: + t[i][j] = [i-1, j]; + break; + case Traceback.UP: + t[i][j] = [i, j-1]; + break; } - } - // traceback - i = iMax; - j = jMax; - let diag; - let up; - let left; - let score = max; - const outRef = []; - const outSeq = []; - while (i > 0 && j > 0) { - score = a[i][j]; - if (score === 0) { - break; + if (m[i][j] > max) { + max = m[i][j]; + maxCell = [i, j]; } - diag = a[i - 1][j - 1]; - up = a[i][j - 1]; - left = a[i - 1][j]; - if (score === diag + computeScore( - ref[i - 1], seq[j - 1], options.accessor, options.scores, options.ignore, - )) { - outRef.unshift(ref[i - 1]); - outSeq.unshift(seq[j - 1]); - i -= 1; - j -= 1; - } else if (score === left + options.scores.gap) { - outRef.unshift(ref[i - 1]); - outSeq.unshift(null); - i -= 1; - } else if (score === up + options.scores.gap) { - outRef.unshift(null); - outSeq.unshift(seq[j - 1]); - j -= 1; - } else { - break; + } + } + + // construct alignment via traceback + const alignment = {coordinates: [], scores: []}; + let insertion = 0; + let [i, j] = maxCell; + while (m[i][j] !== 0) { + let [i2, j2] = t[i][j]; + // insertion + if (j2 === j && j !== 0) { + insertion += 1; + // (mis)match + } else if (j2 === j-1 && i2 === i-1) { + // backfill insertion + if (insertion > 0) { + let step = 1/(insertion+1); + for (let k = insertion-1; k >= 0; k--) { + const x = j + (k+1)*step; + alignment.coordinates.unshift(x); + alignment.scores.unshift(m[i][j+k+1]-m[i][j+k]); + } + insertion = 0; } + alignment.coordinates.unshift(j) + alignment.scores.unshift(m[i][j]-m[i2][j2]); } - while (j > 0) { - outRef.unshift(null); - outSeq.unshift(seq[j - 1]); - j -= 1; + // end alignment + if (m[i2][j2] === 0) { + if (alignment.coordinates.filter((x) => x !== null).length > 2) { + const fill = Array(ref.length-alignment.coordinates.length).fill(null); + alignment.coordinates.unshift(...fill); + alignment.scores.unshift(...fill); + } } - return {sequence: outSeq, reference: outRef, score: max}; + [i, j] = [i2, j2]; } + return alignment; +} + + +/** + * The Smith-Waterman algorithm. + * @param {Array} reference - The sequence to be aligned to. + * @param {Array} sequence - The sequence to align to the reference. + * @param {object} options - Optional parameters. + * @return {Alignment[]} - Local alignments in the form of coordinates + * describing which element in reference each element in sequence aligned to. In + * the case of an insertion, the value will be a decimal between the flanking + * match states. Non-aligned genes get a null coordinate value. + * Ex: [0, 1, 2, 3, 4, 5] // a perfect alignment + * [0, 0.5, 1, 2, 3, 4] // an insertion + * [0, 0.3, 0.6, 1, 2 ] // sequential insertions + * [0, 1, 2, 5, 6, 7] // deletions + * [null, null, 3, 4, 5, null] // ends weren't aligned + */ +export function smithWaterman( + reference: T[], + sequence: T[], + options: any={}): Alignment[] +{ + // parse optional parameters - options = options || {}; - options.accessor = options.accessor || ((e) => e); - options.scores = options.scores || {}; - options.scores.match = options.scores.match || 5; - options.scores.mismatch = options.scores.mismatch || 0; - options.scores.gap = options.scores.gap || -1; - options.ignore = options.ignore || []; - options.reverse = options.reverse || ((s) => { - const r = s.slice(); - r.reverse(); - return r; - }); + options = Object.assign({}, options); + options.scores = Object.assign({}, options.scores); + setOption(options.scores, "match", 5); + setOption(options.scores, "mismatch", 0); + setOption(options.scores, "gap", -1); + setOption(options, "omit", new Set()); + setOption(options, "reverse", true); + setOption(options, "inversions", true); // perform forward and reverse alignments - const forward = align(sequence, reference); - const reverseReference = options.reverse(reference); - const reverse = align(sequence, reverseReference); + const forward = sequence; + const forwardAlignment = + align(reference, forward, options.scores, options.omit); + const reverse = (options.reverse || options.inversions) ? + [...forward].reverse() : []; + const reverseAlignment = + align(reference, reverse, options.scores, options.omit); + // reverse the reverse alignment + reverseAlignment.coordinates.reverse(); + reverseAlignment.scores.reverse(); + + // merge alignments + const alignments = mergeAlignments( + [forwardAlignment], + [reverseAlignment], + options.reversals, + options.inversions) + .map((a) => ({alignment: a.coordinates, score: sum(a.scores)})); - // output the highest scoring alignment - const output = forward.score >= reverse.score ? forward : reverse; - return [output]; + return alignments; } diff --git a/client/src/assets/js/gcv/alignment/trackify.ts b/client/src/assets/js/gcv/alignment/trackify.ts deleted file mode 100644 index 10adac29..00000000 --- a/client/src/assets/js/gcv/alignment/trackify.ts +++ /dev/null @@ -1,92 +0,0 @@ -/** - * Converts the inverted portion of an alignment into a track portion. - * @param {number} i - The current position in the reference alignment. - * @param {number} insertionCount - The genes to be inserted. - * @param {number} queryCount - The number of query genes that have been - * inserted. - * @param {object} alignment - The alignment object to use. - * @param {Array} track - The current track being constructed. - */ -function insertion(i, insertionCount, queryCount, alignment, track) { - const start = i - insertionCount; - const step = 1.0 / (insertionCount + 1); - for (let j = start; j < i; j++) { - if (alignment.reference[j] != null) { - alignment.reference[j].x = queryCount + (step * (j - start + 1)) - 1; - alignment.reference[j].y = 0; - track.genes.push(alignment.reference[j]); - } - } -} - -/** - * Converts alignments into micro-synteny viewer tracks. - * @param {object} query - The reference tracck that all others were aligned to. - * @param {object} tracks - The original viewer tracks. - * @param {Array} alignments - The alignments of the tracks to be trackified. - * @return {object} - A new set of MicroTracks derived from the alignments. - */ -export default function trackify(query, tracks, alignments) { - // make a copy of the data (tracks) and only save the first group (query) - //const aligned = JSON.parse(JSON.stringify(data)); - const aligned = {families: tracks.families, groups: []}; - //const query = aligned.groups[0]; - if (query !== undefined) { - //aligned.groups = [query]; - // initialize letiables - const length = query.genes.length; - // update the context data with the alignment - for (const alignment of alignments) { - let queryCount = 0; - let preQuery = 0; - let insertionCount = 0; - const track = alignment.track; - track.score = alignment.score; - track.genes = []; - for (let i = 0; i < alignment.sequence.length; i++) { - // keep track of how many selected genes come before the query genes - if (alignment.sequence[i] === null && queryCount === 0) { - preQuery += 1; - // an insertion - } else if (alignment.sequence[i] == null) { - // position the genes that come after the query genes - if (queryCount >= length) { - alignment.reference[i].x = queryCount++; - alignment.reference[i].y = 0; - track.genes.push(alignment.reference[i]); - // track how many genes were inserted - } else { - insertionCount++; - } - // a deletion - } else if (alignment.reference[i] == null) { - if (insertionCount > 0) { - insertion(i, insertionCount, queryCount, alignment, track); - insertionCount = 0; - } - queryCount++; - // a (mis)match - } else { - // position the genes that came before the query - if (preQuery > 0) { - for (let j = 0; j < preQuery; j++) { - alignment.reference[j].x = -(preQuery - (j + 1)); - alignment.reference[j].y = 0; - track.genes.push(alignment.reference[j]); - } - preQuery = 0; - // position the genes that go between query genes - } else if (insertionCount > 0) { - insertion(i, insertionCount, queryCount, alignment, track); - insertionCount = 0; - } - alignment.reference[i].x = queryCount++; - alignment.reference[i].y = 0; - track.genes.push(alignment.reference[i]); - } - } - aligned.groups.push(track); - } - } - return aligned; -} diff --git a/client/src/assets/js/gcv/alignment/utils/alignment-interval.ts b/client/src/assets/js/gcv/alignment/utils/alignment-interval.ts new file mode 100644 index 00000000..7ce0d5cf --- /dev/null +++ b/client/src/assets/js/gcv/alignment/utils/alignment-interval.ts @@ -0,0 +1,22 @@ +/** + * Takes a (local) alignment and converts it into an interval describing the + * non-null portion of the alignment. + * @param {Array} a - The alignment to be converted. + * @return {[number, number]} - A tuple describing the non-null interval of the + * alignment. + */ +export function alignmentInterval(a: number[]): [number, number] { + let begin: number; + let end: number; + let i = 0; + while (i < a.length && a[i] === null) { + i += 1; + } + begin = i; + i = a.length-1; + while (i >= 0 && a[i] === null) { + i -= 1; + } + end = i; + return [begin, end]; +} diff --git a/client/src/assets/js/gcv/alignment/utils/compute-score.ts b/client/src/assets/js/gcv/alignment/utils/compute-score.ts new file mode 100644 index 00000000..a8a19fcc --- /dev/null +++ b/client/src/assets/js/gcv/alignment/utils/compute-score.ts @@ -0,0 +1,13 @@ +/** + * Uses an accessor to compute a score by comparing the given elements. + * @param {T} a - The first element. + * @param {T} b - The second element. + * @param {object} scores - The contains the match and mismatch scores. + * @return {int} - The computed score. + */ +export function computeScore(a: T, b: T, scores, omit=new Set()): number { + if (a === b && !omit.has(a)) { + return scores.match; + } + return scores.mismatch; +} diff --git a/client/src/assets/js/gcv/alignment/utils/index.ts b/client/src/assets/js/gcv/alignment/utils/index.ts new file mode 100644 index 00000000..d4f156a1 --- /dev/null +++ b/client/src/assets/js/gcv/alignment/utils/index.ts @@ -0,0 +1,4 @@ +export { alignmentInterval } from "./alignment-interval"; +export { computeScore } from "./compute-score"; +export { intervalsToSets } from "./intervals-to-sets"; +export { mergeAlignments } from "./merge-alignments"; diff --git a/client/src/assets/js/gcv/alignment/utils/intervals-to-sets.ts b/client/src/assets/js/gcv/alignment/utils/intervals-to-sets.ts new file mode 100644 index 00000000..423e9061 --- /dev/null +++ b/client/src/assets/js/gcv/alignment/utils/intervals-to-sets.ts @@ -0,0 +1,62 @@ +// app +import { Interval } from "../models"; +// dependencies +import { StaticDisjointSet } from "mnemonist"; + + +/** + * Given forward and reverse interval arrays (assumes both are sorted by start + * position and elements don't overlap), the function creates sets of indices in + * the forward and reverse interval arrays that correspond to overlapping + * intervals. + * @param {Array<[number, number]>} forward - The forward intervals. + * @param {Array<[number, number]>} reverse - The reverse intervals. + * @return {Array} - An array containing an object for each set of + * overlapping intervals. Each object has a "forward" attribute containing an + * array of the forward intervals in the set and a "reverse" attribute + * containing the reverse intervals in the set. The intervals in each array are + * ordered the same as in the input arrays. + */ +export function intervalsToSets( + forward: Array, + reverse: Array, + inversions: boolean +): Array<{forward: Array, reverse: Array}> { + // use a disjoint set and union-find to track interval sets + const sets = new StaticDisjointSet(forward.length+reverse.length); + // dovetail iterate intervals to find overlaps + let f = 0; + let r = 0; + while (f < forward.length && r < reverse.length) { + let [fBegin, fEnd] = forward[f]; + let [rBegin, rEnd] = reverse[r]; + // consider neighboring intervals as overlapping for inversions + if (inversions) { + rBegin -= 1; + rEnd += 1; + } + if (fEnd < rBegin) { + f += 1; + } else if (rEnd < fBegin) { + r += 1; + } else if (rBegin <= fEnd && fEnd <= rEnd) { + sets.union(f, r+forward.length); + f += 1; + } else if (fBegin <= rEnd && rEnd <= fEnd) { + sets.union(f, r+forward.length); + r += 1; + } + } + // convert disjoint sets to interval objects + const reducer = (accumulator, i) => { + if (i < forward.length) { + accumulator.forward.push(i); + } else { + accumulator.reverse.push(i-forward.length); + } + return accumulator; + }; + const overlaps = sets.compile() + .map((s) => s.reduce(reducer, {forward: [], reverse: []})); + return overlaps; +} diff --git a/client/src/assets/js/gcv/alignment/utils/merge-alignments.ts b/client/src/assets/js/gcv/alignment/utils/merge-alignments.ts new file mode 100644 index 00000000..b04a3d8b --- /dev/null +++ b/client/src/assets/js/gcv/alignment/utils/merge-alignments.ts @@ -0,0 +1,167 @@ +import { filterByIndex, sum } from "../../common"; +import { InternalAlignment, Interval } from "../models"; +import { alignmentInterval } from "./alignment-interval"; +import { intervalsToSets } from "./intervals-to-sets"; + + +function combineAlignments( + alignments: InternalAlignment[], + intervals: Interval[], + indexes: number[]): InternalAlignment +{ + const l = alignments[0].coordinates.length; + const alignment = { + coordinates: Array(l).fill(null), + scores: Array(l).fill(null) + }; + indexes.forEach((i) => { + let [begin, end] = intervals[i]; + const coordinates = alignments[i].coordinates.slice(begin, end+1); + const scores = alignments[i].scores.slice(begin, end+1); + alignment.coordinates.splice(begin, end+1-begin, ...coordinates); + alignment.scores.splice(begin, end+1-begin, ...scores); + }); + return alignment; +} + + +function reversalsAndInversions( + flatForward: InternalAlignment, + flatReverse: InternalAlignment): InternalAlignment +{ + const resolveOverlap = (begin, end) => { + const fScores = flatForward.scores.slice(begin, end); + const rScores = flatReverse.scores.slice(begin, end); + // add the inversion + if (sum(fScores) < sum(rScores)) { + const rCoordinates = flatReverse.coordinates.slice(begin, end); + return [rCoordinates, rScores]; + } + // add the forward orientation + const fCoordinates = flatForward.coordinates.slice(begin, end); + return [fCoordinates, fScores]; + }; + // iterate flattened alignments simultaneously and resolves inversions + const coordinates = []; + const scores = []; + let overlap = 0; + for (let i = 0; i < flatForward.coordinates.length; i++) { + const fCoordinate = flatForward.coordinates[i]; + const rCoordinate = flatReverse.coordinates[i]; + if (fCoordinate === null || rCoordinate === null) { + // back add overlap + if (overlap > 0) { + let [c, s] = resolveOverlap(i-overlap, i); + coordinates.push(...c); + scores.push(...s); + overlap = 0; + } + // add forward + if (fCoordinate !== null) { + coordinates.push(fCoordinate); + scores.push(flatForward.scores[i]); + // add reverse + } else if (rCoordinate !== null) { + coordinates.push(rCoordinate); + scores.push(flatReverse.scores[i]); + } else { + coordinates.push(null); + scores.push(null); + } + } else { + overlap += 1; + } + } + // back add overlap + if (overlap > 0) { + let [c, s] = resolveOverlap(-overlap, undefined); + coordinates.push(...c); + scores.push(...s); + } + return {coordinates, scores}; +} + + +function reversalsOnly( + flatForward: InternalAlignment, + flatReverse: InternalAlignment): InternalAlignment +{ + // TODO: implement w/ weighted interval scheduling dynamic program + return reversalsAndInversions(flatForward, flatReverse); +} + + +function inversionsOnly( + flatForward: InternalAlignment, + flatReverse: InternalAlignment): InternalAlignment +{ + // TODO: implement + return reversalsAndInversions(flatForward, flatReverse); +} + + +/** + * Given a set of forward and reverse alignment segments for the same sequence + * and reference, the function merges overlapping segments in a manner dictated + * by whether or not reversals or inversions are allowed. + * @param {Array} forwardAlignments - The forward alignments. + * @param {Array} reverseAlignments - The reverse alignments. + * @param {boolean} reversals - Whether or not reversals are allowed. + * @param {boolean} inversions - Whether or not inversions are allowed. + * @return{Array} - The merged alignments. + */ +export function mergeAlignments( + forwardAlignments: InternalAlignment[], + reverseAlignments: InternalAlignment[], + reversals: boolean, + inversions: boolean): InternalAlignment[] +{ + + // make a sequence coordinate interval for each alignment + const forwardIntervals = forwardAlignments.map((a): Interval => { + return alignmentInterval(a.coordinates); + }); + const reverseIntervals = reverseAlignments.map((a): Interval => { + return alignmentInterval(a.coordinates); + }); + + // find sets of overlapping intervals + const overlaps = + intervalsToSets(forwardIntervals, reverseIntervals, inversions); + + // convert the sets of overlapping intervals into alignments + const alignments = []; + overlaps.forEach(({forward, reverse}) => { + // just save the forward alignments + if (forward.length !== 0 && (reverse.length === 0 || + (!reversals && !inversions))) { + alignments.push(...filterByIndex(forwardAlignments, forward)); + // just save the reverse alignments + } else if (forward.length === 0 && reverse.length !== 0) { + if (reversals) { + alignments.push(...filterByIndex(reverseAlignments, reverse)); + } + } else { + // combine each orientation's alignments and scores + const flatForward = + combineAlignments(forwardAlignments, forwardIntervals, forward); + const flatReverse = + combineAlignments(reverseAlignments, reverseIntervals, reverse); + let alignment; + // combine all alignments such that the score is maximized + if (reversals && inversions) { + alignment = reversalsAndInversions(flatForward, flatReverse); + // keep a set of non-overlapping alignments with maximized score + } else if (reversals && !inversions) { + alignment = reversalsOnly(flatForward, flatReverse); + // only sub-intervals of forward segments can be inverted + } else if (!reversals && inversions) { + alignment = inversionsOnly(flatForward, flatReverse); + } + // else handled by "just save the forward alignments" + alignments.push(alignment); + } + }); + + return alignments; +} diff --git a/client/src/assets/js/gcv/common/filter-by-index.ts b/client/src/assets/js/gcv/common/filter-by-index.ts new file mode 100644 index 00000000..d469c291 --- /dev/null +++ b/client/src/assets/js/gcv/common/filter-by-index.ts @@ -0,0 +1,11 @@ +/** + * Given an array of elements and an array of indexes, the function returns an + * array containing elements only from the given indexes. + * @param {Array} a - The array to filter. + * @param {Array} idx - The indexes to filter the array with. + * @return {Array} - The filtered array. Note, the elements are ordered + * according to the index array. + */ +export function filterByIndex(a: T[], idx: number[]): T[] { + return idx.map((i) => a[i]); +} diff --git a/client/src/assets/js/gcv/common/index.ts b/client/src/assets/js/gcv/common/index.ts index 6fab9824..7c133646 100644 --- a/client/src/assets/js/gcv/common/index.ts +++ b/client/src/assets/js/gcv/common/index.ts @@ -1,4 +1,8 @@ export { clone } from "./clone"; export { colors } from "./colors"; export { eventBus } from "./eventbus"; +export { filterByIndex } from "./filter-by-index"; export { getFamilySizeMap } from "./get-family-size-map"; +export { matrix } from "./matrix"; +export { setOption } from "./set-option"; +export { sum } from "./sum"; diff --git a/client/src/assets/js/gcv/common/matrix.ts b/client/src/assets/js/gcv/common/matrix.ts new file mode 100644 index 00000000..9a9d94cd --- /dev/null +++ b/client/src/assets/js/gcv/common/matrix.ts @@ -0,0 +1,21 @@ +/** + * Creates a matrix with the given initial value. + * @param {number} m - Number of columns. + * @param {number} n - Number of rows. + * @param {number} initial - The initial value to put in each cell. + * @return {number[][]} - An m x n array with initial value in each cell. + */ +export function matrix(m, n, initial) { + let a; + let i; + let j; + const mat = []; + for (i = 0; i < m; i += 1) { + a = []; + for (j = 0; j < n; j += 1) { + a[j] = initial; + } + mat[i] = a; + } + return mat; +} diff --git a/client/src/assets/js/gcv/common/set-option.ts b/client/src/assets/js/gcv/common/set-option.ts new file mode 100644 index 00000000..5a35402a --- /dev/null +++ b/client/src/assets/js/gcv/common/set-option.ts @@ -0,0 +1,12 @@ +/** + * A help for assigning default values to an optional parameter object. + * @param {object} options - The optional parameter object. + * @param {string} option - The name of the option to set. + * @param {any} value - The default value for the given option. + */ +export function setOption(options: object, option: string, value: any): void { + if (options[option] === undefined || + typeof options[option] !== typeof value) { + options[option] = value; + } +} diff --git a/client/src/assets/js/gcv/common/sum.ts b/client/src/assets/js/gcv/common/sum.ts new file mode 100644 index 00000000..63ad8b51 --- /dev/null +++ b/client/src/assets/js/gcv/common/sum.ts @@ -0,0 +1,8 @@ +/** + * Sums the numbers in the given arrays. + * @param {Array} - The array to sum. + * @return {number} - The sum. + */ +export function sum(l) { + return l.reduce((a, b) => a + b); +} diff --git a/client/src/assets/js/gcv/graph/msa-hmm.ts b/client/src/assets/js/gcv/graph/msa-hmm.ts index 6bec2dcf..051e24e8 100644 --- a/client/src/assets/js/gcv/graph/msa-hmm.ts +++ b/client/src/assets/js/gcv/graph/msa-hmm.ts @@ -711,6 +711,7 @@ export class MSAHMM extends Directed { // parse options this._setOption(options, "reverse", true); this._setOption(options, "inversions", true); + // compute Viterbi paths and their emission probabilities let forward = sequence; const forwardPath = this.viterbi(forward); @@ -719,6 +720,7 @@ export class MSAHMM extends Directed { [...forward].reverse() : []; const reversePath = this.viterbi(reverse); let reverseEmissions = this.sequenceEmissions(reverse, reversePath); + // swap orientations depending on score const reversed = options.reverse && forwardPath.probability < reversePath.probability; @@ -726,6 +728,7 @@ export class MSAHMM extends Directed { [forward, forwardEmissions, reverseEmissions] = [reverse, reverseEmissions, forwardEmissions]; } + // compute inversions let alignmentPath = forwardPath; if (options.inversions) { @@ -733,12 +736,14 @@ export class MSAHMM extends Directed { this._emissionsToOrientation(forwardEmissions, reverseEmissions); alignmentPath = this._mergePaths(forwardPath, reversePath, orientations); } + // convert path to hmm state independent alignment let alignment = this.pathToAlignment(alignmentPath); if (reversed) { const l = sequence.length-1; alignment = alignment.map((x) => l-x); } + return alignment; }