<|start_filename|>test/transcript-test.js<|end_filename|>
var assert = require('assert');
let data = require("./data.js")
let transcript = require("../src/seqcoverpkg/transcript.js")

let tr = new transcript.Transcript(data.tr_data)

describe('Transcript', function () {
  describe('cdsstart', function () {
    it('should equal cdsstart', function () {
        assert.equal(tr.cdsstart, data.tr_data.cdsstart)
    });
  });

  describe('transcript tests', function () {
      it('inferred genomic coords should match actual', function () {
          let parts = tr.parts()
          parts.forEach((p, i) => {
              //  u:Transcript(cdsstart: 15321505, cdsend: 15335554, chr: "X", position: @[[15319450, 15321772], [15324664, 15324871], [15325019, 15325152], [15325913, 15326046], [15331215, 15331992], [15335500, 15335554]], strand: -1, transcript: "union", txstart: 15319450, txend: 15335554)
              let s = p.hoverinfo(data.xs, data.gs)
              if(i == 0){
                  assert.equal(s.start, 15319450)
                  assert.equal(s.stop, 15321505)
              } else if (i == parts.length) {
                  assert.equal(s.start, 15335554)
                  assert.equal(s.stop, 15335554)

              } else {

                  if(s.type == transcript.FeatureType.CDS) {
                      assert(s.start >= 15321505, s.start.toString() + " " + i.toString())
                      assert(s.stop <= 15335554, s.stop.toString() + " " + i.toString())
                  }
                  if(s.type == transcript.FeatureType.EXON) {
                      assert(s.start >= 15319450, s.start.toString() + " " + i.toString())
                      assert(s.stop <= 15335554, s.stop.toString() + " " + i.toString())
                  }
              }

          })

      });
      it('traces should work', function () {
          let traces = tr.traces(0, data.xs, data.gs)
          assert.equal(traces.length, 3)
          assert.equal(traces[0].x.length, 2, "only need start and end")
          assert.equal(traces[0].y[0], 0, "y offset should be 0")
          assert.equal(traces[0].genome_x.length, 2, "only need start and end")

          assert.equal(traces[0].name, tr.name)


          var trace = traces[traces.length - 1]; // get a CDS trace
          assert(trace.text[0].startsWith("exon 6 / 6"))
          assert(trace.text[trace.text.length - 1].startsWith("exon 1 / 6"))


      })

      it('stats returns expected values', function() {
          var depths = {
              "sampleA": Array(60).fill(10),
              "sampleB": Array(60).fill(0),
          }
          var background_depths = {
              "p10": Array(60).fill(5),
              "p90": Array(60).fill(90),
          }

          background_depths.p10[4] = 500;
          //{ sampleA: { lt_background: 0, low: 0, mean: 10, median: 10 },
          //  sampleB: { lt_background: 20, low: 20, mean: 0, median: 0 } }

          let stats = tr.stats([{start:20, stop:40}], depths, background_depths, 5);
          assert.equal(stats.sampleA.median, 10)
          assert.equal(stats.sampleB.median, 0)

          assert.equal(stats.sampleA.low, 0)
          assert.equal(stats.sampleB.low, 20)

          assert.equal(stats.sampleA["lt-background"], 0)
          assert.equal(stats.sampleB["lt-background"], 20)



      })
  });
  describe('zscore tests', function () {
      it('should return correct zscore', function (){
          var vals = [2, 4, 6,7,8,9]
          var expected = [-1.6803361 , -0.84016805,  0.        ,  0.42008403,  0.84016805,
        1.26025208]

          let zs = transcript.z_transform(vals)
          for(var i=0; i <vals.length; i++){
              assert(Math.abs(zs[i] - expected[i]) < 0.0001, "expected:" + expected[i] + " obs:" + zs[i])
          }

      })

  })
});


<|start_filename|>src/seqcoverpkg/transcript.js<|end_filename|>
"use strict";

function binary_search(A, v) {
    var result = 0;
    var j = A.length
    while (j != 0) {
        let step = j >> 1
        let pos = result + step;
        if (A[pos] < v) {
            result = pos + 1
            j -= step + 1
        } else {
            j = step
        }

    }
    return result
}

// enum
const FeatureType = Object.freeze({
    EXON: "exon",
    CDS: "CDS",
    UTR: "UTR",
    TRANSCRIPT: "transcript"
})

const aesthetics = {
    TRANSCRIPT_COLOR: "rgba(65, 65, 65, 0.6)",
    TRANSCRIPT_WIDTH: 4,
    EXON_COLOR: "rgba(105,105,105, 0.6)",
    EXON_WIDTH: 19,
    CDS_COLOR: 'rgb(195, 155, 155)',
    CDS_WIDTH: 12,
}

