0a01f6a29963ba42555b02782532a5702b649f29 braney Fri Mar 10 10:18:36 2023 -0800 first cut at drawing base probability logos in wiggle tracks calculated at run time using an associated MAF related diff --git src/hg/hgTracks/wigTrack.c src/hg/hgTracks/wigTrack.c index 0df427a..ce30adb 100644 --- src/hg/hgTracks/wigTrack.c +++ src/hg/hgTracks/wigTrack.c @@ -13,52 +13,60 @@ #include "hdb.h" #include "hgTracks.h" #include "wiggle.h" #include "hmmstats.h" #include "scoredRef.h" #ifndef GBROWSE #include "customTrack.h" #endif /* GBROWSE */ #include "wigCommon.h" #include "imageV2.h" #include "memgfx.h" #include "udc.h" #include "trashDir.h" #include "jsonWrite.h" #include "dnaMotif.h" +#include "maf.h" +#include "hgMaf.h" struct wigItem /* A wig track item. */ { struct wigItem *next; int start, end; /* Start/end in chrom (aka browser) coordinates. */ char *db; /* Database */ int ix; /* Position in list. */ int height; /* Pixel height of item. */ unsigned span; /* each value spans this many bases */ unsigned count; /* number of values to use */ unsigned offset; /* offset in File to fetch data */ char *file; /* path name to data file, one byte per value */ double lowerLimit; /* lowest data value in this block */ double dataRange; /* lowerLimit + dataRange = upperLimit */ unsigned validCount; /* number of valid data values in this block */ double sumData; /* sum of the data points, for average and stddev calc */ double sumSquares; /* sum of data points squared, for stddev calc */ double graphUpperLimit; /* filled in by DrawItems */ double graphLowerLimit; /* filled in by DrawItems */ }; +static boolean doLogo(struct trackDb *tdb) +/* Are we going to draw a logo? */ +{ +return trackDbSettingOn(tdb, "logo") || (trackDbSetting(tdb, "logoMaf") != NULL); +} + static void wigFillInColorLfArray(struct track *wigTrack, Color *colArray, int colSize, struct track *colorTrack) /* Fill in a color array with the linkedFeatures based colorTrack's color where it would normally have an exon. */ { struct linkedFeatures *lf = NULL, *lfList = colorTrack->items; struct simpleFeature *sf = NULL; double scale = scaleForPixels(colSize); int x1 = 0, x2 = 0; int i = 0; for(lf = lfList; lf != NULL; lf = lf->next) { for (sf = lf->components; sf != NULL; sf = sf->next) { x1 = round((double)((int)sf->start-winStart)*scale); @@ -1173,56 +1181,112 @@ /* honor the viewLimits, data below is white, data above is black */ grayValue = max(dataValue,graphLowerLimit); grayValue = min(grayValue,graphUpperLimit); grayIndex = ((grayValue-graphLowerLimit)/graphRange)*MAX_WIG_VALUE; drawColor = tg->colorShades[grayInRange(grayIndex, 0, MAX_WIG_VALUE)]; doLine(image, x, yOff, tg->lineHeight, drawColor); } /* vis == tvDense || vis == tvSquish */ } /* if (preDraw[].count) */ } /* for (x1 = 0; x1 < width; ++x1) */ return(mouseOverData); } /* graphPreDraw() */ +struct probVal // the probability of a certain nucleotide being seen in a column + { + struct probVal *next; + char *nuc; + double prob; + unsigned long color; + }; + +struct probVal probVals[4]; // one per nucleotide +struct probVal *probList = probVals; // a sortable list + +static int probListCmp(const void *va, const void *vb) +/* Compare probabilities for sorting. */ +{ +const struct probVal *a = *((struct probVal **)va); +const struct probVal *b = *((struct probVal **)vb); + +if (a->prob < b->prob) + return 1; +return -1; +} + + static struct wigMouseOver *logoPreDrawContainer(struct preDrawContainer *preDrawContainer, int preDrawZero, int width, struct track *tg, struct hvGfx *hvg, int xOff, int yOff, double graphUpperLimit, double graphLowerLimit, double graphRange, enum trackVisibility vis, struct wigCartOptions *wigCart, int seqStart, int seqEnd) { boolean baseCmpl = cartUsualBooleanDb(cart, database, COMPLEMENT_BASES_VAR, FALSE); struct preDrawElement *preDraw = preDrawContainer->preDraw; struct wigGraphOutput *wgo = tg->wigGraphOutput; //struct wigMouseOver *mouseOverData = NULL; unsigned numBases = seqEnd - seqStart; struct dnaSeq *seq = hChromSeq(database, chromName, seqStart, seqEnd); if (baseCmpl) complement(seq->dna, seq->size); struct pixelCountBin *pixelBins = wgo->pixelBins; double *yOffsets = wgo->yOffsets; int numTrack = wgo->numTrack; Color clipColor = MG_MAGENTA; WigVerticalLineVirtual vLine = wgo->vLine; void *image = wgo->image; #define doLine(image, x, y, height, color) {vLine(image, x, y, height, color); } int h = tg->lineHeight; /* the height of our drawing window */ double scaleFactor = h/graphRange; struct wigMouseOver *mouseOverData = getMouseOverData(tg, preDraw, width, xOff, preDrawZero); +struct mafBaseProbs *baseProbs = wigCart->baseProbs; +if (baseProbs != NULL) + { + // initialize our nucleotide probability data structure + probVals[0].next = &probVals[1]; + probVals[1].next = &probVals[2]; + probVals[2].next = &probVals[3]; + + boolean baseCmpl = cartUsualBooleanDb(cart, database, COMPLEMENT_BASES_VAR, FALSE); + if (baseCmpl) + { + probVals[3].nuc = "A"; + probVals[3].color = MG_RED; + probVals[2].nuc = "C"; + probVals[2].color = MG_BROWN; + probVals[1].nuc = "G"; + probVals[1].color = MG_BLUE; + probVals[0].nuc = "T"; + probVals[0].color = MG_GREEN; + } + else + { + probVals[0].nuc = "A"; + probVals[0].color = MG_RED; + probVals[1].nuc = "C"; + probVals[1].color = MG_BROWN; + probVals[2].nuc = "G"; + probVals[2].color = MG_BLUE; + probVals[3].nuc = "T"; + probVals[3].color = MG_GREEN; + } + } + double xIncr = (double)width / numBases; unsigned baseNum; for(baseNum = 0; baseNum < numBases; baseNum++) { int x1 = ceil(baseNum * xIncr); int x = x1 + xOff; int baseWidth = xIncr; int base = seq->dna[baseNum]; int preDrawIndex = x1 + preDrawZero; struct preDrawElement *p = &preDraw[preDrawIndex]; assert(x1/pixelBins->binSize < pixelBins->binCount); unsigned long *bitCount = &pixelBins->bins[x1/pixelBins->binSize]; /* count is non-zero meaning valid data exists here */ @@ -1296,41 +1360,66 @@ string[0] = toupper(base); string[1] = 0; MgFont *font = tl.font; int height = dataValue * scaleFactor; unsigned color = MG_BLACK; if (base == 'a') color = MG_RED; else if (base == 't') color = MG_GREEN; else if (base == 'c') color = MG_BROWN; else if (base == 'g') color = MG_BLUE; if (height != 0) { + if (baseProbs) + { + int thisHeight; + + // we want a sorted list so the most probable gets drawn on top + probVals[0].prob = baseProbs[baseNum].aProb; + probVals[1].prob = baseProbs[baseNum].cProb; + probVals[2].prob = baseProbs[baseNum].gProb; + probVals[3].prob = baseProbs[baseNum].tProb; + slSort(&probList,probListCmp); + + int y = yOff+boxTop; + struct probVal *pl = probList; + for(; pl; pl = pl->next) + { + thisHeight = pl->prob * boxHeight; + if (thisHeight) + hvGfxTextInBox(hvg, x, y, baseWidth - 1, thisHeight, + pl->color, font, pl->nuc); + y += thisHeight; + } + } + else + { if (dataValue < 0) { hvGfxTextInBox(hvg, x, yOff+boxTop, baseWidth - 1, -boxHeight, color, font, string); } else { hvGfxTextInBox(hvg, x, yOff+boxTop, baseWidth - 1, boxHeight, color, font, string); } } + } if (((boxTop+boxHeight) == 0) && !isnan(dataValue)) boxHeight += 1; } double stackValue = dataValue; if ((yOffsets != NULL) && (numTrack > 0)) stackValue += yOffsets[(numTrack-1) * width + x1]; if (stackValue > graphUpperLimit) { doLine(image, x, yOff, 2, clipColor); } else if (stackValue < graphLowerLimit) { doLine(image, x, yOff + h - 1, 2, clipColor); } @@ -1598,32 +1687,31 @@ } /* if we're autoscaling and the range is 0 this implies that all values * in the given range are the same. We create a bottom of the scale * by subtracting one from the only value. * This results in drawing a box that fills the range. */ if (graphUpperLimit == graphLowerLimit) { graphLowerLimit = graphUpperLimit - 1; } graphRange = graphUpperLimit - graphLowerLimit; wigTrackSetGraphOutputDefault(tg, xOff, yOff, width, hvg); struct wigMouseOver *mouseOverData = NULL; -if (zoomedToCodonLevel && trackDbSettingOn(tg->tdb, "logo") && vis != tvDense) -//if (zoomedToBaseLevel && trackDbSettingOn(tg->tdb, "logo")) +if (zoomedToCodonLevel && doLogo(tg->tdb) && vis != tvDense) mouseOverData = logoPreDrawContainer(preContainer, preDrawZero, width, tg, hvg, xOff, yOff, graphUpperLimit, graphLowerLimit, graphRange, vis, wigCart, seqStart, seqEnd); else mouseOverData = graphPreDrawContainer(preContainer, preDrawZero, width, tg, hvg, xOff, yOff, graphUpperLimit, graphLowerLimit, graphRange, vis, wigCart); drawZeroLine(vis, wigCart->horizontalGrid, graphUpperLimit, graphLowerLimit, hvg, xOff, yOff, width, tg->lineHeight); drawArbitraryYLine(vis, (enum wiggleGridOptEnum)wigCart->yLineOnOff, graphUpperLimit, graphLowerLimit, hvg, xOff, yOff, width, tg->lineHeight, wigCart->yLineMark, graphRange, @@ -1838,30 +1926,38 @@ } /* for (wi = tg->items; wi != NULL; wi = wi->next) */ if (wibFH > 0) { udcFileClose(&wibFH); wibFH = 0; freeMem(currentFile); } return pre; } static void wigPreDrawItems(struct track *tg, int seqStart, int seqEnd, struct hvGfx *hvg, int xOff, int yOff, int width, MgFont *font, Color color, enum trackVisibility vis) /* Draw wiggle items that resolve to doing a box for each pixel. */ { +char *logoMaf = trackDbSetting(tg->tdb, "logoMaf"); + +if (logoMaf != NULL) + { + struct wigCartOptions *wigCart = tg->wigCartData; + wigCart->baseProbs = hgMafProbs(database, logoMaf, chromName, seqStart, seqEnd, '+'); + } + struct preDrawContainer *pre = wigLoadPreDraw(tg, seqStart, seqEnd, width); if (pre != NULL) { wigPreDrawPredraw(tg, seqStart, seqEnd, hvg, xOff, yOff, width, font, color, vis, pre, pre->preDrawZero, pre->preDrawSize, &tg->graphUpperLimit, &tg->graphLowerLimit); } } void wigMultiRegionGraphLimits(struct track *tg) /* Set common graphLimits across all windows */ { double graphUpperLimit = -BIGDOUBLE; double graphLowerLimit = BIGDOUBLE; struct track *tgSave = tg;