b20d369b166f18509c1fd18682ffc4fc5a5aac09 braney Thu Aug 21 17:38:37 2025 -0700 working on detail page for quickLift Chain. diff --git src/hg/lib/quickLift.c src/hg/lib/quickLift.c index 4389e8ce638..f333560c9db 100644 --- src/hg/lib/quickLift.c +++ src/hg/lib/quickLift.c @@ -1,260 +1,442 @@ /* 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 "bigChain.h" +#include "bigLink.h" +#include "chromAlias.h" + +#define QUICKTYPE_INSERT 0 +#define QUICKTYPE_DEL 1 +#define QUICKTYPE_DOUBLE 2 +#define QUICKTYPE_MISMATCH 3 +#define QUICKTYPE_NOTHING 4 + +char *quickTypeStrings[] = +{ + "insert", + "deletion", + "double", + "mismatch" +}; 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); } // 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; 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; } 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. */ { 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 ((error = remapBlockedBed(chainHash, bed, 0.0, 0.1, TRUE, TRUE, NULL, NULL)) == NULL) return bed; //else //printf("bed %s error:%s
", bed->name, error); return NULL; } unsigned quickLiftGetChain(char *fromDb, char *toDb) /* Return the id from the quickLiftChain table for given assemblies. */ { if (!quickLiftEnabled()) return 0; unsigned ret = 0; struct sqlConnection *conn = hConnectCentral(); char query[2048]; 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; } 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. { // need to add some padding to these coordinates int padStart = start - 100000; if (padStart < 0) padStart = 0; char *linkFileName = bigChainGetLinkFile(quickLiftFile); struct chain *chain, *chainList = chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, padStart, end+100000, -1); struct slList *item, *itemList = NULL; int rowOffset = 0; struct sqlResult *sr = NULL; char **row = NULL; 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; } if (chain->qStrand == '-') qStart = chain->qSize - qStart; // 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 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]); } else error = remapBlockedBed(chainHash, bed, 0.0, 0.1, TRUE, TRUE, NULL, NULL); if (error == NULL) { slAddHead(&liftedBedList, bed); } } return liftedBedList; } boolean quickLiftEnabled() /* Return TRUE if feature is available */ { char *cfgEnabled = cfgOption("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 *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); + // crop the chain block if it's bigger than the window + int tMin, tMax; + int qMin, qMax; + tMin = bl->chromStart; + tMax = bl->chromEnd; + qMin = bl->qStart; + if (seqStart > bl->chromStart) + { + tMin = seqStart; + qMin = qStart + (seqStart - bl->chromStart); + } + if (seqEnd < bl->chromEnd) + { + tMax = seqEnd; + } + qMax = qMin + (tMax - tMin); + + if (bc->strand[0] == '-') + { + qMin = bc->qSize - qMax; + qMax = qMin + (tMax - tMin); + } + + struct quickLiftRegions *hr; + AllocVar(hr); + slAddHead(&hrList, hr); + hr->chromStart = tMin; + hr->chromEnd = tMax; + hr->oChromStart = qMin; + hr->oChromEnd = qMax; + //jhr->type = highlightColors[NOTHING]; + hr->type = QUICKTYPE_NOTHING; + hr->id = bc->name; + + struct dnaSeq *tSeq = hDnaFromSeq(ourDb, chrom, tMin, tMax, dnaUpper); + struct dnaSeq *qSeq = hDnaFromSeq(liftDb, bc->qName, qMin, qMax, dnaUpper); + if (bc->strand[0] == '-') + reverseComplement(qSeq->dna, qSeq->size); + if ((previousTEnd != -1) && (previousTEnd == tStart)) + { + AllocVar(hr); + slAddHead(&hrList, hr); + hr->chromStart = previousTEnd; + hr->chromEnd = tStart; + hr->oChromStart = previousQEnd; + hr->oChromEnd = qStart; + hr->type = QUICKTYPE_DEL; + hr->otherBases = &qSeq->dna[qStart - qMin]; + hr->otherCount = hr->oChromEnd - hr->oChromStart; + hr->id = bc->name; + } + else if ( (previousQEnd != -1) && (previousQEnd == qStart)) + { + AllocVar(hr); + slAddHead(&hrList, hr); + hr->chromStart = previousTEnd; + hr->chromEnd = tStart; + hr->oChromStart = previousQEnd; + hr->oChromEnd = qStart; + hr->type = QUICKTYPE_INSERT; + hr->id = bc->name; + } + else if ( ((previousQEnd != -1) && (previousQEnd != qStart)) + && ((previousTEnd != -1) && (previousTEnd != tStart))) + { + AllocVar(hr); + slAddHead(&hrList, hr); + hr->chromStart = previousTEnd; + hr->chromEnd = tStart; + hr->oChromStart = previousQEnd; + hr->oChromEnd = qStart; + hr->type = QUICKTYPE_DOUBLE; + hr->id = bc->name; + } + + previousQEnd = qEnd; + previousTEnd = tEnd; + + + unsigned tAddr = tMin; + unsigned qAddr = qMin; + int count = 0; + for(; tAddr < tMax; tAddr++, qAddr++, count++) + { + if (tSeq->dna[count] != qSeq->dna[count]) + { + AllocVar(hr); + slAddHead(&hrList, hr); + hr->chromStart = tAddr; + hr->chromEnd = tAddr + 1; + hr->oChromStart = qAddr; + hr->oChromEnd = qAddr + 1; + hr->otherBases = &qSeq->dna[count]; + hr->otherCount = 1; + hr->type = QUICKTYPE_MISMATCH; + hr->id = bc->name; + } + } + } + } + +slSort(&hrList, hrCmp); +hashAdd(highLightsHash, quickLiftFile, hrList); + +return hrList; +}