class Feature {
    constructor(start, stop, type, transcript) {
        this.start = start
        this.stop = stop
        this.type = type
        this.transcript = transcript
    }
    get coding() {
        return this.type == FeatureType.CDS
    }


    hoverinfo(xs, gs) {
        // xs is plot coordinates, gs is genome coordinates
        // get genomic coordinate by finding index of plot-coord
        // and then looking it up in genomic array
        var start = gs[binary_search(xs, this.start)]
        var stop = gs[binary_search(xs, this.stop)]
        return {
            "transcript": this.transcript.data.transcript,
            "strand": this.transcript.strand,
            "start": start,
            "stop": stop,
            "type": this.type.toString()
        }
    }
}

class Transcript {
    constructor(data) {
        this.data = data
    }
    get cdsstart() {
        return this.data.cdsstart
    }
    get cdsend() {
        return this.data.cdsend
    }
    get chr() {
        return this.data.chr
    }

    get exons() {
        return this.data.position
    }

    get position() {
        return this.data.position
    }
    get strand() {
        return this.data.strand == -1 ? "-" : "+"
    }

    get name() {
        return this.data.transcript
    }

    get txstart() {
        return this.data.txstart
    }
    get txend() {
        return this.data.txend
    }

    parts() {
        // return CDS,exon,UTR in an array. exon and CDS are potentially
        // (often) duplicated.
        var that = this
        var result = []
        result.push(new Feature(this.data.txstart, this.data.cdsstart, FeatureType.UTR, that))
        this.data.position.forEach((exon, i) => {
            result.push(new Feature(exon[0], exon[1], FeatureType.EXON, that))
            if (exon[1] < this.data.cdsstart || exon[0] > this.data.cdsend) {
                // continue
            } else {
                result.push(new Feature(Math.max(this.data.cdsstart, exon[0]), Math.min(this.data.cdsend, exon[1]), FeatureType.CDS, that))
            }

        })
        result.push(new Feature(this.data.cdsend, this.data.txend, FeatureType.UTR, that))
        return result.filter(f => f.stop - f.start > 0)

    }

    overlaps(position) {
        // return parts() that overlap with this position
        let that = this
        var result = []
        if (position < this.txstart || position > this.txend) { return result }
        if (position < this.cdsstart) { result.push(new Feature(this.data.txstart, this.data.cdsstart, FeatureType.UTR, that)) }
        this.data.position.forEach((exon, i) => {
            if (exon[0] > position || exon[1] < position) {
                return
            }
            result.push(new Feature(exon[0], exon[1], FeatureType.EXON, that))
            if (exon[1] < that.data.cdsstart || exon[0] > that.data.cdsend) {
                // continue
            } else {
                var f = new Feature(Math.max(this.data.cdsstart, exon[0]), Math.min(this.data.cdsend, exon[1]), FeatureType.CDS, that)
                if (f.stop - f.start > 0) {
                    result.push(f)
                }
            }
        })
        if (position >= this.cdsend) {
            result.push(new Feature(this.data.cdsend, this.data.txend, FeatureType.UTR, that))
        }
        return result;
    }

    traces(y_offset, xs, gs) {

        function get_genomic_coord(x) {
            return isNaN(x) ? NaN : gs[binary_search(xs, x)]
        }

        var transcript_trace = {
            name: this.data.transcript, x: [this.data.txstart, this.data.txend], y: [y_offset, y_offset],
            type: "scatter", mode: "lines", showlegend: false,
            hoverinfo: "none",
            line: { color: aesthetics.TRANSCRIPT_COLOR, width: aesthetics.TRANSCRIPT_WIDTH }
        }

        let parts = this.parts()

        var exon_trace = {
            name: this.data.transcript + " exons", x: [], y: [], text:[],
            type: "scatter", mode: "lines", showlegend: false,
            hoverinfo: "text",
            line: { color: aesthetics.EXON_COLOR, width: aesthetics.EXON_WIDTH }
        }
        let exons = parts.filter(p => p.type == FeatureType.EXON)
        exons.forEach((e, i) => {
            if ((exon_trace.x.length) > 0) {
                exon_trace.x.push(NaN)
                exon_trace.y.push(y_offset)
                exon_trace.text.push(undefined)
            }
            var iex = this.data.strand == -1 ? (exons.length - i) : (i + 1)

            let txt = `exon ${iex} / ${exons.length}`;
            exon_trace.text.push(txt, txt)
            exon_trace.x.push(e.start, e.stop)
            exon_trace.y.push(y_offset, y_offset)
        })

        var cds_trace = {
            name: this.data.transcript + " CDS", x: [], y: [], text:[],
            type: "scatter", mode: "lines", showlegend: false,
            hoverinfo: "text",
            hovermode:"closest-x",
            line: { color: aesthetics.CDS_COLOR, width: aesthetics.CDS_WIDTH }
        }
        parts.filter(p => p.type == FeatureType.CDS).forEach(c => {
            if ((cds_trace.x.length) > 0) {
                cds_trace.x.push(NaN)
                cds_trace.y.push(y_offset)
                cds_trace.text.push(undefined)
            }
            // index back into exon array even for CDS hover
            var ei = 0
            for(var e of exons) {
                ei += 1
                if(e.start <= c.start && e.stop >= c.stop) {
                    break;
                }
            }

            var iex = this.data.strand == -1 ? (exons.length - ei + 1) : (ei + 1)
            let txt = `exon ${iex} / ${exons.length} (CDS)`
            cds_trace.text.push(txt, txt)
            cds_trace.x.push(c.start, c.stop)
            cds_trace.y.push(y_offset, y_offset)
        })

        var result = [transcript_trace, exon_trace, cds_trace]
        result.forEach(trace => {
            trace.genome_x = trace.x.map(x => get_genomic_coord(x))
        })
        return result

    }

    stats(xranges, depths, background_depths, low_depth_cutoff) {
        // NOTE xranges is in reduced, plot coords, not genomic coords and it
        // has form [{start: 23, stop: 44}, ...]
        var result = {}
        var background_low;
        // handle missing or undefined backgrounds
        if (background_depths != undefined && background_depths != {} && background_depths != []) {
            // background depths might hvae string keys like p5, p95. so we sort to
            // get the lowest value.
            var lo_key = Object.keys(background_depths).sort((a, b) => {
                return parseInt(a.replace(/^\D+/g, "")) - parseInt(b.replace(/^\D+/g, ""))
            })[0]
            background_low = background_depths[lo_key] //.slice(xstart, xstop)
        }
        var H = Array(16384);
        for (var sample in depths) {
            var S = 0; var N = 0; var lo = 0; var bg_lo = 0;
            H.fill(0);
            var dps = depths[sample];

            for (var rng of xranges) {
                for (var i = rng.start; i < rng.stop; i++) {
                    let d = dps[i]
                    if (d < 0 || isNaN(d)) { continue; }
                    lo += (d < low_depth_cutoff);
                    if (background_low != undefined) {
                        bg_lo += d < background_low[i];
                    }
                    H[Math.min(d, H.length - 1)] += 1;
                    S += d;
                    N += 1
                }
            }
            result[sample] = { "low": lo, "lt-background": bg_lo }
            var mid = N * 0.5; // 50% of data below this number of samples.
            var j = 0
            var nc = H[j]
            while (nc < mid) {
                j += 1
                nc += H[j]
            }
            result[sample]["mean"] = S / N
            result[sample]["median"] = j
        }
        return result;
    }
}

function z_transform(data) {
    // TODO: use 1-pass for mean, sd
    var s = 0
    for(let d of data) s += d;
    let mean = s / data.length;

    var sd = 0;
    for(let d of data) { sd += Math.pow(d - mean, 2); }
    sd = Math.sqrt(sd / data.length);

    return data.map(d => (d - mean) / sd)
}

try {
    // node.js stuff to allow testing
    exports.Transcript = Transcript
    exports.FeatureType = FeatureType
    exports.z_transform = z_transform


    if (require.main === module) {
        // xs and gs data for testing.
        let data = require("./test/data.js")
        let tr = new Transcript(data.tr_data)

        tr.parts().forEach(p => console.log(p.hoverinfo(data.xs, data.gs)))

        console.log(tr.traces(0))

    }
} catch (e) {
    // browser
}


<|start_filename|>src/seqcoverpkg/seqcover.js<|end_filename|>
var nan = NaN; // hack to support json dumped with NaN values.
const gene_plot_height = 500
const color_list = [
    "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b",
    "#e377c2", "#7f7f7f", "#bcbd22", "#17becf", "#7CB5EC", "#434348",
    "#90ED7D", "#F7A35C", "#8085E9", "#F15C80", "#E4D354", "#2B908F",
    "#F45B5B", "#91E8E1", "#4E79A7", "#F28E2C", "#E15759", "#76B7B2",
    "#59A14F", "#EDC949", "#AF7AA1", "#FF9DA7", "#9C755F", "#BAB0AB",
]
const HOVER_TEMPLATE = '<b>position</b>:%{text}<br><b>depth</b>:%{y}<br>(debug) x: %{x}';


