How does htslib/samtools access optional BAM fields?

Question

I am confused regarding how various software libraries deal with optional fields in a BAM:

Based upon the BAM specification, there are 11 mandatory fields to a BAM:

QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL

The BAM files of today normally have many optional fields as well, which must be of the format TAG:TYPE:VALUE, whereby TAG is a two-character string that matches /[A-Za-z][A-Za-z0-9]/ and TYPE is a single case-sensitive letter which defines the format of VALUE. Each TAG can only appear once in one alignment line.

Certain BAMs may have no optional fields, while others have several.

My understanding is that optional tags may be in a different order from row to row, and then sometimes they simply do not exist.

My questions:

1. How does htslib/samtools (e.g. samtools view) access these optional fields, which differ in number from BAM to BAM (possible row to row)?

   In htslib/sam.h, a BAM alignment struct is defined as follows:

   /*! @typedef
    @abstract Structure for one alignment.
    @field  core       core information about the alignment
    @field  l_data     current length of bam1_t::data
    @field  m_data     maximum length of bam1_t::data
    @field  data       all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux
   typedef struct {
       bam1_core_t core;
       int l_data;
       uint32_t m_data;
       uint8_t *data;
   #ifndef BAM_NO_ID
       uint64_t id;
   #endif
   } bam1_t;
   

   and followed by several methods which access mandatory BAM fields, e.g.

   /*! @function
    @abstract  Get the name of the query
    @param  b  pointer to an alignment
    @return    pointer to the name string, null terminated
    */
   #define bam_get_qname(b) ((char*)(b)->data)
   

2. How does htslib access each of the optional fields, given this structure? Does it keep the order of the optional fields, or is this information unavailable?

3. When a user runs htslib or even samtools view, how does the BAM row "end", i.e. what data structure exists to show that the BAM row has no more optional fields, and that the algorithm should "move on" to the next BAM row.

Any clarification for how the source code works in this case is much appreciated.

Answer 1

How does htslib access each of the optional fields, given this structure?

htslib uses a macro, bam_get_aux, for this purpose. You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data. But you shouldn’t care about this, and just use the relevant exported functions for reading optional (“auxiliary”) data:

• bam_aux_get returns a pointer to an optional record, given by its two-character tag name.
• bam_aux2‹C› converts the pointer returned by the previous function to the actual type of the optional tag; here, ‹C› is a type identifier: i for an integer, Z for a zero-terminated character string, etc. Refer to the documentation and the htslib/sam.h file for a complete list.

Does it keep the order of the optional fields, or is this information unavailable?

As long as you don’t modify optional fields, the order is kept, but unspecified. Do not rely on the order!

When a user runs htslib or even samtools view, how does the BAM row "end"

BAM records store a header that contains the record size.

Answer 2

1. In a BAM file, all of the auxiliary tags (the optional fields you refer to) are held what is effectively a void * field, one after the other. Each entry within this field is structured as a two chars storing the tag name (e.g., NM, followed by a char storing the type (e.g., Z or i), followed by the data in whatever format is appropriate (char, float, etc). Strings in the aux tags are NULL terminated, so their lengths aren't explicitly stored.
2. To access a given field, htslib has to iterate over all of them, in whatever order they happen to be stored, until it finds the correct one. The output order with samtools view is the order in which they're stored.
3. The size of *data is stored as l_data, so once it hits that then it knows that the record is finished.

Regarding point 1, this is a standard C trick of casting between data types. So, while *data is actually a uint8_t * (aka, a byte array), the first two bytes of each tag are chars and whatever is after them (well, after the byte specifying the type) is then cast to the appropriate type (be it float or char* or what not). This is all probably clearer if you look at the code.