src/hg/hgTracks/bamTrack.c 1.5

1.5 2009/08/03 22:00:24 angie
Libified more commonly used BAM code from hgc and hgTracks.
Index: src/hg/hgTracks/bamTrack.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/hgTracks/bamTrack.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -b -B -U 4 -r1.4 -r1.5
--- src/hg/hgTracks/bamTrack.c	27 Jul 2009 21:52:09 -0000	1.4
+++ src/hg/hgTracks/bamTrack.c	3 Aug 2009 22:00:24 -0000	1.5
@@ -22,21 +22,20 @@
     struct hash *pairHash;
     };
 
 struct simpleFeature *sfFromNumericCigar(const bam1_t *bam, int *retLength)
-/* Translate BAM's numeric CIGAR encoding into a list of simpleFeatures  */
+/* Translate BAM's numeric CIGAR encoding into a list of simpleFeatures, 
+ * and tally up length on reference sequence while we're at it. */
 {
 const bam1_core_t *core = &bam->core;
 struct simpleFeature *sf, *sfList = NULL;
 int tLength=0, tPos = core->pos, qPos = 0;
 unsigned int *cigar = bam1_cigar(bam);
 int i;
 for (i = 0;  i < core->n_cigar;  i++)
     {
-    // decoding lifted from bam.c bam_format1(), long may it remain stable:
-    int n = cigar[i]>>BAM_CIGAR_SHIFT;
-    int opcode = cigar[i] & BAM_CIGAR_MASK;
-    char op = "MIDNSHP"[opcode];
+    char op;
+    int n = bamUnpackCigarElement(cigar[i], &op);
     switch (op)
 	{
 	case 'M': // match or mismatch (gapless aligned block)
 	    AllocVar(sf);
@@ -59,9 +58,9 @@
 	case 'H': // skipped query bases not stored in record's query sequence ("hard clipping")
 	case 'P': // P="silent deletion from padded reference sequence" -- ignore these.
 	    break;
 	default:
-	    errAbort("Unrecognized CIGAR op index %d -- has bam.c bam_format1() changed?", opcode);
+	    errAbort("sfFromNumericCigar: unrecognized CIGAR op %c -- update me", op);
 	}
     }
 if (retLength != NULL)
     *retLength = tLength;
@@ -71,9 +70,8 @@
 struct linkedFeatures *bamToLf(const bam1_t *bam, void *data)
 /* Translate a BAM record into a linkedFeatures item. */
 {
 const bam1_core_t *core = &bam->core;
-uint8_t *s = bam1_seq(bam);
 struct linkedFeatures *lf;
 AllocVar(lf);
 lf->score = core->qual;
 safef(lf->name, sizeof(lf->name), bam1_qname(bam));
@@ -81,15 +79,9 @@
 int length;
 lf->components = sfFromNumericCigar(bam, &length);
 lf->start = lf->tallStart = core->pos;
 lf->end = lf->tallEnd = core->pos + length;
-char *qSeq = needMem(core->l_qseq + 1);
-int i;
-for (i = 0; i < core->l_qseq; i++)
-    qSeq[i] = bam_nt16_rev_table[bam1_seqi(s, i)];
-if (lf->orientation == -1)
-    reverseComplement(qSeq, core->l_qseq);
-lf->extra = qSeq;
+lf->extra = bamGetQuerySequence(bam);
 return lf;
 }
 
 int addBam(const bam1_t *bam, void *data)