// https://stackoverflow.com/a/4198132
function getHashParams(h=false) {

    var hashParams = {};
    var e,
        a = /\+/g,  // Regex for replacing addition symbol with a space
        r = /([^&;=]+)=?([^&;]*)/g,
        d = function (s) { return decodeURIComponent(s.replace(a, " ")); },
        q = h ? h : window.location.hash.substring(1)

    while (e = r.exec(q))
       hashParams[d(e[1])] = d(e[2]);

    return hashParams;
}

function setHashParams(obj) {
    var old = getHashParams();
    for(var k in obj){
        if(obj[k] == null){
            delete old[k];
        } else {
          old[k] = obj[k]
        }
    }
    var str = "#";
    for(var k in old){
        str += `${k}=${old[k]}&`
    }

    window.location.hash = str.slice(0, -1); // drop final '&'
    return old;

}

window.onhashchange = function(event) {
    let old_hp = getHashParams(new URL(event.oldURL).hash.substring(1))
    let new_hp = getHashParams()
}

function update_selected_gene_header(gene) {
    let g = plot_data.filter(g => g.symbol == gene)[0]
    if(g == undefined) g = plot_data[0];
    let chrom = g.transcripts[0].chr
    let start = g.plot_coords.g[0]
    let stop = g.plot_coords.g[g.plot_coords.g.length - 1]
    document.getElementById("header-selected-gene").innerHTML = `${gene} <i>${chrom}:${start}-${stop}</i>`
}

function update_header_metadata(n_samples, n_genes) {
    document.getElementById("sample-count").innerHTML = n_samples
    document.getElementById("gene-count").innerHTML = n_genes
}

// https://stackoverflow.com/questions/8069315/create-array-of-all-integers-between-two-numbers-inclusive-in-javascript-jquer/8069367
function range(start, end, step = 1) {
    const len = Math.floor((end - start) / step) + 1
    return Array(len).fill().map((_, idx) => start + (idx * step))
}

function get_gene_plot_layout(gene) {
    let step = Math.round(gene.plot_coords.x.length / 10)
    let layout = {
        grid: { rows: 3, columns: 1, },
        autosize: true,
        height: gene_plot_height,
        margin: { t: 10, r: 0, b: 30 },
        hovermode: 'x',
        xaxis: {
            tickvals: range(gene.plot_coords.x[0], gene.plot_coords.x[gene.plot_coords.x.length - 10], step),
            ticktext: range(gene.plot_coords.g[0], gene.plot_coords.g[gene.plot_coords.x.length - 10], step),
            showgrid: false,
            ticklen: 5,
        },
        yaxis: { title: "Depth", domain: [0.28, 1]},
        yaxis2: {
            range: [0.5, -0.5],
            showlegend: false,
            zeroline: false,
            tickvals: [0],
            ticktext: ["Merged<br>Transcript"],
            tickangle: 40,
            domain: [0.25, 0.28],
        },
        yaxis3: {
            range: [0, 2],
            zeroline: false,
            showlegend: false,
            domain: [0.0, 0.25],
            tickangle: 40,
            ticktext: gene.transcripts.map(t => t.data.transcript),
        },
        hoverdistance: 100000,
        hoverinfo: "none",
        showlegend: false,
    }

    return (layout)
}

//Get a by position trace per sample of depth
function get_depth_trace(gene) {

    var traces = [];
    var low_lvl = 10000;
    var low_dp;
    for (var slvl in gene.plot_coords.background_depths) {
        let lvl = slvl.replace(/^\D+/g, "")

        let dp = gene.plot_coords.background_depths[slvl].map(function (v) { return v < -1000 ? NaN : v })
        if (lvl < low_lvl) { low_dp = dp; low_lvl = lvl }

        var trace = {
            x: gene.plot_coords.x,
            text: gene.plot_coords.g,
            y: dp,
            type: "scatter",
            mode: "lines",
            name: `Percentile: ${lvl}`,
            tracktype: "background",
            line: { width: 1, /* dash: "dot", */ color: 'rgb(200, 200, 200)' },
            yaxis: "y",
            hoverinfo: "skip",
        }
        traces.push(trace)
    }

    var i = -1;
    var samples_ = []

    for (var sample in gene.plot_coords.depths) {
        samples_.push(sample)
    }
    for (var sample of samples_.sort()) {


        i += 1
        var dp = gene.plot_coords.depths[sample]
        dp = dp.map(function (v) { return v < -1000 ? NaN : v })
        var color = color_list[i % color_list.length];
        var trace = {
            x: gene.plot_coords.x, text: gene.plot_coords.g, y: dp,
            type: 'scattergl', mode: 'lines', name: sample, line: { width: 0.36, color: color_list[i % color_list.length] },
            hovertemplate: HOVER_TEMPLATE,
            hoverdistance: 1000,
            hoverinfo: "none",
            yaxis: "y",
            tracktype: "sample",
        };

        traces.push(trace);

        // extract parts of this sample that are below the threshold and plot
        // in a similar trace with wider line.
        if (low_dp) {
            var low_trace = {
                x: [], y: [], text: [], type: 'scatter', mode: "lines", name: sample, line: { width: 3.0, color: color_list[i % color_list.length] },
                hoverinfo: "none",
                connectgaps: false,
                yaxis: "y",
            };

            dp.forEach((d, i) => {
                if (d < low_dp[i]) {
                    low_trace.x.push(gene.plot_coords.x[i])
                    low_trace.y.push(d)
                } else {
                    if (low_trace.x.length > 0 && !isNaN(low_trace.x[low_trace.x.length - 1])) {
                        low_trace.x.push(gene.plot_coords.x[i])
                        low_trace.y.push(null)
                    }
                }

            })
            traces.push(low_trace)

        }
    };

    return (traces)

};

