c28867b853a769949049ae459cee687bcad25aca
kate
  Fri Jul 24 13:07:01 2015 -0700
Add support for comparison graphs (2 graphs, one top and one below inverted, based on splitting samples by sex or age).  First cut of display only. refs #15645

diff --git src/hg/hgTracks/gtexTracks.c src/hg/hgTracks/gtexTracks.c
index ba91a5a..c4d9df3 100644
--- src/hg/hgTracks/gtexTracks.c
+++ src/hg/hgTracks/gtexTracks.c
@@ -1,407 +1,503 @@
 /* GTEX tracks  */
 
 /* Copyright (C) 2015 The Regents of the University of California 
  * See README in this or parent directory for licensing information. */
 
 #include "common.h"
 #include "hgTracks.h"
 #include "hvGfx.h"
 #include "rainbow.h"
 #include "gtexGeneBed.h"
 #include "gtexTissue.h"
 #include "gtexUi.h"
 
 #define WIN_MAX_GRAPH 20000
 #define MAX_GRAPH_HEIGHT 100
 #define MAX_BAR_WIDTH 5
 #define MAX_GRAPH_PADDING 2
 
 #define WIN_MED_GRAPH 500000
 #define MED_GRAPH_HEIGHT 60
 #define MED_BAR_WIDTH 3
 #define MED_GRAPH_PADDING 1
 
 #define MIN_GRAPH_HEIGHT 20
 #define MIN_BAR_WIDTH 1
 #define MIN_GRAPH_PADDING 0
 
 #define MARGIN_WIDTH 1
 
 //#define GENCODE_CODING_COLOR 0x0c0c78    // rgb(12,12,120)
 //#define GENCODE_NONCODING_COLOR 0x006400 // rgb(0,100,0)
 //#define GENCODE_PROBLEM_COLOR 0xfe000  //rgb(254,0,0)  
 
 #define GENCODE_CODING_COLOR "12,12,120"
 #define GENCODE_NONCODING_COLOR "0,100,0"
 #define GENCODE_PROBLEM_COLOR "254,0,0"
 
 #define UNKNOWN_COLOR "1,1,1"
 
 static struct statusColors
     {
     Color coding;
     Color nonCoding;
     Color problem;
     Color unknown;
     } statusColors = {0,0,0};
 
 static int findColorIx(struct hvGfx *hvg, char *rgb)
 {
 unsigned char red, green, blue;
 parseColor(rgb, &red, &green, &blue);
 return (hvGfxFindColorIx(hvg, red, green, blue));
 //return MAKECOLOR_32(red, green, blue);
 }
 
 static void initGeneColors(struct hvGfx *hvg)
 {
 if (statusColors.coding != 0)
     return;
 statusColors.coding = findColorIx(hvg, GENCODE_CODING_COLOR);
 statusColors.nonCoding = findColorIx(hvg, GENCODE_NONCODING_COLOR);
 statusColors.problem = findColorIx(hvg, GENCODE_PROBLEM_COLOR);
 statusColors.unknown = findColorIx(hvg, UNKNOWN_COLOR);
 }
 
 static int gtexBarWidth()
 {
 int winSize = winEnd - winStart;
 if (winSize < WIN_MAX_GRAPH)
     return MAX_BAR_WIDTH;
 else if (winSize < WIN_MED_GRAPH)
     return MED_BAR_WIDTH;
 else
     return MIN_BAR_WIDTH;
 }
 
 static int gtexGraphPadding()
 {
 int winSize = winEnd - winStart;
 if (winSize < WIN_MAX_GRAPH)
     return MAX_GRAPH_PADDING;
 else if (winSize < WIN_MED_GRAPH)
     return MED_GRAPH_PADDING;
 else
     return MIN_GRAPH_PADDING;
 }
 
 static int gtexGraphHeight()
 {
 int winSize = winEnd - winStart;
 if (winSize < WIN_MAX_GRAPH)
     return MAX_GRAPH_HEIGHT;
 else if (winSize < WIN_MED_GRAPH)
     return MED_GRAPH_HEIGHT;
 else
     return MIN_GRAPH_HEIGHT;
 }
 
 static int gtexGraphWidth(struct gtexGeneBed *gtex)
 /* Width of GTEx graph in pixels */
 {
 int barWidth = gtexBarWidth();
 int padding = gtexGraphPadding();
 int count = gtex->expCount;
 return (barWidth * count) + (padding * (count-1));
 }
 
