bc9f54cc38a4f72f3a53c3e7bedccec41a1870be max Mon Apr 29 06:26:58 2024 -0700 adding space and skipping long cigar strings on hgc page, refs #33292 diff --git src/hg/hgc/bamClick.c src/hg/hgc/bamClick.c index 8c058e4..e582f1c 100644 --- src/hg/hgc/bamClick.c +++ src/hg/hgc/bamClick.c @@ -1,260 +1,265 @@ /* 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" 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); +if (core->n_cigar > 50) + printf("CIGAR string: Cannot show long CIGAR string, more than 50 operations. Contact us if you need to see the full CIGAR string here.
\n"); +else + { 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, 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 (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); 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); } } }