//Get the transcript shapes for the current region
function get_transcript_traces(gene, transcripts, plotRef) {
    var traces = []
    // use negative values for y_offset so we can tell from y-value if we
    // are in depth plot (which must be positive) or in transcript plot

    transcripts.forEach((transcript, y_offset) => {
        var trs = transcript.traces(-(y_offset), gene.plot_coords.x, gene.plot_coords.g)
        trs.forEach(t => {
            t.yaxis = plotRef
            traces.push(t)
        })
    })
    return traces
};


function plot_per_base_depth(gene) {
    setHashParams({'gene': gene.symbol})
    update_selected_gene_header(gene.symbol)
    let gene_layout = get_gene_plot_layout(gene)
    let depth_traces = get_depth_trace(gene)
    let unioned_transcript_traces = get_transcript_traces(gene, [gene.unioned_transcript], "y2")
    unioned_transcript_traces.forEach(t => depth_traces.push(t))

    let transcript_traces = get_transcript_traces(gene, gene.transcripts, "y3")
    transcript_traces.forEach(t => depth_traces.push(t))
    // each transcript is centered on the negative integers.
    gene_layout.yaxis3.range = [0.5, -(1 + transcript_traces.length / 3)]
    gene_layout.yaxis3.tickvals = [...Array(gene.transcripts.length).keys()].map(f => -f)
    gene_layout.yaxis3.ticktext = gene.transcripts.map(t => t.data.transcript)

    let d = document.getElementById("gene_plot")
    let p = Plotly.newPlot(d, depth_traces, gene_layout, {responsive: true})
    d.on("plotly_hover", data => {
        // handle_hover(data, depth_traces, gene, gene_layout)
        Plotly.react("gene_plot", depth_traces, gene_layout)
        data.event.stopPropagation()
    }).on("plotly_unhover", data => {
        if (gene_layout.shapes != undefined && gene_layout.shapes.length > 0) {
            gene_layout.shapes.pop()
        }
        Plotly.react("gene_plot", depth_traces, gene_layout)
        data.event.stopPropagation()
    })
    return p
}

function mean(data) {
    var result = 0
    var n = 0
    for (var i = 0; i < data.length; i++) {
        let d = data[i]
        if (d >= 0) {
            n += 1
            result += d
        }
    }
    return result / n
}

