f863c012ea39c498eaf4cff84f85f92fef210c5f
jcasper
  Mon May 13 00:51:57 2019 -0700
Initial commit of Hi-C display support via .hic files, refs #18842

diff --git src/hg/hgTracks/hicTrack.c src/hg/hgTracks/hicTrack.c
new file mode 100644
index 0000000..3d8d02f
--- /dev/null
+++ src/hg/hgTracks/hicTrack.c
@@ -0,0 +1,513 @@
+/* hic -- draw Hi-C maps */
+
+/* Copyright (C) 2018 The Regents of the University of California 
+ * See README in this or parent directory for licensing information. */
+
+#include "common.h"
+#include "obscure.h"
+#include "hgTracks.h"
+#include "bedCart.h"
+#include "bigWarn.h"
+#include "interact.h" // for interact structures
+#include "hicUi.h"
+#include "Cstraw.h"
+#include "hic.h"
+#include "htmlColor.h"
+
+
+static int hicTotalHeight(struct track *tg, enum trackVisibility vis)
+/* Calculate height of the track across all windows.  For arc and triangle, this
+ * depends on the actual interactions being drawn.  For square, it just depends
+ * on the height of the highest square (yes, this currently means smaller
+ * squares are stretched vertically in multi-region - it's not ideal). */
+{
+int maxHeight = 0;
+
+struct track *trackScan = tg;
+while (trackScan != NULL) // Canvas max height for this and previous windows ...
+    {
+    if ((int)trackScan->maxRange > maxHeight)
+        maxHeight = (int)trackScan->maxRange;
+    trackScan = trackScan->prevWindow;
+    }
+trackScan = tg->nextWindow;
+while (trackScan != NULL) // ... and look at following windows ...
+    {
+    if ((int)trackScan->maxRange > maxHeight)
+        maxHeight = (int)trackScan->maxRange;
+    trackScan = trackScan->nextWindow;
+    }
+
+// Later, try making pack 0.5, squish 0.25, and dense 0.125.  The various
+// draw functions will also have to be modified accordingly.
+if ( tg->visibility == tvDense)
+    maxHeight *= 0.5;
+return maxHeight;
+}
+
+struct hicMeta *grabHeader(struct track *tg)
+/* Fetch a hicMeta structure that describes the Hi-C data associated with
+ * the track. */
+{
+char *filename = trackDbSettingOrDefault(tg->tdb, "bigDataUrl", NULL);
+struct hicMeta *metaResult = NULL;
+if (filename == NULL)
+    {
+    warn("Missing bigDataUrl setting for track %s", tg->track);
+    return NULL;
+    }
+char *errMsg = hicLoadHeader(filename, &metaResult);
+if (errMsg != NULL)
+    {
+    tg->networkErrMsg = errMsg;
+    return NULL;
+    }
+return metaResult;
+}
+
+
+int fetchResolution(struct track *tg)
+/* Retrieve the binsize to draw the track in.  If the track UI has the resolution
+ * set as "auto", find the widest resolution that still splits the window into 5000
+ * pieces (if possible, otherwise select the smallest resolution). */
+{
+if (tg->customPt == NULL)
+    tg->customPt = grabHeader(tg);
+if (tg->customPt == 0)
+    return -1;
+
+struct hicMeta *meta = (struct hicMeta*) tg->customPt;
+char *binSizeString = hicUiFetchResolution(cart, tg->track, meta);
+if (binSizeString == NULL)
+    {
+    warn("Empty binSize string");
+    return 0;
+    }
+int result;
+if (sameString(binSizeString, "Auto"))
+    {
+    int idealRes = (winEnd-winStart)/5000;
+    char *autoRes = meta->resolutions[meta->nRes-1];
+    int i;
+    for (i=meta->nRes-1; i>= 0; i--)
+        {
+        if (atoi(meta->resolutions[i]) < idealRes)
+            {
+            autoRes = meta->resolutions[i];
+            break;
+            }
+        }
+    result = atoi(autoRes);
+    }
+else
+    result = atoi(binSizeString);
+return result;
+}
+
+char *fetchNormalization(struct track *tg)
+/* Fetch the normalization to use for this track. */
+{
+if (tg->customPt == NULL)
+    tg->customPt = grabHeader(tg);
+if (tg->customPt == NULL)
+    return NULL;
+struct hicMeta *meta = (struct hicMeta*) tg->customPt;
+return hicUiFetchNormalization(cart, tg->track, meta);
+}
+
+
+static void loadAndFilterItems(struct track *tg)
+/* Load all Hi-C items in the current region and identify the window height
+ * and median value for this region. */
+{
+int *x, *y, numRecords;
+double *counts;
+int binSize = fetchResolution(tg);
+if (binSize == 0)
+    return;
+struct hicMeta *meta = (struct hicMeta*)tg->customPt;
+char *normalization = fetchNormalization(tg);
+
+// Later, we should validate that this file is for the current assembly (see the hicMeta structure)
+
+// Note: This is giving it a 0-based, full-closed window. Straw seems to use 1-based coordinates
+// by default, but accepts 0 as the start of a window without complaint.
+struct dyString *windowPos = dyStringNew(0);
+// Pad the start because we want to display partial interactions if the end of an interacting block is
+// in view but not the start.  Straw won't report interactions if the start of the block isn't in the
+// supplied position range.
+int strawStart = winStart - binSize + 1;
+if (strawStart < 0) strawStart = 0;
+
+char *hicChromName = chromName;
+if (meta->ucscToAlias != NULL)
+    {
+    hicChromName = (char*) hashFindVal(meta->ucscToAlias, chromName);
+    if (hicChromName == NULL)
+        {
+        hicChromName = chromName;
+        }
+    }
+dyStringPrintf(windowPos, "%s:%d:%d", hicChromName, strawStart, winEnd-1);
+
+char *filename = trackDbSettingOrDefault(tg->tdb, "bigDataUrl", NULL);
+if (filename == NULL)
+    return;
+
+tg->networkErrMsg = Cstraw(normalization, filename, binSize, dyStringContents(windowPos), dyStringContents(windowPos), "BP", &x, &y, &counts, &numRecords);
+
+// Using the interact structure because it has convenient fields, but this is not interact data and
+// shouldn't be passed to those functions.
+struct interact *hicItems = NULL;
+int i = 0, filtNumRecords = 0;
+tg->maxRange = 0.0; // the max height of an interaction in this window
+double *countsCopy = (double*) calloc(numRecords, sizeof(double));
+for (i=0; i<numRecords; i++)
+    {
+    char *drawMode = hicUiFetchDrawMode(cart, tg->track);
+    if (isnan(counts[i]))
+        {
+        // Yes, apparently NAN is possible with normalized values in some methods.  Ignore those.
+        continue;
+        }
+    if (sameString(drawMode, HIC_DRAW_MODE_ARC))
+        {
+        // we omit self-interactions in arc mode (they'd just be weird vertical lines)
+        if (x[i] == y[i])
+            continue;
+        }
+    countsCopy[filtNumRecords++] = counts[i];
+    struct interact *new = NULL;
+    AllocVar(new);
+    new->chrom = cloneString(chromName);
+    // x is always <= y, so x is always the start
+    new->chromStart = x[i];
+    new->chromEnd = y[i]+binSize;
+    new->name = NULL;
+    new->score = 0; // ignored
+    new->value = counts[i];
+    new->exp = NULL;
+    new->color = 0;
+    new->sourceChrom = cloneString(chromName);
+    new->sourceStart = x[i];
+    new->sourceEnd = x[i]+binSize;
+    new->targetChrom = cloneString(chromName);
+    new->targetStart = y[i];
+    new->targetEnd = y[i]+binSize;
+    new->targetName = new->targetStrand = NULL;
+
+    // Calculate the track draw height required to see this item
+    int leftx = new->chromStart > winStart ? new->chromStart : winStart;
+    int rightx = new->chromEnd < winEnd ? new->chromEnd : winEnd;
+    double thisHeight = scaleForWindow(insideWidth, winStart, winEnd)*(rightx - leftx)/2.0; // triangle or arc
+    if (sameString(drawMode,HIC_DRAW_MODE_SQUARE))
+        thisHeight = scaleForWindow(insideWidth, winStart, winEnd)*(winEnd-winStart); // square - always draw the full square
+
+    if (thisHeight > tg->maxRange)
+        tg->maxRange = thisHeight;
+
+    slAddHead(&hicItems, new);
+    }
+
+// Heuristic for auto-scaling the color gradient based on the scores in view - draw the max color value
+// at or above 2*median score.
+tg->graphUpperLimit = 2.0*doubleMedian(filtNumRecords, countsCopy);
+free(countsCopy);
+tg->items = hicItems;
+}
+
+void hicLoadItems(struct track *tg)
+/* Load Hi-C items in (mostly) interact format */
+{
+char *filename = trackDbSettingOrDefault(tg->tdb, "bigDataUrl", NULL);
+if (filename == NULL)
+    return;
+tg->customPt = grabHeader(tg);
+if (tg->customPt == NULL)
+    return;
+loadAndFilterItems(tg);
+}
+
+
+Color *ColorSetForHic(struct hvGfx *hvg, struct track *tg, int bucketCount)
+/* Create the gradient color array for drawing a Hi-C heatmap */
+{
+struct rgbColor rgbLow;
+char *lowColorText = hicUiFetchBgColor(cart, tg->track); // This is an HTML color like #ffed02
+unsigned lowRgbVal = 0;
+if (!htmlColorForCode(lowColorText, &lowRgbVal))
+{
+    warn("Bad RGB background color value %s for track %s", lowColorText, tg->track);
+    return NULL;
+}
+int r, g, b;
+htmlColorToRGB(lowRgbVal, &r, &g, &b);
+
+rgbLow.r=(unsigned char)r;
+rgbLow.g=(unsigned char)g;
+rgbLow.b=(unsigned char)b;
+
+struct rgbColor rgbHigh;
+char *highColorText = hicUiFetchDrawColor(cart, tg->track); // This is an HTML color like #ffed02
+unsigned highRgbVal = 0;
+if (!htmlColorForCode(highColorText, &highRgbVal))
+{
+    warn("Bad RGB color value %s for track %s", highColorText, tg->track);
+    return NULL;
+}
+htmlColorToRGB(highRgbVal, &r, &g, &b);
+rgbHigh.r=(unsigned char)r;
+rgbHigh.g=(unsigned char)g;
+rgbHigh.b=(unsigned char)b;
+Color *colorIxs = (Color*) calloc (bucketCount, sizeof(Color));
+hvGfxMakeColorGradient(hvg, &rgbLow, &rgbHigh, bucketCount, colorIxs);
+return colorIxs;
+}
+
+
+double getHicMaxScore(struct track *tg)
+/* Return the score at which we should reach the maximum intensity color. */
+{
+if (hicUiFetchAutoScale(cart, tg->track))
+    return tg->graphUpperLimit;
+else
+    return hicUiFetchMaxValue(cart, tg->track);
+}
+
+
+static void drawHicTriangle(struct track *tg, int seqStart, int seqEnd,
+        struct hvGfx *hvg, int xOff, int yOff, int width, 
+        MgFont *font, Color color, enum trackVisibility vis)
+/* Draw a list of Hi-C interactions in a triangle */
+{
+double xScale = scaleForWindow(width, seqStart, seqEnd);
+double yScale = xScale;
+int maxHeight = tg->height;
+struct interact *hicItem = NULL;
+int binSize = fetchResolution(tg);
+if (binSize == 0)
+    return;
+
+if (vis == tvDense)
+    {
+    yScale *= 0.5;
+    }
+
+double maxScore = getHicMaxScore(tg);
+Color *colorIxs = ColorSetForHic(hvg, tg, HIC_SCORE_BINS+1);
+if (colorIxs == NULL)
+    return; // something went wrong with colors
+
+for (hicItem = (struct interact *)tg->items; hicItem; hicItem = hicItem->next)
+    {
+    int leftStart = hicItem->sourceStart - winStart;
+    int leftEnd = leftStart + binSize;
+    int rightStart = hicItem->targetStart - winStart;
+    int rightEnd = rightStart + binSize;
+    if (leftStart < 0)
+        leftStart = 0;
+    if (leftEnd > winEnd-winStart)
+        leftEnd = winEnd-winStart;
+    if (rightStart < 0)
+        rightStart = 0;
+    if (rightEnd > winEnd-winStart)
+        rightEnd = winEnd-winStart;
+
+    int colorIx;
+    if (hicItem->value > maxScore)
+        colorIx = colorIxs[HIC_SCORE_BINS];
+    else
+        colorIx = colorIxs[(int)(HIC_SCORE_BINS * hicItem->value/maxScore)];
+    // adding four polygon points for a diamond, based on the starts and ends of the source and target coordinate ranges.
+    struct gfxPoly *diamond = gfxPolyNew();
+
+    // left point of the diamond
+    double x = xScale * (leftStart+rightStart)/2.0;
+    double y = yScale * (rightStart-leftStart)/2.0;
+    gfxPolyAddPoint(diamond, (int)x+xOff, maxHeight-(int)y+yOff);
+
+    // top point of the diamond
+    x = xScale * (leftStart+rightEnd)/2.0;
+    y = yScale * (rightEnd-leftStart)/2.0;
+    gfxPolyAddPoint(diamond, (int)x+xOff, maxHeight-(int)y+yOff);
+
+    // right point of the diamond
+    x = xScale * ((leftEnd + rightEnd)/2.0);
+    y = yScale * (rightEnd-leftEnd)/2.0;
+    gfxPolyAddPoint(diamond, (int)x+xOff, maxHeight-(int)y+yOff);
+
+    // bottom point of the diamond
+    x = xScale * (leftEnd+rightStart)/2.0;
+    y = yScale * (rightStart-leftEnd)/2.0;
+    gfxPolyAddPoint(diamond, (int)x+xOff, maxHeight-(int)y+yOff);
+    hvGfxDrawPoly(hvg, diamond, colorIx, TRUE);
+    gfxPolyFree(&diamond);
+    }
+}
+
+
+
+static void drawHicSquare(struct track *tg, int seqStart, int seqEnd,
+        struct hvGfx *hvg, int xOff, int yOff, int width, 
+        MgFont *font, Color color, enum trackVisibility vis)
+/* Draw a Hi-C heatmap as a square */
+{
+double xScale = scaleForWindow(width, seqStart, seqEnd);
+double yScale = xScale;
+int maxHeight = tg->height;
+struct interact *hicItem = NULL;
+int binSize = fetchResolution(tg);
+if (binSize == 0)
+    return;
+
+if (vis == tvDense)
+    {
+    yScale *= 0.5;
+    }
+
+double maxScore = getHicMaxScore(tg);
+Color *colorIxs = ColorSetForHic(hvg, tg, HIC_SCORE_BINS+1);
+if (colorIxs == NULL)
+    return; // something went wrong with colors
+
+for (hicItem = (struct interact *)tg->items; hicItem; hicItem = hicItem->next)
+    {
+    int leftStart = hicItem->sourceStart - winStart;
+    int leftEnd = leftStart + binSize;
+    int rightStart = hicItem->targetStart - winStart;
+    int rightEnd = rightStart + binSize;
+    if (leftStart < 0)
+        leftStart = 0;
+    if (leftEnd > winEnd-winStart)
+        leftEnd = winEnd-winStart;
+    if (rightStart < 0)
+        rightStart = 0;
+    if (rightEnd > winEnd-winStart)
+        rightEnd = winEnd-winStart;
+
+    int colorIx;
+    if (hicItem->value > maxScore)
+        colorIx = colorIxs[HIC_SCORE_BINS];
+    else
+        colorIx = colorIxs[(int)(HIC_SCORE_BINS * hicItem->value/maxScore)];
+    double x = xScale * leftStart;
+    double y = yScale * ((winEnd-winStart)-rightStart);
+    hvGfxBox(hvg, (int)x+xOff, maxHeight-(int)y+yOff, (int)(xScale*(leftEnd-leftStart))+1, (int)(yScale*(rightEnd-rightStart))+1, colorIx);
+    if (leftStart != rightStart)
+        {
+        x = xScale * rightStart;
+        y = yScale * ((winEnd-winStart)-leftStart);
+        hvGfxBox(hvg, (int)x+xOff, maxHeight-(int)y+yOff, (int)(xScale*(rightEnd-rightStart))+1, (int)(yScale*(leftEnd-leftStart))+1, colorIx);
+        }
+    }
+    // Draw top-left to bottom-right diagonal axis line in black
+    int colorIx = hvGfxFindColorIx(hvg, 0, 0, 0);
+    hvGfxLine(hvg, xOff, yOff, ((winEnd-winStart)*xScale)+xOff, maxHeight+yOff, colorIx);
+}
+
+
+int cmpHicItems(const void *elem1, const void *elem2)
+/* Comparison function for sorting Hi-C interactions by score */
+{
+struct interact *item1 = *((struct interact**)elem1);
+struct interact *item2 = *((struct interact**)elem2);
+if (item1->value > item2->value)
+    return 1;
+if (item1->value < item2->value)
+    return -1;
+return 0;
+}
+
+static void drawHicArc(struct track *tg, int seqStart, int seqEnd,
+        struct hvGfx *hvg, int xOff, int yOff, int width, 
+        MgFont *font, Color color, enum trackVisibility vis)
+/* Draw Hi-C interactions in arc mode */
+{
+double xScale = scaleForWindow(width, seqStart, seqEnd);
+double yScale = xScale;
+int maxHeight = tg->height;
+struct interact *hicItem = NULL;
+int binSize = fetchResolution(tg);
+if (binSize == 0)
+    return;
+
+if (vis == tvDense)
+    {
+    yScale *= 0.5;
+    }
+
+double maxScore = getHicMaxScore(tg);
+Color *colorIxs = ColorSetForHic(hvg, tg, HIC_SCORE_BINS+1);
+if (colorIxs == NULL)
+    return; // something went wrong with colors
+
+slSort(&tg->items, cmpHicItems); // So that the darkest arcs are drawn on top and not lost
+for (hicItem = (struct interact *)tg->items; hicItem; hicItem = hicItem->next)
+    {
+    int leftStart = hicItem->sourceStart - winStart;
+    int leftMidpoint = leftStart + binSize/2;
+    int rightStart = hicItem->targetStart - winStart;
+    int rightMidpoint = rightStart + binSize/2;
+    if ((leftMidpoint < 0) || (leftMidpoint > winEnd-winStart))
+        continue;  // skip this item - we'd be drawing to a point off the screen
+    if ((rightMidpoint < 0) || (rightMidpoint > winEnd-winStart))
+        continue;  // skip this item - we'd be drawing to a point off the screen
+
+    int colorIx;
+    if (hicItem->value > maxScore)
+        colorIx = colorIxs[HIC_SCORE_BINS];
+    else
+        colorIx = colorIxs[(int)(HIC_SCORE_BINS * hicItem->value/maxScore)];
+
+    double leftx = xScale * leftMidpoint;
+    double rightx = xScale * rightMidpoint;
+    double midx = xScale * (rightMidpoint+leftMidpoint)/2.0;
+    double midy = yScale * (rightMidpoint-leftMidpoint)/2.0;
+    midy *= 1.5; // Heuristic scaling for better use of vertical space
+    hvGfxCurve(hvg, (int)leftx+xOff, maxHeight+yOff, (int)midx+xOff, maxHeight-(int)midy+yOff,
+            (int)rightx+xOff, maxHeight+yOff, colorIx, FALSE);
+    }
+}
+
+void hicDrawItems(struct track *tg, int seqStart, int seqEnd,
+        struct hvGfx *hvg, int xOff, int yOff, int width, 
+        MgFont *font, Color color, enum trackVisibility vis)
+/* Draw a set of Hi-C interactions with the current user settings. */
+{
+char *drawMode = hicUiFetchDrawMode(cart, tg->track);
+if (sameString(drawMode,HIC_DRAW_MODE_SQUARE))
+    drawHicSquare(tg, seqStart, seqEnd, hvg, xOff, yOff, width, font, color, vis);
+else if (sameString(drawMode,HIC_DRAW_MODE_TRIANGLE))
+    drawHicTriangle(tg, seqStart, seqEnd, hvg, xOff, yOff, width, font, color, vis);
+else if (sameString(drawMode,HIC_DRAW_MODE_ARC))
+    drawHicArc(tg, seqStart, seqEnd, hvg, xOff, yOff, width, font, color, vis);
+else
+    warn ("Unknown draw mode %s for track %s", drawMode, tg->track);
+}
+
+
+void doLeftLabelsExceptNot(struct track *tg, int seqStart, int seqEnd, struct hvGfx *hvg,
+                               int xOff, int yOff, int width, int height, boolean withCenterLabels,
+                               MgFont *font, Color color, enum trackVisibility vis)
+/* A no-op function. There are no left labels associated with this track type. */
+{
+return;
+}
+
+void hicMethods(struct track *tg)
+/* Interact track type methods */
+{
+tg->bedSize = 12;
+tg->loadItems = hicLoadItems;
+tg->drawItems = hicDrawItems;
+tg->totalHeight = hicTotalHeight;
+tg->drawLeftLabels = doLeftLabelsExceptNot; // If we don't pretend to do it ourselves, hgTracks tries and fails badly
+tg->mapsSelf=1;
+}
+
+void hicCtMethods(struct track *tg)
+/* Interact track methods for custom track */
+{
+hicMethods(tg);
+}
+