-static int valToY(double val, double maxVal, int height)
+
+static int valToHeight(double val, double maxVal, int maxHeight)
+/* Convert a value from 0 to maxVal to 0 to maxHeight-1 */
+{
+double scaled = val/maxVal;
+return (scaled * (maxHeight-1));
+}
+
+static int valToY(double val, double maxVal, int maxHeight)
 /* Convert a value from 0 to maxVal to 0 to height-1 */
 {
 double scaled = val/maxVal;
-int y = scaled * (height);
-return (height-1) - y;
+int y = scaled * (maxHeight);
+return (maxHeight-1) - y;
+
 //int y = scaled * (height-1);
 //return (height - 1) - y;
 }
 
 struct gtexGeneExtras 
     {
-    double maxExp;
+    double maxExp; // TODO: remove this local normalization factor
+    char *graphType;
+    boolean isComparison;
     };
 
-
 /* Cache tissue metadata */
 
 struct gtexTissue *getGtexTissues()
 /* Get tissue metadata from database */
 {
 static struct gtexTissue *gtexTissues = NULL;
 if (gtexTissues == NULL)
     {
     char query[1024];
     struct sqlConnection *conn = sqlConnect("hgFixed");
     sqlSafef(query, sizeof(query), "select * from gtexTissue order by id");
     gtexTissues = gtexTissueLoadByQuery(conn, query);
     }
 return gtexTissues;
 }
 
 struct rgbColor *getGtexTissueColors()
 /* Get RGB colors from tissue table */
 {
 struct gtexTissue *tissues = getGtexTissues();
 struct gtexTissue *tissue = NULL;
 int count = slCount(tissues);
 struct rgbColor *colors;
 AllocArray(colors, count);
 int i = 0;
 for (tissue = tissues; tissue != NULL; tissue = tissue->next)
     {
     // TODO: reconcile 
     colors[i] = (struct rgbColor){.r=COLOR_32_BLUE(tissue->color), .g=COLOR_32_GREEN(tissue->color), .b=COLOR_32_RED(tissue->color)};
     //colors[i] = mgColorIxToRgb(NULL, tissue->color);
     i++;
     }
 return colors;
 }
 
 static int gtexGraphX(struct gtexGeneBed *gtex)
 /* Locate graph on X, relative to viewport. Return -1 if it won't fit */
 {
 int start = max(gtex->chromStart, winStart);
 //int end = min(gtex->chromEnd, winEnd);
 //if (start > end)
     //return -1;
 double scale = scaleForWindow(insideWidth, winStart, winEnd);
 int x1 = round((start - winStart) * scale);
 //int x2 = round((end - winStart) * scale);
 //int width = gtexGraphWidth(gtex);
 //if (x1 + width > x2)
     //return -1;
 return x1;
 }
 
 static int gtexGeneHeight()
 {
     return 8; 
 }
 
 static int gtexGeneMargin()
 {
     return 1;
 }
 
