e70c7f68f61a3477631bd364c31f99c157528fc8 braney Fri Sep 4 12:48:02 2026 -0700 quickLift: lift chains, and let chain and bigChain tracks into the hub A chain is an alignment between the assembly the track came from and some other species, so lifting one composes two alignments and leaves the user with chains between the assembly on screen and that species. quickLiftChain does it by the same route the psl work uses: chainToPsl, then quickLiftPsl, then chainFromPsl, which is the inverse chainToPsl never had. It does not modify the chain handed to it. That matters more than it sounds: every chain loader leaves the header describing the whole chain while loading only the blocks in range, so the header has to be corrected for the conversion, but the callers still want the original. The details page reports it, and the track sorts on it, which is what puts a lifted chain in the same row as the native one. quickLiftSourceRanges hands back the ranges in the other assembly that the window maps to, for callers whose items cannot be had from a query quickLiftSql knows how to make. It works those ranges out by intersecting the window with each chain block rather than taking the whole block, which matters here because one block can be enormous: hg19 and hg38 run identical for 12.8Mb on chr7, and asking for all of it turned a 39 chain window into a 1892 chain one. quickLiftIsOwnChainTrack keeps quickLift's own chain track out of all this. That stanza carries quickLiftUrl and quickLiftDb like any lifted track and is loaded by bigChainLoadItems like any bigChain, but its data is already in reference coordinates. The giveaway is that its bigDataUrl is the quickLift chain file. validateOneTdb now accepts chain and bigChain. netAlign is still refused: a net is a hierarchy of gaps rather than a plain alignment and wants its own thought. refs #38249 diff --git src/hg/lib/quickLift.c src/hg/lib/quickLift.c index 0275143cf1b..57d99647552 100644 --- src/hg/lib/quickLift.c +++ src/hg/lib/quickLift.c @@ -1,918 +1,1104 @@ /* 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" 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. */ { // need to add some padding to these coordinates -int padStart = start - 100000; +int padStart = start - QUICKLIFT_RANGE_PAD; if (padStart < 0) padStart = 0; char *linkFileName = bigChainGetLinkFile(quickLiftFile); -return chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, padStart, end+100000, -1); +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; } -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. +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 hash *chainHash = newHash(8); 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); 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)) return FALSE; for (i = 0; i < lifted->blockCount; i++) if ((lifted->blockSizes[i] % 3) || (lifted->qStarts[i] % 3)) return FALSE; lifted->qStart /= 3; lifted->qEnd /= 3; lifted->qSize /= 3; for (i = 0; i < lifted->blockCount; i++) { lifted->blockSizes[i] /= 3; lifted->qStarts[i] /= 3; } // 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; return TRUE; } 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); 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; +} + +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")) ; } 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); }