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)