ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GBam.cpp
(Generate patch)
# Line 248 | Line 248
248    this->add_aux(tag, atype, alen, adata);
249    }//add_aux()
250  
251 + void interpret_CIGAR(char cop, int cl, int aln_start) {
252 + // NB: Closed ranges
253 + // gpos = current genomic position (will end up as right coordinate on the genome)
254 + // rpos = read position (will end up as the length of the read)
255 + // cop = CIGAR operation, cl = operation length
256 + int rpos = 0;
257 + int gpos = aln_start;
258 + int num_mismatches=0; //NM tag value = edit distance
259 + switch (cop) {
260 + case BAM_CDIFF:
261 +      num_mismatches+=cl;
262 + case BAM_CMATCH:
263 +      //have to actually check for mismatches: num_mismatches+=count_mismatches;
264 + case BAM_CEQUAL:
265 +      printf("[%d-%d]", gpos, gpos + cl - 1);
266 +      gpos+=cl;
267 +      rpos+=cl;
268 +      break;
269 + case BAM_CPAD:
270 +      // printf("[%d-%d]", pos, pos + cl - 1); // Spans positions, No Coverage
271 +      gpos+=cl;
272 +      break;
273 +
274 + case BAM_CHARD_CLIP:
275 +      // printf("[%d]", pos);  // No coverage
276 +      // gpos is not advanced by this operation
277 +      break;
278 +
279 + case BAM_CSOFT_CLIP:
280 +      //soft clipped bases, present in SEQ
281 +      rpos+=cl;
282 +      break;
283 + case BAM_CINS:
284 +      // No Coverage
285 +      // adds cl bases "throughput" but not genomic position "coverage" (gap in genomic seq)
286 +      // should also add cl to the number of "mismatches" (unaligned read bases)
287 +      num_mismatches+=cl;
288 +      // How you handle this is application dependent
289 +      // gpos is not advanced by this operation
290 +      rpos+=cl;
291 +      break;
292 +
293 + case BAM_CDEL:
294 +      printf("D"); //deletion in reference sequence relative to the read (gap in read sequence)
295 +      // printf("[%d-%d]", pos, pos + cl - 1);  
296 +      // Spans positions
297 +      num_mismatches+=cl;
298 +      gpos += cl;
299 +      break;
300 +
301 + case BAM_CREF_SKIP:
302 +      // intron
303 +      printf("S");
304 +      // printf("[%d-%d]", pos, pos + cl - 1); // Spans positions, No Coverage
305 +      //special skip operation, not contributing to "edit distance",
306 +      //   so num_mismatches is not updated
307 +      gpos+=cl;
308 +      break;
309 +
310 + default:
311 +      fprintf(stderr, "Unhandled cigar_op %d:%d\n", cop, cl);
312 +      //printf("?");
313 +  }
314 + } // interpret_CIGAR(), just a reference of CIGAR operations interpretation
315  
316   void GBamRecord::setupCoordinates() {
317     uint32_t *cigar = bam1_cigar(b);
318     const bam1_core_t *c = &b->core;
319     if (c->flag & BAM_FUNMAP) return; /* skip unmapped reads */
320     int l=0;
321 <   start=c->pos+1;
321 >   start=c->pos+1; //genomic start coordinate
322     int exstart=c->pos;
323     for (int i = 0; i < c->n_cigar; ++i) {
324          int op = cigar[i]&0xf;
# Line 273 | Line 337
337          }
338     GSeg exon(exstart+1,c->pos+l);
339     exons.Add(exon);
340 <   end=c->pos+l;
340 >   end=c->pos+l; //genomic start coordinate
341   }
342   /*
343   int GBamRecord::find_tag(const char tag[2], uint8_t* & s, char& tag_type) {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines