1084bf14211e7e4ace072d7a213672f3a8f98be1 kate Fri Jul 24 17:30:21 2015 -0700 Scale wrt all genes, not just those in window. We now save the max tissue median score in the gtexInfo table. refs #15645 diff --git src/hg/hgTracks/gtexTracks.c src/hg/hgTracks/gtexTracks.c index c4d9df3..9cfb0a7 100644 --- src/hg/hgTracks/gtexTracks.c +++ src/hg/hgTracks/gtexTracks.c @@ -1,503 +1,505 @@ /* 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 "gtexInfo.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 valToHeight(double val, double maxVal, int maxHeight) /* Convert a value from 0 to maxVal to 0 to maxHeight-1 */ { -double scaled = val/maxVal; +// FIXME: This sort of works. Seems to be dropping some tho. e.g. check MRAP lung +if (val == 0.0) + return 0; +// smallest counts are 1x10e-3, translate to counter negativity +double scaled = (log10(val)+3.001)/(log10(maxVal)+3.001); +if (scaled < 0) + warn("scaled=%f\n", scaled); +//uglyf("%.2f -> %.2f height %d", val, scaled, (int)scaled * (maxHeight-1)); return (scaled * (maxHeight-1)); } +// TODO: whack this 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; +if (val == 0.0) + return 0; +double scaled = (log10(val)+3.001)/(log10(maxVal)+3.001); +if (scaled < 0) + warn("scaled=%f\n", scaled); int y = scaled * (maxHeight); +//uglyf("%.2f -> %.2f height %d", val, scaled, (int)scaled * (maxHeight-1)); return (maxHeight-1) - y; - -//int y = scaled * (height-1); -//return (height - 1) - y; } struct gtexGeneExtras { - double maxExp; // TODO: remove this local normalization factor + double maxMedian; 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; +double maxMedian = ((struct gtexGeneExtras *)tg->extraUiData)->maxMedian; 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 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); + int height = valToHeight(expScore, maxMedian, gtexGraphHeight()); + // TODO: adjust yGene to get gene track distance as desired + //if (i ==0) uglyf("DRAW: expScore=%.2f, maxMedian=%.2f, graphHeight=%d, y=%d<br>", expScore, maxMedian, 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); + hvGfxBox(hvg, x1, yZero-height, barWidth, height, fillColorIx); else - hvGfxOutlinedBox(hvg, x1, yMedian+1, barWidth, height, fillColorIx, lineColorIx); + hvGfxOutlinedBox(hvg, x1, yZero-height, 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()); + int height = valToHeight(expScore, maxMedian, 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: expScore=%.2f, maxMedian=%.2f, graphHeight=%d, y=%d<br>", expScore, maxMedian, 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; +double maxMedian = ((struct gtexGeneExtras *)tg->extraUiData)->maxMedian; 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; + //TODO: use valToHeight + int yMedian = valToY(expScore, maxMedian, 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: expScore=%.2f, maxMedian=%.2f, graphHeight=%d, y=%d<br>", expScore, maxMedian, 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 *tg) { 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; -// 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); -extras->maxExp = maxExp; +extras->maxMedian = gtexMaxMedianScore(NULL); } static int gtexGeneItemHeight(struct track *tg, void *item) { if ((item == NULL) || (tg->visibility == tvSquish) || (tg->visibility == tvDense)) return 0; 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 *tg, enum trackVisibility vis) /* Figure out total height of track */ { int height; 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() + extra; //uglyf("GTEX totalHeight = %d<br>", height); return tgFixedTotalHeightOptionalOverflow(tg, vis, height, height, FALSE); } 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 *tg) { tg->drawItemAt = gtexGeneDrawAt; tg->loadItems = gtexGeneLoadItems; tg->mapItem = gtexGeneMapItem; tg->itemHeight = gtexGeneItemHeight; tg->itemEnd = gtexGeneItemEnd; tg->totalHeight = gtexTotalHeight; }