Read count in a specific region [samtools]

[serhatg]

Dear all,

I'd like to get the number of reads in a specified region with samtools. Is there a function which does that ? The method shouldn't be OS dependent.

Thanks in advance,
Serhat

[GenoMax]

http://seqanswers.com/forums/showthread.php?t=33357

Using an indexed BAM file

$ samtools view input.bam chr1:1234-2345 | wc -l

[serhatg]

http://seqanswers.com/forums/showthread.php?t=33357

Using an indexed BAM file

$ samtools view input.bam chr1:1234-2345 | wc -l

Thanks for your answer GenoMax;however, I am looking for a method which is not depended on the OS.

P.S. I am using samtools libraries in a cpp project.

[GenoMax]

I missed the programmatic part of your original question. Someone else should be along with the right answer.

[dpryan]

I guess that someone is me.

With the older samtools 0.1.19 API, you can just use the bam_fetch() function and give it a function to just increment a counter with each call.

That function doesn't exist in htslib, so you either need to roll your own version:

//N.B., there are some functions here you'll need to replace
int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func)
{
    int ret;
    bam_iter_t iter;
    bam1_t *b;
    b = bam_init1();
    iter = bam_iter_query(idx, tid, beg, end);
    while ((ret = bam_iter_read(fp, iter, b)) >= 0) func(b, data);
    bam_iter_destroy(iter);
    bam_destroy1(b);
    return (ret == -1)? 0 : ret;
}

Alternatively, you can statically compile in libbam.a from the newer version of samtools and the API should be more or less the same as the previous version (annoyingly, there's little documentation at the moment).

[serhatg]

I guess that someone is me.

With the older samtools 0.1.19 API, you can just use the bam_fetch() function and give it a function to just increment a counter with each call.

That function doesn't exist in htslib, so you either need to roll your own version:

//N.B., there are some functions here you'll need to replace
int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func)
{
    int ret;
    bam_iter_t iter;
    bam1_t *b;
    b = bam_init1();
    iter = bam_iter_query(idx, tid, beg, end);
    while ((ret = bam_iter_read(fp, iter, b)) >= 0) func(b, data);
    bam_iter_destroy(iter);
    bam_destroy1(b);
    return (ret == -1)? 0 : ret;
}

Alternatively, you can statically compile in libbam.a from the newer version of samtools and the API should be more or less the same as the previous version (annoyingly, there's little documentation at the moment).

Thank you for the answer. I was planning to use exactly the same method; however, then I thought there might be a smarter way or a method which I am not aware of. Apparently, I will need to stick to this method.

Thanks again

[dpryan]

Well there might be a smarter way, but this is at least the simplest way

[rnaeye]

You can call samtools command from a Perl script that also counts the lines.

[syfo]

samtools view input.bam chr1:1234-2345 | wc -l

You can call samtools command from a Perl script that also counts the lines.

On a side note, you can use the "samtools view -c" option instead of counting "yourself" the number of lines.