1
$\begingroup$

I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the following samtools command in C++:

samtools view name.bam -o name.sam scaffold:pos-pos

Since I must extract reads from thousands of regions, I do not want to iterate through the whole bam file each time I extract reads from a region. I am wondering if HTSlib can extract reads occurring in a specified scaffold and range without iterating through the whole bam. I found the function bam_itr_querys() in the htslib/sam.h documentation on GitHub, but don't know how to format the inputs for this function, or what the "idx" input is. I also saw somewhere that I would need to pair this function with bam_itr_next(), but don't know the format to do so.

EDIT: My bam files are sorted and indexed, if that matters.

$\endgroup$
3
  • $\begingroup$ Have you looked at how samtools view works with a bed file of coordinates as input? $\endgroup$
    – swbarnes2
    Jun 30 at 20:14
  • $\begingroup$ @swbarnes2 I have not, but my issue is that I cannot use 'samtools view'. I should provide some context (I will edit the original question to include this info): in my first draft of my C++ program, I simply used the system() command to run 'samtools view' where needed. I am now wrapping the original program in an R package, and simply running the "system" command does not work. So, I need to run something with the same functionality as "samtools view" from within C++ itself. $\endgroup$ Jul 1 at 19:40
  • $\begingroup$ @annabelperry I can strongly recommend looking around at the internals of samtools if you need to work with htslib. Samtools is very well put together and I've found scrolling through its code to be a great way to gain insight into how to structure my own programs using htslib. $\endgroup$
    – Devon Ryan
    Jul 3 at 20:27
2
$\begingroup$
#include <stdio.h>
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  hts_idx_t *idx = sam_index_load(fp, argv[1]);
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]);
  while (sam_itr_next(fp, itr, b) >= 0) {
    if (b->core.tid < 0) continue;
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos,
      b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)));
  }
  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
}

The usage would then be ./a.out alignments.bam chr1:123-456, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via sam_itr_querys() so it's not iterating over the entire file.

$\endgroup$
3
  • $\begingroup$ So I'm looking for a solution which doesn't require iterating through the whole bam file. Unless I am missing something, this solution requires that you iterate through the whole bam file and add "if-else" statements within the "for" loop to extract out your reads of interest. I am trying to avoid this, as it takes so much time. I'm looking for a solution which searches the binary bam file for reads in specific scaffolds/regions WITHOUT iterating through every read $\endgroup$ Jul 3 at 19:19
  • 2
    $\begingroup$ @annabelperry That isn't reading through the entire file. sam_itr_querys is accepting a string like chr1:123-456 and then making an iterator over just those using the index (idx). This is all that samtools is doing internally. I've added a bit onto the end of this message to hopefully clarify that. $\endgroup$
    – Devon Ryan
    Jul 3 at 20:19
  • 1
    $\begingroup$ Thank you, @DevonRyan, that clears it up! $\endgroup$ Jul 3 at 23:18
0
$\begingroup$

Have you considered using a sorted and indexed bam file? e.g.

samtools sort in.bam > out.sorted.bam
samtools index out.sorted.bam

This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the AlignmentFile.fetch() functionality.

Surely there is an analog in HTSlib API, though I don't know what to search to get it. It is in samtools as demonstrated here.

This answer seems to put the components together that you are interested in, among them sam_index_load(), sam_itr_querys().

$\endgroup$
1
  • $\begingroup$ My bam files are already sorted and indexed - I will check out this other answer $\endgroup$ Jul 3 at 19:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.