40a345b04f483a7e67494e62afc8528675eb587e braney Mon Mar 6 12:06:43 2023 -0800 make sure alignment blockse are single coverage over the query. diff --git src/hg/hgTracks/snake.c src/hg/hgTracks/snake.c index 61e8425..215f29e 100644 --- src/hg/hgTracks/snake.c +++ src/hg/hgTracks/snake.c @@ -1,1210 +1,1256 @@ /* Copyright (C) 2014 The Regents of the University of California * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ /* snakeTrack - stuff to load and display snake type tracks in browser. */ #include "common.h" #include "hash.h" #include "localmem.h" #include "linefile.h" #include "jksql.h" #include "hdb.h" #include "hgTracks.h" #include "chainBlock.h" #include "chainLink.h" #include "chainDb.h" #include "chainCart.h" #include "errCatch.h" #include "twoBit.h" #include "bigWarn.h" #include <pthread.h> #include "trackHub.h" #include "limits.h" #include "snakeUi.h" #include "bits.h" #include "trix.h" #include "chromAlias.h" #include "bigPsl.h" #include "snake.h" +static int snakeFeatureCmpScore(const void *va, const void *vb) +/* sort by score of the alignment. */ +{ +const struct snakeFeature *a = *((struct snakeFeature **)va); +const struct snakeFeature *b = *((struct snakeFeature **)vb); +int diff = b->score - a->score; + +return diff; +} + static int snakeFeatureCmpQStart(const void *va, const void *vb) /* sort by start position on the query sequence */ { const struct snakeFeature *a = *((struct snakeFeature **)va); const struct snakeFeature *b = *((struct snakeFeature **)vb); int diff = a->qStart - b->qStart; if (diff == 0) { diff = a->start - b->start; } return diff; } // static machinery to calculate full and pack snake #define NUM_LEVELS 1000 struct level { boolean init; /* has this level been initialized */ int orientation; /* strand.. see above */ unsigned long edge; /* the leading edge of this level */ int adjustLevel; /* used to compress out the unused levels */ boolean hasBlock; /* are there any blocks in this level */ // these fields are used to compact full snakes unsigned *connectsToLevel; /* what levels does this level connect to */ Bits *pixels; /* keep track of what pixels are used */ struct snakeFeature *blocks; /* the ungapped blocks attached to this level */ }; static struct level Levels[NUM_LEVELS]; /* for packing the snake, not re-entrant! */ static int maxLevel = 0; /* deepest level */ /* blocks that make it through the min size filter */ static struct snakeFeature *newList = NULL; static void clearLevels() /* clear out the data structure that we use to calculate snakes */ { int ii; for(ii=0; ii < sizeof(Levels) / sizeof(Levels[0]); ii++) Levels[ii].init = FALSE; maxLevel = 0; } static void calcFullSnakeHelper(struct snakeFeature *list, int level) // calculate a full snake // updates global newList with unsorted blocks that pass the min // size filter. { struct snakeFeature *cb = list; struct snakeFeature *proposedList = NULL; if (level > maxLevel) maxLevel = level; if (level > ArraySize(Levels)) errAbort("too many levels"); if (Levels[level].init == FALSE) { // initialize the level if this is the first time we've seen it Levels[level].init = TRUE; Levels[level].orientation = list->orientation; if (list->orientation == -1) Levels[level].edge = LONG_MAX; // bigger than the biggest chrom else Levels[level].edge = 0; } // now we step through the blocks and assign them to levels struct snakeFeature *next; struct snakeFeature *insertHead = NULL; for(; cb; cb = next) { // we're going to add this block to a different list // so keep track of the next pointer next = cb->next; // if this block can't fit on this level add to the insert // list and move on to the next block if ((Levels[level].orientation != cb->orientation) || ((cb->orientation == 1) && (Levels[level].edge > cb->start)) || ((cb->orientation == -1) && (Levels[level].edge < cb->end))) { // add this to the list of an insert slAddHead(&insertHead, cb); continue; } // the current block will fit on this level // go ahead and deal with the blocks that wouldn't if (insertHead) { // we had an insert, go ahead and calculate where that goes slReverse(&insertHead); calcFullSnakeHelper(insertHead, level + 1); insertHead = NULL; } // assign the current block to this level cb->level = level; if (cb->orientation == 1) Levels[level].edge = cb->end; else Levels[level].edge = cb->start; // add this block to the list of proposed blocks for this level slAddHead(&proposedList, cb); } // we're at the end of the list. Deal with any blocks // that didn't fit on this level if (insertHead) { slReverse(&insertHead); int nextLevel = level + 1; calcFullSnakeHelper(insertHead, nextLevel); insertHead = NULL; } // do we have any proposed blocks if (proposedList == NULL) return; // we parsed all the blocks in the list to see if they fit in our level // now let's see if the list of blocks for this level is big enough // to actually add struct snakeFeature *temp; double scale = scaleForWindow(insideWidth, winStart, winEnd); int start, end; // order the blocks so lowest start position is first if (Levels[level].orientation == 1) slReverse(&proposedList); start=proposedList->start; for(temp=proposedList; temp->next; temp = temp->next) ; end=temp->end; // calculate how big the string of blocks is in screen coordinates double apparentSize = scale * (end - start); if (apparentSize < 0) errAbort("size of a list of blocks should not be less than zero"); // check to see if the apparent size is big enough if (apparentSize < 1) return; // transfer proposedList to new block list for(temp=proposedList; temp; temp = next) { next = temp->next; temp->next = NULL; slAddHead(&newList, temp); } } static void freeFullLevels() // free the connection levels and bitmaps for each level { int ii; for(ii=0; ii <= maxLevel; ii++) { freez(&Levels[ii].connectsToLevel); bitFree(&Levels[ii].pixels); } } static void initializeFullLevels(struct snakeFeature *sfList) /* initialize levels for compact snakes */ { int ii; for(ii=0; ii <= maxLevel; ii++) { Levels[ii].connectsToLevel = needMem((maxLevel+1) * sizeof(unsigned)); Levels[ii].pixels = bitAlloc(insideWidth); Levels[ii].blocks = NULL; } struct snakeFeature *sf, *next; int prevLevel = -1; double scale = scaleForWindow(insideWidth, winStart, winEnd); for(sf=sfList; sf; sf = next) { next = sf->next; if (positiveRangeIntersection(sf->start, sf->end, winStart, winEnd)) { int sClp = (sf->start < winStart) ? winStart : sf->start; sf->pixX1 = round((sClp - winStart)*scale); int eClp = (sf->end > winEnd) ? winEnd : sf->end; sf->pixX2 = round((eClp - winStart)*scale); } else { sf->pixX1 = sf->pixX2 = 0; } bitSetRange(Levels[sf->level].pixels, sf->pixX1, sf->pixX2 - sf->pixX1); if (prevLevel != -1) Levels[sf->level].connectsToLevel[prevLevel] = 1; if(next != NULL) Levels[sf->level].connectsToLevel[next->level] = 1; prevLevel = sf->level; slAddHead(&Levels[sf->level].blocks, sf); } } static int highestLevelBelowUs(int level) // is there a level below us that doesn't connect to us { int newLevel; // find the highest level below us that we connect to for(newLevel = level - 1; newLevel >= 0; newLevel--) if (Levels[level].connectsToLevel[newLevel]) break; // if the next level is more than one below us, we may // be able to push our blocks to it. if ((newLevel < 0) || (newLevel + 1 == level)) return -1; return newLevel + 1; } // the number of pixels we need clear on either side of the blocks we push #define EXTRASPACE 10 static boolean checkBlocksOnLevel(int level, struct snakeFeature *sfList) // is there enough pixels on this level to hold our blocks? { struct snakeFeature *sf; for(sf=sfList; sf; sf = sf->next) { if (bitCountRange(Levels[level].pixels, sf->pixX1 - EXTRASPACE, (sf->pixX2 - sf->pixX1) + 2*EXTRASPACE)) return FALSE; } return TRUE; } static void collapseLevel(int oldLevel, int newLevel) // push level up to higher level { struct snakeFeature *sf; struct snakeFeature *next; for(sf=Levels[oldLevel].blocks; sf; sf = next) { next = sf->next; sf->level = newLevel; slAddHead(&Levels[newLevel].blocks, sf); bitSetRange(Levels[sf->level].pixels, sf->pixX1, sf->pixX2 - sf->pixX1); } Levels[oldLevel].blocks = NULL; Levels[oldLevel].pixels = bitAlloc(insideWidth); } static void remapConnections(int oldLevel, int newLevel) // if we've moved a level, we need to remap all the connections { int ii, jj; for(ii=0; ii <= maxLevel; ii++) { for(jj=0; jj <= maxLevel; jj++) { if (Levels[ii].connectsToLevel[oldLevel]) { Levels[ii].connectsToLevel[oldLevel] = 0; Levels[ii].connectsToLevel[newLevel] = 1; } } } } static struct snakeFeature *compactSnakes(struct snakeFeature *sfList) // collapse levels that will fit on a higer level { int ii; initializeFullLevels(sfList); // start with third level to see what can be compacted for(ii=2; ii <= maxLevel; ii++) { int newLevel; // is there a level below us that doesn't connect to us if ((newLevel = highestLevelBelowUs(ii)) > 0) { int jj; for (jj=newLevel; jj < ii; jj++) { if (checkBlocksOnLevel(jj, Levels[ii].blocks)) { collapseLevel(ii, jj); remapConnections(ii, jj); break; } } } } // reattach blocks in the levels to linkedFeature struct snakeFeature *newList = NULL; for(ii=0; ii <= maxLevel; ii++) { newList = slCat(Levels[ii].blocks, newList); } slSort(&newList, snakeFeatureCmpQStart); // figure out the new max int newMax = 0; for(ii=0; ii <= maxLevel; ii++) if (Levels[ii].blocks != NULL) newMax = ii; freeFullLevels(); maxLevel = newMax; return newList; } static void calcFullSnake(struct track *tg, void *item) // calculate a full snake { struct linkedFeatures *lf = (struct linkedFeatures *)item; // if there aren't any blocks, don't bother if (lf->components == NULL) return; // we use the snakeInfo field to keep track of whether we already // calculated the height of this snake if (lf->snakeInfo == NULL) { clearLevels(); struct snakeFeature *sf; // this will destroy lf->components, and add to newList calcFullSnakeHelper((struct snakeFeature *)lf->components, 0); lf->components = (struct simpleFeature *)newList; newList = NULL; slSort(&lf->components, snakeFeatureCmpQStart); // now we're going to compress the levels that aren't used // to do that, we need to see which blocks are on the screen, // or connected to something on the screen int oldMax = maxLevel; clearLevels(); struct snakeFeature *prev = NULL; for(sf=(struct snakeFeature *)lf->components; sf; prev = sf, sf = sf->next) { if (Levels[sf->level].init == FALSE) { Levels[sf->level].init = TRUE; Levels[sf->level].orientation = sf->orientation; Levels[sf->level].hasBlock = FALSE; } else { if (Levels[sf->level].orientation != sf->orientation) errAbort("found snakeFeature with wrong orientation"); if ((sf->orientation == 1) && (Levels[sf->level].edge > sf->start)) errAbort("found snakeFeature that violates edge"); if ((sf->orientation == -1) && (Levels[sf->level].edge < sf->end)) errAbort("found snakeFeature that violates edge"); } if (sf->orientation == 1) Levels[sf->level].edge = sf->end; else Levels[sf->level].edge = sf->start; if (sf->level > maxLevel) maxLevel = sf->level; if(positiveRangeIntersection(winStart, winEnd, sf->start, sf->end)) { Levels[sf->level].hasBlock = TRUE; if (sf->next != NULL) Levels[sf->next->level].hasBlock = TRUE; if (prev != NULL) Levels[prev->level].hasBlock = TRUE; } } // now figure out how to remap the blocks int ii; int count = 0; for(ii=0; ii < oldMax + 1; ii++) { Levels[ii].adjustLevel = count; if ((Levels[ii].init) && (Levels[ii].hasBlock)) count++; } maxLevel = count; // remap blocks for(sf=(struct snakeFeature *)lf->components; sf; sf = sf->next) sf->level = Levels[sf->level].adjustLevel; // now compact the snakes lf->components = (void *)compactSnakes((struct snakeFeature *)lf->components); struct snakeInfo *si; AllocVar(si); si->maxLevel = maxLevel; lf->snakeInfo = si; } } static void calcPackSnakeHelper(struct snakeFeature *list, int level) // calculate a packed snake // updates global newList with unsorted blocks that pass the min // size filter. { struct snakeFeature *cb = list; struct snakeFeature *proposedList = NULL; if (level > maxLevel) maxLevel = level; if (level > ArraySize(Levels)) errAbort("too many levels"); if (Levels[level].init == FALSE) { // initialize the level if this is the first time we've seen it Levels[level].init = TRUE; Levels[level].edge = 0; } // now we step through the blocks and assign them to levels struct snakeFeature *next; struct snakeFeature *insertHead = NULL; for(; cb; cb = next) { // we're going to add this block to a different list // so keep track of the next pointer next = cb->next; // if this block can't fit on this level add to the insert // list and move on to the next block if ( Levels[level].edge > cb->start) { // add this to the list of an insert slAddHead(&insertHead, cb); continue; } // the current block will fit on this level // go ahead and deal with the blocks that wouldn't if (insertHead) { // we had an insert, go ahead and calculate where that goes slReverse(&insertHead); calcPackSnakeHelper(insertHead, level + 1); insertHead = NULL; } // assign the current block to this level cb->level = level; Levels[level].edge = cb->end; // add this block to the list of proposed blocks for this level slAddHead(&proposedList, cb); } // we're at the end of the list. Deal with any blocks // that didn't fit on this level if (insertHead) { slReverse(&insertHead); calcPackSnakeHelper(insertHead, level + 1); insertHead = NULL; } // do we have any proposed blocks if (proposedList == NULL) return; // we parsed all the blocks in the list to see if they fit in our level // now let's see if the list of blocks for this level is big enough // to actually add // order the blocks so lowest start position is first slReverse(&proposedList); struct snakeFeature *temp; double scale = scaleForWindow(insideWidth, winStart, winEnd); int start, end; start=proposedList->start; for(temp=proposedList; temp->next; temp = temp->next) ; end=temp->end; // calculate how big the string of blocks is in screen coordinates double apparentSize = scale * (end - start); if (apparentSize < 0) errAbort("size of a list of blocks should not be less than zero"); // check to see if the apparent size is big enough if (apparentSize < 1) return; // transfer proposedList to new block list for(temp=proposedList; temp; temp = next) { next = temp->next; temp->next = NULL; slAddHead(&newList, temp); } } static void calcPackSnake(struct track *tg, void *item) { struct linkedFeatures *lf = (struct linkedFeatures *)item; if (lf->components == NULL) return; // we use the snakeIfno field to keep track of whether we already // calculated the height of this snake if (lf->snakeInfo == NULL) { clearLevels(); // this will destroy lf->components, and add to newList calcPackSnakeHelper((struct snakeFeature *)lf->components, 0); lf->components = (struct simpleFeature *)newList; newList = NULL; //slSort(&lf->components, snakeFeatureCmpQStart); struct snakeInfo *si; AllocVar(si); si->maxLevel = maxLevel; lf->snakeInfo = si; } } static int snakeItemHeight(struct track *tg, void *item) // return height of a single packed snake { if ((item == NULL) || (tg->visibility == tvSquish) || (tg->visibility == tvDense)) return 0; struct linkedFeatures *lf = (struct linkedFeatures *)item; if (lf->components == NULL) return 0; if (tg->visibility == tvFull) calcFullSnake(tg, item); else if (tg->visibility == tvPack) calcPackSnake(tg, item); struct snakeInfo *si = (struct snakeInfo *)lf->snakeInfo; int lineHeight = tg->lineHeight ; int multiplier = 1; if (tg->visibility == tvFull) multiplier = 2; return (si->maxLevel + 1) * (multiplier * lineHeight); } static int snakeHeight(struct track *tg, enum trackVisibility vis) /* calculate height of all the snakes being displayed */ { if (tg->networkErrMsg != NULL) { // we had a parallel load failure tg->drawItems = bigDrawWarning; tg->totalHeight = bigWarnTotalHeight; return bigWarnTotalHeight(tg, vis); } if (vis == tvDense) return tg->lineHeight; if (vis == tvSquish) return tg->lineHeight/2; int height = 0; struct linkedFeatures *item = tg->items; for (;item; item = item->next) { height += tg->itemHeight(tg, item); } return height; } // this is a 16 color palette with every other color being a lighter version of // the color before it static int snakePalette2[] = { 0x1f77b4, 0xaec7e8, 0xff7f0e, 0xffbb78, 0x2ca02c, 0x98df8a, 0xd62728, 0xff9896, 0x9467bd, 0xc5b0d5, 0x8c564b, 0xc49c94, 0xe377c2, 0xf7b6d2, 0x7f7f7f, 0xc7c7c7, 0xbcbd22, 0xdbdb8d, 0x17becf, 0x9edae5 }; static Color hashColor(char *name) { bits32 hashVal = hashString(name); unsigned int colorInt = snakePalette2[hashVal % (sizeof(snakePalette2)/sizeof(Color))]; return MAKECOLOR_32(((colorInt >> 16) & 0xff),((colorInt >> 8) & 0xff),((colorInt >> 0) & 0xff)); } static void boundMapBox(struct hvGfx *hvg, int start, int end, int x, int y, int width, int height, char *track, char *item, char *statusLine, char *directUrl, boolean withHgsid, char *extra) // make sure start x and end x position are on the screen // otherwise the tracking box code gets confused { if (x > insideWidth) return; if (x < 0) { width -= -x; x = 0; } if (x + width > insideWidth) { width -= x + width - insideWidth; } mapBoxHgcOrHgGene(hvg, start, end, x, y, width, height, track, item, statusLine, directUrl, withHgsid, extra); } static void snakeDraw(struct track *tg, int seqStart, int seqEnd, struct hvGfx *hvg, int xOff, int yOff, int width, MgFont *font, Color color, enum trackVisibility vis) /* Draw snake items. */ { struct slList *item; int y; struct linkedFeatures *lf; double scale = scaleForWindow(width, seqStart, seqEnd); int height = snakeHeight(tg, vis); hvGfxSetClip(hvg, xOff, yOff, width, height); if ((tg->visibility == tvFull) || (tg->visibility == tvPack)) { // score snakes by how many bases they cover for (item = tg->items; item != NULL; item = item->next) { lf = (struct linkedFeatures *)item; struct snakeFeature *sf; lf->score = 0; for (sf = (struct snakeFeature *)lf->components; sf != NULL; sf = sf->next) { lf->score += sf->end - sf->start; } } slSort(&tg->items, linkedFeaturesCmpScore); } y = yOff; for (item = tg->items; item != NULL; item = item->next) { if(tg->itemColor != NULL) color = tg->itemColor(tg, item, hvg); tg->drawItemAt(tg, item, hvg, xOff, y, scale, font, color, vis); if (vis == tvFull) y += tg->itemHeight(tg, item); } } static void snakeDrawAt(struct track *tg, void *item, struct hvGfx *hvg, int xOff, int y, double scale, MgFont *font, Color color, enum trackVisibility vis) /* Draw a single simple bed item at position. */ { unsigned showSnpWidth = cartOrTdbInt(cart, tg->tdb, SNAKE_SHOW_SNP_WIDTH, SNAKE_DEFAULT_SHOW_SNP_WIDTH); struct linkedFeatures *lf = (struct linkedFeatures *)item; if (tg->visibility == tvFull) calcFullSnake(tg, item); else if (tg->visibility == tvPack) calcPackSnake(tg, item); if (lf->components == NULL) return; struct snakeFeature *sf = (struct snakeFeature *)lf->components, *prevSf = NULL; int s = tg->itemStart(tg, item); int sClp = (s < winStart) ? winStart : s; int x1 = round((sClp - winStart)*scale) + xOff; int textX = x1; int yOff = y; boolean withLabels = (withLeftLabels && (vis == tvFull) && !tg->drawName); unsigned labelColor = MG_BLACK; // draw the labels if (withLabels) { char *name = tg->itemName(tg, item); int nameWidth = mgFontStringWidth(font, name); int dotWidth = tl.nWidth/2; boolean snapLeft = FALSE; boolean drawNameInverted = FALSE; textX -= nameWidth + dotWidth; snapLeft = (textX < fullInsideX); /* Special tweak for expRatio in pack mode: force all labels * left to prevent only a subset from being placed right: */ snapLeft |= (startsWith("expRatio", tg->tdb->type)); #ifdef IMAGEv2_NO_LEFTLABEL_ON_FULL if (theImgBox == NULL && snapLeft) #else///ifndef IMAGEv2_NO_LEFTLABEL_ON_FULL if (snapLeft) /* Snap label to the left. */ #endif ///ndef IMAGEv2_NO_LEFTLABEL_ON_FULL { textX = leftLabelX; assert(hvgSide != NULL); hvGfxUnclip(hvgSide); hvGfxSetClip(hvgSide, leftLabelX, yOff, fullInsideX - leftLabelX, tg->height); if(drawNameInverted) { int boxStart = leftLabelX + leftLabelWidth - 2 - nameWidth; hvGfxBox(hvgSide, boxStart, y, nameWidth+1, tg->heightPer - 1, color); hvGfxTextRight(hvgSide, leftLabelX, y, leftLabelWidth-1, tg->heightPer, MG_WHITE, font, name); } else hvGfxTextRight(hvgSide, leftLabelX, y, leftLabelWidth-1, tg->heightPer, labelColor, font, name); hvGfxUnclip(hvgSide); hvGfxSetClip(hvgSide, insideX, yOff, insideWidth, tg->height); } else { int pdfSlop=nameWidth/5; hvGfxUnclip(hvg); hvGfxSetClip(hvg, textX-1-pdfSlop, y, nameWidth+1+pdfSlop, tg->heightPer); if(drawNameInverted) { hvGfxBox(hvg, textX - 1, y, nameWidth+1, tg->heightPer-1, color); hvGfxTextRight(hvg, textX, y, nameWidth, tg->heightPer, MG_WHITE, font, name); } else hvGfxTextRight(hvg, textX, y, nameWidth, tg->heightPer, labelColor, font, name); hvGfxUnclip(hvg); hvGfxSetClip(hvg, insideX, yOff, insideWidth, tg->height); } } // now we're going to draw the boxes s = sf->start; int lastE = -1; int lastS = -1; int offY = y; int lineHeight = tg->lineHeight ; int qs, qe; int heightPer = tg->heightPer; int lastX = -1; int lastQEnd = 0; int lastLevel = -1; int e; qe = lastQEnd = 0; for (sf = (struct snakeFeature *)lf->components; sf != NULL; lastQEnd = qe, prevSf = sf, sf = sf->next) { qs = sf->qStart; qe = sf->qEnd; if (vis == tvDense) y = offY; else if ((vis == tvPack) || (vis == tvSquish)) y = offY + (sf->level * 1) * lineHeight; else if (vis == tvFull) y = offY + (sf->level * 2) * lineHeight; s = sf->start; e = sf->end; int sx, ex; if (!positiveRangeIntersection(winStart, winEnd, s, e)) continue; sx = round((double)((int)s-winStart)*scale) + xOff; ex = round((double)((int)e-winStart)*scale) + xOff; // color by strand static Color darkBlueColor = 0; static Color darkRedColor = 0; if (darkRedColor == 0) { //the light blue: rgb(149, 204, 252) //the light red: rgb(232, 156, 156) darkRedColor = hvGfxFindColorIx(hvg, 232,156,156); darkBlueColor = hvGfxFindColorIx(hvg, 149,204,252); } char *colorBy = cartOrTdbString(cart, tg->tdb, SNAKE_COLOR_BY, SNAKE_DEFAULT_COLOR_BY); extern Color getChromColor(char *name, struct hvGfx *hvg); if (sameString(colorBy, SNAKE_COLOR_BY_STRAND_VALUE)) color = (sf->orientation == -1) ? darkRedColor : darkBlueColor; else if (sameString(colorBy, SNAKE_COLOR_BY_CHROM_VALUE)) color = hashColor(sf->qName); else color = darkBlueColor; int w = ex - sx; if (w == 0) w = 1; assert(w > 0); char buffer[1024]; if (vis == tvFull) safef(buffer, sizeof buffer, "%d-%d",sf->qStart,sf->qEnd); else safef(buffer, sizeof buffer, "%s:%d-%d",sf->qName,sf->qStart,sf->qEnd); if (sx < insideX) { int olap = insideX - sx; sx = insideX; w -= olap; } char qAddress[4096]; if ((vis == tvFull) || (vis == tvPack) ) { safef(qAddress, sizeof qAddress, "qName=%s&qs=%d&qe=%d&qWidth=%d",sf->qName, qs, qe, winEnd - winStart); boundMapBox(hvg, s, e, sx+1, y, w-2, heightPer, tg->track, sf->qName, sf->qName, NULL, TRUE, qAddress); } hvGfxBox(hvg, sx, y, w, heightPer, color); // now draw the mismatches if we're at high enough resolution and have the sequence char *twoBitString = trackDbSetting(tg->tdb, "twoBit"); if (twoBitString && ((winBaseCount < showSnpWidth) && ((vis == tvFull) || (vis == tvPack)))) { static struct twoBitFile *tbf = NULL; static char *lastTwoBitString = NULL; static struct dnaSeq *seq = NULL; static char *lastQName = NULL; // sequence for chain snakes is in 2bit files which we cache if ((lastTwoBitString == NULL) || differentString(lastTwoBitString, twoBitString)) { if (tbf != NULL) { lastQName = NULL; twoBitClose(&tbf); } tbf = twoBitOpen(twoBitString); } // we're reading in the whole chrom if ((lastQName == NULL) || differentString(sf->qName, lastQName)) seq = twoBitReadSeqFrag(tbf, sf->qName, 0, 0); lastQName = sf->qName; lastTwoBitString = twoBitString; char *ourDna; ourDna = &seq->dna[sf->qStart]; int seqLen = sf->qEnd - sf->qStart; toUpperN(ourDna, seqLen); if (sf->orientation == -1) reverseComplement(ourDna,seqLen); // get the reference sequence char *refDna; struct dnaSeq *extraSeq = hDnaFromSeq(database, chromName, sf->start, sf->end, dnaUpper); refDna = extraSeq->dna; int si = s; char *ptr1 = refDna; char *ptr2 = ourDna; for(; si < e; si++,ptr1++,ptr2++) { // if mismatch! If reference is N ignore, if query is N, paint yellow if ( (*ptr1 != *ptr2) && !((*ptr1 == 'N') || (*ptr1 == 'n'))) { int misX1 = round((double)((int)si-winStart)*scale) + xOff; int misX2 = round((double)((int)(si+1)-winStart)*scale) + xOff; int w1 = misX2 - misX1; if (w1 < 1) w1 = 1; Color boxColor = MG_RED; if ((*ptr2 == 'N') || (*ptr2 == 'n')) boxColor = hvGfxFindRgb(hvg, &undefinedYellowColor); hvGfxBox(hvg, misX1, y, w1, heightPer, boxColor); } } // if we're zoomed to base level, draw sequence of mismatch if (zoomedToBaseLevel) { int mysx = round((double)((int)s-winStart)*scale) + xOff; int myex = round((double)((int)e-winStart)*scale) + xOff; int myw = myex - mysx; spreadAlignString(hvg, mysx, y, myw, heightPer, MG_WHITE, font, ourDna, refDna, seqLen, TRUE, FALSE); } } sf->drawn = TRUE; lastLevel = sf->level; //lastX = x; } if (vis != tvFull) return; // now we're going to draw the lines between the blocks lastX = -1; lastQEnd = 0; lastLevel = 0; qe = lastQEnd = 0; prevSf = NULL; for (sf = (struct snakeFeature *)lf->components; sf != NULL; lastQEnd = qe, prevSf = sf, sf = sf->next) { int y1, y2; int sx, ex; qs = sf->qStart; qe = sf->qEnd; if (lastLevel == sf->level) { y1 = offY + (lastLevel * 2) * lineHeight + lineHeight/2; y2 = offY + (sf->level * 2) * lineHeight + lineHeight/2; } else if (lastLevel > sf->level) { y1 = offY + (lastLevel * 2 ) * lineHeight; y2 = offY + (sf->level * 2 + 1) * lineHeight + lineHeight/3; } else { y1 = offY + (lastLevel * 2 + 1) * lineHeight - 1; y2 = offY + (sf->level * 2 ) * lineHeight - lineHeight/3;; } s = sf->start; e = sf->end; sx = round((double)((int)s-winStart)*scale) + xOff; ex = round((double)((int)e-winStart)*scale) + xOff; color = (sf->orientation == -1) ? MG_RED : MG_BLUE; if (lastX != -1) { char buffer[1024]; #define MG_ORANGE 0xff0082E6 int color = MG_GRAY; if (lastQEnd != qs) { long long queryInsertSize = llabs(lastQEnd - qs); if (queryInsertSize > 100) color = MG_ORANGE; } double queryGapNFrac = 0.0; double queryGapMaskedFrac = 0.0; // draw the vertical orange bars if there is an insert in the other sequence if ((winBaseCount < showSnpWidth) ) { if ((sf->orientation == 1) && (qs != lastQEnd) && (lastE == s)) { hvGfxLine(hvg, sx, y2 - lineHeight/2 , sx, y2 + lineHeight/2, MG_ORANGE); safef(buffer, sizeof buffer, "%dbp (%.1lf%% N, %.1lf%% masked)", qs - lastQEnd, queryGapNFrac*100, queryGapMaskedFrac*100); boundMapBox(hvg, s, e, sx, y2 - lineHeight/2, 1, lineHeight, tg->track, "foo", buffer, NULL, TRUE, NULL); safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //hvGfxTextCentered(hvg, sx - 10, y2 + lineHeight/2, 20, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } else if ((sf->orientation == -1) && (qs != lastQEnd) && (lastS == e)) { hvGfxLine(hvg, ex, y2 - lineHeight/2 , ex, y2 + lineHeight/2, MG_ORANGE); safef(buffer, sizeof buffer, "%dbp (%.1lf%% N, %.1lf%% masked)", qs - lastQEnd, queryGapNFrac*100, queryGapMaskedFrac*100); boundMapBox(hvg, s, e, ex, y2 - lineHeight/2, 1, lineHeight, tg->track, "foo", buffer, NULL, TRUE, NULL); safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //hvGfxTextCentered(hvg, ex - 10, y2 + lineHeight/2, 20, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } } // now draw the lines between blocks if ((!((lastX == sx) && (y1 == y2))) && (sf->drawn || ((prevSf != NULL) && (prevSf->drawn))) && (((lastE >= winStart) && (lastE <= winEnd)) || ((s > winStart) && (s < winEnd)))) { if (lastLevel == sf->level) { if (sf->orientation == 1) safef(buffer, sizeof buffer, "%dbp (%dbp in ref) (%.1lf%% N, %.1lf%% masked)", qs - lastQEnd, s - lastE, queryGapNFrac*100, queryGapMaskedFrac*100); else safef(buffer, sizeof buffer, "%dbp (%dbp in ref) (%.1lf%% N, %.1lf%% masked)", qs - lastQEnd, lastS - e, queryGapNFrac*100, queryGapMaskedFrac*100); if (sf->orientation == -1) { if (lastX != ex) { hvGfxLine(hvg, ex, y1, lastX, y2, color); boundMapBox(hvg, s, e, ex, y1, lastX-ex, 1, tg->track, "", buffer, NULL, TRUE, NULL); if (lastQEnd != qs) { safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //hvGfxTextCentered(hvg, ex, y2 + lineHeight/2, lastX-ex, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } } } else { if (lastX != sx) { hvGfxLine(hvg, lastX, y1, sx, y2, color); boundMapBox(hvg, s, e, lastX, y1, sx-lastX, 1, tg->track, "", buffer, NULL, TRUE, NULL); if (lastQEnd != qs) { safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //hvGfxTextCentered(hvg, lastX, y2 + lineHeight/2, sx-lastX, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } } } } else if (lastLevel > sf->level) { hvGfxLine(hvg, lastX, y1, sx, y2, color); hvGfxLine(hvg, sx, y2, sx, y2 - lineHeight - lineHeight/3, color); char buffer[1024]; safef(buffer, sizeof buffer, "%d-%d %dbp gap (%.1lf%% N, %.1lf%% masked)",prevSf->qStart,prevSf->qEnd, qs - lastQEnd, queryGapNFrac*100, queryGapMaskedFrac*100); boundMapBox(hvg, s, e, sx, y2 - lineHeight - lineHeight/3, 2, lineHeight + lineHeight/3, tg->track, "", buffer, NULL, TRUE, NULL); safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //if (lastQEnd != qs) //hvGfxTextCentered(hvg, sx - 10, y2 + lineHeight/2, 20, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } else { char buffer[1024]; safef(buffer, sizeof buffer, "%d-%d %dbp gap (%.1lf%% N, %.1lf%% masked)",prevSf->qStart,prevSf->qEnd, qs - lastQEnd, queryGapNFrac*100, queryGapMaskedFrac*100); if (sf->orientation == -1) { hvGfxLine(hvg, lastX-1, y1, ex, y2, color); hvGfxLine(hvg, ex, y2, ex, y2 + lineHeight , color); boundMapBox(hvg, s, e, ex-1, y2, 2, lineHeight , tg->track, "", buffer, NULL, TRUE, NULL); safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //if (lastQEnd != qs) //hvGfxTextCentered(hvg, ex - 10, y2 + lineHeight, 20, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } else { hvGfxLine(hvg, lastX-1, y1, sx, y2, color); hvGfxLine(hvg, sx, y2, sx, y2 + lineHeight , color); boundMapBox(hvg, s, e, sx-1, y2, 2, lineHeight , tg->track, "", buffer, NULL, TRUE, NULL); safef(buffer, sizeof buffer, "%dbp", qs - lastQEnd); //if (lastQEnd != qs) //hvGfxTextCentered(hvg, sx - 10, y2 + lineHeight, 20, INSERT_TEXT_HEIGHT, MG_ORANGE, font, buffer); } } } } if (sf->orientation == -1) lastX = sx; else lastX = ex; lastS = s; lastE = e; lastLevel = sf->level; lastQEnd = qe; } } struct cartOptions { enum chainColorEnum chainColor; /* ChromColors, ScoreColors, NoColors */ int scoreFilter ; /* filter chains by score if > 0 */ }; static void fixItems(struct linkedFeatures *lf) // put all chain blocks from a single query chromosome into one // linkedFeatures structure { struct linkedFeatures *firstLf, *next; struct snakeFeature *sf, *nextSf; firstLf = lf; for (;lf; lf = next) { next = lf->next; if (!sameString(firstLf->name, lf->name) && (lf->components != NULL)) { slSort(&firstLf->components, snakeFeatureCmpQStart); firstLf->next = lf; firstLf = lf; } for (sf = (struct snakeFeature *)lf->components; sf != NULL; sf = nextSf) { sf->qName = lf->name; sf->orientation = lf->orientation; + sf->score = lf->score; nextSf = sf->next; if (firstLf != lf) { lf->components = NULL; slAddHead(&firstLf->components, sf); } } } if (firstLf != NULL) { slSort(&firstLf->components, snakeFeatureCmpQStart); firstLf->next = 0; } } void snakeDrawLeftLabels() { } static void makeSnakeFeatures(struct linkedFeatures *lf) +/* allocate an info structure for every block in the alignment. */ { for(; lf; lf = lf->next) { struct simpleFeature *sf = lf->components; struct snakeFeature *sfList = NULL; for(;sf; sf = sf->next) { struct snakeFeature *this; AllocVar(this); *(struct simpleFeature *)this = *sf; this->orientation = lf->orientation; slAddHead(&sfList, this); } lf->components = (struct simpleFeature *)sfList; } } +static void makeSingleCoverage(struct linkedFeatures *lf) +/* Because snakes are a display of O+O of the query sequence projected onto the reference sequence, + * we need to only have one copy of each piece of the query sequence. */ +{ +for(; lf; lf = lf->next) // for each query sequence + { + slSort(&lf->components, snakeFeatureCmpScore); // we want to choose the higher scoring alignements + Bits *qBits = bitAlloc(lf->qSize); // the bit array we use to keep track if we've seen this query sequence + struct simpleFeature *sf = lf->components; + struct simpleFeature *prevSf = NULL; + + for(; sf ; sf = sf->next) + { + // have we seen this query range? + if (bitCountRange(qBits, sf->qStart, sf->qEnd - sf->qStart)) + { + // if we've seen this query sequence, then delete it + if (prevSf == NULL) + lf->components = sf->next; + else + prevSf->next = sf->next; + } + else + { + bitSetRange(qBits, sf->qStart, sf->qEnd - sf->qStart); + prevSf = sf; + } + } + + slSort(&lf->components, snakeFeatureCmpQStart); + } +} + void maybeLoadSnake(struct track *track) /* check to see if we're doing snakes and if so load the correct methods. */ { boolean doSnake = cartOrTdbBoolean(cart, track->tdb, "doSnake", FALSE); if (doSnake) { track->drawLeftLabels = snakeDrawLeftLabels; track->itemHeight = snakeItemHeight; track->totalHeight = snakeHeight; track->drawItemAt = snakeDrawAt; track->drawItems = snakeDraw; slSort(&track->items, linkedFeaturesCmpName); makeSnakeFeatures(track->items); fixItems(track->items); + makeSingleCoverage(track->items); track->visibility = track->limitedVis = tvFull; track->canPack = FALSE; } }