function handle_hover(data, depth_traces, gene, gene_layout) {
    //  this function handles hover events in the transcript plots. it
    //  finds the correct transcript and feature (exon/transcript/UTR, etc)
    //  and calculates stats, and updates the gene_layout to draw the
    //  region shape.
    //
    let ax = data.yaxes[0]._attr;
    // don't handle hover in depth plot
    if (ax == "yaxis") { return false; }
    // this is the x in plot coordinates.
    var x = Math.round(data.xvals[0])
    var y = Math.round(Math.abs(data.yvals[0])) // can now use this as an index to get the transcript.
    var transcript = ax == "yaxis2" ? gene.unioned_transcript : gene.transcripts[y]


    if (transcript == undefined) { return false }

    let overlaps = transcript.overlaps(x)
    if (overlaps.length > 1) {
        if (overlaps.some(p => p.type == FeatureType.CDS)) {
            overlaps = overlaps.filter(p => p.type == FeatureType.CDS)
        } else if (overlaps.some(p => p.type == FeatureType.EXON)) {
            overlaps = overlaps.filter(p => p.type == FeatureType.EXON)
        }
    }
    if (overlaps.length == 0) {
        overlaps.push(new Feature(transcript.txstart, transcript.txend, FeatureType.TRANSCRIPT))
    }
    if (gene_layout.shapes == undefined) { gene_layout.shapes = [] }
    gene_layout.shapes.push({
        type: "rect",
        xref: "x",
        yref: "paper",
        y0: 0,
        y1: 1,
        fillcolor: "#cccccc",
        opacity: 0.2,
        x0: overlaps[0].start,
        x1: overlaps[0].stop,
    })

    var start_idx = binary_search(gene.plot_coords.x, overlaps[0].start)
    var stop_idx = binary_search(gene.plot_coords.x, overlaps[0].stop)

    let low_depth_cutoff = 7; // TODO: get this from user-form input.
    var selection_stats = transcript.stats([{ start: start_idx, stop: stop_idx }], gene.plot_coords.depths, gene.plot_coords.background_depths, low_depth_cutoff)

    var tx_stats = transcript.stats([{ start: 0, stop: gene.plot_coords.x.length }], gene.plot_coords.depths, gene.plot_coords.background_depths, low_depth_cutoff)

    // example of how to get stats for all CDSs. NOTE: this should only be done
    // once per gene and result cached.
    var cds = transcript.parts().filter(p => p.type == FeatureType.CDS)
    var cds_stats = transcript.stats(cds, gene.plot_coords.depths, gene.plot_coords.background_depths, low_depth_cutoff)

    var gstart = gene.plot_coords.g[start_idx]
    var gstop = gene.plot_coords.g[stop_idx]

    var columns = [{ title: "sample" }, { title: "transcript mean" }, { title: "CDS mean" }, { title: "selection mean" },
    { title: "transcript bases < lower" }, { title: "CDS bases < lower" }, { title: "selection bases < lower" }];
    var table = []
    var samples_ = [];

    //var means = {}
    //hoverInfo.innerHTML = `${gene.unioned_transcript.chr}:${gstart}-${gstop}<br><ul>`
    for (var sample in gene.plot_coords.depths) { samples_.push(sample) }
    for (var sample of samples_.sort()) {

        var row = [sample, tx_stats[sample].mean.toFixed(2), cds_stats[sample].mean.toFixed(2), selection_stats[sample].mean.toFixed(2), tx_stats[sample].low, cds_stats[sample].low, selection_stats[sample].low]
        table.push(row)
        //let depths = gene.plot_coords.depths[sample];
        //means[sample] = mean(depths.slice(start_idx, stop_idx)).toPrecision(4)
    }
    if (datatable) {
        datatable.clear()
        datatable.rows.add(table)
        datatable.draw()
    } else {
        datatable = jQuery('table#stats_table').DataTable({ data: table, columns: columns })
    }
}

function highlight_sample(sample) {
    setHashParams({'sample': sample})
    let d = document.getElementById("gene_plot")
    let vals = d.data.map((t, i) => {
        if (t.tracktype == "background") {
            return [1, 1.5]
        }
        if (t.tracktype != 'sample') {
            return [undefined, undefined]
        } else {
            if (t.name == sample) {
                return [1, 1.8]
            } else {
                return [0.15, 0.36]
            }
        }
    })

    Plotly.restyle(d, {'opacity': vals.map(i => i[0]), 'line.width': vals.map(i => i[1]), 'hovertemplate': vals.map(i => i[1] > 0.8 ? HOVER_TEMPLATE: null)})
}

function tie_heatmap_to_line_plot() {
    // clicks on heatmap x-axis labels highlights sample in line plot
    // https://issue.life/questions/44297012
    var xticks = jQuery("#heatmap_plot .xaxislayer-above > .xtick > text")
    xticks.css({ "cursor": "pointer" })
    xticks.attr('class', 'unselected_label')
    var d = document.getElementById("gene_plot")

    xticks.each(function (i) {
        var item = jQuery(this);
        item.css({ "fill": color_list[i % color_list.length] })
        item.attr('pointer-events', 'all')
        item.on("click", function (e) {
            var n = jQuery(this)
            var undo = n.css("font-weight") == "800";
            xticks.attr('class', 'unselected_label')
            var sample = n.text()
            var vals = null
            if (undo) {
                vals = d.data.map((t, i) => [1, (t.tracktype == 'sample' || t.tracktype == 'background') ? 0.36 : undefined])
                Plotly.restyle(d, {'opacity': vals.map(i => i[0]), 'line.width': vals.map(i => i[1]), 'hovertemplate': d.data.map((t) => t.tracktype == "sample" ? HOVER_TEMPLATE: null)})
                setHashParams({'sample': null})
            } else {
                n.attr('class', 'selected_label')
                highlight_sample(sample)
            }
            // see: https://codepen.io/etpinard/pen/RQQqzq?editors=0010
        }).on("unhover", function (e) {
        })
    })
}

