896eb991fff50f44f18811a05171e36fa111163d braney Thu Sep 24 17:13:20 2026 -0700 quickLift: leave empty blocks out of a protein alignment before lifting it, refs #38249 The UniProt bigPsl files store block sizes in bases, and pslFromBigPsl divides them by three, so a block shorter than a codon loads with size 0. pslTransMap rejects such an alignment and aborts, which replaced the whole lifted SwissProt track with a BUG message. About one alignment in eight has a block like this. quickLiftPsl now lifts a copy with the empty blocks removed. quickLiftTester has a new case in the shape of Q96ME1-2. diff --git src/hg/lib/quickLift.c src/hg/lib/quickLift.c index c2769612503..5516bfe6e79 100644 --- src/hg/lib/quickLift.c +++ src/hg/lib/quickLift.c @@ -1,1244 +1,1295 @@ /* Copyright (C) 2023 The Regents of the University of California * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "obscure.h" #include "limits.h" #include "float.h" #include "asParse.h" #include "chain.h" #include "binRange.h" #include "basicBed.h" #include "liftOver.h" #include "hash.h" #include "bigBed.h" #include "bbiFile.h" #include "chainNetDbLoad.h" #include "hdb.h" #include "jksql.h" #include "hgConfig.h" #include "quickLift.h" #include "genePredReader.h" #include "bigChain.h" #include "bigLink.h" #include "chromAlias.h" #include "customTrack.h" #include "encode/encodePeak.h" #include "psl.h" #include "chainToPsl.h" #include "pslTransMap.h" #include "maf.h" struct bigBedInterval *quickLiftGetIntervals(char *quickLiftFile, struct bbiFile *bbi, char *chrom, int start, int end, struct hash **pChainHash) /* Return intervals from "other" species that will map to the current window. * These intervals are NOT YET MAPPED to the current assembly. */ { char *linkFileName = bigChainGetLinkFile(quickLiftFile); int maxGapBefore = 0; int maxGapAfter = 0; struct chain *chain, *chainList = chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, start, end, -1); struct lm *lm = lmInit(0); struct bigBedInterval *bbList = NULL, *bb; for(chain = chainList; chain; chain = chain->next) { struct cBlock *cb; cb = chain->blockList; if (cb == NULL) continue; int qStart = cb->qStart; int qEnd = cb->qEnd; // get the range for the links on the "other" species for(; cb; cb = cb->next) { if (cb->qStart < qStart) qStart = cb->qStart; if (cb->qEnd > qEnd) qEnd = cb->qEnd; } // now grab the items , probably we should parameterize the max number of items, but to what? struct bigBedInterval *thisInterval = NULL; if (chain->qStrand == '-') thisInterval = bigBedIntervalQuery(bbi, chain->qName, chain->qSize - qEnd, chain->qSize - qStart, 1000000, lm); else thisInterval = bigBedIntervalQuery(bbi, chain->qName, qStart, qEnd, 1000000, lm); // find how much of the items are beyond the viewport for(bb=thisInterval; bb; bb = bb->next) { if (bb->start < qStart) { int gap = qStart - bb->start; if (gap > maxGapBefore) maxGapBefore = gap; } if (bb->end > qEnd) { int gap = bb->end - qEnd; if (gap > maxGapAfter) maxGapAfter = gap; } } bbList = slCat(thisInterval, bbList); } // We are done with the chains we used to bound the data query; // release them before loading the wider set for the lift map below. // Without this, every quickLifted track leaked the cBlocks of every // chain overlapping the window. chainFreeList(&chainList); // now we need to grab the links outside of our viewport so we can map long items // probably we could reuse the chains from above but for the moment this is easier // For the moment we use the same padding on both sides so we don't have to worry about strand if (maxGapBefore > maxGapAfter) maxGapAfter = maxGapBefore; else maxGapBefore = maxGapAfter; // Cap the padding so a single oversized item doesn't drag in chains // (and their dense cBlock lists) covering many megabases. #define QUICKLIFT_MAX_GAP_PAD 1000000 if (maxGapBefore > QUICKLIFT_MAX_GAP_PAD) maxGapBefore = maxGapAfter = QUICKLIFT_MAX_GAP_PAD; int newStart = start - maxGapBefore * 2; if (newStart < 0) newStart = 0; int newEnd = end + maxGapAfter * 2; chainList = chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, newStart, newEnd, -1); for(chain = chainList; chain; chain = chain->next) { chainSwap(chain); if (*pChainHash == NULL) *pChainHash = newHash(0); liftOverAddChainHash(*pChainHash, chain); } return bbList; } static void make12(struct bed *bed) /* Make a bed12 out of something less than that. */ { bed->blockCount = 1; bed->blockSizes = needMem(sizeof(int)); bed->blockSizes[0] = bed->chromEnd - bed->chromStart; bed->chromStarts = needMem(sizeof(int)); bed->chromStarts[0] = 0; } static int snapToBlock(struct chain *chain, int pos, boolean forward) /* remapRangeList will only place a coordinate that falls inside an aligned block, so a * clipped end has to land on real alignment, not merely inside the chain. Return pos if * it is already aligned, otherwise the nearest aligned coordinate looking forward (for a * start) or backward (for an end). Return -1 if there is no such coordinate. * The asymmetry matches remapRangeList: a start needs b->tStart <= start < b->tEnd, so * b->tStart is legal; an end needs b->tStart < end <= b->tEnd, so b->tEnd is legal. */ { struct cBlock *b, *prev = NULL; for (b = chain->blockList; b != NULL; b = b->next) { if ((b->tStart <= pos) && (pos < b->tEnd)) return pos; if (forward && (b->tStart > pos)) return b->tStart; if (b->tEnd <= pos) prev = b; } return (!forward && (prev != NULL)) ? prev->tEnd : -1; } static boolean clipBedToChains(struct hash *chainHash, struct bed *bed) /* An item can be far bigger than the region we loaded chains for -- ClinVar has copy * number variants spanning most of a chromosome. Its ends then sit where no chain * reaches, remapRangeList can place neither of them, and the whole item is dropped even * though the part on screen maps perfectly well. Pull the ends in to the nearest aligned * base so the visible part can lift. Return TRUE if the item was clipped. */ { struct chain *chain = liftOverChainForRange(chainHash, bed->chrom, bed->chromStart, bed->chromEnd); if (chain == NULL) return FALSE; // nothing covers it, let it fail the way it used to int newStart = snapToBlock(chain, bed->chromStart, TRUE); int newEnd = snapToBlock(chain, bed->chromEnd, FALSE); if ((newStart < 0) || (newEnd < 0) || (newStart >= newEnd)) return FALSE; if ((newStart == bed->chromStart) && (newEnd == bed->chromEnd)) return FALSE; // both ends already sit on alignment, nothing to do /* Trim the blocks to the new range. Blocks are in ascending order, so walk them and * keep the part that survives; a block entirely outside the range is dropped. */ if (bed->blockCount > 0) { int i, keep = 0; for (i = 0; i < bed->blockCount; ++i) { int bStart = bed->chromStart + bed->chromStarts[i]; int bEnd = bStart + bed->blockSizes[i]; if (bStart < newStart) bStart = newStart; if (bEnd > newEnd) bEnd = newEnd; if (bStart >= bEnd) continue; bed->chromStarts[keep] = bStart - newStart; bed->blockSizes[keep] = bEnd - bStart; keep++; } if (keep == 0) return FALSE; bed->blockCount = keep; } bed->chromStart = newStart; bed->chromEnd = newEnd; if (bed->thickStart < newStart) bed->thickStart = newStart; if (bed->thickEnd > newEnd) bed->thickEnd = newEnd; if (bed->thickStart > bed->thickEnd) bed->thickStart = bed->thickEnd; return TRUE; } static struct bed *quickLiftBed(struct bbiFile *bbi, struct hash *chainHash, struct bigBedInterval *bb, boolean clip); struct bed *quickLiftIntervalsToBed(struct bbiFile *bbi, struct hash *chainHash, struct bigBedInterval *bb) /* Using chains stored in chainHash, port a bigBedInterval from another assembly to a bed * on the reference. */ { return quickLiftBed(bbi, chainHash, bb, FALSE); } struct bed *quickLiftIntervalsToBedClip(struct bbiFile *bbi, struct hash *chainHash, struct bigBedInterval *bb) /* Like quickLiftIntervalsToBed, but an item too big for the chains we loaded is pulled in * to what they cover rather than dropped. Callers that need the item's true extent (the * details page) should use quickLiftIntervalsToBed instead. */ { // quickLiftClipToChains=off restores the old behavior, where an item whose ends fall // outside the chains we loaded is dropped instead of being pulled in. boolean clip = cfgOptionBooleanDefault("quickLiftClipToChains", TRUE); return quickLiftBed(bbi, chainHash, bb, clip); } static struct bed *quickLiftBed(struct bbiFile *bbi, struct hash *chainHash, struct bigBedInterval *bb, boolean clip) /* Port a bigBedInterval to a bed on the reference. If clip, an item too big for the * chains we loaded is pulled in to what they cover rather than dropped. */ { char startBuf[16], endBuf[16]; char *bedRow[bbi->fieldCount]; char chromName[256]; static int lastChromId = -1; bbiCachedChromLookup(bbi, bb->chromId, lastChromId, chromName, sizeof(chromName)); bigBedIntervalToRow(bb, chromName, startBuf, endBuf, bedRow, ArraySize(bedRow)); struct bed *bed = bedLoadN(bedRow, bbi->definedFieldCount); char *error; if (bbi->definedFieldCount < 12) make12(bed); if (clip) clipBedToChains(chainHash, bed); if ((error = remapBlockedBed(chainHash, bed, 0.0, 0.1, TRUE, TRUE, NULL, NULL)) == NULL) return bed; //else //printf("bed %s error:%s<BR>", bed->name, error); return NULL; } char *quickLiftGetChainPath(struct cart *cart, char *fromDb, char *toDb) /* Return the path from the quickLiftChain table for given assemblies. */ { if (!quickLiftEnabled(cart)) return 0; struct sqlConnection *conn = hConnectCentral(); char query[2048]; sqlSafef(query, sizeof(query), "select q.path from %s q where q.fromDb='%s' and q.toDb='%s'", quickLiftChainTable(), fromDb, toDb); char *path = sqlQuickString(conn, query); hDisconnectCentral(&conn); return path; } unsigned quickLiftGetChainId(struct cart *cart, char *fromDb, char *toDb) /* Return the id from the quickLiftChain table for given assemblies. */ { if (!quickLiftEnabled(cart)) return 0; unsigned ret = 0; struct sqlConnection *conn = hConnectCentral(); char query[2048]; // this needs to use the hg.conf setting sqlSafef(query, sizeof(query), "select q.id from quickLiftChain q where q.fromDb='%s' and q.toDb='%s'", fromDb, toDb); char *geneId = sqlQuickString(conn, query); hDisconnectCentral(&conn); if (geneId) ret = atoi(geneId); return ret; } #define QUICKLIFT_RANGE_PAD 100000 static struct chain *quickLiftLoadChains(char *quickLiftFile, char *chrom, int start, int end) /* Load the chains from quickLiftFile that overlap a padded window around the * destination range. */ { // A track can name the assembly it came from without naming a chain file, since nothing // stops a hub from setting one of the pair and not the other. With no chains there is // nothing to lift, and every caller copes with an empty answer. if (quickLiftFile == NULL) return NULL; // need to add some padding to these coordinates int padStart = start - QUICKLIFT_RANGE_PAD; if (padStart < 0) padStart = 0; char *linkFileName = bigChainGetLinkFile(quickLiftFile); return chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, padStart, end + QUICKLIFT_RANGE_PAD, -1); } static void quickLiftChainQueryRange(struct chain *chain, int *retQStart, int *retQEnd) /* Return the query-side ("other" species) coordinate range spanned by the * aligned blocks of chain, corrected for query strand. chain->blockList must * not be NULL. */ { struct cBlock *cb = chain->blockList; int qStart = cb->qStart; int qEnd = cb->qEnd; for(; cb; cb = cb->next) { if (cb->qStart < qStart) qStart = cb->qStart; if (cb->qEnd > qEnd) qEnd = cb->qEnd; } // correct for strand if (chain->qStrand == '-') { int saveStart = qStart; qStart = chain->qSize - qEnd; qEnd = chain->qSize - saveStart; } *retQStart = qStart; *retQEnd = qEnd; } static boolean quickLiftChainRangeIn(struct chain *chain, int tStart, int tEnd, int *retQStart, int *retQEnd) /* The query side range matching tStart..tEnd on the target, rather than the whole extent * of the chain's blocks. One block can be enormous: hg19 and hg38 run identical for * 12.8Mb on chr7, so the whole-block answer would ask the other assembly for millions of * bases either side of the window. Within a block the two sides are colinear, so the * part that matters can be worked out exactly. Returns FALSE if no block overlaps. */ { struct cBlock *cb; int qStart = 0, qEnd = 0; boolean any = FALSE; for (cb = chain->blockList; cb != NULL; cb = cb->next) { int s = max(cb->tStart, tStart); int e = min(cb->tEnd, tEnd); if (s >= e) continue; int qLo = cb->qStart + (s - cb->tStart); int qHi = cb->qStart + (e - cb->tStart); if (!any || (qLo < qStart)) qStart = qLo; if (!any || (qHi > qEnd)) qEnd = qHi; any = TRUE; } if (!any) return FALSE; // correct for strand if (chain->qStrand == '-') { int saveStart = qStart; qStart = chain->qSize - qEnd; qEnd = chain->qSize - saveStart; } *retQStart = qStart; *retQEnd = qEnd; return TRUE; } struct quickLiftRange *quickLiftSourceRanges(char *quickLiftFile, char *chrom, int start, int end, struct hash *chainHash) // The ranges in the other assembly that map into chrom:start-end on the reference. The // chains that do the mapping are added to chainHash, which is the form the lift functions // read. Use this when the items cannot be had from a query quickLiftSql knows how to make. { struct chain *chain, *chainList = quickLiftLoadChains(quickLiftFile, chrom, start, end); struct quickLiftRange *rangeList = NULL; for(chain = chainList; chain; chain = chain->next) { if (chain->blockList == NULL) continue; // pad the window the same way quickLiftLoadChains does, so an item that reaches into // the window from just outside it is still found int qStart, qEnd; int padStart = start - QUICKLIFT_RANGE_PAD; if (padStart < 0) padStart = 0; if (quickLiftChainRangeIn(chain, padStart, end + QUICKLIFT_RANGE_PAD, &qStart, &qEnd)) { struct quickLiftRange *range; AllocVar(range); range->chrom = cloneString(chain->qName); range->start = qStart; range->end = qEnd; slAddHead(&rangeList, range); } // the query range was read off the chain as it came, so swap only afterwards chainSwap(chain); liftOverAddChainHash(chainHash, chain); } slReverse(&rangeList); return rangeList; } struct hash *quickLiftChainHash(char *quickLiftFile, char *chrom, int start, int end) // Load the quickLift chains covering chrom:start-end on the reference and return them in a // hash keyed on the other assembly's sequence names, which is the shape the lift functions // want. Use this when the items were fetched some other way, so quickLiftSql was not the // thing that collected the chains. { struct hash *chainHash = newHash(8); quickLiftSourceRanges(quickLiftFile, chrom, start, end, chainHash); return chainHash; } struct slList *quickLiftSql(struct sqlConnection *conn, char *quickLiftFile, char *table, char *chrom, int start, int end, char *query, char *extraWhere, ItemLoader2 loader, int numFields,struct hash *chainHash) // retrieve items for which we have a loader from a SQL database for which we have a set quickLift chains. // Save the chains we used to map the item back to the current reference. { struct chain *chain, *chainList = quickLiftLoadChains(quickLiftFile, chrom, start, end); struct slList *item, *itemList = NULL; int rowOffset = 0; struct sqlResult *sr = NULL; char **row = NULL; for(chain = chainList; chain; chain = chain->next) { if (chain->blockList == NULL) continue; int qStart, qEnd; quickLiftChainQueryRange(chain, &qStart, &qEnd); // now grab the items if (query == NULL) sr = hRangeQuery(conn, table, chain->qName, qStart, qEnd, extraWhere, &rowOffset); else sr = sqlGetResult(conn, query); // numFields is what the loader will read, so it is also the least the row can have. // The native loaders check this; without it a table of the wrong type walks off the // end of the row. if ((numFields > 0) && (sqlCountColumns(sr) < numFields + rowOffset)) errAbort("table %s in %s has %d columns, need at least %d", table, sqlGetDatabase(conn), sqlCountColumns(sr), numFields + rowOffset); while ((row = sqlNextRow(sr)) != NULL) { item = loader(row + rowOffset, numFields); slAddHead(&itemList, item); } // now squirrel the swapped chains we used to use to make the retrieved items back to us chainSwap(chain); liftOverAddChainHash(chainHash, chain); } return itemList; } struct genePred *quickLiftGenePreds(struct sqlConnection *conn, char *quickLiftFile, char *table, char *chrom, int start, int end, char *extraWhere, struct hash *chainHash) // Like quickLiftSql, but load genePreds with a genePredReader so the actual set // of (extended) genePred columns in the table is honored. A fixed 15-column // loader misreads classic knownGene-style tables, whose trailing proteinID and // alignID columns are not extended genePred fields. { struct chain *chain, *chainList = quickLiftLoadChains(quickLiftFile, chrom, start, end); struct genePred *gpList = NULL; for(chain = chainList; chain; chain = chain->next) { if (chain->blockList == NULL) continue; int qStart, qEnd; quickLiftChainQueryRange(chain, &qStart, &qEnd); struct genePredReader *gpr = genePredReaderRangeQuery(conn, table, chain->qName, qStart, qEnd, extraWhere); struct genePred *gp; while ((gp = genePredReaderNext(gpr)) != NULL) slAddHead(&gpList, gp); genePredReaderFree(&gpr); // now squirrel the swapped chains we used to use to map the retrieved items back to us chainSwap(chain); liftOverAddChainHash(chainHash, chain); } return gpList; } struct bed *quickLiftBeds(struct bed *bedList, struct hash *chainHash, boolean blocked) // Map a list of bedd in query coordinates to our current reference { struct bed *liftedBedList = NULL; struct bed *nextBed; struct bed *bed; for(bed = bedList; bed; bed = nextBed) { // remapBlockedBed may want to add new beds after this bed if the region maps to more than one location nextBed = bed->next; bed->next = NULL; char *error; if (!blocked) { error = liftOverRemapRange(chainHash, 0.0, bed->chrom, bed->chromStart, bed->chromEnd, bed->strand[0], 0.001, &bed->chrom, (int *)&bed->chromStart, (int *)&bed->chromEnd, &bed->strand[0]); // probably this should keep track of cases where the input does NOT have thickStart == chromStart bed->thickStart = bed->chromStart; bed->thickEnd = bed->chromEnd; } else error = remapBlockedBed(chainHash, bed, 0.0, 0.1, TRUE, TRUE, NULL, NULL); if (error == NULL) { slAddHead(&liftedBedList, bed); } } return liftedBedList; } static long pslAlignedBases(struct psl *psl) /* Total size of the alignment's blocks, in whatever units the blocks are in. */ { long total = 0; int i; for (i = 0; i < psl->blockCount; i++) total += psl->blockSizes[i]; return total; } static void quickLiftPslCounts(struct psl *psl, struct psl *lifted) /* Put the original match, mismatch, repeat and N counts back on a lifted alignment, * scaled by how much of it survived the lift. pslTransMap recounts them off the blocks, * which reads every lifted alignment as a perfect match: the details page then claims * 100% identity and the browser draws every item at full shade. */ { // A protein alignment comes back from the lift in nucleotide space, so its block sizes, // and therefore its counts, are in different units than the ones we started with. double protMul = (pslIsProtein(psl) && !pslIsProtein(lifted)) ? 3.0 : 1.0; long origBases = pslAlignedBases(psl); long newBases = pslAlignedBases(lifted); if ((origBases <= 0) || (newBases <= 0)) return; double survived = newBases / (origBases * protMul); if (survived > 1.0) survived = 1.0; lifted->match = round(psl->match * protMul * survived); lifted->misMatch = round(psl->misMatch * protMul * survived); lifted->repMatch = round(psl->repMatch * protMul * survived); lifted->nCount = round(psl->nCount * protMul * survived); } static struct psl *mapPslForChain(struct hash **pMapPsls, struct chain *chain) /* The mapping alignment for one chain, made once and kept. The chains in a quickLift * chainHash have the other assembly on the target side, which is what remapBlockedBed * wants. pslTransMap wants it the other way round: query on the other assembly, target * on the reference. * Building this per item costs nothing at gene zoom, where a chain covers a handful of * blocks, and a great deal zoomed out, where the chain covering the window carries * thousands of blocks and an alignment track can have hundreds of thousands of items. */ { char key[32]; if (*pMapPsls == NULL) *pMapPsls = newHash(8); safef(key, sizeof key, "%p", chain); struct psl *mapPsl = hashFindVal(*pMapPsls, key); if (mapPsl == NULL) { mapPsl = chainToPsl(chain); pslSwap(mapPsl, FALSE); hashAdd(*pMapPsls, key, mapPsl); } return mapPsl; } static boolean quickLiftPslBackToProtein(struct psl *lifted) /* pslTransMap puts a protein alignment into nucleotide space to do the mapping and leaves * it there, so the query start, end and size come back three times too large and the base * alignment view refuses the alignment ("size of rna X is 604, has changed since alignment * was performed when it was 1812"). Put the query side back into protein units. Returns * FALSE, leaving the alignment alone, when the lift split a codon so the query side no * longer divides evenly. */ { int i; if ((lifted->qStart % 3) || (lifted->qEnd % 3) || (lifted->qSize % 3) || (lifted->qBaseInsert % 3)) return FALSE; for (i = 0; i < lifted->blockCount; i++) if ((lifted->blockSizes[i] % 3) || (lifted->qStarts[i] % 3)) return FALSE; // A protein psl always has its query on the forward strand, "++" or "+-". pslTransMap can // hand back strand[0] == '-' (it reverse complements the input when the two alignments // disagree about the shared sequence's strand), and "-+" would tell pslShow to reverse // complement the protein as though it were DNA. Turn it over so the minus lands on the // target side, where the protein display expects it. // pslRc makes the target strand explicit as it turns the alignment over, so the // assignment below is the other half of this test, not something to do as well. if (lifted->strand[0] == '-') pslRc(lifted); else { // A protein psl carries the target strand explicitly, and pslTransMap normalized the // target onto the forward strand on the way out. lifted->strand[1] = '+'; lifted->strand[2] = 0; } lifted->qStart /= 3; lifted->qEnd /= 3; lifted->qSize /= 3; lifted->qBaseInsert /= 3; for (i = 0; i < lifted->blockCount; i++) { lifted->blockSizes[i] /= 3; lifted->qStarts[i] /= 3; } return TRUE; } +static struct psl *pslWithoutEmptyBlocks(struct psl *psl) +/* A copy of a protein alignment with its zero-length blocks taken out, or NULL if it has + * none. The copy can be left with no blocks at all. The UniProt bigPsl files store block sizes in bases, and a block shorter than a + * codon comes out of pslFromBigPsl with size 0. Drawing does not mind, but pslTransMap + * checks the alignment after converting it to nucleotides and aborts on the empty block, + * which takes down every item in the track. About one alignment in eight has one. */ +{ +int i, j; + +for (i = 0; i < psl->blockCount; i++) + if (psl->blockSizes[i] == 0) + break; +if (i == psl->blockCount) + return NULL; + +struct psl *copy = pslClone(psl); +for (i = 0, j = 0; i < copy->blockCount; i++) + { + if (copy->blockSizes[i] == 0) + continue; + copy->blockSizes[j] = copy->blockSizes[i]; + copy->qStarts[j] = copy->qStarts[i]; + copy->tStarts[j] = copy->tStarts[i]; + j++; + } +copy->blockCount = j; +if (j == 0) + return copy; + +// Removing an end block moves the bounds, and pslIsProtein compares tEnd with the last +// block, so put them back in step. pslRecalcBounds does not know about protein units. +int last = j - 1; +int tStart = copy->tStarts[0]; +int tEnd = copy->tStarts[last] + 3 * copy->blockSizes[last]; +if (copy->strand[1] == '-') + reverseIntRange(&tStart, &tEnd, copy->tSize); +copy->tStart = tStart; +copy->tEnd = tEnd; +copy->qStart = copy->qStarts[0]; +copy->qEnd = copy->qStarts[last] + copy->blockSizes[last]; +pslComputeInsertCounts(copy); +return copy; +} + struct psl *quickLiftPsl(struct hash *chainHash, struct hash **pMapPsls, struct psl *psl) // Map the target side of an alignment from the other assembly onto our current reference. // The query side (the mRNA, EST or protein the alignment is to) is left alone. Returns // NULL if the alignment doesn't map. pMapPsls points at a hash of mapping alignments the // caller keeps across a run of items; point it at a NULL hash to start. { struct chain *chain = liftOverChainForRange(chainHash, psl->tName, psl->tStart, psl->tEnd); if (chain == NULL) return NULL; struct psl *mapPsl = mapPslForChain(pMapPsls, chain); // pslTransMap aborts when the two alignments disagree about the size of the sequence they // share. That means the chain and the track were built against different versions of the // other assembly, so drop the item rather than taking the CGI down with it. if (psl->tSize != mapPsl->qSize) return NULL; -struct psl *lifted = pslTransMap(pslTransMapNoOpts, psl, pslTypeUnspecified, - mapPsl, pslTypeUnspecified); +struct psl *trimmed = pslIsProtein(psl) ? pslWithoutEmptyBlocks(psl) : NULL; +if (trimmed && (trimmed->blockCount == 0)) + { + pslFree(&trimmed); + return NULL; + } +struct psl *lifted = pslTransMap(pslTransMapNoOpts, trimmed ? trimmed : psl, + pslTypeUnspecified, mapPsl, pslTypeUnspecified); +pslFree(&trimmed); if (lifted != NULL) { // before counting, so quickLiftPslCounts sees both sides in the same units if (pslIsProtein(psl)) quickLiftPslBackToProtein(lifted); quickLiftPslCounts(psl, lifted); } return lifted; } static struct chain *chainFromPsl(struct psl *psl) /* The inverse of chainToPsl. Score and id are the caller's to fill in, since an alignment * does not carry them. */ { struct chain *chain; struct cBlock *blockList = NULL, *b; int i; AllocVar(chain); chain->tName = cloneString(psl->tName); chain->tSize = psl->tSize; chain->tStart = psl->tStart; chain->tEnd = psl->tEnd; chain->qName = cloneString(psl->qName); chain->qSize = psl->qSize; chain->qStrand = psl->strand[0]; // chainToPsl turns the chain's query bounds the right way up for a psl, so turn them back if (chain->qStrand == '-') { chain->qStart = psl->qSize - psl->qEnd; chain->qEnd = psl->qSize - psl->qStart; } else { chain->qStart = psl->qStart; chain->qEnd = psl->qEnd; } for (i = psl->blockCount - 1; i >= 0; i--) { AllocVar(b); b->tStart = psl->tStarts[i]; b->tEnd = b->tStart + psl->blockSizes[i]; b->qStart = psl->qStarts[i]; b->qEnd = b->qStart + psl->blockSizes[i]; slAddHead(&blockList, b); } chain->blockList = blockList; return chain; } struct mafAli *quickLiftMafs(struct hash *chainHash, struct mafAli *mafList, char *sourceDb, char *refSrc, int refSrcSize) // Map MAF blocks from the other assembly onto our current reference. // // A MAF block has to be one contiguous run on its first row, and the lift does not keep // the reference contiguous: where the reference assembly has lost bases the columns for // them go away, and where it has gained bases the alignment says nothing about them. So a // block is cut at every chain block boundary. Inside one chain block the two assemblies // run in step, which is what lets the columns be carried over untouched: only the first // row's coordinates change, and mafSubset does the rest of the arithmetic. // // refSrc is the name the browser expects on the reference row, "<db>.<chrom>", with no hub // prefix. Blocks whose reference does not map are dropped. { struct mafAli *outList = NULL; struct mafAli *maf, *nextMaf; for (maf = mafList; maf != NULL; maf = nextMaf) { nextMaf = maf->next; maf->next = NULL; // The first row of a MAF is its reference, and a reference row is always forward. struct mafComp *ref = maf->components; if ((ref == NULL) || (ref->strand != '+') || (ref->size <= 0)) { mafAliFree(&maf); continue; } // the chains are keyed on the sequence name in the other assembly // mafSplitSrcGetChrom writes into what it is given, so it needs a copy, and the copy // has to be allocated: a maf component name comes from a hub and safecpy into a // fixed buffer would abort on a long one rather than truncate. char *srcBuf = cloneString(ref->src); char *srcChrom = mafSplitSrcGetChrom(srcBuf, sourceDb); int refStart = ref->start; int refEnd = refStart + ref->size; struct chain *chain = liftOverChainForRange(chainHash, srcChrom, refStart, refEnd); if (chain == NULL) { freeMem(srcBuf); mafAliFree(&maf); continue; } struct cBlock *cb; for (cb = chain->blockList; cb != NULL; cb = cb->next) { int runStart = max(cb->tStart, refStart); int runEnd = min(cb->tEnd, refEnd); if (runStart >= runEnd) continue; struct mafAli *sub = mafSubset(maf, ref->src, runStart, runEnd); if (sub == NULL) continue; int destStart = cb->qStart + (runStart - cb->tStart); if (chain->qStrand == '-') { // The lift turns the block over, so turn every row over with it. A chain // keeps its query side reverse complemented, so the forward start of the run // comes from the far end of it. mafFlipStrand(sub); destStart = chain->qSize - (cb->qStart + (runEnd - cb->tStart)); } struct mafComp *subRef = sub->components; freeMem(subRef->src); subRef->src = cloneString(refSrc); subRef->srcSize = refSrcSize; subRef->strand = '+'; subRef->start = destStart; slAddHead(&outList, sub); } freeMem(srcBuf); mafAliFree(&maf); } slReverse(&outList); return outList; } boolean quickLiftIsLifted(struct trackDb *tdb) // TRUE when this track's data comes from another assembly and there is enough to lift it. // Both halves have to be there: the chain file that does the lifting and the assembly the // data came from. A hub can set either one on its own, and half the pair is no use. { return (tdb != NULL) && (trackDbSetting(tdb, "quickLiftUrl") != NULL) && (trackDbSetting(tdb, "quickLiftDb") != NULL); } boolean quickLiftIsOwnChainTrack(struct trackDb *tdb) // TRUE when this is the chain track quickLift builds to show the lift itself. That stanza // carries quickLiftUrl and quickLiftDb like any lifted track, but its data is already in // reference coordinates and must not be lifted a second time. The giveaway is that its // bigDataUrl IS the quickLift chain file. { char *quickLiftFile = trackDbSetting(tdb, "quickLiftUrl"); if (quickLiftFile == NULL) return FALSE; if (startsWithNoCase("bigQuickLiftChain", tdb->type)) return TRUE; char *bigDataUrl = trackDbSetting(tdb, "bigDataUrl"); return (bigDataUrl != NULL) && sameString(bigDataUrl, quickLiftFile); } struct chain *quickLiftChain(struct hash *chainHash, struct hash **pMapPsls, struct chain *chain) // Map a chain's target side from the other assembly onto our current reference. A chain is // an alignment between that assembly and some other species, so this composes the two and // leaves a chain between the reference and that species. The query side is left alone. // Returns NULL if the chain doesn't map. The chain handed in is not modified. { // chainToPsl copies the header, and every chain loader leaves the header describing the // whole chain while loading only the blocks that overlap the range asked for. Correct it // for the conversion, then put it back: callers still want the whole-chain header, which // is what the native details page reports. int saveTStart = chain->tStart, saveTEnd = chain->tEnd; int saveQStart = chain->qStart, saveQEnd = chain->qEnd; struct cBlock *b = chain->blockList; if (b == NULL) return NULL; int tStart = b->tStart, tEnd = b->tEnd, qStart = b->qStart, qEnd = b->qEnd; for (; b != NULL; b = b->next) { if (b->tStart < tStart) tStart = b->tStart; if (b->tEnd > tEnd) tEnd = b->tEnd; if (b->qStart < qStart) qStart = b->qStart; if (b->qEnd > qEnd) qEnd = b->qEnd; } chain->tStart = tStart; chain->tEnd = tEnd; chain->qStart = qStart; chain->qEnd = qEnd; struct psl *psl = chainToPsl(chain); chain->tStart = saveTStart; chain->tEnd = saveTEnd; chain->qStart = saveQStart; chain->qEnd = saveQEnd; struct psl *lifted = quickLiftPsl(chainHash, pMapPsls, psl); pslFree(&psl); if (lifted == NULL) return NULL; struct chain *out = chainFromPsl(lifted); out->score = chain->score; out->id = chain->id; pslFree(&lifted); return out; } struct psl *quickLiftPsls(struct hash *chainHash, struct psl *pslList) // Map a list of alignments in the other assembly's coordinates onto our current reference. // Alignments that don't map are dropped. { struct psl *liftedList = NULL; struct psl *psl, *nextPsl; struct hash *mapPsls = NULL; for(psl = pslList; psl; psl = nextPsl) { nextPsl = psl->next; psl->next = NULL; struct psl *lifted = quickLiftPsl(chainHash, &mapPsls, psl); if (lifted != NULL) slAddHead(&liftedList, lifted); pslFree(&psl); } slReverse(&liftedList); return liftedList; } struct encodePeak *quickLiftPeaks(struct encodePeak *peakList, struct hash *chainHash) // Map a list of encodePeaks in query coordinates to our current reference. These can't go // through quickLiftBeds: the thickStart and thickEnd it assigns overlay signalValue and // pValue in struct encodePeak. { struct encodePeak *liftedList = NULL; struct encodePeak *nextPeak; struct encodePeak *peak; for(peak = peakList; peak; peak = nextPeak) { nextPeak = peak->next; peak->next = NULL; char *error = liftOverRemapRange(chainHash, 0.0, peak->chrom, peak->chromStart, peak->chromEnd, peak->strand[0], 0.001, &peak->chrom, (int *)&peak->chromStart, (int *)&peak->chromEnd, &peak->strand[0]); if (error == NULL) slAddHead(&liftedList, peak); } return liftedList; } boolean quickLiftEnabled(struct cart *cart) /* Return TRUE if feature is available */ { char *cfgEnabled = cartOrCfgOption(cart, "browser.quickLift"); return cfgEnabled && (sameString(cfgEnabled, "on") || sameString(cfgEnabled, "true")) ; } boolean quickLiftAlignmentsEnabled(struct cart *cart) /* Return TRUE if quickLift is allowed to lift alignment tracks: psl, bigPsl, chain, * bigChain, maf, bigMaf and wigMaf. Off unless hg.conf says * browser.quickLiftAlignments=on, and a cart variable of the same name overrides that so * one machine can show both answers. The hg.conf half is read with a literal * cfgOptionBooleanDefault rather than cartOrCfgOption because harvestHgConf.py only sees * the cfgOption* accessors, which is why browser.quickLift itself is missing from the * hg.conf catalog. */ { char *cartEnabled = cartOptionalString(cart, "browser.quickLiftAlignments"); if (cartEnabled != NULL) return sameString(cartEnabled, "on") || sameString(cartEnabled, "true") || sameString(cartEnabled, "yes"); return cfgOptionBooleanDefault("browser.quickLiftAlignments", FALSE); } static int hrCmp(const void *va, const void *vb) /* Compare to sort based on chromStart. */ { const struct quickLiftRegions *a = *((struct quickLiftRegions **)va); const struct quickLiftRegions *b = *((struct quickLiftRegions **)vb); return a->chromStart - b->chromStart; } struct quickLiftRegions *getMismatches(char *ourDb, char strand, char *chrom, char *liftDb, char *liftChrom, struct bigLink *bl, int querySize, int seqStart, int seqEnd, char * chainId) // Helper function to calculate mismatches in a bigLink block { struct quickLiftRegions *hrList = NULL, *hr; int tStart = bl->chromStart; int tEnd = bl->chromEnd; int width = tEnd - tStart; int qStart = bl->qStart; int qEnd = qStart + width; if (strand == '-') { int saveStart = qStart; qStart = querySize - qEnd; qEnd = querySize - saveStart; } // grab that DNA struct dnaSeq *tSeq = hDnaFromSeq(ourDb, chrom, tStart, tEnd, dnaUpper); struct dnaSeq *qSeq = hDnaFromSeq(liftDb, liftChrom, qStart, qEnd, dnaUpper); if (strand == '-') reverseComplement(qSeq->dna, qSeq->size); // now step through looking for mismatches char *tDna = tSeq->dna; char *qDna = qSeq->dna; unsigned tAddr = tStart; unsigned qAddr = qStart; for(; tAddr < tEnd; tAddr++, qAddr++, tDna++, qDna++) { if (tAddr < seqStart) continue; if (tAddr > seqEnd) break; if (*tDna != *qDna) { AllocVar(hr); slAddHead(&hrList, hr); hr->chrom = cloneString(chrom); hr->oChrom = cloneString(liftChrom); hr->chromStart = tAddr; hr->chromEnd = tAddr + 1; hr->oChromStart = qAddr; hr->oChromEnd = qAddr + 1; hr->bases = tDna; hr->otherBases = qDna; hr->baseCount = 1; hr->otherBaseCount = 1; hr->type = QUICKTYPE_MISMATCH; hr->id = chainId; } } return hrList; } struct quickLiftRegions *fillWithGap(struct bigChain *bc, unsigned previousTEnd, unsigned tStart, unsigned previousQEnd, unsigned qStart) { struct quickLiftRegions *hr; AllocVar(hr); hr->id = bc->name; hr->chrom = cloneString(bc->chrom); hr->oChrom = cloneString(bc->qName); hr->chromStart = previousTEnd; hr->chromEnd = tStart; if (bc->strand[0] == '-') { hr->oChromStart = bc->qSize - qStart; hr->oChromEnd = bc->qSize - previousQEnd; } else { hr->oChromStart = previousQEnd; hr->oChromEnd = qStart; } return hr; } struct quickLiftRegions *quickLiftGetRegions(char *ourDb, char *liftDb, char *quickLiftFile, char *chrom, int seqStart, int seqEnd) /* Figure out the highlight regions and cache them. */ { static struct hash *highLightsHash = NULL; struct quickLiftRegions *hrList = NULL; unsigned lengthLimit = atoi(cfgOptionDefault("quickLift.lengthLimit", "10000")); if (seqEnd - seqStart > lengthLimit) return hrList; if (highLightsHash != NULL) { if ((hrList = (struct quickLiftRegions *)hashFindVal(highLightsHash, quickLiftFile)) != NULL) return hrList; } else { highLightsHash = newHash(0); } struct bbiFile *bbiChain = bigBedFileOpenAlias(quickLiftFile, chromAliasFindAliases); struct lm *lm = lmInit(0); struct bigBedInterval *bbChain, *bbChainList = bigBedIntervalQuery(bbiChain, chrom, seqStart, seqEnd, 0, lm); char *links = bigChainGetLinkFile(quickLiftFile); struct bbiFile *bbiLink = bigBedFileOpenAlias(links, chromAliasFindAliases); struct bigBedInterval *bbLink, *bbLinkList = bigBedIntervalQuery(bbiLink, chrom, seqStart, seqEnd, 0, lm); char *chainRow[1024]; char *linkRow[1024]; char startBuf[16], endBuf[16]; for (bbChain = bbChainList; bbChain != NULL; bbChain = bbChain->next) { bigBedIntervalToRow(bbChain, chrom, startBuf, endBuf, chainRow, ArraySize(chainRow)); struct bigChain *bc = bigChainLoad(chainRow); int previousTEnd = -1; int previousQEnd = -1; for (bbLink = bbLinkList; bbLink != NULL; bbLink = bbLink->next) { bigBedIntervalToRow(bbLink, chrom, startBuf, endBuf, linkRow, ArraySize(linkRow)); struct bigLink *bl = bigLinkLoad(linkRow); if (!sameString(bl->name, bc->name)) continue; int tStart = bl->chromStart; int tEnd = bl->chromEnd; int qStart = bl->qStart; int qEnd = qStart + (tEnd - tStart); struct quickLiftRegions *hr; if ((previousTEnd != -1) && (previousTEnd == tStart)) { hr = fillWithGap(bc, previousTEnd, tStart, previousQEnd, qStart); slAddHead(&hrList, hr); hr->type = QUICKTYPE_DEL; struct dnaSeq *qSeq = NULL; if (bc->strand[0] == '-') { qSeq = hDnaFromSeq(liftDb, bc->qName, bc->qSize - hr->oChromEnd, bc->qSize - hr->oChromStart, dnaUpper); reverseComplement(qSeq->dna, qSeq->size); } else qSeq = hDnaFromSeq(liftDb, bc->qName, hr->oChromStart, hr->oChromEnd, dnaUpper); hr->otherBases = qSeq->dna; hr->otherBaseCount = hr->oChromEnd - hr->oChromStart; } else if ( (previousQEnd != -1) && (previousQEnd == qStart)) { hr = fillWithGap(bc, previousTEnd, tStart, previousQEnd, qStart); slAddHead(&hrList, hr); hr->type = QUICKTYPE_INSERT; struct dnaSeq *tSeq = hDnaFromSeq(ourDb, chrom, hr->chromStart, hr->chromEnd, dnaUpper); hr->bases = tSeq->dna; hr->baseCount = hr->chromEnd - hr->chromStart; } else if ( ((previousQEnd != -1) && (previousQEnd != qStart)) && ((previousTEnd != -1) && (previousTEnd != tStart))) { hr = fillWithGap(bc, previousTEnd, tStart, previousQEnd, qStart); hr->type = QUICKTYPE_DOUBLE; hr->baseCount = hr->chromEnd - hr->chromStart; hr->otherBaseCount = hr->oChromEnd - hr->oChromStart; slAddHead(&hrList, hr); } previousQEnd = qEnd; previousTEnd = tEnd; // now find the mismatches in this block struct quickLiftRegions *mismatches = getMismatches(ourDb, bc->strand[0], chrom, liftDb, bc->qName, bl, bc->qSize, seqStart, seqEnd, bc->name); hrList = slCat(mismatches, hrList); } } slSort(&hrList, hrCmp); hashAdd(highLightsHash, quickLiftFile, hrList); return hrList; } void quickLiftResolveTable(struct trackDb *tdb, char *trackTable, char **retTable, char **retLiftDb) /* Resolve the table name and liftDb for a quickLift track. For custom tracks, * sets *retLiftDb to CUSTOM_TRASH and *retTable to the dbTableName setting; * otherwise sets *retTable to trackTable. Caller should have already set * *retLiftDb to trackDbSetting(tdb, "quickLiftDb"). */ { if (isCustomTrack(trackTable)) { *retLiftDb = CUSTOM_TRASH; *retTable = trackDbSetting(tdb, "dbTableName"); } else *retTable = trackTable; } struct bed *quickLiftSqlLoadBeds(struct trackDb *tdb, char *trackTable, char *liftDb, char *chrom, int start, int end, char *extraWhere, ItemLoader2 loader, int numFields, boolean blocked) /* Load items from another assembly via quickLift SQL, map them back to the reference, * and return the lifted beds. Handles custom track table resolution internally. * Caller provides liftDb from trackDbSetting(tdb, "quickLiftDb"). */ { char *table; quickLiftResolveTable(tdb, trackTable, &table, &liftDb); struct hash *chainHash = newHash(8); struct sqlConnection *conn = hAllocConn(liftDb); char *quickLiftFile = cloneString(trackDbSetting(tdb, "quickLiftUrl")); struct bed *bed = (struct bed *)quickLiftSql(conn, quickLiftFile, table, chrom, start, end, NULL, extraWhere, loader, numFields, chainHash); struct bed *liftedBeds = quickLiftBeds(bed, chainHash, blocked); hFreeConn(&conn); return liftedBeds; } char *quickLiftChainTable() /* Return the name of the quickLiftChain table. */ { static char *quickLiftChainTable = NULL; if (quickLiftChainTable == NULL) quickLiftChainTable = cfgOptionEnvDefault("QUICKLIFTCHAINNAME", quickLiftChainTableConfVariable, defaultQuickLiftChainTableName); return quickLiftChainTable; } boolean quickLiftLiftPos(char *sourceDb, char *destDb, char *chrom, int start, int end, char **retChrom, int *retStart, int *retEnd) /* Map a position from source (sourceDb) coords to destination (destDb) coords * using the liftOver chain for sourceDb -> destDb. This is used to remap * hgFind results from quickLifted bigBed tracks (which return hits in the * source assembly's coordinates) back to the destination assembly the user * is viewing. Returns TRUE on success. */ { static struct hash *fileToChainHash = NULL; if (fileToChainHash == NULL) fileToChainHash = newHash(0); char key[1024]; safef(key, sizeof(key), "%s->%s", sourceDb, destDb); struct hash *chainHash = hashFindVal(fileToChainHash, key); if (chainHash == NULL) { char *chainFile = liftOverChainFile(sourceDb, destDb); if (chainFile == NULL) return FALSE; chainHash = newHash(0); // This reads every chain in the file up front. A bigChain-format // liftOver chain indexed on the source (fromDb) side would let us // load just the chains overlapping the hit; worth revisiting if the // upfront cost becomes an issue. readLiftOverMap(chainFile, chainHash); hashAdd(fileToChainHash, key, chainHash); } char strand = '+'; char *error = liftOverRemapRange(chainHash, 0.0, chrom, start, end, strand, 0.001, retChrom, retStart, retEnd, &strand); return (error == NULL); }