What's the correct way to store bam record in a vector and free them?

Here's my demo code

#include <stdio.h>
#include "htslib/sam.h"
#include <vector>
#include <fstream>
#include <string>
#include <iostream>

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_queryi(idx, 1,2125034,2125394);
  std::vector<bam1_t> record_list(0);
  record_list.reserve(100000);
  while (sam_itr_next(fp, itr, b) >= 0) {
    record_list.resize(record_list.size()+1);
    bam_copy1(&record_list[record_list.size()-1],b);
  }

  for(auto& bb : record_list)
  {
    printf("%d\t%d\t%d\t%s\n", bam_name2id(h,h->target_name[bb.core.tid]), bb.core.pos,
      bb.core.pos + bam_cigar2rlen(bb.core.n_cigar, bam_get_cigar(&bb)), bam_get_qname(&bb));
    //bam_destroy1(&bb);  //segment fault
  }

  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
  return 0;
}

For this code everything looks fine, but if we uncomment the line 'bam_destroy1(&bb)' would lead to segment fault. Here's two questions:

Is the resize-copy a correct way to store bam record in a vector?
Why bam_destroy1() lead to segment fault? what's the correct way to free those memory?

htslib

• 1,049 views • link

updated 3.2 years ago by Pierre Lindenbaum

• written 3.2 years ago by PeterWu •

1 answer

You should allocate a new pointer for each SAM record. I would do something like this (not tested):

  std::vector<bam1_t*> record_list;
  record_list.reserve(100000);
  for(;;) {
    bam1_t *b = bam_init1();
   if(sam_itr_next(fp, itr, b) >= 0) {
        record_list.push_back(b);
        }
  else {
       bam_destroy1(b); 
       break;
       }
  }
(...)
for(std::vector<bam1_t*>::size_type i=0;i< record_list.size();i++) {
   bam_destroy1(record_list[i]); 
}

• 0 views • link

written 3.2 years ago by Pierre Lindenbaum