function draw_heatmap() {

    let sample_z_scale = jQuery('#z_scale_samples').is(":checked")
    let gene_z_scale = jQuery('#z_scale_genes').is(":checked")

    let stat_metric = jQuery("#metric_select").val()

    var low_depth_cutoff = 7;
    var z = []
    var x = [] // samples
    var y = [] // genes

    var is_cds = stat_metric.startsWith("cds")
    var metric = stat_metric.split("_")[1]
    for (var i = 0; i < plot_data.length; i++) {
        var g = plot_data[i]
        y.push(g.symbol)

        var transcript = g.unioned_transcript
        var stats = null
        if (is_cds) {
            var cds = transcript.parts().filter(p => p.type == FeatureType.CDS)
            stats = transcript.stats(cds, g.plot_coords.depths, g.plot_coords.background_depths, low_depth_cutoff)
        } else {
            stats = transcript.stats([{ start: 0, stop: g.plot_coords.x.length }], g.plot_coords.depths, g.plot_coords.background_depths, low_depth_cutoff)
        }

        var row = []

        // get samples in alphabetical order
        if (i == 0) {
            for(var s in stats) { x.push(s) }
            x = x.sort()
        }

        for(var s of x) {
            row.push(stats[s][metric])
        }
        // row is a set of samples for a given metric
        if(gene_z_scale) {
            row = z_transform(row)
        }

        z.push(row)
    }

    if(sample_z_scale) {
        z = z_transform_columns(z)
    }

    // heatmap draws from bottom up by default
    z = z.reverse()
    y = y.reverse()

    let hlayout = {
        title: "",
        margin: { t: 10 },
        height: 20 * y.length + 60,
        width: 20 * x.length + 30,
        paper_bgcolor: '#f8f9fa',
        uniformtext: { mode: 'show', minsize: 2 },
        xaxis: { tickfont: { size: 12 }, fixedrange: true, automargin: true, type: "category"},
        yaxis: { fixedrange: true, automargin: true, },
        font: { size: 14 },
    }
    // TODO: this should be based on window width.
    if(hlayout.width < jQuery('#heatmap_plot').width()) {
        hlayout.width = jQuery('#heatmap_plot').width();
    }
    if(hlayout.height < 300) {
        hlayout.height = 300;
    }
    let heatmap_config = {
        displaylogo: false,
        toImageButtonOptions: {
            format: 'svg', // one of png, svg, jpeg, webp
            filename: 'seqcover_summary',
            height: 18 * y.length + 60,
            width: 1200,
            scale: 1 // Multiply title/legend/axis/canvas sizes by this factor
        },
        // responsive: true,
    }
    let trace = [{ x: x, y: y, z: z, type: 'heatmap', hoverongaps: false, colorscale: "Cividis" }]

    //https://plotly.com/javascript/reference/heatmap/
    let p = document.getElementById('heatmap_plot')

    let yticks = jQuery("#heatmap_plot .yaxislayer-above > .ytick > text")
    let selected = false
    yticks.each(function (i) {
        let elem = jQuery(this)
        if (elem.css("font-weight") == "800") {
            selected = i
        }
    })

    Plotly.newPlot(p, trace, hlayout, heatmap_config)
    p.removeAllListeners("plotly_click")
    p.on('plotly_click', function (click_data) {
        let selected_gene_idx = click_data.points[0].pointIndex[0]
        let gene_idx = plot_data.length - selected_gene_idx - 1
        let selected_sample_idx = click_data.points[0].pointIndex[1]
        let sample = click_data.points[0].x
        let gene = click_data.points[0].y

        // y-axis ticks
        let yticks = jQuery("#heatmap_plot .yaxislayer-above > .ytick > text")
        yticks.attr('class', 'unselected_label')
        yticks.each(function (i) {
            if (i == selected_gene_idx) {
                jQuery(this).attr('class', 'selected_label')
            }
        })
        // x-axis ticks
        let xticks = jQuery("#heatmap_plot .xaxislayer-above > .xtick > text")
        xticks.attr("class", "unselected_label")
        xticks.each(function (i) {
            if (i == selected_sample_idx) {
                jQuery(this).attr("class", "selected_label")
            }
        })

        if(selected_gene_idx != selected) {
            // only redraw if it's a new gene
            plot_per_base_depth(plot_data[gene_idx]).then(highlight_sample(sample))
            selected = selected_gene_idx;
            setHashParams({'gene': gene, 'sample': sample});
        } else {
            highlight_sample(sample)
        }

    })
    p.on('plotly_afterplot', function () {
        // initialize first gene as selected; or re-use selected (when changing metric)
        let yticks = jQuery("#heatmap_plot .yaxislayer-above > .ytick > text")
        yticks.attr('class', 'unselected_label')
        if (!selected) {
            jQuery(yticks[plot_data.length - 1]).attr("class", "selected_label")
        } else {
            jQuery(yticks[selected]).attr("class", "selected_label")
        }

        // x-axis ticks
        tie_heatmap_to_line_plot()
        // y-axis ticks
        Plotly.d3.selectAll(".yaxislayer-above").selectAll('text')
            .on("click", function (d) {
                yticks.attr("class", "unselected_label")
                jQuery(this).attr('class', 'selected_label')
                var hash = setHashParams({'gene': d.text})
                plot_per_base_depth(plot_data[plot_data.length - d.x - 1])
                if("sample" in hash){ highlight_sample(hash.sample) }
            })
    })
}

