src/hg/hgc/bamClick.c 1.13

1.13 2009/12/10 15:02:12 angie
Got rid of bogus bamIgnoreStrand() -- make useStrand explicit. Added auto-detect of missing 'chr' in sequence names, so stripPrefix setting is no longer necessary. Better message for samopen failures.
Index: src/hg/hgc/bamClick.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/hgc/bamClick.c,v
retrieving revision 1.12
retrieving revision 1.13
diff -b -B -U 4 -r1.12 -r1.13
--- src/hg/hgc/bamClick.c	26 Nov 2009 00:29:11 -0000	1.12
+++ src/hg/hgc/bamClick.c	10 Dec 2009 15:02:12 -0000	1.13
@@ -18,16 +18,21 @@
     char *itemName;
     struct hash *pairHash;
     };
 
+/* 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 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 = bamIsRc(bam);
+boolean isRc = useStrand && bamIsRc(bam);
 printPosOnChrom(seqName, tStart, tEnd, NULL, FALSE, itemName);
 printf("<B>Alignment Quality: </B>%d<BR>\n", core->qual);
 printf("<B>CIGAR string: </B><tt>%s</tt> (", bamGetCigar(bam));
 bamShowCigarEnglish(bam);
@@ -40,16 +45,16 @@
 puts("<BR><BR>");
 char nibName[HDB_MAX_PATH_STRING];
 hNibForChrom(database, seqName, nibName);
 struct dnaSeq *genoSeq = hFetchSeq(nibName, seqName, tStart, tEnd);
-struct ffAli *ffa = bamToFfAli(bam, genoSeq, tStart);
+struct ffAli *ffa = bamToFfAli(bam, genoSeq, tStart, useStrand);
 char *qSeq = ffa->nStart;
 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);
 printf("<B>Sequence quality scores:</B><BR>\n<TT><TABLE><TR>\n");
-UBYTE *quals = bamGetQueryQuals(bam);
+UBYTE *quals = bamGetQueryQuals(bam, useStrand);
 int i;
 for (i = 0;  i < core->l_qseq;  i++)
     {
     if (i > 0 && (i % 24) == 0)
@@ -64,12 +69,13 @@
 {
 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), *rightSeq = bamGetQuerySequence(rightBam);
-if (bamIsRc(leftBam))
+char *leftSeq = bamGetQuerySequence(leftBam, useStrand);
+char *rightSeq = bamGetQuerySequence(rightBam, useStrand);
+if (useStrand && bamIsRc(leftBam))
     reverseComplement(leftSeq, strlen(leftSeq));
-if (bamIsRc(rightBam))
+if (useStrand && bamIsRc(rightBam))
     reverseComplement(rightSeq, strlen(rightSeq));
 if ((rightStart > leftStart && leftStart + leftLen > rightStart) ||
     (leftStart > rightStart && rightStart+rightLen > leftStart))
     {
@@ -137,16 +143,10 @@
 char varName[1024];
 safef(varName, sizeof(varName), "%s_pairEndsByName", tdb->tableName);
 boolean isPaired = cartUsualBoolean(cart, varName,
 				    (trackDbSetting(tdb, "pairEndsByName") != NULL));
-char *seqNameForBam = seqName;
-char *stripPrefix = trackDbSetting(tdb, "stripPrefix");
-if (stripPrefix && startsWith(stripPrefix, seqName))
-    seqNameForBam = seqName + strlen(stripPrefix);
-char posForBam[512];
-safef(posForBam, sizeof(posForBam), "%s:%d-%d", seqNameForBam, winStart, winEnd);
-
-bamIgnoreStrand();
+char position[512];
+safef(position, sizeof(position), "%s:%d-%d", seqName, winStart, winEnd);
 struct hash *pairHash = isPaired ? hashNew(0) : NULL;
 struct bamTrackData btd = {start, item, pairHash};
 char *fileName;
 if (isCustomTrack(tdb->tableName))
@@ -155,10 +155,10 @@
     if (fileName == NULL)
 	errAbort("doBamDetails: can't find bigDataUrl for custom track %s", tdb->tableName);
     }
 else
-    fileName = bamFileNameFromTable(database, tdb->tableName, seqNameForBam);
-bamFetch(fileName, posForBam, oneBam, &btd);
+    fileName = bamFileNameFromTable(database, tdb->tableName, seqName);
+bamFetch(fileName, position, oneBam, &btd);
 if (isPaired)
     {
     char *setting = trackDbSettingOrDefault(tdb, "pairSearchRange", "20000");
     int pairSearchRange = atoi(setting);
@@ -166,11 +166,11 @@
 	{
 	// Repeat the search for item in a larger window:
 	struct hash *newPairHash = hashNew(0);
 	btd.pairHash = newPairHash;
-	safef(posForBam, sizeof(posForBam), "%s:%d-%d", seqNameForBam,
+	safef(position, sizeof(position), "%s:%d-%d", seqName,
 	      max(0, winStart-pairSearchRange), winEnd+pairSearchRange);
-	bamFetch(fileName, posForBam, oneBam, &btd);
+	bamFetch(fileName, position, oneBam, &btd);
 	}
     struct hashEl *hel;
     struct hashCookie cookie = hashFirst(btd.pairHash);
     while ((hel = hashNext(&cookie)) != NULL)