+static struct gtexGeneBed *loadComputedMedians(struct gtexGeneBed *geneBed, char *graphType)
+/* Compute medians based on graph type.  Returns a list of 2 for comparison graph types */
+/* TODO: add support for filter function */
+{
+/* FIXME: dummy load of two for display implementation */
+struct gtexGeneBed *medians = NULL, *medians2 = NULL;
+AllocVar(medians);
+/* TODO: move to lib */
+medians->expCount = geneBed->expCount;
+AllocArray(medians->expScores, medians->expCount);
+int i;
+for (i = 0; i < medians->expCount; ++i)
+    medians->expScores[i] = geneBed->expScores[i];
+
+AllocVar(medians2);
+medians2->expCount = geneBed->expCount;
+AllocArray(medians2->expScores, medians2->expCount);
+for (i = 0; i < medians2->expCount; ++i)
+    medians2->expScores[i] = geneBed->expScores[i];
+
+medians->next = medians2;
+return medians;
+}
+
 static void gtexGeneDrawAt(struct track *tg, void *item, struct hvGfx *hvg, int xOff, int y, 
                 double scale, MgFont *font, Color color, enum trackVisibility vis)
 {
 struct gtexGeneBed *geneBed = item;
 initGeneColors(hvg);
 //warn("item: %s, xOff=%d\n", geneBed->name, xOff);
 // Color using transcriptClass
 Color statusColor;
 if (geneBed->transcriptClass == NULL)
     statusColor = statusColors.unknown;
 else if (sameString(geneBed->transcriptClass, "coding"))
     statusColor = statusColors.coding;
 else if (sameString(geneBed->transcriptClass, "nonCoding"))
     statusColor = statusColors.nonCoding;
 else if (sameString(geneBed->transcriptClass, "problem"))
     statusColor = statusColors.problem;
 else 
     statusColor = statusColors.unknown;
 
 if (vis != tvFull && vis != tvPack)
     {
     bedDrawSimpleAt(tg, item, hvg, xOff, y, scale, font, statusColor, vis);
     return;
     }
+struct gtexGeneBed *computedMedians = NULL;     // 1 or 2 (if comparison) 
+                                                // with medians computed for sample subsets
+
+struct gtexGeneExtras *extras = (struct gtexGeneExtras *)tg->extraUiData;
+if ((extras->isComparison) &&
+        (tg->visibility == tvFull || tg->visibility == tvPack) && 
+        gtexGraphHeight() != MIN_GRAPH_HEIGHT)
+    {
+    // compute medians based on configuration (comparisons, and later, filters)
+    computedMedians = loadComputedMedians(geneBed, extras->graphType);
+    }
 int i;
 int expCount = geneBed->expCount;
 double maxExp = ((struct gtexGeneExtras *)tg->extraUiData)->maxExp;
 struct rgbColor lineColor = {.r=0};
 int lineColorIx = hvGfxFindColorIx(hvg, lineColor.r, lineColor.g, lineColor.b);
 int heightPer = tg->heightPer;
 
 int graphX = gtexGraphX(geneBed);
 if (graphX < 0)
     return;
 // x1 is at left of graph
 int x1 = xOff + graphX;
+int startX = x1;
 
 // yZero is at bottom of graph
 int yZero = gtexGraphHeight() + y - 1;
 
 // draw faint line under graph to delineate extent when bars are missing (tissue w/ 0 expression)
+// TODO: skip missing bars -- then we can lose the gray line (at least for non-comparison mode)
 Color lightGray = MAKECOLOR_32(0xD1, 0xD1, 0xD1);
 int graphWidth = gtexGraphWidth(geneBed);
 hvGfxBox(hvg, x1, yZero+1, graphWidth, 1, lightGray);
 
 //uglyf("DRAW: xOff=%d, x1=%d, y=%d, yZero=%d<br>", xOff, x1, y, yZero);
 
 int barWidth = gtexBarWidth();
 int graphPadding = gtexGraphPadding();
 char *colorScheme = cartUsualStringClosestToHome(cart, tg->tdb, FALSE, GTEX_COLORS, 
                         GTEX_COLORS_DEFAULT);
 struct rgbColor *colors;
 if (sameString(colorScheme, GTEX_COLORS_GTEX))
     {
     // retrieve from table
     // TODO: cache this
     colors = getGtexTissueColors();
     }
 else
     {
     // currently the only other choice
     // TODO: cache this
     colors = getRainbow(&saturatedRainbowAtPos, expCount);
     //colors = getRainbow(&lightRainbowAtPos, expCount);
     }
 for (i=0; i<expCount; i++)
     {
     struct rgbColor fillColor = colors[i];
     if (barWidth == 1 && sameString(colorScheme, GTEX_COLORS_GTEX))
         {
         // brighten colors a bit so they'll be more visible at this scale
         struct hslColor hsl = mgRgbToHsl(fillColor);
         hsl.s = min(1000, hsl.s + 300);
         fillColor = mgHslToRgb(hsl);
         }
     int fillColorIx = hvGfxFindColorIx(hvg, fillColor.r, fillColor.g, fillColor.b);
     double expScore = geneBed->expScores[i];
     int yMedian = valToY(expScore, maxExp, gtexGraphHeight()) + y;
-    //int height = yZero - yMedian + 1;
+    //int foo = valToY(expScore, maxExp, gtexGraphHeight());
+    //int yMedian = yZero - height;
     int height = yZero - yMedian;
     // TODO: adjust yGene instead of yMedian+1 to get gene track distance as desired
     //if (i ==0) uglyf("DRAW: expScore=%.2f, maxExp=%.2f, graphHeight=%d, y=%d<br>", expScore, maxExp, gtexGraphHeight(), y);
     //if (i ==0) uglyf("DRAW: yZero=%d, yMedian=%d, height=%d<br>", yZero, yMedian, height);
     if (graphPadding == 0 || sameString(colorScheme, GTEX_COLORS_GTEX))
         hvGfxBox(hvg, x1, yMedian+1, barWidth, height, fillColorIx);
     else
         hvGfxOutlinedBox(hvg, x1, yMedian+1, barWidth, height, fillColorIx, lineColorIx);
     x1 = x1 + barWidth + graphPadding;
     }
 
 // mark gene extent
 int yGene = yZero + gtexGeneMargin() - 1;
 
 // load & draw gene model
 
 char query[1024];
 char **row;
 sqlSafef(query, sizeof query, "select * from gtexGeneModel where name='%s'", geneBed->geneId);
 struct sqlConnection *conn = hAllocConn(database);
 if (conn == NULL)
     return;
 //uglyf("query: %s<br>", query);
 struct sqlResult *sr = sqlGetResult(conn, query);
 struct genePred *geneModel = NULL;
 if (sr != NULL)
     {
     if ((row = sqlNextRow(sr)) != NULL)
         geneModel = genePredLoad(row);
     sqlFreeResult(&sr);
 }
 if (geneModel == NULL)
 return;
 struct linkedFeatures *lf = linkedFeaturesFromGenePred(tg, geneModel, FALSE);
 tg->heightPer = gtexGeneHeight()+1;
 lf->filterColor = statusColor;
 linkedFeaturesDrawAt(tg, lf, hvg, xOff, yGene, scale, font, color, tvSquish);
 tg->heightPer = heightPer;
 hFreeConn(&conn);
+
+if (!extras->isComparison || slCount(computedMedians) != 2)
+    return;
+
+// draw comparison graph (upside down)
+
+x1 = startX;
+// yZero is at top of graph
+yZero = yGene + gtexGeneHeight();
+for (i=0; i<expCount; i++)
+    {
+    struct rgbColor fillColor = colors[i];
+    if (barWidth == 1 && sameString(colorScheme, GTEX_COLORS_GTEX))
+        {
+        // brighten colors a bit so they'll be more visible at this scale
+        struct hslColor hsl = mgRgbToHsl(fillColor);
+        hsl.s = min(1000, hsl.s + 300);
+        fillColor = mgHslToRgb(hsl);
+        }
+    int fillColorIx = hvGfxFindColorIx(hvg, fillColor.r, fillColor.g, fillColor.b);
+    double expScore = geneBed->expScores[i];
+    int height = valToHeight(expScore, maxExp, gtexGraphHeight());
+    // TODO: adjust yGene instead of yMedian+1 to get gene track distance as desired
+    //if (i ==0) uglyf("DRAW2: expScore=%.2f, maxExp=%.2f, graphHeight=%d, y=%d<br>", expScore, maxExp, gtexGraphHeight(), y);
+    //if (i ==0) uglyf("DRAW2: yZero=%d, height=%d<br>", yZero, height);
+    if (graphPadding == 0 || sameString(colorScheme, GTEX_COLORS_GTEX))
+        hvGfxBox(hvg, x1, yZero, barWidth, height, fillColorIx);
+    else
+        hvGfxOutlinedBox(hvg, x1, yZero, barWidth, height, fillColorIx, lineColorIx);
+    x1 = x1 + barWidth + graphPadding;
+    }
 }
 
 static void gtexGeneMapItem(struct track *tg, struct hvGfx *hvg, void *item, char *itemName, 
                         char *mapItemName, int start, int end, int x, int y, int width, int height)
 /* Create a map box for each tissue (bar in the graph) or a single map for squish/dense modes */
 {
 //uglyf("map item: itemName=%s, mapItemName=%s, start=%d, end=%d, x=%d, y=%d, width=%d, height=%d, insideX=%d<br>",
         //itemName, mapItemName, start, end, x, y, width, height, insideX);
 
 if (tg->visibility == tvDense || tg->visibility == tvSquish)
     {
     genericMapItem(tg, hvg, item, itemName, itemName, start, end, x, y, width, height);
     }
 
 struct gtexTissue *tissues = getGtexTissues();
 struct gtexTissue *tissue = NULL;
 struct gtexGeneBed *gtex = item;
 int barWidth = gtexBarWidth();
 int padding = gtexGraphPadding();
 double maxExp = ((struct gtexGeneExtras *)tg->extraUiData)->maxExp;
 
 int graphX = gtexGraphX((struct gtexGeneBed *)item);
 if (graphX < 0)
     return;
 // x1 is at left of graph
 int x1 = insideX + graphX;
 
 int i = 0;
 int yZero = gtexGraphHeight() + y - 1;
 for (tissue = tissues; tissue != NULL; tissue = tissue->next, i++)
     {
     double expScore = gtex->expScores[i];
     int yMedian = valToY(expScore, maxExp, gtexGraphHeight()) + y;
     int height = yZero - yMedian;
     // TODO: call genericMapItem
     //genericMapItem(tg, hvg, item, itemName, tissue->description, start, end, x1, y, barWidth, height);
     mapBoxHc(hvg, start, end, x1, yMedian+1, barWidth, height, tg->track, mapItemName, tissue->description);
     //if (i==0) uglyf("MAP: expScore=%.2f, maxExp=%.2f, graphHeight=%d, y=%d<br>", expScore, maxExp, gtexGraphHeight(), y);
     //if (i==0) uglyf("MAP: x=%d, x1=%d, y=%d, yZero=%d<br>", x, x1, y, yZero); 
     //if (i==0) uglyf("MAP: yZero=%d, yMedian=%d, height=%d<br>", yZero, yMedian, height); 
     x1 = x1 + barWidth + padding;
     }
 }
 
 static struct gtexGeneBed *loadGtexGeneBed(char **row)
 {
 return gtexGeneBedLoad(row);
 }
 
