How to build a BAM header file with htslib in C++?

Question

I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly.

https://github.com/samtools/htslib/blob/develop/test/test-bcf-translate.c

The two key functions are:

• bcf_hdr_append
• bcf_hdr_add_sample

Now when I navigate to the header file for SAM/BAM functions:

https://github.com/samtools/htslib/blob/develop/htslib/sam.h

The only "append" function is bam_aux_append. But that's for optional fields.

My questions:

• Given a bam_hdr_t header (potentially empty), how to add a new chromosome reference? (For example, add a new chrAB to the BAM file)
• Given a non-empty bam_hdr_t * header, how to remove a chromosome? (For example, remove chr4 from the BAM header)

My sample code:

// Empty BAM header, nothing inside
bam_hdr_t * hdr = bam_hdr_init();

// I'd like to add a new synthetic `chrAB` into the header, how?  

Answer 1

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there).

With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

1. Increment hdr->ntargets
2. realloc(hdr->target_len, ...) and hdr->target_name.
3. You likely want to append the name/length to the end of hdr->target_name and hdr->target_len, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).

Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.

WARNING, doing this is prone to error!

I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into hdr->target_name, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify all of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with samtools reheader).