function update_popovers() {
    $('[data-toggle="tooltip"]').tooltip()

    $("#navigation-popover").popover({
        container: "body",
        html: true,
        title: "Navigating the report",
        content: `
            <ul class="m-0 ml-1 p-2">
                <li>X and Y axes labels are clickable</li>
                <li>Clicking gene names (y-axis labels) updates 'Depth per base' plot</li>
                <li>Clicking sample names (x-axis labels) highlights the sample trace</li>
                <li>Clicking heatmap values updates the selected gene and highlights the sample</li>
            </ul>
        `,
    })
    $("#depth-plot-popover").popover({
        container: "body",
        html: true,
        title: "Depths by position across gene",
        content: `
            <ul class="m-0 ml-1 p-2">
                <li>Depths per base per sample are displayed over gene length</li>
                <li>Individual transcripts are merged into the Merged transcript</li>
                <li>CDS segments are colored light red</li>
                <li>UTR segments are light gray</li>
            </ul>
        `,
    })
}

jQuery(document).ready(function () {
    update_header_metadata(Object.keys(plot_data[0].plot_coords.depths).length, plot_data.length)
    update_popovers()

    if (Object.keys(plot_data[0].plot_coords.background_depths).length == 0) {
        jQuery('#cds_lt_bg').remove()
    }

    for(var g of plot_data) {
        if (g.unioned_transcript.constructor.name == "Object") {
            g.unioned_transcript = new Transcript(g.unioned_transcript)
            g.transcripts = g.transcripts.map(t => new Transcript(t))
            for(var s in g.plot_coords.depths) {
                g.plot_coords.depths[s] = Float32Array.from(g.plot_coords.depths[s])
            }

            g.plot_coords.x = Int32Array.from(g.plot_coords.x)
            g.plot_coords.g = Int32Array.from(g.plot_coords.g)
        }
    }

    jQuery('#metric_select,.z_scale').on('change', function() { draw_heatmap() })
    jQuery('#metric_select').trigger('change')
    var p = getHashParams();
    var i = 0; // TODO: rename `i`
    if ("gene" in p) {
        var gene = p.gene;
        for(var j = 0; j < plot_data.length; j++){
            if(plot_data[j].symbol == gene){
                i = j
                break
            }
        }

    }
    plot_per_base_depth(plot_data[i]).then(function() {
        // NOTE we should put these events directly in the plot fn so we don't
        // have code duplication, but here for now...
        let yticks = jQuery("#heatmap_plot .yaxislayer-above > .ytick > text")
        yticks.attr('class', 'unselected_label')
        yticks.each(function (k) {
            if (i == plot_data.length - k - 1) {
                jQuery(this).attr('class', 'selected_label')
            }
        })
        if("sample" in p) {
            highlight_sample(p.sample)
            jQuery("#heatmap_plot .xaxislayer-above > .xtick > text").each(function(i, v) {
                if(v.innerHTML == p.sample){
                    jQuery(v).attr('class', 'selected_label')
                }
            })

        }
    })
})

function z_transform_columns(z) {
    var ztr = new Array(z.length);
    for(var rowi=0; rowi < z.length; rowi++){
        ztr[rowi] = new Array(z[0].length)
    }

    for(var coli=0; coli < z[0].length; coli++){
        var col = []
        for(var rowi=0; rowi < z.length; rowi++) {
            col.push(z[rowi][coli])
        }
        var tr = z_transform(col)
        for(var rowi=0; rowi < z.length; rowi++) {
            ztr[rowi][coli] = tr[rowi];
        }

    }
    return ztr
}
