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)