-static void gtexGeneLoadItems(struct track *track)
+static void gtexGeneLoadItems(struct track *tg)
 {
-bedLoadItem(track, track->table, (ItemLoader)loadGtexGeneBed);
+bedLoadItem(tg, tg->table, (ItemLoader)loadGtexGeneBed);
+
+struct gtexGeneExtras *extras;
+AllocVar(extras);
+tg->extraUiData = extras;
+
+// TODO: move test to lib
+char *graphType = cartUsualStringClosestToHome(cart, tg->tdb, FALSE, GTEX_GRAPH, 
+                                                GTEX_GRAPH_DEFAULT);
+extras->graphType = cloneString(graphType);
+if (sameString(graphType, GTEX_GRAPH_AGE) || sameString(graphType, GTEX_GRAPH_SEX))
+    extras->isComparison = TRUE;
+
 double maxExp = 0.0;
 int i;
 struct gtexGeneBed *geneBed;
-for (geneBed = track->items; geneBed != NULL; geneBed = geneBed->next)
+// TODO: Remove this fake (window-based) normalization
+for (geneBed = tg->items; geneBed != NULL; geneBed = geneBed->next)
     for (i=0; i<geneBed->expCount; i++)
         maxExp = (geneBed->expScores[i] > maxExp ? geneBed->expScores[i] : maxExp);
-struct gtexGeneExtras *extras;
-AllocVar(extras);
 extras->maxExp = maxExp;
-track->extraUiData = extras;
-
 }
 
