f29b65452a1cca8a5bc62f007100b313e589d036 braney Sun Sep 6 15:11:36 2026 -0700 quickLift: lift a bigNet track, refs #20824 A net is not an ordinary track to lift. It is an alignment of two assemblies, so only its target side moves to the new reference; the query side names a third assembly and is carried across untouched. And the browser draws a net by recursion, so what really has to survive the lift is the tree: a row's level only means anything relative to the row above it. chainNetLoadRangeQuickLift() maps each row's target range with quickLiftIntervalsToBedClip, which is the same code every other quickLift track uses, so a net row lands where a bed of the same span would. The surviving rows then go to the same helpToNet() the unlifted path uses, so there is one tree builder and not two. The row collection either path does is now cnlHelperNew and cnlHelperAddBigNet. validateOneTdb lets bigNet into a quickLift hub. Its chain track is not offered, so a lifted net's details page has no chain to follow and says so rather than looking for a track this assembly does not have. bigNetLoadOne lifts the same way for that page, unclipped, so it reports the item's whole extent. bigNetFromInterval passes -1 as the cached chromId. bbiCachedChromLookup leaves the buffer alone when the id matches the one before it, so a cache that outlives the buffer hands back stale bytes. Measured against the standalone liftOver tool, on hg19 chr22's mouse net lifted to hg38: 36 of 36 source rows in the sampled window land on the same coordinates, counting the one liftOver will not take whole, whose two ends it does place exactly where the browser puts them. Against hg38's own mouse net, computed independently, 1.02% of the drawn pixels differ. diff --git src/hg/lib/chainNetDbLoad.c src/hg/lib/chainNetDbLoad.c index c349b6a89d9..f5d5dba18b4 100644 --- src/hg/lib/chainNetDbLoad.c +++ src/hg/lib/chainNetDbLoad.c @@ -1,435 +1,515 @@ /* chainNetDbLoad - This will load a database representation of * a net into a chainNet representation. Also helps database * representation of chain into chain. */ /* 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 "hash.h" #include "linefile.h" #include "jksql.h" #include "hdb.h" #include "chain.h" #include "chainDb.h" #include "chainLink.h" #include "chainNet.h" #include "netAlign.h" #include "bigNet.h" #include "bigBed.h" +#include "bed.h" +#include "quickLift.h" #include "chainNetDbLoad.h" #include "chromAlias.h" struct cnFill *cnFillFromNetAlign(struct netAlign *na, struct hash *nameHash) /* Convert netAlign to cnFill. Name hash is a place to store * the strings. */ { struct cnFill *fill; AllocVar(fill); fill->tStart = na->tStart; fill->tSize = na->tEnd - na->tStart; fill->qName = hashStoreName(nameHash, na->qName); fill->qStrand = na->strand[0]; fill->qStart = na->qStart; fill->qSize = na->qEnd - na->qStart; fill->chainId = na->chainId; fill->score = na->score; fill->ali = na->ali; fill->qOver = na->qOver; fill->qFar = na->qFar; fill->qDup = na->qDup; if (!sameString(na->type, "gap")) fill->type = hashStoreName(nameHash, na->type); fill->tN = na->tN; fill->qN = na->qN; fill->tR = na->tR; fill->qR = na->qR; fill->tNewR = na->tNewR; fill->qNewR = na->qNewR; fill->tOldR = na->tOldR; fill->qOldR = na->qOldR; fill->tTrf = na->tTrf; fill->qTrf = na->qTrf; return fill; } struct cnlHelper /* Structure to help us make tree from flattened database * representation. */ { struct cnlHelper *next; /* Next in list. */ char *tName; /* Name of target. Not allocated here. */ struct hash *nameHash; /* Place to keep names. */ struct cnFill **levels; /* Array to sort leaves by depth. */ int maxDepth; /* Maximum depth allowed. */ }; struct chainNet *helpToNet(struct cnlHelper **pHelp) /* Make a chainNet from cnlHelper. This will destroy * *pHelp in the process. */ { struct chainNet *net; struct cnlHelper *help = *pHelp; struct cnFill *parentList, *childList, *parent, *child, *nextChild; struct cnFill **levels = help->levels; int depth, maxDepth = 0; /* Note that level 0 is always empty. */ /* Sort everybody by target start. */ for (depth=1; depth < help->maxDepth; ++depth) { if (levels[depth] == NULL && depth != 0) break; maxDepth = depth; slSort(&levels[depth], cnFillCmpTarget); } /* Assign children to parents. */ for (depth=maxDepth; depth >= 2; --depth) { childList = levels[depth]; parentList = levels[depth-1]; child = childList; for (parent = parentList; parent != NULL; parent = parent->next) { while (child != NULL && child->tStart + child->tSize <= parent->tStart + parent->tSize) { nextChild = child->next; slAddHead(&parent->children, child); child = nextChild; } slReverse(&parent->children); } } /* Make up net structure and fill it in. */ AllocVar(net); net->name = cloneString(help->tName); net->fillList = levels[1]; net->nameHash = help->nameHash; /* Cleanup what's left of help and go home. */ freeMem(help->levels); freez(pHelp); return net; } struct chainNet *chainNetLoadResult(struct sqlResult *sr, int rowOffset) /* Given a query result that returns a bunch netAligns, make up * a list of chainNets that has the equivalent information. * Note the net->size field is not filled in. */ { char **row; struct netAlign na; struct cnFill *fill; struct cnlHelper *helpList = NULL, *help, *nextHelp; struct chainNet *netList = NULL, *net; struct hash *helpHash = hashNew(0); /* Convert database rows to cnFills, and store * them in help->levels. */ while ((row = sqlNextRow(sr)) != NULL) { netAlignStaticLoad(row + rowOffset, &na); help = hashFindVal(helpHash, na.tName); if (help == NULL) { AllocVar(help); hashAddSaveName(helpHash, na.tName, help, &help->tName); help->nameHash = hashNew(8); help->maxDepth = 40; slAddHead(&helpList, help); AllocArray(help->levels, help->maxDepth); } fill = cnFillFromNetAlign(&na, help->nameHash); if (na.level >= help->maxDepth) errAbort("can't handle level %d, only go up to %d", na.level, help->maxDepth-1); slAddHead(&help->levels[na.level], fill); } /* Convert helps to nets. */ for (help = helpList; help != NULL; help = nextHelp) { nextHelp = help->next; net = helpToNet(&help); slAddHead(&netList, net); } helpList = NULL; /* Clean up and go home. */ hashFree(&helpHash); return netList; } struct chainNet *chainNetLoadRange(char *database, char *track, char *chrom, int start, int end, char *extraWhere) /* Load parts of a net track that intersect range. */ { int rowOffset; struct sqlConnection *conn; struct sqlResult *sr; struct chainNet *net; conn = sqlConnect(database); sr = hRangeQuery(conn, track, chrom, start, end, extraWhere, &rowOffset); net = chainNetLoadResult(sr, rowOffset); sqlFreeResult(&sr); if (net != NULL) net->size = hChromSize(database, chrom); sqlDisconnect(&conn); return net; } static struct cnFill *cnFillFromBigNet(struct bigNet *bn, struct hash *nameHash) /* Convert a bigNet row to cnFill. Name hash is a place to store * the strings. */ { struct cnFill *fill; AllocVar(fill); fill->tStart = bn->chromStart; fill->tSize = bn->chromEnd - bn->chromStart; fill->qName = hashStoreName(nameHash, bn->name); fill->qStrand = bn->strand[0]; fill->qStart = bn->qStart; fill->qSize = bn->qEnd - bn->qStart; fill->chainId = bn->chainId; fill->score = bn->chainScore; fill->ali = bn->ali; fill->qOver = bn->qOver; fill->qFar = bn->qFar; fill->qDup = bn->qDup; if (!sameString(bn->type, "gap")) fill->type = hashStoreName(nameHash, bn->type); fill->tN = bn->tN; fill->qN = bn->qN; fill->tR = bn->tR; fill->qR = bn->qR; fill->tNewR = bn->tNewR; fill->qNewR = bn->qNewR; fill->tOldR = bn->tOldR; fill->qOldR = bn->qOldR; fill->tTrf = bn->tTrf; fill->qTrf = bn->qTrf; return fill; } +static struct cnlHelper *cnlHelperNew(char *tName) +/* A place to collect net rows by level, before helpToNet rebuilds the tree from + * them. tName is not allocated here; helpToNet clones it. */ +{ +struct cnlHelper *help; +AllocVar(help); +help->tName = tName; +help->nameHash = hashNew(8); +help->maxDepth = 40; +AllocArray(help->levels, help->maxDepth); +return help; +} + +static void cnlHelperAddBigNet(struct cnlHelper *help, char *fileName, struct bigNet *bn, + int tStart, int tEnd) +/* Add one bigNet row to the helper at the target range given. That is the row's own + * range for a plain bigNet, and where the row landed for a quickLifted one. */ +{ +struct cnFill *fill; + +if (bn->level < 1 || bn->level >= help->maxDepth) + errAbort("%s has level %d, net levels run from 1 to %d", + fileName, bn->level, help->maxDepth-1); +fill = cnFillFromBigNet(bn, help->nameHash); +fill->tStart = tStart; +fill->tSize = tEnd - tStart; +slAddHead(&help->levels[bn->level], fill); +} + struct chainNet *chainNetLoadRangeHub(char *fileName, char *chrom, int start, int end) /* Load the parts of a bigNet file that intersect range into a chainNet. * Note the net->size field is not filled in. */ { struct lm *lm = lmInit(0); struct bbiFile *bbi = bigBedFileOpenAlias(fileName, chromAliasFindAliases); struct bigBedInterval *bb, *bbList = bigBedIntervalQuery(bbi, chrom, start, end, 0, lm); char *bedRow[BIGNET_NUM_COLS]; char startBuf[16], endBuf[16]; struct cnlHelper *help = NULL; struct chainNet *net; if (bbList == NULL) { bbiFileClose(&bbi); lmCleanup(&lm); return NULL; } -AllocVar(help); -help->tName = chrom; /* helpToNet clones this. */ -help->nameHash = hashNew(8); -help->maxDepth = 40; -AllocArray(help->levels, help->maxDepth); - +help = cnlHelperNew(chrom); for (bb = bbList; bb != NULL; bb = bb->next) { struct bigNet bn; int fieldCount = bigBedIntervalToRow(bb, chrom, startBuf, endBuf, bedRow, ArraySize(bedRow)); if (fieldCount != BIGNET_NUM_COLS) errAbort("%s has %d fields, bigNet needs %d", fileName, fieldCount, BIGNET_NUM_COLS); bigNetStaticLoad(bedRow, &bn); - if (bn.level < 1 || bn.level >= help->maxDepth) - errAbort("%s has level %d, net levels run from 1 to %d", - fileName, bn.level, help->maxDepth-1); - slAddHead(&help->levels[bn.level], cnFillFromBigNet(&bn, help->nameHash)); + cnlHelperAddBigNet(help, fileName, &bn, bn.chromStart, bn.chromEnd); } net = helpToNet(&help); bbiFileClose(&bbi); lmCleanup(&lm); return net; } +struct bigNet *bigNetFromInterval(struct bbiFile *bbi, struct bigBedInterval *bb, + char *fileName, struct bigNet *bn) +/* Fill in bn from one interval of a bigNet file. The chrom name is the one the file + * carries, which for a quickLifted net is in the source assembly. */ +{ +char *bedRow[BIGNET_NUM_COLS]; +char startBuf[16], endBuf[16]; +char chromName[256]; +int fieldCount; + +/* -1 rather than a remembered chromId: bbiCachedChromLookup leaves the buffer alone + * when the id matches, so a cache that outlives the buffer returns stale bytes. */ +bbiCachedChromLookup(bbi, bb->chromId, -1, chromName, sizeof(chromName)); +fieldCount = bigBedIntervalToRow(bb, chromName, startBuf, endBuf, bedRow, ArraySize(bedRow)); +if (fieldCount != BIGNET_NUM_COLS) + errAbort("%s has %d fields, bigNet needs %d", fileName, fieldCount, BIGNET_NUM_COLS); +bigNetStaticLoad(bedRow, bn); +return bn; +} + +struct chainNet *chainNetLoadRangeQuickLift(char *quickLiftFile, char *fileName, + char *chrom, int start, int end) +/* Load the part of a bigNet file that quickLifts into chrom:start-end, and build a + * chainNet in the destination assembly's coordinates. Only the target side of the net + * moves; the query side describes a third assembly and is carried across untouched. + * Note the net->size field is not filled in. */ +{ +struct bbiFile *bbi = bigBedFileOpenAlias(fileName, chromAliasFindAliases); +struct hash *chainHash = NULL; +struct bigBedInterval *bb, *bbList = quickLiftGetIntervals(quickLiftFile, bbi, chrom, + start, end, &chainHash); +struct cnlHelper *help = NULL; +struct chainNet *net; + +for (bb = bbList; bb != NULL; bb = bb->next) + { + /* Lift through the same code every other quickLift track uses, so a net row lands + * where a bed of the same span would. Only the target range comes from the lifted + * bed; the rest of the row is read from the interval it came from. */ + struct bed *bed = quickLiftIntervalsToBedClip(bbi, chainHash, bb); + struct bigNet bn; + + if ((bed == NULL) || !sameString(bed->chrom, chrom)) + continue; + bigNetFromInterval(bbi, bb, fileName, &bn); + if (help == NULL) + help = cnlHelperNew(chrom); + cnlHelperAddBigNet(help, fileName, &bn, bed->chromStart, bed->chromEnd); + } + +bbiFileClose(&bbi); +if (help == NULL) + return NULL; +net = helpToNet(&help); +return net; +} + struct chainNet *chainNetLoadChrom(char *database, char *track, char *chrom, char *extraWhere) /* Load net on whole chromosome. */ { int rowOffset; struct sqlConnection *conn; struct sqlResult *sr; struct chainNet *net; conn = sqlConnect(database); sr = hChromQuery(conn, track, chrom, extraWhere, &rowOffset); net = chainNetLoadResult(sr, rowOffset); sqlFreeResult(&sr); net->size = hChromSize(database, chrom); sqlDisconnect(&conn); return net; } static void chainLinkAddResult(struct sqlResult *sr, int rowOffset, struct chain *chain) /* Add links that we get from sr to chain. */ { char **row; struct cBlock *list = NULL, *b; struct chainLink link; while ((row = sqlNextRow(sr)) != NULL) { chainLinkStaticLoad(row+rowOffset, &link); AllocVar(b); b->tStart = link.tStart; b->tEnd = link.tEnd; b->qStart = link.qStart; b->qEnd = link.qStart + (link.tEnd - link.tStart); slAddHead(&list, b); } slSort(&list, cBlockCmpTarget); chain->blockList = list; } struct chain *chainLoadIdRangeHub(char *db, char *fileName, char *linkFileName, char *chrom, int start, int end, int id) /* Load parts of chain of given ID from bigChain file. Note the chain header * including score, tStart, tEnd, will still reflect the whole chain, * not just the part in range. However only the blocks of the chain * overlapping the range will be loaded. */ { struct lm *lm = lmInit(0); struct bbiFile *bbi = bigBedFileOpenAlias(fileName, chromAliasFindAliases); struct bigBedInterval *bb, *bbList = bigBedIntervalQuery(bbi, chrom, start, end, 0, lm); struct bbiFile *linkBbi = bigBedFileOpenAlias(linkFileName, chromAliasFindAliases); struct bigBedInterval *linkBb, *linkBbList = bigBedIntervalQuery(linkBbi, chrom, start, end, 0, lm); char *bedRow[12]; char startBuf[16], endBuf[16]; struct chain *chainList = NULL; struct chain *chain; for (bb = bbList; bb != NULL; bb = bb->next) { bigBedIntervalToRow(bb, chrom, startBuf, endBuf, bedRow, ArraySize(bedRow)); unsigned chainId = sqlUnsigned(bedRow[3]); if ((chainId == id) || (id == -1)) { AllocVar(chain); chain->tName = cloneString(chrom); chain->tSize = sqlUnsigned(bedRow[6]); chain->tStart = sqlUnsigned(bedRow[1]); chain->tEnd = sqlUnsigned(bedRow[2]); chain->qName = cloneString(bedRow[7]); chain->qSize = sqlUnsigned(bedRow[8]); chain->qStrand = *bedRow[5]; chain->qStart = sqlUnsigned(bedRow[9]); chain->qEnd = sqlUnsigned(bedRow[10]); chain->score = atof(bedRow[11]); chain->id = chainId; for (linkBb = linkBbList; linkBb != NULL; linkBb = linkBb->next) { bigBedIntervalToRow(linkBb, chrom, startBuf, endBuf, bedRow, ArraySize(bedRow)); unsigned chainId = sqlUnsigned(bedRow[3]); if (chainId == chain->id) { struct cBlock *cBlock; AllocVar(cBlock); slAddHead(&chain->blockList, cBlock); cBlock->tStart = sqlUnsigned(bedRow[1]); cBlock->tEnd = sqlUnsigned(bedRow[2]); unsigned size = cBlock->tEnd - cBlock->tStart; cBlock->qStart = sqlUnsigned(bedRow[4]); cBlock->qEnd = cBlock->qStart + size; } } slReverse(&chain->blockList); slAddHead(&chainList, chain); if (id > 0) // bail out if we found the id we were looking for break; } } bbiFileClose(&bbi); bbiFileClose(&linkBbi); lmCleanup(&lm); if ((id > 0) && (bb == NULL)) errAbort("chain %d is not in %s", id, fileName); return chainList; } static struct chain *chainLoadIdSome(char *database, char *track, char *chrom, int start, int end, int id, boolean loadAll) /* Load some or all of chain. */ { struct sqlConnection *conn; struct sqlResult *sr; char **row; char table[HDB_MAX_TABLE_STRING]; boolean hasBin; struct chain *chain; char query[256]; struct dyString *dy = dyStringNew(128); /* Load chain header. */ if (!hFindSplitTable(database, chrom, track, table, sizeof table, &hasBin)) errAbort("%s table is not in %s", track, database); conn = sqlConnect(database); sqlSafef(query, sizeof(query), "select * from %s where id = %d", table, id); sr = sqlGetResult(conn, query); if ((row = sqlNextRow(sr)) == NULL) errAbort("chain %d is not in %s", id, table); chain = chainHeadLoad(row+hasBin); sqlFreeResult(&sr); /* Load links. */ if (loadAll) { sqlDyStringPrintf(dy, "select * from %sLink where chainId = %d", table, id); } else { sqlDyStringPrintf(dy, "select * from %sLink where ",table ); hAddBinToQuery(start, end, dy); sqlDyStringPrintf(dy," chainId = %d and tStart < %d and tEnd > %d", id, end, start); } sr = sqlGetResult(conn, dy->string); chainLinkAddResult(sr, hasBin, chain); sqlFreeResult(&sr); sqlDisconnect(&conn); dyStringFree(&dy); return chain; } struct chain *chainLoadIdRange(char *database, char *track, char *chrom, int start, int end, int id) /* Load parts of chain of given ID from database. Note the chain header * including score, tStart, tEnd, will still reflect the whole chain, * not just the part in range. However only the blocks of the chain * overlapping the range will be loaded. */ { return chainLoadIdSome(database, track, chrom, start, end, id, FALSE); } struct chain *chainLoadId(char *database, char *track, char *chrom, int id) /* Load chain of given ID from database. */ { return chainLoadIdSome(database, track, chrom, 0, 0, id, TRUE); }