<p>Hi, this question is somehow related to <a href="http://biostar.stackexchange.com/questions/1681">this previous question</a>. I'm playing with the <strong>C API for BAM</strong> and I wrote the following code: <a href="https://gist.github.com/736059">https://gist.github.com/736059</a></p>

<p>I'm now trying to find the indexes of the genomic bases covered by the <strong>CIGAR</strong> string (my final goal is to create a <strong>WIG</strong> file containing the coverage of the genome)</p>

<p>What would be the correct way to change my code to get the genomic indexes of the bases covered by each CIGAR element ?</p>

<pre><code>(...)
for( k=0;k&lt; b-&gt;core.n_cigar;++k)
{
int cop =cigar[k] &amp; BAM_CIGAR_MASK; // operation
int cl = cigar[k] &gt;&gt; BAM_CIGAR_SHIFT; // length
switch(cop)
 {
 case BAM_CMATCH: printf("M");break;
 case BAM_CINS: printf("I");break;
 case BAM_CDEL: printf("D");break;
 case BAM_CREF_SKIP: printf("N"); break;
 case BAM_CSOFT_CLIP: printf("S");break;
 case BAM_CHARD_CLIP: printf("R");break;
 case BAM_CPAD: printf("P");break;
 default:printf("?");break;
 }
printf("%d",cl);
}
(...)
</code></pre>

<p>Thanks</p>

For example

```
int depth[REF_LEN], x, j, k;
uint32_t *cigar = bam1_cigar(b);
for (k = 0, x = b->core.pos; k < b->core.n_cigar; ++k) {
  int op = cigar[k]&16;
  int l = cigar[k]>>4;
  if (op == BAM_CMATCH) {
    for (j = x; j < x + l; ++j) ++depth[j];
    x += l;
  } else if (op == BAM_CREF_SKIP || op == BAM_CDEL) x += l;
}
```

Something like this...