-static int gtexGeneItemHeight(struct track *track, void *item)
+static int gtexGeneItemHeight(struct track *tg, void *item)
 {
-if ((item == NULL) || (track->visibility == tvSquish) || (track->visibility == tvDense))
+if ((item == NULL) || (tg->visibility == tvSquish) || (tg->visibility == tvDense))
     return 0;
-return gtexGraphHeight() + + gtexGeneMargin() + gtexGeneHeight();
+int extra = 0;
+if (((struct gtexGeneExtras *)tg->extraUiData)->isComparison)
+    extra = gtexGraphHeight() + 2;
+//uglyf("GTEX itemHeight extra = %d<br>", extra);
+return gtexGraphHeight() + gtexGeneMargin() + gtexGeneHeight() + extra;
 }
 
-static int gtexTotalHeight(struct track *track, enum trackVisibility vis)
+static int gtexTotalHeight(struct track *tg, enum trackVisibility vis)
 /* Figure out total height of track */
 {
 int height;
-if (track->visibility == tvSquish || track->visibility == tvDense)
+int extra = 0;
+if (((struct gtexGeneExtras *)tg->extraUiData)->isComparison)
+    extra = gtexGraphHeight() + 2;
+if (tg->visibility == tvSquish || tg->visibility == tvDense)
     height = 10;
 else 
-    height = gtexGraphHeight() + gtexGeneMargin() + gtexGeneHeight();
-return tgFixedTotalHeightOptionalOverflow(track, vis, height, height, FALSE);
+    height = gtexGraphHeight() + gtexGeneMargin() + gtexGeneHeight() + extra;
+//uglyf("GTEX totalHeight = %d<br>", height);
+return tgFixedTotalHeightOptionalOverflow(tg, vis, height, height, FALSE);
 }
 
