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");
}
}
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");
UBYTE *quals = bamGetQueryQuals(bam, useStrand);
int i;
for (i = 0; i < core->l_qseq; i++)
{
if (i > 0 && (i % 24) == 0)
printf(" \n");
printf(" %c ", qSeq[i], quals[i]);
}
printf("
%d
\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
"); } showOverlap(leftBam, rightBam); printf("
Left end read\n"); singleBamDetails(leftBam); printf(" | Right end read\n"); singleBamDetails(rightBam); printf(" |