Calculate insert size stats and color abnormally small and large

GMOD · Oct 31, 2018 · ec1fcd3 · ec1fcd3
1 parent d7fe480
commit ec1fcd3
Show file tree

Hide file tree

Showing 3 changed files with 34 additions and 6 deletions.
diff --git a/src/JBrowse/Store/SeqFeature/BAM.js b/src/JBrowse/Store/SeqFeature/BAM.js
@@ -6,6 +6,22 @@ const { Buffer } = cjsRequire('buffer')
 const bamIndexedFilesCache = LRU(5)
 
 const BlobFilehandleWrapper = cjsRequire('../../Model/BlobFilehandleWrapper')
+const MAX_INSERT_SIZE_FOR_STATS = 30000
+
+function percentile(arr, p) {
+    if (arr.length === 0) return 0;
+    if (typeof p !== 'number') throw new TypeError('p must be a number');
+    if (p <= 0) return arr[0];
+    if (p >= 1) return arr[arr.length - 1];
+
+    var index = arr.length * p,
+            lower = Math.floor(index),
+            upper = lower + 1,
+            weight = index % 1;
+
+    if (upper >= arr.length) return arr[lower];
+    return arr[lower] * (1 - weight) + arr[upper] * weight;
+}
 
 class PairedBamRead {
     id() {
@@ -23,24 +39,21 @@ class PairedBamRead {
             return this.read1._get('name')
         } else if(field === 'pair_orientation') {
             return this.read1._get('pair_orientation')
+        } else if(field === 'template_length') {
+            return this.read1._get('template_length')
         }
     }
     pairedFeature() { return true }
     children() {}
 }
 
 class BamSlightlyLazyFeature {
-    _get_name() { return this.record._get('name') }
-    _get_start() { return this.record._get('start') }
-    _get_end() { return this.record._get('end') }
     _get_type() { return 'match'}
     _get_score() { return this.record._get('mq')}
     _get_mapping_quality() { return this.record.mappingQuality}
     _get_flags() { return `0x${this.record.flags.toString(16)}`}
     _get_strand() { return this.record.isReverseComplemented() ? -1 : 1 }
     _get_read_group_id() { return this.record.readGroupId }
-    _get_qual() { return this.record._get('qual')}
-    _get_cigar() { return this.record._get('cigar')}
     _get_seq_id() { return this._store._refIdToName(this.record._refID)}
     _get_qc_failed() { return this.record.isFailedQc()}
     _get_duplicate() { return this.record.isDuplicate()}
@@ -60,7 +73,6 @@ class BamSlightlyLazyFeature {
         ? ( this._store._refIdToName(this.record._next_refid())+':'+(this.record._next_pos()+1)) : undefined}
     _get_tags() { return this.record._tags() }
     _get_seq() { return this.record.getReadBases() }
-    _get_md() { return this.record._get('md') }
 
     constructor(record, store) {
         this.record = record
@@ -312,6 +324,14 @@ return declare( [ SeqFeatureStore, DeferredStatsMixin, DeferredFeaturesMixin, In
         })
     },
 
+    getStatsForPairCache() {
+        var total = 0
+        var stddev = 0
+        var tlens = Object.entries(this.featureCache).map(([k, v]) => Math.abs(v.get('template_length'))).filter(x => x < MAX_INSERT_SIZE_FOR_STATS).sort((a, b) => a - b)
+        var upper = percentile(tlens, 0.995)
+        var lower = percentile(tlens, 0.005)
+        return { upper, lower }
+    },
 
     _bamRecordToFeature(record) {
         return new BamSlightlyLazyFeature(record, this)

diff --git a/src/JBrowse/View/FeatureGlyph/Alignment.js b/src/JBrowse/View/FeatureGlyph/Alignment.js
@@ -80,6 +80,11 @@ return declare( [BoxGlyph,MismatchesMixin], {
                       if(Math.abs(strand) != 1 && strand != '+' && strand != '-') {
                         return track.colorForBase('reference');
                       }
+                      else if(track.upperPercentile < Math.abs(feature.get('template_length'))) {
+                        return 'red'
+                      } else if(track.lowerPercentile > Math.abs(feature.get('template_length'))) {
+                        return 'pink'
+                      }
                       else if(track.config.colorByOrientation) {
                         const type = orientationTypes[track.config.orientationType]
                         const orientation = type[feature.get('pair_orientation')]

diff --git a/src/JBrowse/View/Track/Alignments2.js b/src/JBrowse/View/Track/Alignments2.js
@@ -186,6 +186,9 @@ return declare( [ CanvasFeatureTrack, AlignmentsMixin ], {
             const min = Math.max(0, region.start - len*4)
             const max = region.end + len*4
             this.store.getFeatures({ ref: this.refSeq.name, start: min, end: max, viewAsPairs: true }, () => { /* do nothing */}, () => {
+                var stats = this.store.getStatsForPairCache()
+                this.upperPercentile = stats.upper
+                this.lowerPercentile = stats.lower
                 if(this.removeFeaturesFromCacheAfterDelay) {
                     let f = args.finishCallback
                     args.finishCallback = () => {