2de160a334786c2168e379c5f02838f4aed40a3b chmalee Thu Jul 27 12:20:35 2023 -0700 Fix two separate hgPcr related bugs. One, the primers entered by the user do not have to match the genome/transcript fully, and so when looking up the primers for showing mismatches we need to ensure we find them by relying on the psl qName rather than the primer file that holds the users input. Second, when clicking through to hgc, ensure that we check that the tStart and tEnd of the psl match what was clicked on so the right amplicon sequence can be shown, in the case that two separate searches that hit the same transcript are visible, refs #31784 diff --git src/hg/hgTracks/cds.c src/hg/hgTracks/cds.c index 47a2aed..79bfd4f 100644 --- src/hg/hgTracks/cds.c +++ src/hg/hgTracks/cds.c @@ -1,2280 +1,2296 @@ /* cds.c - code for coloring of bases, codons, or alignment differences. */ /* Copyright (C) 2014 The Regents of the University of California * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "hCommon.h" #include "hash.h" #include "jksql.h" #include "memalloc.h" #include "dystring.h" #include "memgfx.h" #include "hvGfx.h" #include "dnaseq.h" #include "dnautil.h" #include "hdb.h" #include "psl.h" #include "fa.h" #include "genePred.h" #include "cds.h" #include "genbank.h" #include "twoBit.h" #include "cacheTwoBit.h" #include "hgTracks.h" #include "cdsSpec.h" #include "axt.h" #include "lrg.h" /* * WARNING: this code is incomprehensible: * - variables named codon often contain amino acids, not condons * - it does two passes, one undocumented function encodes both a * color and an amino acid into the struct simpleFeature grayIx * field and this is decoded in the second pass. * - baseColorDrawItem doesn't draw a item, it draws a single codon, * or a maybe even a single base. * - there are many assumptions and state shared between this module * and the simpleTracks.c. * may The Force be with you.. */ #ifndef GBROWSE #include "pcrResult.h" #endif /* GBROWSE */ Color lighterShade(struct hvGfx *hvg, Color color, double percentLess); /* Find a color which is a percentless 'lighter' shade of color */ /* Forward declaration */ /* Array of colors used in drawing codons/bases/differences: */ Color cdsColor[CDS_NUM_COLORS]; boolean cdsColorsMade = FALSE; Color getCdsColor(int index) { assert(index < CDS_NUM_COLORS); return cdsColor[index]; } boolean pslTargetToQueryRangeMap(struct psl *psl, int tStart, int tEnd, int *retQStart, int *retQEnd) /* Comes up with qStart/qEnd that corresponds to tStart/tEnd in psl. Returns * FALSE if there's no corresponding part for target in query */ { if (rangeIntersection(tStart, tEnd, psl->tStart, psl->tEnd) <= 0) return FALSE; int qStart = psl->qStart, qEnd = psl->qEnd; boolean foundStart = FALSE, foundEnd = FALSE; int lastBqEnd = psl->qStart; int i; for (i=0; i<psl->blockCount; ++i) { int bSize = psl->blockSizes[i]; int bqStart = psl->qStarts[i]; int bqEnd = bqStart + bSize; int btStart = psl->tStarts[i]; int btEnd = btStart + bSize; if (!foundStart) { if (btStart <= tStart && tStart < btEnd) { qStart = bqStart + (tStart - btStart); foundStart = TRUE; } else if (btStart >= tStart) { qStart = bqStart; foundStart = TRUE; } } if (!foundEnd) { if (btStart < tEnd && tEnd <= btEnd) { qEnd = bqStart + (tEnd - btStart); foundEnd = TRUE; break; } else if (btStart >= tEnd) { qEnd = lastBqEnd; foundEnd = TRUE; break; } } lastBqEnd = bqEnd; } *retQStart = qStart; *retQEnd = qEnd; return TRUE; } static void drawScaledBoxWithText(struct hvGfx *hvg, int chromStart, int chromEnd, double scale, int xOff, int y, int height, Color color, int score, MgFont *font, char *text, bool zoomed, int winStart, int maxPixels, boolean isCoding, boolean justifyString) /* Draw a box scaled from chromosome to window coordinates with a codon or set of 3 or less bases drawn in the box. */ { /*first draw the box itself*/ drawScaledBox(hvg, chromStart, chromEnd, scale, xOff, y, height, color); /*draw text in box if space, and align properly for codons or DNA*/ if (zoomed) { int i; Color textColor = hvGfxContrastingColor(hvg, color); int x1, x2, w; x1 = round((double)(chromStart-winStart)*scale) + xOff; x2 = round((double)(chromEnd-winStart)*scale) + xOff; if (x2 >= maxPixels) x2 = maxPixels - 1; w = x2-x1; if (w < 1) w = 1; if (chromEnd - chromStart == 3 && isCoding) { if (justifyString) spreadBasesString(hvg, x1, y, w, height, textColor, font, text, strlen(text), TRUE); else hvGfxTextCentered(hvg, x1, y, w, height, textColor, font, text); } else if (chromEnd - chromStart < 3 && isCoding) { if (justifyString) spreadBasesString(hvg, x1, y, w, height, cdsColor[CDS_PARTIAL_CODON], font, text, strlen(text), TRUE); else hvGfxTextCentered(hvg, x1, y, w, height, cdsColor[CDS_PARTIAL_CODON], font, text); } else { int thisX,thisX2; char c[2]; c[1] = '\0'; int iMin = max(0, (winStart-chromStart)); int iMax = min((chromEnd-chromStart), (winEnd-chromStart)); for (i=iMin; i<iMax; i++) { if (text[i] == ' ') continue; c[0] = text[i]; thisX = round((double)(chromStart+i-winStart)*scale) + xOff; thisX2 = round((double)(chromStart+1+i-winStart)*scale) + xOff; hvGfxTextCentered(hvg, thisX, y, thisX2-thisX, height, textColor,font,c); } } } } static int convertCoordUsingPsl( int s, struct psl *psl ) /*return query position corresponding to target position in one * of the coding blocks of a psl file, or return -1 if target position * is not in an exon of this psl entry*/ { int i; int idx = -1; unsigned *qStarts = psl->qStarts; unsigned tStart = 0; unsigned thisQStart = 0; if (psl == NULL) errAbort("NULL passed to convertCoordUsingPsl<br>\n"); if (s < psl->tStart || s >= psl->tEnd) { warn("Position %d is outside of this psl entry q(%s:%d-%d), t(%s:%d-%d)." "<br>\n", s, psl->qName, psl->qStart, psl->qEnd, psl->tName, psl->tStart, psl->tEnd); return(-1); } for ( i=0; i<psl->blockCount; i++ ) { unsigned tEnd; tStart = psl->tStarts[i]; tEnd = tStart + psl->blockSizes[i]; if (psl->strand[1] == '-') reverseUnsignedRange(&tStart, &tEnd, psl->tSize); if (s >= tStart && s < tEnd) { idx = i; break; } } if (idx < 0) return(-1); if(psl->strand[1] == '-') thisQStart = psl->qSize - (qStarts[idx]+psl->blockSizes[idx]); else thisQStart = qStarts[idx]; return(thisQStart + (s - tStart)); } /* Calls to hDnaFromSeq are rather expensive for 2bit, so cache genomic sequence */ struct genoCacheWindow { struct genoCacheWindow *next; struct window *window; char *initedTrack; struct dnaSeq *cachedGenoDna; int cachedGenoStart; int cachedGenoEnd; }; static struct genoCacheWindow *gcWindows = NULL; static struct genoCacheWindow *gcWindow = NULL; static struct genoCacheWindow *gcWindowOld = NULL; static bool setGcWindow() /* scan genoCache windows. create new one if not found */ { if (gcWindow && gcWindow->window == currentWindow) return FALSE; gcWindowOld = gcWindow; for (gcWindow = gcWindows; gcWindow; gcWindow = gcWindow->next) { if (gcWindow->window == currentWindow) { return TRUE; } } AllocVar(gcWindow); gcWindow->window = currentWindow; slAddTail(&gcWindows, gcWindow); return TRUE; } static char *initedTrack = NULL; static struct dnaSeq *cachedGenoDna = NULL; static int cachedGenoStart = 0; static int cachedGenoEnd = 0; static void setGc() /* set up globals for the current window */ { if (setGcWindow()) { if (gcWindowOld) { gcWindowOld->initedTrack = initedTrack; gcWindowOld->cachedGenoDna = cachedGenoDna; gcWindowOld->cachedGenoStart = cachedGenoStart; gcWindowOld->cachedGenoEnd = cachedGenoEnd; } initedTrack = gcWindow->initedTrack; cachedGenoDna = gcWindow->cachedGenoDna; cachedGenoStart = gcWindow->cachedGenoStart; cachedGenoEnd = gcWindow->cachedGenoEnd; } } static void getLinkedFeaturesSpan(struct linkedFeatures *lfList, int *retStart, int *retEnd, boolean isSeries) /* Find the overall lowest and highest coords in lfList. If any items hang off the * edge of the window, we will end up with coords <winStart and/or >winEnd which is * what we want. */ { int start = winStart, end = winEnd; if (isSeries) { struct linkedFeaturesSeries *lfs, *lfsList = (struct linkedFeaturesSeries *)lfList; for (lfs = lfsList; lfs != NULL; lfs = lfs->next) { if (lfs->start < start) start = lfs->start; if (lfs->end > end) end = lfs->end; } } else { struct linkedFeatures *lf; for (lf = lfList; lf != NULL; lf = lf->next) { if (lf->start < start) start = lf->start; if (lf->end > end) end = lf->end; } } if (retStart) *retStart = start; if (retEnd) *retEnd = end; } static char *getCachedDna(int chromStart, int chromEnd) /* Return a pointer into our cached genomic dna. chromEnd is just for * bounds-checking (honor system). Do not change or free the return value! */ { setGc(); if (!initedTrack || !cachedGenoDna) errAbort("getCachedDnaAt called before baseColorInitTrack?!"); if (chromStart < cachedGenoStart || chromEnd > cachedGenoEnd) errAbort("getCachedDnaAt: coords %d,%d are out of cached range %d,%d", chromStart, chromEnd, cachedGenoStart, cachedGenoEnd); return &(cachedGenoDna->dna[chromStart - cachedGenoStart]); } static void getNextCodonDna(char *retStr, int n, struct genePred *gp, int startI, boolean posStrand) /* Get at most n bases from coding exons following exon startI. */ { int i, j, thisN; int cdsExonStart, cdsExonEnd; char *codonDna; for (i = 0; i < n; i++) retStr[i] = 'N'; retStr[n] = '\0'; if (posStrand) { for (thisN = 0, i = startI+1; thisN < n && i < gp->exonCount; i++) { // get dna for exon to the right codonDna = getCachedDna(gp->exonStarts[i], gp->exonEnds[i]); cdsExonStart = gp->exonStarts[i]; cdsExonEnd = gp->exonEnds[i]; if (gp->cdsStart < gp->txEnd && cdsExonStart < gp->cdsStart) cdsExonStart = gp->cdsStart; if (gp->cdsEnd > gp->txStart && cdsExonEnd > gp->cdsEnd) cdsExonEnd = gp->cdsEnd; for (j=0; j < (cdsExonEnd - cdsExonStart); j++) { retStr[thisN] = codonDna[j+cdsExonStart-gp->exonStarts[i]]; thisN++; if (thisN >= n) break; } } } else { for (thisN = n-1, i = startI-1; thisN >= 0 && i >= 0; i--) { // get dna for exon to the left codonDna = getCachedDna(gp->exonStarts[i], gp->exonEnds[i]); cdsExonStart = gp->exonStarts[i]; cdsExonEnd = gp->exonEnds[i]; if (gp->cdsStart < gp->txEnd && cdsExonStart < gp->cdsStart) cdsExonStart = gp->cdsStart; if (gp->cdsEnd > gp->txStart && cdsExonEnd > gp->cdsEnd) cdsExonEnd = gp->cdsEnd; for (j=0; j < (cdsExonEnd - cdsExonStart); j++) { retStr[thisN] = codonDna[cdsExonEnd-j-1-gp->exonStarts[i]]; thisN--; if (thisN < 0) break; } } } } #ifdef UNUSED #endif /* UNUSED */ static void drawVertLine(struct linkedFeatures *lf, struct hvGfx *hvg, int chromStart, int xOff, int y, int height, double scale, Color color) /* Draw a 1-pixel wide vertical line at the given chromosomal coord. * The line is 0 bases wide (chromStart==chromEnd) but that doesn't * matter if we're zoomed out to >1base/pixel, so this is OK for diffs * when zoomed way out and for insertion points at any scale. */ { int thisX = round((double)(chromStart-winStart)*scale) + xOff; int thisY = y; height -= 1; int thisHeight = height; if ((chromStart < lf->tallStart) || (chromStart > lf->tallEnd)) { /* adjust for UTR. WARNING: this duplicates shortOff & shortHeight * calculations in linkedFeaturesDrawAt */ thisY += height/4; thisHeight = height - height/2; } hvGfxBox(hvg, thisX-1, thisY, 2, thisHeight, color); } static void drawMidNumber(struct linkedFeatures *lf, struct hvGfx *hvg, int chromStart, int xOff, int y, int height, double scale, Color color, MgFont *font, int size) /* Draw a short string encoding size around chromStart */ /* Draw a 1-pixel wide vertical line at the given chromosomal coord. * The line is 0 bases wide (chromStart==chromEnd) but that doesn't * matter if we're zoomed out to >1base/pixel, so this is OK for diffs * when zoomed way out and for insertion points at any scale. */ { char sizeString[32]; safef(sizeString, sizeof(sizeString), "%d", size); drawScaledBoxLabel(hvg, chromStart-1, chromStart+1, scale, xOff, y, height, color, font, sizeString); } #ifdef SOON static void drawLeftNumber(struct linkedFeatures *lf, struct hvGfx *hvg, int chromStart, int xOff, int y, int height, double scale, Color color, MgFont *font, int size) /* Draw a short string encoding size around chromStart */ /* Draw a 1-pixel wide vertical line at the given chromosomal coord. * The line is 0 bases wide (chromStart==chromEnd) but that doesn't * matter if we're zoomed out to >1base/pixel, so this is OK for diffs * when zoomed way out and for insertion points at any scale. */ { drawMidNumber(lf, hvg, chromStart+1, xOff, y, height, scale, color, font, size); } #endif /* SOON */ #ifdef SOON static void drawRightNumber(struct linkedFeatures *lf, struct hvGfx *hvg, int chromStart, int xOff, int y, int height, double scale, Color color, MgFont *font, int size) /* Draw a short string encoding size around chromStart */ /* Draw a 1-pixel wide vertical line at the given chromosomal coord. * The line is 0 bases wide (chromStart==chromEnd) but that doesn't * matter if we're zoomed out to >1base/pixel, so this is OK for diffs * when zoomed way out and for insertion points at any scale. */ { drawMidNumber(lf, hvg, chromStart-1, xOff, y, height, scale, color, font, size); } #endif /* SOON */ static void drawCdsDiffBaseTickmarksOnly(struct track *tg, struct linkedFeatures *lf, struct hvGfx *hvg, int xOff, int y, double scale, int heightPer, struct dnaSeq *qSeq, int qOffset, struct psl *psl, int winStart) /* Draw thin vertical red lines only where mRNA bases differ from genomic. * This assumes that lf has been drawn already, we're zoomed out past * zoomedToBaseLevel, we're not in dense mode etc. */ { struct simpleFeature *sf = NULL; char *winDna = getCachedDna(winStart, winEnd); Color c = cdsColor[CDS_STOP]; // check if we need a contrasting color instead of default 'red' (CDS_STOP) char *tickColor = NULL; if ( tg->itemColor && (tickColor = trackDbSetting(tg->tdb, "baseColorTickColor"))) { Color ci = tg->itemColor(tg, lf, hvg); if (sameString(tickColor, "contrastingColor")) c = hvGfxContrastingColor(hvg, ci); else if (sameString(tickColor, "lighterShade")) c = lighterShade(hvg, ci, 1.5); } for (sf = lf->components; sf != NULL; sf = sf->next) { int s = max(winStart, sf->start); int e = min(winEnd, sf->end); if (s > winEnd || e < winStart) continue; if (e > s) { int mrnaS = -1, mrnaE = 0; if (psl) { // mrnaS = convertCoordUsingPsl(s, psl); pslTargetToQueryRangeMap(psl, s, e, &mrnaS, &mrnaE); } else { mrnaS = sf->qStart + (s - sf->start); } if(mrnaS >= 0) { int i; for (i=0; i < (e - s); i++) { if (qSeq->dna[mrnaS+i-qOffset] != winDna[s-winStart+i]) drawScaledBox(hvg, s+i, s+i+1, scale, xOff, y+1, heightPer-2, c); } } } } } static void maskDiffString( char *retStr, char *s1, char *s2, char mask, int size) /*copies s1, masking off similar characters, and returns result into retStr. *if strings are of different size it stops after s1 is done.*/ { memset(retStr, mask, size); int i; for (i=0; i < size; i++) { char c = s1[i]; if (c != s2[i]) { retStr[i] = c; } } retStr[i] = '\0'; } Color lighterShade(struct hvGfx *hvg, Color color, double percentLess) /* Get lighter shade of a color, with a variable level */ { struct rgbColor rgbColor = hvGfxColorIxToRgb(hvg, color); rgbColor.r = (int)((rgbColor.r+127)/percentLess); rgbColor.g = (int)((rgbColor.g+127)/percentLess); rgbColor.b = (int)((rgbColor.b+127)/percentLess); return hvGfxFindColorIx(hvg, rgbColor.r, rgbColor.g, rgbColor.b); } // We need an uppercase alpha character to represent the stop codon in // the alpha-offset grayIx scheme. J is in neither dnautil.c's // codonTable nor axt.c's blosum62 peptide scoring matrix. X was // formerly used to represent stop codon, but unfortunately it also is // the error return from dnautil.c's lookupCodon (and it is in // blosum62, so maybe it is a valid peptide in some contexts?). #define GRAYIX_STOP_CODON_ALPHA 'J' // In addition to the alpha-offset encoding of peptide and alternating // shade, grayIx can take on several negative values to represent colors // for special cases: #define GRAYIX_CDS_START -1 #define GRAYIX_CDS_ERROR -2 #define GRAYIX_CDS_STOP -3 #define GRAYIX_CDS_SYN_PROT -4 static Color colorAndCodonFromGrayIx(struct hvGfx *hvg, char *codon, int grayIx, Color ixColor) /*convert grayIx value to color and codon which * are both encoded in the grayIx*/ { Color color; if (grayIx == GRAYIX_CDS_ERROR) { color = cdsColor[CDS_ERROR]; sprintf(codon,"X"); } else if (grayIx == GRAYIX_CDS_START) { color = cdsColor[CDS_START]; sprintf(codon,"M"); } else if (grayIx == GRAYIX_CDS_STOP) { color = cdsColor[CDS_STOP]; sprintf(codon,"*"); } else if (grayIx == GRAYIX_CDS_SYN_PROT) { color = cdsColor[CDS_SYN_PROT]; sprintf(codon,"*"); } #ifdef LOWELAB else if(grayIx == - 'V') { color = cdsColor[CDS_ALT_START]; sprintf(codon,"%c",'V'); } else if(grayIx == - 'L') { color = cdsColor[CDS_ALT_START]; sprintf(codon,"%c",'L'); } #endif else if (grayIx <= 26) { color = ixColor; sprintf(codon,"%c",grayIx + 'A' - 1); if (codon[0] == GRAYIX_STOP_CODON_ALPHA) codon[0] = '*'; } else if (grayIx > 78) // ribosomal slip begin { color = cdsColor[CDS_RIBO_SLIP1]; sprintf(codon,"%c",grayIx - 78 + 'A' - 1); } else if (grayIx > 52) // ribosomal slip end { color = cdsColor[CDS_RIBO_SLIP2]; sprintf(codon,"%c",grayIx - 52 + 'A' - 1); } else if (grayIx > 26) { color = lighterShade(hvg, ixColor,1.5); sprintf(codon,"%c",grayIx - 26 + 'A' - 1); if (codon[0] == GRAYIX_STOP_CODON_ALPHA) codon[0] = '*'; } else { errAbort("colorAndCodonFromGrayIx: invalid grayIx %d", grayIx); color = cdsColor[CDS_ERROR]; sprintf(codon,"X"); } return color; } static char baseColorLookupCodon(DNA *dna) /* Call dnautil's lookupCodon, but translate stop codon '\0' to '*' for display. */ { char peptide; if (isMito(chromName)) peptide = lookupMitoCodon(dna); else peptide = lookupCodon(dna); if (peptide == '\0') peptide = '*'; return peptide; } static int peptideToGrayIx(char peptide, boolean codonFirstColor) /* Encode peptide (a letter or '*') and alternating gray shade into our alpha-offset scheme. */ { if (peptide == '*') peptide = GRAYIX_STOP_CODON_ALPHA; if (codonFirstColor) return(peptide - 'A' + 1); else return(peptide - 'A' + 1 + 26); } static int codonToGrayIx(DNA *dna, bool codonFirstColor, boolean *foundStart, boolean reverse, boolean colorStopStart) /* Return grayIx encoding the codon and color (or alternating shades). */ { if (reverse) reverseComplement(dna,strlen(dna)); char codonChar = baseColorLookupCodon(dna); if (codonChar == 'M' && foundStart != NULL && !(*foundStart)) *foundStart = TRUE; #ifdef LOWELAB if(sameString(dna,"GTG")) { return -'V'; } if(sameString(dna,"TTG")) { return -'L'; } #endif if (codonChar == '*') { if (colorStopStart) return(GRAYIX_CDS_STOP); else return peptideToGrayIx(codonChar, codonFirstColor); } else if (codonChar == 'X') return(GRAYIX_CDS_ERROR); // bad input to lookupCodon, e.g. 'n' base else if (colorStopStart && codonChar == 'M') return(GRAYIX_CDS_START); else return peptideToGrayIx(codonChar, codonFirstColor); } static boolean protEquivalent(int aa1, int aa2) /* returns TRUE if amino acids have a positive blosum62 score i.e. I, V or R, K else FALSE Note: X is a valid code in the default scheme (blosum62), so don't pass in dnautil.c's X which is an error return... */ { static struct axtScoreScheme *ss = NULL; if (ss == NULL) ss = axtScoreSchemeProteinDefault(); if ((ss->matrix[aa1][aa2]) > 0) return TRUE; else return FALSE; } static int mrnaCodonToGrayIx(DNA *rna, char genomicCodon, bool codonFirstColor) /* Difference ==> red, otherwise keep the alternating shades. */ { char rnaCodon = baseColorLookupCodon(rna); if (genomicCodon != rnaCodon) { if (genomicCodon != 'X' && rnaCodon != 'X' && protEquivalent(genomicCodon, rnaCodon)) return(GRAYIX_CDS_SYN_PROT); // yellow, "synonymous" protein else return(GRAYIX_CDS_STOP); } else return peptideToGrayIx(genomicCodon, codonFirstColor); } static void getGenbankCds(char *acc, struct genbankCds* cds) /* Get cds start and stop from genbank tables, if available. Otherwise it * does nothing */ { static boolean first = TRUE, haveGbCdnaInfo = FALSE; struct sqlConnection *conn = hAllocConn(database); if (first) { haveGbCdnaInfo = sqlTableExists(conn, gbCdnaInfoTable); first = FALSE; } if (haveGbCdnaInfo) { char query[4096], buf[4096], *cdsStr; sqlSafef(query, sizeof query, "select c.name from %s g,%s c where (acc = '%s') and (g.cds = c.id)", gbCdnaInfoTable, cdsTable, acc); cdsStr = sqlQuickQuery(conn, query, buf, sizeof(buf)); if (cdsStr != NULL) genbankCdsParse(cdsStr, cds); } hFreeConn(&conn); } static void getCdsFromTbl(char *acc, char *baseColorSetting, struct genbankCds* cds) /* Get CDS from a specified table, doing nothing if not found */ { char *p = skipToSpaces(baseColorSetting); char *cdsSpecTbl = skipLeadingSpaces(p); if (*cdsSpecTbl == '\0') errAbort("%s table requires a table name as an argument", BASE_COLOR_USE_CDS); struct sqlConnection *conn = hAllocConnDbTbl(cdsSpecTbl, &cdsSpecTbl, database); // allow multiple, but only use the first, since transMapGene table might have // multiple entries for same gene from different source dbs. struct cdsSpec *cdsSpec = sqlQueryObjs(conn, (sqlLoadFunc)cdsSpecLoad, sqlQueryMulti, "SELECT * FROM %s WHERE id=\"%s\"", cdsSpecTbl, acc); hFreeConn(&conn); if (cdsSpec != NULL) genbankCdsParse(cdsSpec->cds, cds); if (cds != NULL) cdsSpecFreeList(&cdsSpec); } static void getPslCds(struct psl *psl, struct track *tg, struct linkedFeatures *lf, struct genbankCds *cds) /* get CDS defintion for a PSL */ { ZeroVar(cds); if (pslIsProtein(psl)) { cds->start=0; cds->end=psl->qSize*3; cds->startComplete = TRUE; cds->endComplete = TRUE; } else if (startsWith("bigPsl", tg->tdb->type)) { if (lf->cds) genbankCdsParse(lf->cds, cds); } else { char *setting = trackDbSetting(tg->tdb, BASE_COLOR_USE_CDS); char *dataName = getItemDataName(tg, psl->qName); if ((setting != NULL) && startsWith("table", setting)) getCdsFromTbl(dataName, setting, cds); else getGenbankCds(dataName, cds); } } static void getHiddenGaps(struct psl *psl, unsigned *gaps) /*return the insertions in the query sequence of a psl that are 'hidden' in the browser. This lets these insertions be indicated with the color orange.*/ { int i; gaps[0] = psl->qStarts[0] - psl->qStart; for (i=1; i<psl->blockCount; i++) { gaps[i] = psl->qStarts[i] - (psl->qStarts[i-1] + psl->blockSizes[i-1]); } gaps[psl->blockCount] = psl->qSize - (psl->qStarts[psl->blockCount-1] + psl->blockSizes[psl->blockCount-1]); } static struct simpleFeature *splitPslByCodon(char *chrom, struct linkedFeatures *lf, struct psl *psl, int sizeMul, boolean isXeno, int maxShade, enum baseColorDrawOpt drawOpt, boolean colorStopStart, struct track *tg) { struct simpleFeature *sfList = NULL; /* if cds not in genbank or not parsable - revert to normal*/ /* if showing bases we don't want to have thin bars ala genePred format */ if (drawOpt == baseColorDrawItemBases || drawOpt == baseColorDrawDiffBases) { int grayIx = pslGrayIx(psl, isXeno, maxShade); sfList = sfFromPslX(psl, grayIx, sizeMul); } else { /* Previous code didn't use exon frames for baseColorDrawGenomicCodons. * This meant simply counting off aligned bases to define frames. It * didn't work very well for TransMap alignments and not clear that its * the right thing to do for any alignment. By using exonFrames for * genomic codons, this is letting the query sequence define the frame. */ struct genbankCds cds; getPslCds(psl, tg, lf, &cds); int insertMergeSize = -1; unsigned opts = genePredCdsStatFld|genePredExonFramesFld; struct genePred *gp = genePredFromPsl2(psl, opts, &cds, insertMergeSize); lf->start = gp->txStart; lf->end = gp->txEnd; lf->tallStart = gp->cdsStart; lf->tallEnd = gp->cdsEnd; sfList = baseColorCodonsFromGenePred(lf, gp, colorStopStart, FALSE); genePredFree(&gp); } return(sfList); } struct simpleFeature *baseColorCodonsFromPsl(struct linkedFeatures *lf, struct psl *psl, int sizeMul, boolean isXeno, int maxShade, enum baseColorDrawOpt drawOpt, struct track *tg) /* Given an lf and the psl from which the lf was constructed, * return a list of simpleFeature elements, one per codon (or partial * codon if the codon falls on a gap boundary. sizeMul, isXeno and maxShade * are for defaulting to one-simpleFeature-per-exon if cds is not found. */ { boolean colorStopStart = (drawOpt != baseColorDrawDiffCodons); struct simpleFeature *sfList = splitPslByCodon(psl->tName, lf, psl, sizeMul, isXeno, maxShade, drawOpt, colorStopStart, tg); slReverse(&sfList); return sfList; } enum baseColorDrawOpt baseColorGetDrawOpt(struct track *tg) /* Determine what base/codon coloring option (if any) has been selected * in trackDb/cart, and gate with zoom level. */ { // N.B. the way the option is chosen has become insanely complex enum baseColorDrawOpt ret = baseColorDrawOff; if (cart == NULL || tg == NULL || tg->tdb == NULL) return ret ; ret = baseColorDrawOptEnabled(cart, tg->tdb); boolean showDiffBasesAllScales = (trackDbSetting(tg->tdb, "showDiffBasesAllScales") != NULL); float showDiffBasesMaxZoom = trackDbFloatSettingOrDefault(tg->tdb, "showDiffBasesMaxZoom", -1.0); if ((showDiffBasesMaxZoom >= 0) && ((basesPerPixel > showDiffBasesMaxZoom) || (showDiffBasesMaxZoom == 0.0))) showDiffBasesAllScales = FALSE; enum baseColorDrawOpt zoomOutDefault = baseColorDrawOff; if (trackDbSetting(tg->tdb, "showCdsAllScales") != NULL) { float showCdsMaxZoom = trackDbFloatSettingOrDefault(tg->tdb, "showCdsMaxZoom", -1.0); if ((showCdsMaxZoom < 0.0) || ((basesPerPixel <= showCdsMaxZoom) && (showCdsMaxZoom != 0.0))) zoomOutDefault = baseColorDrawCds; #if 0 /* FIXME: this doesn't actually work, it enables maxItems, but then * doesn't does into dense at all. need to figure this out, but in the * mean time, just usine showCdsMaxZoom to keep performance sane */ // don't show CDS if in dense or squish if (zoomOutDefault == baseColorDrawCds) { enum trackVisibility vis = limitVisibility(tg); if ((vis == tvDense) || (vis == tvSquish)) zoomOutDefault = baseColorDrawOff; } #endif } /* Gate with zooming constraints: */ if (!zoomedToCdsColorLevel && (ret == baseColorDrawGenomicCodons || ret == baseColorDrawItemCodons)) ret = zoomOutDefault; if (!zoomedToBaseLevel && (ret == baseColorDrawItemBases)) ret = zoomOutDefault; if (!(zoomedToCdsColorLevel || showDiffBasesAllScales) && ret == baseColorDrawDiffCodons) ret = zoomOutDefault; if (!(zoomedToBaseLevel || showDiffBasesAllScales) && ret == baseColorDrawDiffBases) ret = zoomOutDefault; return ret; } static struct dnaSeq *maybeGetUserSeq(char *qName) /* Look in user's blat sequence file for qName. */ { #ifndef GBROWSE static struct hash *userSeqHash = NULL; struct dnaSeq *userSeq = NULL; char *ss = cartOptionalString(cart, "ss"); if (ss == NULL || !ssFilesExist(ss)) return NULL; if (userSeqHash == NULL) { struct dnaSeq *seqList = NULL, *seq; char *faFileName, *pslFileName; parseSs(ss, &pslFileName, &faFileName); seqList = faReadAllMixed(faFileName); userSeqHash = hashNew(0); for (seq = seqList; seq != NULL; seq = seq->next) { hashAdd(userSeqHash, seq->name, seq); } } userSeq = hashFindVal(userSeqHash, qName); if (userSeq == NULL) return NULL; else return cloneDnaSeq(userSeq); #else return NULL; #endif } #ifndef GBROWSE static struct dnaSeq *maybeGetPcrResultSeq(struct linkedFeatures *lf) /* Return (if possible) the primer sequences concatenated with the * target sequence between where they match. */ { struct dnaSeq *seq = NULL; char *pslFileName, *primerFileName; struct targetDb *target; if (! pcrResultParseCart(database, cart, &pslFileName, &primerFileName, &target)) return NULL; -char *fPrimer=NULL, *rPrimer=NULL, *nonCompRPrimer; +char *fPrimer = NULL, *rPrimer = NULL, *nonCompRPrimer = NULL; char *primerKey = NULL; if (lf->original) - primerKey = ((struct psl *)lf->original)->qName; + { + // we can use the qName to extract the primer sequence, which + // may be different from the primers the user pasted in! + struct psl *psl = (struct psl *)lf->original; + fPrimer = cloneString(psl->qName); + char *under = strchr(fPrimer, '_'); + if (under) + *under = 0; + else + { + errAbort("Badly formatted qName ('%s', missing '_' character. " + "Please send an email to genome-www@soe.ucsc.edu with the assembly, " + "forward and reverse primers, and other PCR settings.", psl->qName); + } + rPrimer = under + 1; + } +else pcrResultGetPrimers(primerFileName, &fPrimer, &rPrimer, primerKey); if ((fPrimer == NULL) || (rPrimer == NULL)) return NULL; int fPrimerSize = strlen(fPrimer); int rPrimerSize = strlen(rPrimer); // we need to reverse complement the sequence for the display, but we // don't want to when we do the lookup in the psl file nonCompRPrimer = cloneString(rPrimer); reverseComplement(rPrimer, rPrimerSize); if (lf->name && isNotEmpty(lf->name)) { struct psl *tpsl; char *words[3]; int wordCount = chopByChar(cloneString(lf->extra), '|', words, ArraySize(words)); if (wordCount != 3) errAbort("maybeGetPcrResultSeq: expected 3 |-sep'd words but got '%s'", (char *)lf->extra); char *displayName = words[0]; int ampStart = atoi(words[1]), ampEnd = atoi(words[2]); char *realName = pcrResultItemAccName(lf->name, displayName, NULL); /* isPcr results are so sparse that I think the performance impact * of re-reading the psl file in the draw function is negligible. */ pcrResultGetPsl(pslFileName, target, realName, chromName, ampStart, ampEnd, &tpsl, NULL, fPrimer, nonCompRPrimer); /* Use seq+extFile if specified; otherwise just retrieve from seqFile. */ if (isNotEmpty(target->seqTable) && isNotEmpty(target->extFileTable)) { struct sqlConnection *conn = hAllocConn(database); seq = hDnaSeqGet(database, tpsl->tName, target->seqTable, target->extFileTable); hFreeConn(&conn); } else { struct twoBitFile *tbf = twoBitOpen(target->seqFile); seq = twoBitReadSeqFrag(tbf, tpsl->tName, 0, 0); twoBitClose(&tbf); } int start, end; if (tpsl->strand[0] == '-') { reverseComplement(seq->dna, seq->size); start = seq->size - tpsl->tEnd; end = seq->size - tpsl->tStart; } else { start = tpsl->tStart; end = tpsl->tEnd; } CopyArray(fPrimer, (seq->dna + start), fPrimerSize); CopyArray(rPrimer, (seq->dna + end - rPrimerSize), rPrimerSize); if (tpsl->strand[0] == '-') reverseComplement(seq->dna, seq->size); } else { seq = hChromSeq(database, chromName, lf->start, lf->end); if (lf->orientation < 0) reverseComplement(seq->dna, seq->size); CopyArray(fPrimer, seq->dna, fPrimerSize); CopyArray(rPrimer, (seq->dna + seq->size - rPrimerSize), rPrimerSize); } return seq; } #endif /* GBROWSE */ static struct dnaSeq *maybeGetExtFileSeq(char *seqSource, char *name) /* look up sequence name in seq and extFile tables specified in seqSource */ { /* seqSource is: extFile seqTbl extFileTbl */ static struct dyString *buf = NULL; if (buf == NULL) buf = dyStringNew(0); dyStringClear(buf); dyStringAppend(buf, seqSource); char *words[3]; int nwords = chopByWhite(buf->string, words, ArraySize(words)); if ((nwords != ArraySize(words)) || !sameString(words[0], "extFile")) errAbort("invalid %s track setting: %s", BASE_COLOR_USE_SEQUENCE, seqSource); return hDnaSeqGet(database, name, words[1], words[2]); } struct cacheTwoBitRanges *cdsQueryCache = NULL; static struct dnaSeq *fetchCachedTwoBitSeq(char *url, char *seqName, int seqStart, int seqEnd, boolean doRc, int *retSeqOffset) /* fetch a sequence from a 2bit. Caches open two bit files and sequence in * both forward and reverse strand */ { /* Init static url cache */ if (cdsQueryCache == NULL) cdsQueryCache = cacheTwoBitRangesNew(TRUE); return cacheTwoBitRangesMayFetch(cdsQueryCache, url, seqName, seqStart, seqEnd, doRc, retSeqOffset); } static struct dnaSeq *maybeGetSeqUpper(struct linkedFeatures *lf, char *mrnaName, int mrnaStart, int mrnaEnd, char *tableName, struct track *tg, boolean doRc, int *retMrnaOffset) /* Look up the sequence in genbank tables (hGenBankGetMrna also searches * seq if it can't find it in GB tables) or user's blat sequence, * uppercase and return it if we find it, return NULL if we don't find it. */ { boolean doUpper = TRUE; struct dnaSeq *mrnaSeq = NULL; char *name = getItemDataName(tg, mrnaName); if (sameString(tableName,"refGene") || sameString(tableName,"refSeqAli")) mrnaSeq = hGenBankGetMrna(database, name, "refMrna"); else { char *seqSource = trackDbSetting(tg->tdb, BASE_COLOR_USE_SEQUENCE); if (seqSource != NULL) { if (sameString(seqSource, "ss")) mrnaSeq = maybeGetUserSeq(name); #ifndef GBROWSE else if (sameString(seqSource, PCR_RESULT_TRACK_NAME)) mrnaSeq = maybeGetPcrResultSeq(lf); #endif /* GBROWSE */ else if (startsWith("extFile", seqSource)) mrnaSeq = maybeGetExtFileSeq(seqSource, name); else if (endsWith("ExtFile", seqSource)) mrnaSeq = maybeGetExtFileSeq(seqSource, name); else if (sameString("nameIsSequence", seqSource)) { mrnaSeq = newDnaSeq(cloneString(name), strlen(name), name); if (lf->orientation == -1) reverseComplement(mrnaSeq->dna, mrnaSeq->size); } else if (sameString("seq1Seq2", seqSource)) { mrnaSeq = lf->extra; if (lf->orientation == -1) reverseComplement(mrnaSeq->dna, mrnaSeq->size); } else if (sameString("lfExtra", seqSource)) { if (lf->extra == NULL) errAbort("baseColorDrawSetup: sequence for track '%s' not loaded when sequence option is set in trackDb\n", tg->track); mrnaSeq = newDnaSeq(cloneString(lf->extra), strlen(lf->extra), lf->extra); if (lf->orientation == -1) reverseComplement(mrnaSeq->dna, mrnaSeq->size); } else if (sameString("lrg", seqSource)) { struct lrg *lrg = lf->original; mrnaSeq = lrgReconstructSequence(lrg, database); } else if (sameString("2bit", seqSource)) { char *url = trackDbSetting(tg->tdb, "otherTwoBitUrl"); if (url == NULL) errAbort("missing otherTwoBitUrl in baseColorUseSequence 2bit trackDb setting"); mrnaSeq = fetchCachedTwoBitSeq(url, name, mrnaStart, mrnaEnd, doRc, retMrnaOffset); doRc = FALSE; // Handled it already doUpper = FALSE; // Handled it already } else if (startsWith("table ", seqSource)) { char *table = seqSource; nextWord(&table); mrnaSeq = hGenBankGetMrna(database, name, table); } else if (startsWithWord("db", seqSource)) { char *sourceDb = seqSource; nextWord(&sourceDb); if (isEmpty(sourceDb)) sourceDb = database; mrnaSeq = hChromSeq(sourceDb, name, 0, 0); } else mrnaSeq = hGenBankGetMrna(database, name, NULL); } } if (mrnaSeq != NULL && doUpper) touppers(mrnaSeq->dna); if (mrnaSeq != NULL && doRc) reverseComplement(mrnaSeq->dna, mrnaSeq->size); return mrnaSeq; } static void makeCdsShades(struct hvGfx *hvg, Color *cdsColor) /* setup CDS colors */ { cdsColor[CDS_ERROR] = hvGfxFindColorIx(hvg,0,0,0); cdsColor[CDS_ODD] = hvGfxFindColorIx(hvg,CDS_ODD_R,CDS_ODD_G,CDS_ODD_B); cdsColor[CDS_EVEN] = hvGfxFindColorIx(hvg,CDS_EVEN_R,CDS_EVEN_G,CDS_EVEN_B); cdsColor[CDS_START] = hvGfxFindColorIx(hvg,CDS_START_R,CDS_START_G,CDS_START_B); cdsColor[CDS_STOP] = hvGfxFindColorIx(hvg,CDS_STOP_R,CDS_STOP_G,CDS_STOP_B); cdsColor[CDS_SPLICE] = hvGfxFindColorIx(hvg,CDS_SPLICE_R,CDS_SPLICE_G, CDS_SPLICE_B); cdsColor[CDS_PARTIAL_CODON] = hvGfxFindColorIx(hvg,CDS_PARTIAL_CODON_R, CDS_PARTIAL_CODON_G, CDS_PARTIAL_CODON_B); cdsColor[CDS_QUERY_INSERTION] = hvGfxFindColorIx(hvg,CDS_QUERY_INSERTION_R, CDS_QUERY_INSERTION_G, CDS_QUERY_INSERTION_B); cdsColor[CDS_QUERY_INSERTION_AT_END] = hvGfxFindColorIx(hvg, CDS_QUERY_INSERTION_AT_END_R, CDS_QUERY_INSERTION_AT_END_G, CDS_QUERY_INSERTION_AT_END_B); cdsColor[CDS_POLY_A] = hvGfxFindColorIx(hvg,CDS_POLY_A_R, CDS_POLY_A_G, CDS_POLY_A_B); cdsColor[CDS_ALT_START] = hvGfxFindColorIx(hvg,CDS_ALT_START_R, CDS_ALT_START_G, CDS_ALT_START_B); cdsColor[CDS_SYN_PROT] = hvGfxFindColorIx(hvg,CDS_SYN_PROT_R, CDS_SYN_PROT_G, CDS_SYN_PROT_B); cdsColor[CDS_SYN_BLEND] = hvGfxFindColorIx(hvg,CDS_SYN_BLEND_R, CDS_SYN_BLEND_G, CDS_SYN_BLEND_B); cdsColor[CDS_RIBO_SLIP1] = hvGfxFindColorIx(hvg,CDS_RIBO_SLIP1_R, CDS_RIBO_SLIP1_G, CDS_RIBO_SLIP1_B); cdsColor[CDS_RIBO_SLIP2] = hvGfxFindColorIx(hvg,CDS_RIBO_SLIP2_R, CDS_RIBO_SLIP2_G, CDS_RIBO_SLIP2_B); } static void updatePartialCodon(char *retStr, int exonStart, int exonEnd, boolean posStrand) /* Add bases to the appropriate end of retStr, from exonStart until we have 3 bases or * run into exonEnd*/ { char tmpStr[4]; char tmpDna[4]; char *codonDna = getCachedDna(exonStart, exonEnd); int baseCount = min(3-strlen(retStr), abs(exonEnd-exonStart)); memcpy(tmpDna, codonDna, baseCount); tmpDna[baseCount] = '\0'; if (posStrand) safef(tmpStr, sizeof(tmpStr), "%s%s", retStr, tmpDna); else safef(tmpStr, sizeof(tmpStr), "%s%s", tmpDna, retStr); memcpy(retStr, tmpStr, 3); retStr[3] = '\0'; } struct simpleFeature *baseColorCodonsFromDna(int frame, int chromStart, int chromEnd, struct dnaSeq *seq, bool reverse) /* Create list of codons from a DNA sequence. The DNA sequence passed in must include 3 bases extra at the start and the end to allow for creating partial codons */ { struct simpleFeature *sfList = NULL, *sf = NULL; char codon[4]; int chromPos = chromStart+1; int i; DNA *start; int seqOffset; if (frame < 0 || frame > 2 || seq == NULL) return NULL; /* pick up bases for partial codon */ seqOffset = frame; chromPos = chromPos - 3 + frame; zeroBytes(codon, 4); for (i = 0, start = seq->dna + seqOffset; i < seq->size; i++, chromPos++) { int offset = i % 3; codon[offset] = *start++; if (offset != 2) continue; /* new codon */ AllocVar(sf); sf->start = chromPos - 3; sf->end = sf->start + 3; if (reverse) { sf->start = winEnd - sf->start + winStart - 3; sf->end = sf->start + 3; } // Base offsets mod 6 for alternating colors: 0,1,2 --> first codon, 3,4,5 --> second codon. bool codonFirstColor = (sf->start % 6 < 3); sf->grayIx = codonToGrayIx(codon, codonFirstColor, NULL, FALSE, TRUE); zeroBytes(codon, 4); slAddHead(&sfList, sf); } slReverse(&sfList); return sfList; } struct simpleFeature *baseColorCodonsFromGenePred(struct linkedFeatures *lf, struct genePred *gp, boolean colorStopStart, boolean codonNumbering) /* Given an lf and the genePred from which the lf was constructed, * return a list of simpleFeature elements, one per codon (or partial * codon if the codon falls on a gap boundary. */ { unsigned *starts = gp->exonStarts; unsigned *ends = gp->exonEnds; int blockCount = gp->exonCount; unsigned cdsStart = gp->cdsStart; unsigned cdsEnd = gp->cdsEnd; int *exonFrames = gp->exonFrames; boolean useExonFrames = (gp->optFields >= genePredExonFramesFld); int frame = 0; int currentStart = 0, currentEnd = 0; char partialCodonSeq[4]; int currentSize; int i; boolean foundStart = FALSE; struct simpleFeature *sfList = NULL, *sf = NULL; int i0, iN, iInc; boolean posStrand; partialCodonSeq[0] = '\0'; if (lf->orientation > 0) //positive strand { i0 = 0; iN = blockCount; iInc = 1; posStrand = TRUE; } else { i0 = blockCount-1; iN=-1; iInc = -1; posStrand = FALSE; } bool altColor = FALSE; unsigned cds5Prime = posStrand ? cdsStart : cdsEnd; int codonIndex = !codonNumbering || !zoomedToCodonNumberLevel ? 0 : 1; for (i=i0; (iInc*i)<(iInc*iN); i=i+iInc) { int exonStart = starts[i]; int exonEnd = ends[i]; if (useExonFrames) { if(exonFrames[i] > 0) frame = 3 - exonFrames[i]; else frame = 0; } if(frame == 0) strcpy(partialCodonSeq,""); // 3' or 5' UTR exon if ((exonStart < cdsStart && exonEnd <= cdsStart) || (exonStart >= cdsEnd && exonEnd > cdsEnd)) { if (exonEnd > winStart && exonStart < winEnd) { AllocVar(sf); sf->start = exonStart; sf->end = exonEnd; slAddHead(&sfList, sf); } continue; } // UTR + coding exon else if (exonEnd > cds5Prime && exonStart < cds5Prime) { int utrStart, utrEnd; if (posStrand) { utrStart = exonStart; utrEnd = cdsStart; currentStart = cdsStart; } else { utrStart = cdsEnd; utrEnd = exonEnd; currentEnd = cdsEnd; } if (utrEnd > winStart && utrStart < winEnd) { AllocVar(sf); sf->start = utrStart; sf->end = utrEnd; slAddHead(&sfList, sf); } } else if(posStrand) currentStart = exonStart; else currentEnd = exonEnd; // If UTR + coding, trim to coding portion: if (exonStart < cdsStart) exonStart = cdsStart; if (exonEnd > cdsEnd) exonEnd = cdsEnd; // break each exon into codons and assign alternating shades. while (TRUE) { int lastEnd = currentEnd; int codonInc = frame; if (frame == 0) { altColor = altColor ? FALSE : TRUE; codonInc = 3; } if(posStrand) currentEnd = currentStart + codonInc; else currentStart = currentEnd - codonInc; //we've gone off the end of the current exon if ((posStrand && currentEnd > exonEnd) || (!posStrand && currentStart < exonStart )) { if (posStrand) { frame = currentEnd - exonEnd; currentEnd = exonEnd; } else { frame = exonStart - currentStart; currentStart = exonStart; } if(frame == 3) { // end of one block and the next block starts less than three bases away frame = 0; if (posStrand && ((iInc*(i + iInc))<(iInc*iN) && (currentEnd + 3 > starts[i+1]))) { int letter = sfList->grayIx; if (sfList->grayIx > 52) letter -= 52; else if (sfList->grayIx > 26) letter -= 52; sfList->grayIx = 78 + letter; } break; } /* accumulate partial codon in case of one base exon, or start a new one. */ updatePartialCodon(partialCodonSeq, currentStart, currentEnd, posStrand); if (currentStart < winEnd && currentEnd > winStart) { // get next 'frame' nt's to see what codon will be (skipping intron sequence) char theRestOfCodon[4]; getNextCodonDna(theRestOfCodon, frame, gp, i, posStrand); /* This code doesn't really work right in all cases of a * one-base blocks. It broke with some TransMap alignments * with indels around the one base. This code is fragile, so * just work around it by truncating the sequence. */ char tempCodonSeq[8]; if(posStrand) safef(tempCodonSeq, sizeof(tempCodonSeq), "%s%s", partialCodonSeq, theRestOfCodon); else safef(tempCodonSeq, sizeof(tempCodonSeq), "%s%s", theRestOfCodon, partialCodonSeq ); tempCodonSeq[4] = '\0'; // no more than 3 bases AllocVar(sf); sf->start = currentStart; sf->end = currentEnd; sf->grayIx = ((posStrand && currentEnd <= cdsEnd) || (!posStrand && currentStart >= cdsStart)) ? codonToGrayIx(tempCodonSeq, altColor, &foundStart, !posStrand, colorStopStart) : GRAYIX_CDS_ERROR; sf->codonIndex = codonIndex; slAddHead(&sfList, sf); } break; } // end if we've gone off the end of the current exon currentSize = currentEnd - currentStart; /*inside a coding block (with 3 bases)*/ if (currentSize == 3) { AllocVar(sf); sf->start = currentStart; sf->end = currentEnd; if ((posStrand && currentEnd <= cdsEnd) || (!posStrand && currentStart >= cdsStart)) { char currentCodon[4]; char *thisDna = getCachedDna(currentStart, currentEnd); memcpy(currentCodon, thisDna, 3); currentCodon[3] = '\0'; sf->grayIx = codonToGrayIx(currentCodon, altColor, &foundStart, !posStrand, colorStopStart); // is this block less than 3 bases away from the previous block (ribo slip) if (posStrand && (lastEnd + 3 > currentEnd)) sf->grayIx = - 'A' + 1 + 52 + baseColorLookupCodon(currentCodon); } else sf->grayIx = GRAYIX_CDS_ERROR; } /*start of a coding block with less than 3 bases*/ else if (currentSize < 3) { updatePartialCodon(partialCodonSeq, currentStart, currentEnd, posStrand); AllocVar(sf); sf->start = currentStart; sf->end = currentEnd; if (strlen(partialCodonSeq) == 3) sf->grayIx = codonToGrayIx(partialCodonSeq, altColor, &foundStart, !posStrand, colorStopStart); else sf->grayIx = GRAYIX_CDS_ERROR; strcpy(partialCodonSeq,"" ); /*update frame based on bases appended*/ frame -= currentSize; } else errAbort("%s: Too much dna (%d - %d = %d)<br>\n", lf->name, currentEnd, currentStart, currentSize); if (codonIndex) sf->codonIndex = codonIndex++; slAddHead(&sfList, sf); if(posStrand) currentStart = currentEnd; else currentEnd = currentStart; } // end loop on codons within exon /* coding + UTR exon */ if (posStrand && (exonEnd < ends[i]) && exonEnd < winEnd && ends[i] > winStart) { AllocVar(sf); sf->start = exonEnd; sf->end = ends[i]; slAddHead(&sfList, sf); } else if (!posStrand && (exonStart > starts[i]) && starts[i] < winEnd && exonStart > winStart) { AllocVar(sf); sf->start = starts[i]; sf->end = exonStart; slAddHead(&sfList, sf); } } // end loop on exons if(posStrand) slReverse(&sfList); return(sfList); } static void getMrnaBases(struct psl *psl, struct dnaSeq *mrnaSeq, int mrnaOffset, int mrnaS, int s, int e, boolean isRc, char retMrnaBases[4], boolean *retQueryInsertion) /* Get mRNA bases for the current mRNA codon triplet. If this is a split * codon, retrieve the adjacent mRNA bases to make a full triplet. */ { int size = e - s; if(size < 3) { if (mrnaSeq != NULL && mrnaS-(3-size) > 0) { int i=0; int idx = -1; int newIdx = -1; unsigned *gaps = NULL; boolean appendAtStart = FALSE; AllocArray(gaps, psl->blockCount+2); for(i=0; i<psl->blockCount; i++) { unsigned tStart = psl->tStarts[i]; unsigned tEnd = tStart + psl->blockSizes[i]; if (psl->strand[1] == '-') reverseUnsignedRange(&tStart, &tEnd, psl->tSize); if (s == tStart) { idx = i; appendAtStart = TRUE; break; } else if (e == tEnd) { idx = i+1; appendAtStart = FALSE; break; } } getHiddenGaps(psl, gaps); if(idx >= 0 && gaps[idx] > 0 && retQueryInsertion != NULL) *retQueryInsertion = TRUE; if (!appendAtStart) { newIdx = mrnaS + size; memcpy(retMrnaBases, &mrnaSeq->dna[mrnaS - mrnaOffset], size); memcpy(retMrnaBases+size, &mrnaSeq->dna[newIdx - mrnaOffset], 3-size); } else { newIdx = mrnaS - (3 - size); memcpy(retMrnaBases, &mrnaSeq->dna[newIdx - mrnaOffset], 3); } } else { memcpy(retMrnaBases, "NNN", 3); } } else memcpy(retMrnaBases, &mrnaSeq->dna[mrnaS - mrnaOffset], 3); retMrnaBases[3] = '\0'; if (isRc) reverseComplement(retMrnaBases, strlen(retMrnaBases)); } static void drawDiffTextBox(struct hvGfx *hvg, int xOff, int y, double scale, int heightPer, MgFont *font, Color color, char *chrom, unsigned s, unsigned e, struct simpleFeature *sf, struct psl *psl, struct dnaSeq *mrnaSeq, int mrnaOffset, struct linkedFeatures *lf, int grayIx, enum baseColorDrawOpt drawOpt, int maxPixels, Color *trackColors, Color ixColor) { int mrnaS = -1, mrnaE = -1; /* Clip s and e to what is actually visible */ if (s < winStart) s = winStart; if (e > winEnd) e = winEnd; if (psl) pslTargetToQueryRangeMap(psl, max(psl->tStart, s), min(psl->tEnd, e), &mrnaS, &mrnaE); else if (sf) mrnaS = sf->qStart; if (mrnaS >= 0) { if (mrnaS < mrnaOffset) { warn("curious mrnaS %d < mrnaOffest %d\n", mrnaS, mrnaOffset); mrnaS = mrnaOffset; } struct dyString *dyMrnaSeq = dyStringNew(256); char mrnaBases[4]; char genomicCodon[2]; char mrnaCodon[2]; boolean queryInsertion = FALSE; boolean isCoding = (drawOpt == baseColorDrawItemCodons || drawOpt == baseColorDrawDiffCodons); mrnaBases[0] = '\0'; if (psl && isCoding) getMrnaBases(psl, mrnaSeq, mrnaOffset, mrnaS, s, e, (lf->orientation == -1), mrnaBases, &queryInsertion); if (queryInsertion && isCoding) color = cdsColor[CDS_QUERY_INSERTION]; dyStringAppendN(dyMrnaSeq, (char*)&mrnaSeq->dna[mrnaS - mrnaOffset], e-s); if (drawOpt == baseColorDrawItemBases) { if (cartUsualBooleanDb(cart, database, COMPLEMENT_BASES_VAR, FALSE)) complement(dyMrnaSeq->string, dyMrnaSeq->stringSize); drawScaledBoxWithText(hvg, s, e, scale, xOff, y, heightPer, color, lf->score, font, dyMrnaSeq->string, zoomedToBaseLevel, winStart, maxPixels, isCoding, TRUE); } else if (drawOpt == baseColorDrawItemCodons) { if (e <= lf->tallEnd) { boolean startColor = FALSE; /* re-set color of this block based on mrna codons rather than * genomic, but keep the odd/even cycle of dark/light shades. */ int mrnaGrayIx = codonToGrayIx(mrnaBases, (grayIx > 26), NULL, FALSE, TRUE); if (color == cdsColor[CDS_START]) startColor = TRUE; color = colorAndCodonFromGrayIx(hvg, mrnaCodon, mrnaGrayIx, ixColor); if (startColor && sameString(mrnaCodon,"M")) color = cdsColor[CDS_START]; drawScaledBoxWithText(hvg, s, e, scale, xOff, y, heightPer, color, lf->score, font, mrnaCodon, zoomedToCodonLevel, winStart, maxPixels, isCoding, TRUE); } else drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); } else if (drawOpt == baseColorDrawDiffBases) { char *diffStr = NULL; char *genoDna = getCachedDna(s, e); diffStr = needMem(sizeof(char) * (e - s + 1)); maskDiffString(diffStr, dyMrnaSeq->string, genoDna, ' ', dyMrnaSeq->stringSize); // fprintf(stderr, "drawOpt =- diffBases. %d %d %d %d\n", (int)strlen(genoDna), (int)strlen(dyMrnaSeq->string), (int)dyMrnaSeq->stringSize, e-s); if (cartUsualBooleanDb(cart, database, COMPLEMENT_BASES_VAR, FALSE)) complement(diffStr, strlen(diffStr)); drawScaledBoxWithText(hvg, s, e, scale, xOff, y, heightPer, color, lf->score, font, diffStr, zoomedToBaseLevel, winStart, maxPixels, isCoding, TRUE); freeMem(diffStr); } else if (drawOpt == baseColorDrawDiffCodons) { if (e <= lf->tallEnd) { /* Color codons red wherever mrna differs from genomic; * keep the odd/even cycle of dark/light shades. */ colorAndCodonFromGrayIx(hvg, genomicCodon, grayIx, ixColor); int mrnaGrayIx = mrnaCodonToGrayIx(mrnaBases, genomicCodon[0], (grayIx > 26)); color = colorAndCodonFromGrayIx(hvg, mrnaCodon, mrnaGrayIx, ixColor); // Look up mrnaCodon again because if mrnaGrayIx is GRAYIX_SYN_PROT, // codon value is lost: safef(mrnaCodon, sizeof(mrnaCodon), "%c", baseColorLookupCodon(mrnaBases)); if (mrnaCodon[0] != genomicCodon[0]) { drawScaledBoxWithText(hvg, s, e, scale, xOff, y, heightPer, color, lf->score, font, mrnaCodon, zoomedToCodonLevel, winStart, maxPixels, isCoding, TRUE); } else drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); } else drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); } else if (drawOpt != baseColorDrawCds) errAbort("Unknown drawOpt: %d<br>\n", drawOpt); dyStringFree(&dyMrnaSeq); } else { if (s < e) { /*show we have an error by coloring entire exon block yellow*/ drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, MG_YELLOW); // FIXME: this shouldn't ever happen, should be an errAbort warn("Bug: drawDiffTextBox: convertCoordUsingPsl failed<br>\n"); } } } void baseColorDrawItem(struct track *tg, struct linkedFeatures *lf, int grayIx, struct hvGfx *hvg, int xOff, int y, double scale, MgFont *font, int s, int e, int heightPer, boolean zoomedToCodonLevel, struct dnaSeq *qSeq, int qOffset, struct simpleFeature *sf, struct psl *psl, enum baseColorDrawOpt drawOpt, int maxPixels, int winStart, Color originalColor) /* Draw codon/base-colored item. */ { char codon[64] = " "; Color color = colorAndCodonFromGrayIx(hvg, codon, grayIx, originalColor); if (sf->codonIndex && ( e - s >= 3)) // don't put exon numbers on split codons because there isn't space. safef(codon, sizeof(codon), "%c %d", codon[0], sf->codonIndex); /* When we are zoomed out far enough so that multiple bases/codons share the * same pixel, we have to draw differences in a separate pass (baseColorOverdrawDiff) * so don't waste time drawing the differences here: */ boolean zoomedOutToPostProcessing = ((drawOpt == baseColorDrawDiffBases && !zoomedToBaseLevel) || (drawOpt == baseColorDrawDiffCodons && !zoomedToCdsColorLevel)); if (drawOpt == baseColorDrawGenomicCodons && (e-s <= 3)) { if (lf->highlightColor) { drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, lf->highlightColor); drawScaledBoxWithText(hvg, s, e, scale, xOff, y+1, heightPer-2, color, lf->score, font, codon, zoomedToCodonLevel, winStart, maxPixels, TRUE, !sf->codonIndex); } else { drawScaledBoxWithText(hvg, s, e, scale, xOff, y, heightPer, color, lf->score, font, codon, zoomedToCodonLevel, winStart, maxPixels, TRUE, !sf->codonIndex); } } else if (qSeq != NULL && (psl != NULL || sf != NULL) && !zoomedOutToPostProcessing && drawOpt != baseColorDrawGenomicCodons && drawOpt != baseColorDrawOff) { if (lf->highlightColor) { drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, lf->highlightColor); drawDiffTextBox(hvg, xOff+1, y+1, scale, heightPer-2, font, color, chromName, s, e, sf, psl, qSeq, qOffset, lf, grayIx, drawOpt, maxPixels, tg->colorShades, originalColor); } else { drawDiffTextBox(hvg, xOff, y, scale, heightPer, font, color, chromName, s, e, sf, psl, qSeq, qOffset, lf, grayIx, drawOpt, maxPixels, tg->colorShades, originalColor); } } else { /* revert to normal coloring */ if (lf->highlightColor) { drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, lf->highlightColor); drawScaledBox(hvg, s, e, scale, xOff+1, y+1, heightPer -2, color); } else { drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); } } } static void drawCdsDiffCodonsOnly(struct track *tg, struct linkedFeatures *lf, struct hvGfx *hvg, int xOff, int y, double scale, int heightPer, struct dnaSeq *qSeq, int qOffset, struct psl *psl, int winStart) /* Draw red boxes only where mRNA codons differ from genomic. This assumes * that lf has been drawn already, we're zoomed out past zoomedToCdsColorLevel, * we're not in dense mode etc. */ { struct simpleFeature *sf = NULL; Color dummyColor = 0; if (lf->codons == NULL) errAbort("drawCdsDiffCodonsOnly: lf->codons is NULL"); for (sf = lf->codons; sf != NULL; sf = sf->next) { int s = sf->start; int e = sf->end; if (s < lf->tallStart) s = lf->tallStart; if (e > lf->tallEnd) e = lf->tallEnd; if (s > winEnd || e < winStart) continue; if (e > s) { int mrnaS = convertCoordUsingPsl( s, psl ); if (mrnaS >= 0) { char mrnaBases[4]; char genomicCodon[2], mrnaCodon; boolean queryInsertion = FALSE; Color color = cdsColor[CDS_STOP]; getMrnaBases(psl, qSeq, qOffset, mrnaS, s, e, (lf->orientation == -1), mrnaBases, &queryInsertion); if (queryInsertion) color = cdsColor[CDS_QUERY_INSERTION]; mrnaCodon = baseColorLookupCodon(mrnaBases); colorAndCodonFromGrayIx(hvg, genomicCodon, sf->grayIx, dummyColor); if (queryInsertion) drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); if (mrnaCodon != genomicCodon[0]) { if (mrnaCodon != genomicCodon[0] && protEquivalent(genomicCodon[0], mrnaCodon)) color = cdsColor[CDS_SYN_PROT]; /* this was a call to drawScaledBoxBlend, but this breaks under * 32-bit color, so for the moment we're going to depend * on the painter's algorithm */ drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, color); } } else { /*show we have an error by coloring entire exon block yellow*/ drawScaledBox(hvg, s, e, scale, xOff, y, heightPer, MG_YELLOW); // FIXME: this shouldn't ever happen, should be an errAbort warn("Bug: drawCdsDiffCodonsOnly: convertCoordUsingPsl failed<br>\n"); } } } } void baseColorOverdrawDiff(struct track *tg, struct linkedFeatures *lf, struct hvGfx *hvg, int xOff, int y, double scale, int heightPer, struct dnaSeq *qSeq, int qOffset, struct psl *psl, int winStart, enum baseColorDrawOpt drawOpt) /* If we're drawing different bases/codons, and zoomed out past base/codon * level, draw 1-pixel wide red lines only where bases/codons differ from * genomic. This tests drawing mode and zoom level but assumes that lf itself * has been drawn already and we're not in dense mode etc. */ { boolean enabled = TRUE; // check max zoom float showDiffBasesMaxZoom = trackDbFloatSettingOrDefault(tg->tdb, "showDiffBasesMaxZoom", -1.0); if ((showDiffBasesMaxZoom >= 0.0) && ((basesPerPixel > showDiffBasesMaxZoom) || (showDiffBasesMaxZoom == 0.0))) enabled = FALSE; if (drawOpt == baseColorDrawDiffCodons && !zoomedToCdsColorLevel && lf->codons && enabled) { drawCdsDiffCodonsOnly(tg, lf, hvg, xOff, y, scale, heightPer, qSeq, qOffset, psl, winStart); } if (drawOpt == baseColorDrawDiffBases && !zoomedToBaseLevel && enabled) { drawCdsDiffBaseTickmarksOnly(tg, lf, hvg, xOff, y, scale, heightPer, qSeq, qOffset, psl, winStart); } } void baseColorOverdrawQInsert(struct track *tg, struct linkedFeatures *lf, struct hvGfx *hvg, int xOff, int y, double scale, int heightPer, struct dnaSeq *qSeq, int qOffset, struct psl *psl, MgFont *font, int winStart, enum baseColorDrawOpt drawOpt, boolean indelShowQInsert, boolean indelShowPolyA) /* If applicable, draw 1-pixel wide orange lines for query insertions in the * middle of the query, 1-pixel wide blue lines for query insertions at the * end of the query, and 1-pixel wide green (instead of blue) when a query * insertion at the end is a valid poly-A tail. */ { assert(psl); int i; int s; int lastBlk = psl->blockCount - 1; boolean gotPolyAStart=FALSE, gotPolyAEnd=FALSE; if (indelShowPolyA && qSeq) { /* Draw green lines for polyA first, so if the entire transcript is * jammed into one pixel and the other end has a blue line, blue is * what the user sees. */ if (psl->qStarts[0] != 0 && psl->strand[0] == '-') { /* Query is -. We reverse-complemented in baseColorDrawSetup, * so test for polyT head: */ int polyTSize = headPolyTSizeLoose(qSeq->dna, qSeq->size); if (polyTSize > 0 && (polyTSize + 3) >= psl->qStarts[0]) { if (psl->strand[1] == '-') s = psl->tSize - psl->tStarts[0] - 1; else s = psl->tStarts[0]; drawScaledBox(hvg, s, s+1, scale, xOff, y+1, heightPer-2, cdsColor[CDS_POLY_A]); gotPolyAStart = TRUE; } } if ((psl->qStarts[lastBlk] + psl->blockSizes[lastBlk] != psl->qSize) && psl->strand[0] == '+') { if (psl->strand[1] == '-') { /* Query is + but target is -. We reverse-complemented in * baseColorDrawSetup, so test for polyT head: */ int polyTSize = headPolyTSizeLoose(qSeq->dna, qSeq->size); int rcQStart = (psl->qSize - (psl->qStarts[lastBlk] + psl->blockSizes[lastBlk])); if (polyTSize > 0 && (polyTSize + 3) >= rcQStart) { s = psl->tStart; drawScaledBox(hvg, s, s+1, scale, xOff, y+1, heightPer-2, cdsColor[CDS_POLY_A]); gotPolyAEnd = TRUE; } } else { /* Both are +. We didn't reverse-complement in * baseColorDrawSetup, so test for polyA tail: */ int polyASize = tailPolyASizeLoose(qSeq->dna, qSeq->size); if (polyASize > 0 && ((polyASize + 3) >= (psl->qSize - (psl->qStarts[lastBlk] + psl->blockSizes[lastBlk])))) { s = psl->tStarts[lastBlk] + psl->blockSizes[lastBlk]; drawScaledBox(hvg, s, s+1, scale, xOff, y+1, heightPer-2, cdsColor[CDS_POLY_A]); gotPolyAEnd = TRUE; gotPolyAEnd = TRUE; } } } } if (indelShowQInsert) { int qStart = psl->qStarts[0]; if (qStart != 0 && !gotPolyAStart) { /* Insert at beginning of query -- draw vertical blue line * unless it's polyA. */ s = (psl->strand[1] == '-') ? (psl->tSize - psl->tStarts[0] - 1) : psl->tStarts[0]; Color color = cdsColor[CDS_QUERY_INSERTION_AT_END]; drawVertLine(lf, hvg, s, xOff, y, heightPer, scale, color); // drawLeftNumber(lf, hvg, s, xOff, y, heightPer, scale, color, font, qStart); } for (i = 1; i < psl->blockCount; i++) { int qBlkStart = psl->qStarts[i]; int qPrevBlkEnd = (psl->qStarts[i-1] + psl->blockSizes[i-1]); if (qBlkStart > qPrevBlkEnd) { int tBlkStart = psl->tStarts[i]; int tPrevBlkEnd = (psl->tStarts[i-1] + psl->blockSizes[i-1]); /* Note: if tBlkStart < tPrevBlkEnd, then we have overlap on * target, possibly indicating a bug in the aligner. */ if (tBlkStart <= tPrevBlkEnd) { /* Insert in query only -- draw vertical orange line. */ s = (psl->strand[1] == '-') ? (psl->tSize - psl->tStarts[i] - 1) : psl->tStarts[i]; Color color = cdsColor[CDS_QUERY_INSERTION]; drawMidNumber(lf, hvg, s, xOff, y, heightPer, scale, color, font, qBlkStart - qPrevBlkEnd); } } /* Note: if qBlkStart < qPrevBlkEnd, then we have overlap on query, * possibly indicating a bug in the aligner. Most likely a gap * should be drawn in that case (really a target insert) but I don't * think this is the place to do it. Should be caught by * pre-screening table data for block coords that overlap. */ } int missingAtEnd = psl->qSize - (psl->qStarts[lastBlk] + psl->blockSizes[lastBlk]); if (missingAtEnd != 0 && !gotPolyAEnd) { /* Insert at end of query -- draw vertical blue line unless it's * all polyA. */ s = (psl->strand[1] == '-') ? (psl->tSize - (psl->tStarts[lastBlk] + psl->blockSizes[lastBlk])) : (psl->tStarts[lastBlk] + psl->blockSizes[lastBlk]); Color color = cdsColor[CDS_QUERY_INSERTION_AT_END]; // drawRightNumber(lf, hvg, s, xOff, y, heightPer, scale, color, font, missingAtEnd); drawVertLine(lf, hvg, s, xOff, y, heightPer, scale, color); } } } void baseColorInitTrack(struct hvGfx *hvg, struct track *tg) /* Set up base coloring state (e.g. cache genomic sequence) for tg. * This must be called by tg->drawItems if baseColorDrawSetup is used * in tg->drawItemAt. Peeks at tg->drawItems method to determine whether * tg is linkedFeatures or linkedFeaturesSeries (currently the only * two supported track types -- bed, psl etc. are subclasses of these). */ { enum baseColorDrawOpt drawOpt = baseColorGetDrawOpt(tg); boolean indelShowDoubleInsert, indelShowQueryInsert, indelShowPolyA; indelEnabled(cart, (tg ? tg->tdb : NULL), basesPerPixel, &indelShowDoubleInsert, &indelShowQueryInsert, &indelShowPolyA); if (drawOpt <= baseColorDrawOff && !(indelShowQueryInsert || indelShowPolyA)) return; setGc(); if (initedTrack == NULL || differentString(tg->track, initedTrack)) { int overallStart, overallEnd; boolean isSeries = FALSE; if (tg->drawItems == linkedFeaturesSeriesDraw || tg->drawItems == bamLinkedFeaturesSeriesDraw) isSeries = TRUE; else if (!baseColorCanDraw(tg)) errAbort("baseColorInitTrack: track %s has a type not recognized by baseColorCanDraw.", tg->track); getLinkedFeaturesSpan((struct linkedFeatures *)tg->items, &overallStart, &overallEnd, isSeries); if (overallStart < cachedGenoStart || overallEnd > cachedGenoEnd) { // leak mem to save time (don't bother freeing old cached dna) cachedGenoStart = overallStart; cachedGenoEnd = overallEnd; cachedGenoDna = hDnaFromSeq(database, chromName, cachedGenoStart, cachedGenoEnd, dnaUpper); } initedTrack = cloneString(tg->track); } /* allocate colors for coding coloring */ if (!cdsColorsMade) { makeCdsShades(hvg, cdsColor); cdsColorsMade = TRUE; } } static void checkTrackInited(struct track *tg, char *what) /* Die if baseColorInitTrack has not been called (most recently) for this track. */ { setGc(); if (initedTrack == NULL || differentString(tg->track, initedTrack)) errAbort("Error: Track %s should have been baseColorInitTrack'd before %s. " "(tg->drawItems may be unrecognized by baseColorCanDraw)", tg->track, what); } struct psl *linkedFeatureToPsl(struct linkedFeatures *lf, char *qName, char *tName, int tSize) /* Fake up a psl from linked features */ { /* Go through link features counting up components and doing min/maxes */ int tStart = BIGNUM, qStart = BIGNUM, tEnd = 0, qEnd = 0; int blockCount = 0; struct simpleFeature *sf; for (sf = lf->components; sf != NULL; sf = sf->next) { if (tStart > sf->start) tStart = sf->start; if (tEnd < sf->end) tEnd = sf->end; if (qStart > sf->qStart) qStart = sf->qStart; if (qEnd < sf->qEnd) qEnd = sf->qEnd; blockCount += 1; } struct psl *psl; AllocVar(psl); psl->strand[0] = (lf->orientation < 0 ? '-' : '+'); psl->qName = cloneString(qName); psl->qSize = lf->qSize; psl->qStart = qStart; psl->qEnd = qEnd; psl->tName = cloneString(tName); psl->tSize = tSize; psl->tStart = tStart; psl->tEnd = tEnd; /* Set block count and allocate block-by-block arrays */ psl->blockCount = blockCount; unsigned *blockSizes = AllocArray(psl->blockSizes, blockCount); unsigned *qStarts = AllocArray(psl->qStarts, blockCount); unsigned *tStarts = AllocArray(psl->tStarts, blockCount); /* Go through link features filling in blockSizes, qStarts, qEnds */ int i; unsigned totalSize = 0; for (i=0, sf = lf->components; sf != NULL; sf = sf->next, ++i) { int size = sf->end - sf->start; totalSize += size; blockSizes[i] = size; qStarts[i] = sf->qStart; tStarts[i] = sf->start; } psl->match = totalSize; /* May want to redo later */ /* Go through linked features one last time filling in gap info */ struct simpleFeature *prev = lf->components; for (sf = prev->next; sf != NULL; sf = sf->next) { int qInsertSize = sf->qStart - prev->qStart; if (qInsertSize > 0) { psl->qNumInsert += 1; psl->qBaseInsert += qInsertSize; } int tInsertSize = sf->start - prev->end; if (tInsertSize > 0) { psl->tNumInsert += 1; psl->tBaseInsert += tInsertSize; } prev = sf; } return psl; } enum baseColorDrawOpt baseColorDrawSetup(struct hvGfx *hvg, struct track *tg, struct linkedFeatures *lf, struct dnaSeq **retMrnaSeq, int *retMrnaOffset, struct psl **retPsl) /* Returns the CDS coloring option, allocates colors if necessary, and * returns the sequence and psl record for the given item if applicable. * Note: even if base coloring is not enabled, this will return psl and * mrna seq if query insert/polyA coloring is enabled. * baseColorInitTrack must be called before this (in tg->drawItems) -- * this is meant to be called by tg->drawItemAt (i.e. linkedFeaturesDrawAt). */ { enum baseColorDrawOpt drawOpt = baseColorGetDrawOpt(tg); boolean indelShowDoubleInsert, indelShowQueryInsert, indelShowPolyA; indelEnabled(cart, (tg ? tg->tdb : NULL), basesPerPixel, &indelShowDoubleInsert, &indelShowQueryInsert, &indelShowPolyA); if (drawOpt <= baseColorDrawOff && !(indelShowQueryInsert || indelShowPolyA)) return drawOpt; checkTrackInited(tg, "calling baseColorDrawSetup"); /* If we are using item sequence, fetch alignment and sequence: */ struct psl *psl = NULL; struct dnaSeq *mrnaSeq = NULL; int mrnaOffset = 0; int mrnaStart = 0, mrnaEnd = 0; if (indelShowQueryInsert || indelShowPolyA || drawOpt > baseColorDrawOff) { char *type = tg->tdb->type; boolean needPsl = FALSE; char *qName = lf->name; if (sameString("lrg", tg->tdb->track)) { psl = lrgToPsl(lf->original, hChromSize(database, chromName)); needPsl = TRUE; } else if (startsWith("psl", type) || sameString("bigPsl", type) || startsWithWord("bam", type)) { psl = (struct psl *)(lf->original); needPsl = TRUE; } else if (startsWithWord("chain", type) || startsWithWord("bigChain", type)) { qName = cloneFirstWord(lf->name); psl = linkedFeatureToPsl(lf, qName, chromName, hChromSize(database, chromName)); needPsl = TRUE; } boolean doRc = FALSE; if (needPsl) { if (psl == NULL) drawOpt = baseColorDrawOff; else { doRc = (psl->strand[0] == '-' || psl->strand[1] == '-'); pslTargetToQueryRangeMap(psl, max(psl->tStart, winStart), min(psl->tEnd, winEnd), &mrnaStart, &mrnaEnd); } } /* Do we need the sequence for display, if so get it */ if (drawOpt == baseColorDrawItemBases || drawOpt == baseColorDrawDiffBases || drawOpt == baseColorDrawItemCodons || drawOpt == baseColorDrawDiffCodons || indelShowPolyA) { mrnaSeq = maybeGetSeqUpper(lf, qName, mrnaStart, mrnaEnd, tg->table, tg, doRc, &mrnaOffset); if (mrnaSeq == NULL) { drawOpt = baseColorDrawOff; } } } *retPsl = psl; *retMrnaSeq = mrnaSeq; *retMrnaOffset = mrnaOffset; return drawOpt; } void baseColorDrawRulerCodons(struct hvGfx *hvg, struct simpleFeature *sfList, double scale, int xOff, int y, int height, MgFont *font, int winStart, int maxPixels, bool zoomedToText) /* Draw amino acid translation of genomic sequence based on a list of codons. Used for browser ruler in full mode*/ { struct simpleFeature *sf; if (!cdsColorsMade) { makeCdsShades(hvg, cdsColor); cdsColorsMade = TRUE; } for (sf = sfList; sf != NULL; sf = sf->next) { char codon[4]; Color color = colorAndCodonFromGrayIx(hvg, codon, sf->grayIx, MG_GRAY); if (zoomedToText) drawScaledBoxWithText(hvg, sf->start, sf->end, scale, insideX, y, height, color, 1.0, font, codon, TRUE, winStart, maxPixels, TRUE, TRUE); else /* zoomed in just enough to see colored boxes */ drawScaledBox(hvg, sf->start, sf->end, scale, xOff, y, height, color); } } void baseColorSetCdsBounds(struct linkedFeatures *lf, struct psl *psl, struct track *tg) /* set CDS bounds in linked features for a PSL. Used when zoomed out too far * for codon or base coloring, but still want to render CDS bounds */ { struct genbankCds cds; getPslCds(psl, tg, lf, &cds); if (cds.start < cds.end) { struct genbankCds genomeCds = genbankCdsToGenome(&cds, psl); if (genomeCds.start < genomeCds.end) { lf->tallStart = genomeCds.start; lf->tallEnd = genomeCds.end; } } }