htslib C API: getting reads from region

Question

[QVINTVS_FABIVS_MAXIMVS]

I want to use the htslib C API to perform this task

my $pos = "1:100-1000";
my $bam_file = "in.bam";
my @reads = get_reads_from_htslib($bam_file,$pos)

foreach(@reads) {
    # do stuff...
}

I wrote it in perl but you get the idea. I can't find an example anywhere online. Does anyone have experience with htslib and getting reads from a region?

Answer 1

[Pierre Lindenbaum]

here is a C code:

Compilation:

gcc -Isamtools -Ihtslib -Lsamtools -Lhtslib biostar151053.c htslib/libhts.a samtools/libbam.a -lz -lpthread

[QVINTVS_FABIVS_MAXIMVS]

Livesaver thanks!

[Danielk]

What function calls here come from samtools.h?

[Pierre Lindenbaum]

may be none :-)

[lghust2011]

Dear Pierre, I have another question. Now I have a sam file without the index file, and I want to iterate all the reads of the sam file, how can I make it? Thanks in advance.

[lghust2011]

It seems that I have to compress the sam file to bam file, and get the index file then I can make it..

[sina.majidian]

Dear Pierre I used the exact code and bash line but I face an error. I've downloaded htslib and samtools and put them in the same folder as the .c file.

.... "_lzma_stream_decoder", referenced from: _cram_uncompress_block in libhts.a(cram_io.o) ld: symbol(s) not found for architecture x86_64 clang: error: linker command failed with exit code 1 (use -v to see invocation)

Should I build samtools and htslib as mentioned here ?

[Pierre Lindenbaum]

I think so, furthermore, the code above is ~5 years old, the htslib has since much changed.

[sina.majidian]

Thanks for fast respond.