ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GBam.cpp
(Generate patch)
# Line 2 | Line 2
2   #include <ctype.h>
3   #include "kstring.h"
4  
5 < //auxiliary functions for BAM record handling
5 > //auxiliary functions for low level BAM record creation
6   uint8_t* realloc_bdata(bam1_t *b, int size) {
7    if (b->m_data < size) {
8          b->m_data = size;
# Line 35 | Line 35
35     novel=true;
36     bam_header=NULL;
37     b=bam_init1();
38 <   b->core.tid=gseq_tid;
39 <   if (pos<=0) {
38 >   if (pos<=0 || gseq_tid<0) {
39                 b->core.pos=-1; //unmapped
41               //if (gseq_tid<0)
40                 b->core.flag |= BAM_FUNMAP;
41 +               gseq_tid=-1;
42                 }
43            else b->core.pos=pos-1; //BAM is 0-based
44 +   b->core.tid=gseq_tid;
45     b->core.qual=255;
46 +   b->core.mtid=-1;
47 +   b->core.mpos=-1;
48     int l_qseq=strlen(qseq);
49     //this may not be accurate, setting CIGAR is the correct way
50     //b->core.bin = bam_reg2bin(b->core.pos, b->core.pos+l_qseq-1);
# Line 168 | Line 170
170  
171   void GBamRecord::add_aux(const char* str) {
172       //requires: being called AFTER add_quals()
171     static char tag[2];
172     static uint8_t abuf[512];
173     //requires: being called AFTER add_quals()
173       int strl=strlen(str);
174       //int doff = b->core.l_qname + b->core.n_cigar*4 + (b->core.l_qseq+1)/2 + b->core.l_qseq + b->l_aux;
175       //int doff0=doff;
# Line 266 | Line 265
265             else if (op == BAM_CREF_SKIP) { //N
266                 //intron starts
267                 //exon ends here
268 <               exons.Add(new GSeg(exstart+1,l));
268 >               GSeg exon(exstart+1,c->pos+l);
269 >               exons.Add(exon);
270                 l += cigar[i]>>4;
271                 exstart=c->pos+l;
272                 }
273          }
274 <   exons.Add(new GSeg(exstart+1, l));
274 >   GSeg exon(exstart+1,c->pos+l);
275 >   exons.Add(exon);
276     end=c->pos+l;
277   }
278 <
278 > /*
279   int GBamRecord::find_tag(const char tag[2], uint8_t* & s, char& tag_type) {
280     //position s at the beginning of tag "data" (after the type char)
281     //returns the length of tag data, and tag type in tag_type
282 <   s=bam1_aux(b);
283 <   while (s < b->data + b->data_len) {
282 >   uint8_t* aux_start=bam1_aux(b);
283 >   s=aux_start;
284 >   while (s < aux_start + b->l_aux - 1) {
285       char key[2];
286       key[0] = (char)s[0]; key[1] = (char)s[1];
287       s += 2; tag_type = (char)*s; ++s;
# Line 317 | Line 319
319                           { m_inc = 2; }
320             else if ('i' == sub_type || 'I' == sub_type || 'f' == sub_type)
321                           { m_inc = 4; }
322 <           if (m_inc==0) GError("Error: invalid 'B' array subtype (%c)!\n",sub_type);
322 >           if (m_inc==0)
323 >              GError("Error: invalid 'B' array subtype (%c)!\n",sub_type);
324             inc += m_inc*n;
325             break;
326         case 'H':
327         case 'Z':
328           while (*(s+inc)) ++inc;
329 +          ++inc; // to skip the terminating \0
330           break;
331         } //switch (tag_type)
332 <     if (tag[0]==key[0] && tag[1]==key[1])
332 >     if (tag[0]==key[0] && tag[1]==key[1]) {
333          return inc;
334 +        }
335 +     s+=inc;
336       }//while aux data
337     return 0;
338     }
339 + */
340 + uint8_t* GBamRecord::find_tag(const char tag[2]) {
341 +   return bam_aux_get(this->b, tag);
342 + }
343  
344   char GBamRecord::tag_char(const char tag[2]) { //retrieve tag data as single char
345 <   uint8_t *s;
346 <   char tag_type=0;
347 <   int vlen=find_tag(tag, s, tag_type);
338 <   if (vlen==0) return 0;
339 <   //if (vlen>1) GWarning("Warning: tag %c%c value has length %d, but char was expected.\n",
340 <   //        tag[0],tag[1],vlen);
341 <   return (char)s[0];
345 >   uint8_t* s=find_tag(tag);
346 >   if (s) return ( bam_aux2A(s) );
347 >   return 0;
348    }
349  
350   int GBamRecord::tag_int(const char tag[2]) { //get the numeric value of tag
351 <   uint8_t *s;
352 <   char tag_type=0;
347 <   int vlen=find_tag(tag, s, tag_type);
348 <   if (vlen==0) return 0;
349 <   if (vlen==1) return (int)(*(int8_t*)s);
350 <    else if (vlen==2) return (int)(*(int16_t*)s);
351 <     else if (vlen==4) return (int)(*(int32_t*)s);
352 <     else GMessage("Warning: tag %c%c value has length %d, but int type was expected.\n",
353 <           tag[0],tag[1], vlen);
351 >   uint8_t *s=find_tag(tag);
352 >   if (s) return ( bam_aux2i(s) );
353     return 0;
354     }
355  
356   char* GBamRecord::tag_str(const char tag[2]) { //return string value for a tag
357 <   uint8_t *s;
358 <   char tag_type=0;
359 <   int vlen=find_tag(tag, s, tag_type);
361 <   if (vlen==0) return NULL; //not found
362 <   //if (vlen>1) GWarning("Warning: tag %c%c value has length %d, but char was expected.\n",
363 <   //        tag[0],tag[1],vlen);
364 <   return (char *)s;
357 >   uint8_t *s=find_tag(tag);
358 >   if (s) return ( bam_aux2Z(s) );
359 >   return NULL;
360     }
361  
362   char GBamRecord::spliceStrand() { // '+', '-' from the XS tag, or 0 if no XS tag
363 <   uint8_t *s;
364 <   char tag_type=0;
365 <   int vlen=find_tag("XS", s, tag_type);
371 <   if (vlen==0) return '.';
372 <   return (char)s[0];
363 >   char c=tag_char("XS");
364 >   if (c) return c;
365 >     else return '.';
366     }
367  
368   char* GBamRecord::sequence() { //user must free this after use

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines