eab3f6024915a08ebe1cce08017271d4a6b142d5 braney Fri Mar 28 11:27:14 2025 -0700 draw quickLift lines for indels, make quickLift bigBed querying smarter about padding diff --git src/hg/lib/quickLift.c src/hg/lib/quickLift.c index beb2035dff1..e1710061e89 100644 --- src/hg/lib/quickLift.c +++ src/hg/lib/quickLift.c @@ -15,67 +15,91 @@ #include "bigBed.h" #include "bbiFile.h" #include "chainNetDbLoad.h" #include "hdb.h" #include "jksql.h" #include "hgConfig.h" #include "quickLift.h" #include "bigChain.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); -// need to add some padding to these coordinates -int padStart = start - 100000; -if (padStart < 0) - padStart = 0; -struct chain *chain, *chainList = chainLoadIdRangeHub(NULL, quickLiftFile, linkFileName, chrom, padStart, end+100000, -1); +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; +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 + // 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, 10000, lm); + thisInterval = bigBedIntervalQuery(bbi, chain->qName, chain->qSize - qEnd, chain->qSize - qStart, 1000000, lm); else - thisInterval = bigBedIntervalQuery(bbi, chain->qName, qStart, qEnd, 10000, lm); + 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); + } - // for the mapping we're going to use the same chain we queried on to map the items, but we need to swap it +// 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 +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)); @@ -83,31 +107,30 @@ 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)); -//lastChromId=bb->chromId; 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<BR>", bed->name, error); return NULL; }