af401081fe03d2bd9eb087487939075d29ff8e03 max Fri May 31 08:43:30 2024 -0700 not showing bam query alignment if seq is too long, refs #33292 diff --git src/hg/hgc/bamClick.c src/hg/hgc/bamClick.c index e582f1c..2b9cb39 100644 --- src/hg/hgc/bamClick.c +++ src/hg/hgc/bamClick.c @@ -33,69 +33,82 @@ 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("<B>Alignment Quality: </B>%d<BR>\n", core->qual); if (core->n_cigar > 50) printf("<B>CIGAR string: </B> Cannot show long CIGAR string, more than 50 operations. Contact us if you need to see the full CIGAR string here.<BR>\n"); else { - printf("<B>CIGAR string: </B><tt>%s</tt> (", bamGetCigar(bam)); - bamShowCigarEnglish(bam); - printf(")<BR>\n"); + printf("<B>CIGAR string: </B><tt>%s</tt>", bamGetCigar(bam)); + //bamShowCigarEnglish(bam); + printf("<BR>\n"); } printf("<B>Tags:</B>"); bamShowTags(bam); puts("<BR>"); printf("<B>Flags: </B><tt>0x%02x:</tt><BR>\n ", core->flag); bamShowFlagsEnglish(bam); puts("<BR>"); if (bamIsRc(bam)) printf("<em>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.</em><BR>\n"); puts("<BR>"); struct dnaSeq *genoSeq = hChromSeq(database, seqName, tStart, tEnd); char *qSeq = bamGetQuerySequence(bam, FALSE); +if (core->l_qseq > 5000) + printf("<B>Alignment not shown, very long sequence</B><BR>\n"); +else + { if (isNotEmpty(qSeq) && !sameString(qSeq, "*")) { char *qSeq = NULL; struct ffAli *ffa = bamToFfAli(bam, genoSeq, tStart, useStrand, &qSeq); printf("<B>Alignment of %s to %s:%d-%d%s:</B><BR>\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) { + if (core->l_qseq > 5000) + { + printf("<B>Sequence too long to show quality scores</B><BR>\n"); + } + else + { printf("<B>Sequence quality scores:</B><BR>\n<TT><TABLE><TR>\n"); UBYTE *quals = bamGetQueryQuals(bam, useStrand); int i; for (i = 0; i < core->l_qseq; i++) { if (i > 0 && (i % 24) == 0) printf("</TR>\n<TR>"); printf("<TD>%c<BR>%d</TD>", qSeq[i], quals[i]); } printf("</TR></TABLE></TT>\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))