Improvements to BAM visualization

lib/underscore.js

@@ -28,5 +28,9 @@ declare module "underscore" {
   declare function sortBy<T>(a: T[], iteratee: (val: T)=>any): T[];
 
   declare function filter<T>(o: {[key:string]: T}, pred: (val: T, k: string)=>boolean): T[];
+
+  declare function isEmpty(o: any): boolean;
+
+  declare function groupBy<T>(a: Array<T>, iteratee: (val: T, index: number)=>any): {[key:string]: T};
 }
 

src/BamDataSource.js

@@ -36,6 +36,9 @@ function createBamSource(remoteSource: BamFile): BamDataSource {
   // Keys are virtualOffset.toString()
   var reads: {[key:string]: SamRead} = {};
 
+  // Mapping from contig name to canonical contig name.
+  var contigNames: {[key:string]: string} = {};
+
   // Ranges for which we have complete information -- no need to hit network.
   var coveredRanges: ContigInterval<string>[] = [];
 
@@ -47,25 +50,48 @@ function createBamSource(remoteSource: BamFile): BamDataSource {
   }
 
   function fetch(range: GenomeRange) {
-    var interval = new ContigInterval(range.contig, range.start, range.stop);
-
-    // Check if this interval is already in the cache.
-    if (interval.isCoveredBy(coveredRanges)) {
-      return Q.when();
+    var refsPromise;
+    if (!_.isEmpty(contigNames)) {
+      refsPromise = Q.when();
+    } else {
+      refsPromise = remoteSource.header.then(header => {
+        header.references.forEach(ref => {
+          var name = ref.name;
+          contigNames[name] = name;
+          contigNames['chr' + name] = name;
+          if (name.slice(0, 3) == 'chr') {
+            contigNames[name.slice(3)] = name;
+          }
+        });
+      });
     }
 
-    interval = expandRange(interval);
-    return remoteSource.getAlignmentsInRange(interval).then(reads => {
-      coveredRanges.push(interval);
-      coveredRanges = ContigInterval.coalesce(coveredRanges);
-      reads.forEach(read => addRead(read));
+    return refsPromise.then(() => {
+      var contigName = contigNames[range.contig];
+      var interval = new ContigInterval(contigName, range.start, range.stop);
+
+      // Check if this interval is already in the cache.
+      if (interval.isCoveredBy(coveredRanges)) {
+        return Q.when();
+      }
+
+      interval = expandRange(interval);
+      return remoteSource.getAlignmentsInRange(interval).then(reads => {
+        coveredRanges.push(interval);
+        coveredRanges = ContigInterval.coalesce(coveredRanges);
+        reads.forEach(read => addRead(read));
+      });
     });
   }
 
   function getAlignmentsInRange(range: ContigInterval<string>): SamRead[] {
     if (!range) return [];
-    // XXX there may be an issue here with adding 'chr' to contig names.
-    return _.filter(reads, read => read.intersects(range));
+    if (_.isEmpty(contigNames)) return [];
+
+    var canonicalRange = new ContigInterval(contigNames[range.contig],
+                                            range.start(), range.stop());
+
+    return _.filter(reads, read => read.intersects(canonicalRange));
   }
 
   var o = {

src/PileupTrack.js

@@ -9,7 +9,9 @@ import type * as SamRead from './SamRead';
 var React = require('react/addons'),
     _ = require('underscore'),
     d3 = require('d3'),
-    types = require('./types');
+    types = require('./types'),
+    Interval = require('./Interval'),
+    {pileup} = require('./pileuputils');
 
 var PileupTrack = React.createClass({
   propTypes: {
@@ -28,7 +30,8 @@ var PileupTrack = React.createClass({
 });
 
 
-var READ_HEIGHT = 10;
+var READ_HEIGHT = 13;
+var READ_SPACING = 2;  // vertical pixels between reads
 
 
 var NonEmptyPileupTrack = React.createClass({
@@ -73,6 +76,9 @@ var NonEmptyPileupTrack = React.createClass({
 
     var scale = this.getScale();
 
+    var rows = pileup(this.props.reads.map(
+        r => new Interval(r.pos, r.pos + r.l_seq)));
+
     svg.attr('width', width)
        .attr('height', height);
 
@@ -87,8 +93,8 @@ var NonEmptyPileupTrack = React.createClass({
     // Update
     reads.attr({
       'x': read => scale(read.pos),
-      'y': (read, i) => i * READ_HEIGHT,
-      'width': read => (scale(read.pos + read.l_seq) - scale(read.pos)),
+      'y': (read, i) => rows[i] * (READ_HEIGHT + READ_SPACING),
+      'width': read => (scale(read.pos + read.l_seq) - scale(read.pos) - 5),
       'height': READ_HEIGHT
     });
 

src/bam.js

@@ -259,7 +259,9 @@ class Bam {
     return this.header.then(header => {
       for (var i = 0; i < header.references.length; i++) {
         var name = header.references[i].name;
-        if (name == contigName || name == 'chr' + contigName) {
+        if (name == contigName ||
+            name == 'chr' + contigName ||
+            'chr' + name == contigName) {
           return {idx: i, name: name};
         }
       }

src/pileuputils.js

@@ -0,0 +1,34 @@
+/** @flow */
+'use strict';
+
+import type * as Interval from './Interval';
+
+/**
+ * Given a list of Intervals, return a parallel list of row numbers for each.
+ * Assuming rows = pileup(reads), then your guarantee is that
+ * rows[i] == rows[j] => !reads[i].intersects(reads[j])
+ */
+function pileup(reads: Interval[]): number[] {
+  var rows = new Array(reads.length),
+      lastReads = [];  // row number --> last Interval in that row
+
+  // For each read, find the first row that it will fit in.
+  // This is potentially O(n^2) in the number of reads; there's probably a
+  // better way.
+  for (var i = 0; i < reads.length; i++) {
+    var r = reads[i];
+    var rowNum = lastReads.length;
+    for (var j = 0; j < lastReads.length; j++) {
+      if (!r.intersects(lastReads[j])) {
+        rowNum = j;
+        break;
+      }
+    }
+    rows[i] = rowNum;
+    lastReads[rowNum] = r;
+  }
+
+  return rows;
+}
+
+module.exports = {pileup};

style/pileup.css

@@ -58,6 +58,6 @@ text.basepair {
   height: 300px;
 }
 .pileup .alignment {
-  fill: red;
+  fill: #c8c8c8;  /* matches IGV */
 }
 

test/BamDataSource-test.js

@@ -54,9 +54,10 @@ describe('BamDataSource', function() {
     this.timeout(5000);
     var source = getTestSource();
 
-    var range       = new ContigInterval('20', 31512050, 31512150),
-        rangeBefore = new ContigInterval('20', 31512000, 31512050),
-        rangeAfter  = new ContigInterval('20', 31512150, 31512199);
+    // Requests are for 'chr20', while the canonical name is just '20'.
+    var range       = new ContigInterval('chr20', 31512050, 31512150),
+        rangeBefore = new ContigInterval('chr20', 31512000, 31512050),
+        rangeAfter  = new ContigInterval('chr20', 31512150, 31512199);
 
     var reads = source.getAlignmentsInRange(range);
     expect(reads).to.deep.equal([]);

test/bam-test.js

@@ -186,7 +186,7 @@ describe('BAM', function() {
       "minBlockIndex": 69454
     });
 
-    var range = new ContigInterval('20', 31511349, 31514172);
+    var range = new ContigInterval('chr20', 31511349, 31514172);
 
     bam.getAlignmentsInRange(range).then(reads => {
       expect(reads).to.have.length(1114);

test/pileuputils-test.js

@@ -0,0 +1,74 @@
+/* @flow */
+'use strict';
+
+var chai = require('chai');
+var expect = chai.expect;
+
+var _ = require('underscore');
+
+var {pileup} = require('../src/pileuputils'),
+    Interval = require('../src/Interval');
+
+describe('pileuputils', function() {
+  // This checks that pileup's guarantee is met.
+  function checkGuarantee(reads: Interval[], rows: number[]) {
+    var readsByRow = _.groupBy(reads, (read, i) => rows[i]);
+    _.each(readsByRow, reads => {
+      // No pair of reads in the same row should intersect.
+      for (var i = 0; i < reads.length - 1; i++) {
+        for (var j = i + 1; j < reads.length; j++) {
+          expect(reads[i].intersects(reads[j])).to.be.false;
+        }
+      }
+    });
+  }
+
+  it('should check the guarantee', function() {
+    var reads = [
+      new Interval(0, 9),
+      new Interval(5, 14),
+      new Interval(10, 19),
+    ];
+    checkGuarantee(reads, [0, 1, 2]);  // ok
+    checkGuarantee(reads, [0, 1, 0]);  // ok
+    expect(() => checkGuarantee(reads, [0, 0, 0])).to.throw();  // not ok
+  });
+
+  it('should pile up a collection of reads', function() {
+    var reads = [
+      new Interval(0, 9),
+      new Interval(5, 14),
+      new Interval(10, 19),
+      new Interval(15, 24)
+    ];
+    var rows = pileup(reads);
+    checkGuarantee(reads, rows);
+    expect(rows).to.deep.equal([0,1,0,1]);
+  });
+
+  it('should pile up a deep collection of reads', function() {
+    var reads = [
+      new Interval(0, 9),
+      new Interval(1, 10),
+      new Interval(2, 11),
+      new Interval(3, 12),
+      new Interval(4, 13)
+    ];
+    var rows = pileup(reads);
+    checkGuarantee(reads, rows);
+    expect(rows).to.deep.equal([0,1,2,3,4]);
+  });
+
+  it('should pile up around a long read', function() {
+    var reads = [
+      new Interval(1, 9),
+      new Interval(0, 100),
+      new Interval(5, 14),
+      new Interval(10, 19),
+      new Interval(15, 24)
+    ];
+    var rows = pileup(reads);
+    checkGuarantee(reads, rows);
+    expect(rows).to.deep.equal([0,1,2,0,2]);
+  });
+});