3c28d2c5a4076d66b330f2c089b2485657fbe9f9 angie Mon Jan 30 18:10:51 2023 -0800 Make sure the loop on aliasList stops when an item is found so that we have the right value of chromName when we need to search for paired items in an expanded range. refs #30566 diff --git src/hg/hgc/bamClick.c src/hg/hgc/bamClick.c index 5696c2c..8c058e4 100644 --- src/hg/hgc/bamClick.c +++ src/hg/hgc/bamClick.c @@ -1,260 +1,260 @@ /* bamClick - handler for alignments in BAM format (produced by MAQ, * BWA and some other short-read alignment tools). */ /* 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 "hdb.h" #include "hgBam.h" #include "hgc.h" #include "knetUdc.h" #include "udc.h" #include "chromAlias.h" #include "hgBam.h" #include "hgConfig.h" -static boolean foundIt = FALSE; // used in chromAlias search - struct bamTrackData { int itemStart; char *itemName; struct hash *pairHash; + boolean foundIt; }; /* Maybe make this an option someday -- for now, I find it too confusing to deal with * CIGAR that is anchored to positive strand while showing rc'd sequence. I think * to do it right, we would need to reverse the CIGAR string for display. */ static boolean useStrand = FALSE; static boolean skipQualityScore = FALSE; static void singleBamDetails(const bam1_t *bam) /* Print out the properties of this alignment. */ { const bam1_core_t *core = &bam->core; char *itemName = bam1_qname(bam); int tLength = bamGetTargetLength(bam); int tStart = core->pos, tEnd = tStart+tLength; boolean isRc = useStrand && bamIsRc(bam); printPosOnChrom(seqName, tStart, tEnd, NULL, FALSE, itemName); if (!skipQualityScore) printf("Alignment Quality: %d
\n", core->qual); printf("CIGAR string: %s (", bamGetCigar(bam)); bamShowCigarEnglish(bam); printf(")
\n"); printf("Tags:"); bamShowTags(bam); puts("
"); printf("Flags: 0x%02x:
\n   ", core->flag); bamShowFlagsEnglish(bam); puts("
"); if (bamIsRc(bam)) printf("Note: although the read was mapped to the reverse strand of the genome, " "the sequence and CIGAR in BAM are relative to the forward strand.
\n"); puts("
"); struct dnaSeq *genoSeq = hChromSeq(database, seqName, tStart, tEnd); char *qSeq = bamGetQuerySequence(bam, FALSE); if (isNotEmpty(qSeq) && !sameString(qSeq, "*")) { char *qSeq = NULL; struct ffAli *ffa = bamToFfAli(bam, genoSeq, tStart, useStrand, &qSeq); printf("Alignment of %s to %s:%d-%d%s:
\n", itemName, seqName, tStart+1, tEnd, (isRc ? " (reverse complemented)" : "")); ffShowSideBySide(stdout, ffa, qSeq, 0, genoSeq->dna, tStart, tLength, 0, tLength, 8, isRc, FALSE); } if (!skipQualityScore && core->l_qseq > 0) { printf("Sequence quality scores:
\n\n"); UBYTE *quals = bamGetQueryQuals(bam, useStrand); int i; for (i = 0; i < core->l_qseq; i++) { if (i > 0 && (i % 24) == 0) printf("\n"); printf("", qSeq[i], quals[i]); } printf("
%c
%d
\n"); } } static void showOverlap(const bam1_t *leftBam, const bam1_t *rightBam) /* If the two reads overlap, show how. */ { const bam1_core_t *leftCore = &(leftBam->core), *rightCore = &(rightBam->core); int leftStart = leftCore->pos, rightStart = rightCore->pos; int leftLen = bamGetTargetLength(leftBam), rightLen = bamGetTargetLength(rightBam); char *leftSeq = bamGetQuerySequence(leftBam, useStrand); char *rightSeq = bamGetQuerySequence(rightBam, useStrand); if (useStrand && bamIsRc(leftBam)) reverseComplement(leftSeq, strlen(leftSeq)); if (useStrand && bamIsRc(rightBam)) reverseComplement(rightSeq, strlen(rightSeq)); if ((rightStart > leftStart && leftStart + leftLen > rightStart) || (leftStart > rightStart && rightStart+rightLen > leftStart)) { int leftClipLow, rightClipLow; bamGetSoftClipping(leftBam, &leftClipLow, NULL, NULL); bamGetSoftClipping(rightBam, &rightClipLow, NULL, NULL); leftStart -= leftClipLow; rightStart -= rightClipLow; printf("Note: End read alignments overlap:
\n
");
     int i = leftStart - rightStart;
     while (i-- > 0)
 	putc(' ', stdout);
     puts(leftSeq);
     i = rightStart - leftStart;
     while (i-- > 0)
 	putc(' ', stdout);
     puts(rightSeq);
     puts("
"); } } static void bamPairDetails(const bam1_t *leftBam, const bam1_t *rightBam) /* Print out details for paired-end reads. */ { if (leftBam && rightBam) { const bam1_core_t *leftCore = &leftBam->core, *rightCore = &rightBam->core; int leftLength = bamGetTargetLength(leftBam), rightLength = bamGetTargetLength(rightBam); int start = min(leftCore->pos, rightCore->pos); int end = max(leftCore->pos+leftLength, rightCore->pos+rightLength); char *itemName = bam1_qname(leftBam); printf("Paired read name: %s
\n", itemName); printPosOnChrom(seqName, start, end, NULL, FALSE, itemName); puts("

"); } showOverlap(leftBam, rightBam); printf("

Left end read

\n"); singleBamDetails(leftBam); printf("

Right end read

\n"); singleBamDetails(rightBam); printf("
\n"); } static int oneBam(const bam1_t *bam, void *data, bam_hdr_t *header) /* This is called on each record retrieved from a .bam file. */ { const bam1_core_t *core = &bam->core; if (core->flag & BAM_FUNMAP) return 0; struct bamTrackData *btd = (struct bamTrackData *)data; if (sameString(bam1_qname(bam), btd->itemName)) { + btd->foundIt = TRUE; if (btd->pairHash == NULL || (core->flag & BAM_FPAIRED) == 0) { if (core->pos == btd->itemStart) { printf("Read name: %s
\n", btd->itemName); singleBamDetails(bam); } } else { bam1_t *firstBam = (bam1_t *)hashFindVal(btd->pairHash, btd->itemName); if (firstBam == NULL) hashAdd(btd->pairHash, btd->itemName, bamClone(bam)); else { bamPairDetails(firstBam, bam); hashRemove(btd->pairHash, btd->itemName); } } } return 0; } void doBamDetails(struct trackDb *tdb, char *item) /* Show details of an alignment from a BAM file. */ { if (item == NULL) errAbort("doBamDetails: NULL item name"); int start = cartInt(cart, "o"); if (!tdb || !trackDbSetting(tdb, "bamSkipPrintQualScore")) skipQualityScore = FALSE; else skipQualityScore = TRUE; // TODO: libify tdb settings table_pairEndsByName, stripPrefix and pairSearchRange knetUdcInstall(); if (udcCacheTimeout() < 300) udcSetCacheTimeout(300); if (sameString(item, "zoom in")) printf("Zoom in to a region with fewer items to enable 'detail page' links for individual items.
"); char varName[1024]; safef(varName, sizeof(varName), "%s_pairEndsByName", tdb->track); boolean isPaired = cartUsualBoolean(cart, varName, (trackDbSetting(tdb, "pairEndsByName") != NULL)); char position[512]; struct hash *pairHash = isPaired ? hashNew(0) : NULL; -struct bamTrackData btd = {start, item, pairHash}; +struct bamTrackData btd = {start, item, pairHash, FALSE}; char *fileName = hReplaceGbdb(trackDbSetting(tdb, "bigDataUrl")); if (fileName == NULL) { if (isCustomTrack(tdb->table)) { errAbort("bamLoadItemsCore: can't find bigDataUrl for custom track %s", tdb->track); } else { struct sqlConnection *conn = hAllocConnTrack(database, tdb); fileName = hReplaceGbdb(bamFileNameFromTable(conn, tdb->table, seqName)); hFreeConn(&conn); } } char *indexName = hReplaceGbdb(trackDbSetting(tdb, "bigDataIndex")); char *cacheDir = cfgOption("cramRef"); char *refUrl = trackDbSetting(tdb, "refUrl"); struct slName *aliasList = chromAliasFindAliases(seqName); struct slName *nativeName = newSlName(seqName); slAddHead(&aliasList, nativeName); char *chromName = NULL; for (; aliasList; aliasList = aliasList->next) { chromName = aliasList->name; safef(position, sizeof(position), "%s:%d-%d", chromName, winStart, winEnd); bamAndIndexFetchPlus(fileName, indexName, position, oneBam, &btd, NULL, refUrl, cacheDir); - if (foundIt) + if (btd.foundIt) break; } if (isPaired) { char *setting = trackDbSettingOrDefault(tdb, "pairSearchRange", "20000"); int pairSearchRange = atoi(setting); if (pairSearchRange > 0 && hashNumEntries(pairHash) > 0) { // Repeat the search for item in a larger window: struct hash *newPairHash = hashNew(0); btd.pairHash = newPairHash; safef(position, sizeof(position), "%s:%d-%d", chromName, max(0, winStart-pairSearchRange), winEnd+pairSearchRange); - bamFetch(fileName, position, oneBam, &btd, NULL); + bamAndIndexFetchPlus(fileName, indexName, position, oneBam, &btd, NULL, refUrl, cacheDir); } struct hashEl *hel; struct hashCookie cookie = hashFirst(btd.pairHash); while ((hel = hashNext(&cookie)) != NULL) { bam1_t *bam = hel->val; const bam1_core_t *core = &bam->core; if (! (core->flag & BAM_FMUNMAP)) printf("Note: unable to find paired end for %s " "within +-%d bases of viewing window %s
\n", item, pairSearchRange, addCommasToPos(database, cartString(cart, "position"))); else printf("Paired read name: %s
\n", item); singleBamDetails(bam); } } }