-static int gtexGeneItemEnd(struct track *track, void *item)
+static int gtexGeneItemEnd(struct track *tg, void *item)
 /* Return end chromosome coordinate of item, including graph */
 {
 struct gtexGeneBed *geneBed = (struct gtexGeneBed *)item;
 double scale = scaleForWindow(insideWidth, winStart, winEnd);
 int graphWidth = gtexGraphWidth(geneBed);
 return max(geneBed->chromEnd, max(winStart, geneBed->chromStart) + graphWidth/scale);
 }
 
 
-void gtexGeneMethods(struct track *track)
+void gtexGeneMethods(struct track *tg)
 {
-track->drawItemAt = gtexGeneDrawAt;
-track->loadItems = gtexGeneLoadItems;
-track->mapItem = gtexGeneMapItem;
-track->itemHeight = gtexGeneItemHeight;
-track->itemEnd = gtexGeneItemEnd;
-track->totalHeight = gtexTotalHeight;
+tg->drawItemAt = gtexGeneDrawAt;
+tg->loadItems = gtexGeneLoadItems;
+tg->mapItem = gtexGeneMapItem;
+tg->itemHeight = gtexGeneItemHeight;
+tg->itemEnd = gtexGeneItemEnd;
+tg->totalHeight = gtexTotalHeight;
 }