bam.h

   1 /* The MIT License
   2
   3    Copyright (c) 2008 Genome Research Ltd (GRL).
   4
   5    Permission is hereby granted, free of charge, to any person obtaining
   6    a copy of this software and associated documentation files (the
   7    "Software"), to deal in the Software without restriction, including
   8    without limitation the rights to use, copy, modify, merge, publish,
   9    distribute, sublicense, and/or sell copies of the Software, and to
  10    permit persons to whom the Software is furnished to do so, subject to
  11    the following conditions:
  12
  13    The above copyright notice and this permission notice shall be
  14    included in all copies or substantial portions of the Software.
  15
  16    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  18    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  20    BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  21    ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  22    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23    SOFTWARE.
  24 */
  25
  26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
  27
  28 #ifndef BAM_BAM_H
  29 #define BAM_BAM_H
  30
  31 /*!
  32   @header
  33
  34   BAM library provides I/O and various operations on manipulating files
  35   in the BAM (Binary Alignment/Mapping) or SAM (Sequence Alignment/Map)
  36   format. It now supports importing from or exporting to TAM, sorting,
  37   merging, generating pileup, and quickly retrieval of reads overlapped
  38   with a specified region.
  39
  40   @copyright Genome Research Ltd.
  41  */
  42
  43 #include <stdint.h>
  44 #include <stdlib.h>
  45 #include <string.h>
  46 #include <stdio.h>
  47
  48 #ifndef BAM_LITE
  49 #define BAM_VIRTUAL_OFFSET16
  50 #include "bgzf.h"
  51 /*! @abstract BAM file handler */
  52 typedef BGZF *bamFile;
  53 #define bam_open(fn, mode) bgzf_open(fn, mode)
  54 #define bam_dopen(fd, mode) bgzf_fdopen(fd, mode)
  55 #define bam_close(fp) bgzf_close(fp)
  56 #define bam_read(fp, buf, size) bgzf_read(fp, buf, size)
  57 #define bam_write(fp, buf, size) bgzf_write(fp, buf, size)
  58 #define bam_tell(fp) bgzf_tell(fp)
  59 #define bam_seek(fp, pos, dir) bgzf_seek(fp, pos, dir)
  60 #else
  61 #define BAM_TRUE_OFFSET
  62 #include <zlib.h>
  63 typedef gzFile bamFile;
  64 #define bam_open(fn, mode) gzopen(fn, mode)
  65 #define bam_dopen(fd, mode) gzdopen(fd, mode)
  66 #define bam_close(fp) gzclose(fp)
  67 #define bam_read(fp, buf, size) gzread(fp, buf, size)
  68 /* no bam_write/bam_tell/bam_seek() here */
  69 #endif
  70
  71 /*! @typedef
  72   @abstract Structure for the alignment header.
  73   @field n_targets   number of reference sequences
  74   @field target_name names of the reference sequences
  75   @field target_len  lengths of the referene sequences
  76   @field dict        header dictionary
  77   @field hash        hash table for fast name lookup
  78   @field rg2lib      hash table for @RG-ID -> LB lookup
  79   @field l_text      length of the plain text in the header
  80   @field text        plain text
  81
  82   @discussion Field hash points to null by default. It is a private
  83   member.
  84  */
  85 typedef struct {
  86         int32_t n_targets;
  87         char **target_name;
  88         uint32_t *target_len;
  89         void *dict, *hash, *rg2lib;
  90         int l_text;
  91         char *text;
  92 } bam_header_t;
  93
  94 /*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
  95 #define BAM_FPAIRED        1
  96 /*! @abstract the read is mapped in a proper pair */
  97 #define BAM_FPROPER_PAIR   2
  98 /*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
  99 #define BAM_FUNMAP         4
 100 /*! @abstract the mate is unmapped */
 101 #define BAM_FMUNMAP        8
 102 /*! @abstract the read is mapped to the reverse strand */
 103 #define BAM_FREVERSE      16
 104 /*! @abstract the mate is mapped to the reverse strand */
 105 #define BAM_FMREVERSE     32
 106 /*! @abstract this is read1 */
 107 #define BAM_FREAD1        64
 108 /*! @abstract this is read2 */
 109 #define BAM_FREAD2       128
 110 /*! @abstract not primary alignment */
 111 #define BAM_FSECONDARY   256
 112 /*! @abstract QC failure */
 113 #define BAM_FQCFAIL      512
 114 /*! @abstract optical or PCR duplicate */
 115 #define BAM_FDUP        1024
 116
 117 #define BAM_OFDEC          0
 118 #define BAM_OFHEX          1
 119 #define BAM_OFSTR          2
 120
 121 /*! @abstract defautl mask for pileup */
 122 #define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP)
 123
 124 #define BAM_CORE_SIZE   sizeof(bam1_core_t)
 125
 126 /**
 127  * Describing how CIGAR operation/length is packed in a 32-bit integer.
 128  */
 129 #define BAM_CIGAR_SHIFT 4
 130 #define BAM_CIGAR_MASK  ((1 << BAM_CIGAR_SHIFT) - 1)
 131
 132 /*
 133   CIGAR operations.
 134  */
 135 /*! @abstract CIGAR: match */
 136 #define BAM_CMATCH      0
 137 /*! @abstract CIGAR: insertion to the reference */
 138 #define BAM_CINS        1
 139 /*! @abstract CIGAR: deletion from the reference */
 140 #define BAM_CDEL        2
 141 /*! @abstract CIGAR: skip on the reference (e.g. spliced alignment) */
 142 #define BAM_CREF_SKIP   3
 143 /*! @abstract CIGAR: clip on the read with clipped sequence present in qseq */
 144 #define BAM_CSOFT_CLIP  4
 145 /*! @abstract CIGAR: clip on the read with clipped sequence trimmed off */
 146 #define BAM_CHARD_CLIP  5
 147 /*! @abstract CIGAR: padding */
 148 #define BAM_CPAD        6
 149
 150 /*! @typedef
 151   @abstract Structure for core alignment information.
 152   @field  tid     chromosome ID, defined by bam_header_t
 153   @field  pos     0-based leftmost coordinate
 154   @field  strand  strand; 0 for forward and 1 otherwise
 155   @field  bin     bin calculated by bam_reg2bin()
 156   @field  qual    mapping quality
 157   @field  l_qname length of the query name
 158   @field  flag    bitwise flag
 159   @field  n_cigar number of CIGAR operations
 160   @field  l_qseq  length of the query sequence (read)
 161  */
 162 typedef struct {
 163         int32_t tid;
 164         int32_t pos;
 165         uint32_t bin:16, qual:8, l_qname:8;
 166         uint32_t flag:16, n_cigar:16;
 167         int32_t l_qseq;
 168         int32_t mtid;
 169         int32_t mpos;
 170         int32_t isize;
 171 } bam1_core_t;
 172
 173 /*! @typedef
 174   @abstract Structure for one alignment.
 175   @field  core       core information about the alignment
 176   @field  l_aux      length of auxiliary data
 177   @field  data_len   current length of bam1_t::data
 178   @field  m_data     maximum length of bam1_t::data
 179   @field  data       all variable-length data, concatenated; structure: cigar-qname-seq-qual-aux
 180
 181   @discussion Notes:
 182
 183    1. qname is zero tailing and core.l_qname includes the tailing '\0'.
 184    2. l_qseq is calculated from the total length of an alignment block
 185       on reading or from CIGAR.
 186  */
 187 typedef struct {
 188         bam1_core_t core;
 189         int l_aux, data_len, m_data;
 190         uint8_t *data;
 191 } bam1_t;
 192
 193 #define bam1_strand(b) (((b)->core.flag&BAM_FREVERSE) != 0)
 194 #define bam1_mstrand(b) (((b)->core.flag&BAM_FMREVERSE) != 0)
 195
 196 /*! @function
 197   @abstract  Get the CIGAR array
 198   @param  b  pointer to an alignment
 199   @return    pointer to the CIGAR array
 200
 201   @discussion In the CIGAR array, each element is a 32-bit integer. The
 202   lower 4 bits gives a CIGAR operation and the higher 28 bits keep the
 203   length of a CIGAR.
 204  */
 205 #define bam1_cigar(b) ((uint32_t*)((b)->data + (b)->core.l_qname))
 206
 207 /*! @function
 208   @abstract  Get the name of the query
 209   @param  b  pointer to an alignment
 210   @return    pointer to the name string, null terminated
 211  */
 212 #define bam1_qname(b) ((char*)((b)->data))
 213
 214 /*! @function
 215   @abstract  Get query sequence
 216   @param  b  pointer to an alignment
 217   @return    pointer to sequence
 218
 219   @discussion Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G,
 220   8 for T and 15 for N. Two bases are packed in one byte with the base
 221   at the higher 4 bits having smaller coordinate on the read. It is
 222   recommended to use bam1_seqi() macro to get the base.
 223  */
 224 #define bam1_seq(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname)
 225
 226 /*! @function
 227   @abstract  Get query quality
 228   @param  b  pointer to an alignment
 229   @return    pointer to quality string
 230  */
 231 #define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + ((b)->core.l_qseq + 1)/2)
 232
 233 /*! @function
 234   @abstract  Get a base on read
 235   @param  s  Query sequence returned by bam1_seq()
 236   @param  i  The i-th position, 0-based
 237   @return    4-bit integer representing the base.
 238  */
 239 #define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
 240
 241 /*! @function
 242   @abstract  Get query sequence and quality
 243   @param  b  pointer to an alignment
 244   @return    pointer to the concatenated auxiliary data
 245  */
 246 #define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2)
 247
 248 #ifndef kroundup32
 249 /*! @function
 250   @abstract  Round an integer to the next closest power-2 integer.
 251   @param  x  integer to be rounded (in place)
 252   @discussion x will be modified.
 253  */
 254 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
 255 #endif
 256
 257 /*!
 258   @abstract Whether the machine is big-endian; modified only in
 259   bam_header_init().
 260  */
 261 extern int bam_is_be;
 262
 263 /*! @abstract Table for converting a nucleotide character to the 4-bit encoding. */
 264 extern unsigned char bam_nt16_table[256];
 265
 266 /*! @abstract Table for converting a 4-bit encoded nucleotide to a letter. */
 267 extern char *bam_nt16_rev_table;
 268
 269 extern char bam_nt16_nt4_table[];
 270
 271 #ifdef __cplusplus
 272 extern "C" {
 273 #endif
 274
 275         /*! @abstract TAM file handler */
 276         typedef struct __tamFile_t *tamFile;
 277
 278         /*!
 279           @abstract   Open a SAM file for reading, either uncompressed or compressed by gzip/zlib.
 280           @param  fn  SAM file name
 281           @return     SAM file handler
 282          */
 283         tamFile sam_open(const char *fn);
 284
 285         /*!
 286           @abstract   Close a SAM file handler
 287           @param  fp  SAM file handler
 288          */
 289         void sam_close(tamFile fp);
 290
 291         /*!
 292           @abstract      Read one alignment from a SAM file handler
 293           @param  fp     SAM file handler
 294           @param  header header information (ordered names of chromosomes)
 295           @param  b      read alignment; all members in b will be updated
 296           @return        0 if successful; otherwise negative
 297          */
 298         int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b);
 299
 300         /*!
 301           @abstract       Read header information from a TAB-delimited list file.
 302           @param  fn_list file name for the list
 303           @return         a pointer to the header structure
 304
 305           @discussion Each line in this file consists of chromosome name and
 306           the length of chromosome.
 307          */
 308         bam_header_t *sam_header_read2(const char *fn_list);
 309
 310         /*!
 311           @abstract       Read header from a SAM file (if present)
 312           @param  fp      SAM file handler
 313           @return         pointer to header struct; 0 if no @SQ lines available
 314          */
 315         bam_header_t *sam_header_read(tamFile fp);
 316
 317         /*!
 318           @abstract       Parse @SQ lines a update a header struct
 319           @param  h       pointer to the header struct to be updated
 320           @return         number of target sequences
 321
 322           @discussion bam_header_t::{n_targets,target_len,target_name} will
 323           be destroyed in the first place.
 324          */
 325         int sam_header_parse(bam_header_t *h);
 326
 327         /*!
 328           @abstract       Parse @RG lines a update a header struct
 329           @param  h       pointer to the header struct to be updated
 330           @return         number of @RG lines
 331
 332           @discussion bam_header_t::rg2lib will be destroyed in the first
 333           place.
 334          */
 335         int sam_header_parse_rg(bam_header_t *h);
 336
 337 #define sam_write1(header, b) bam_view1(header, b)
 338
 339         int bam_strmap_put(void *strmap, const char *rg, const char *lib);
 340         const char *bam_strmap_get(const void *strmap, const char *rg);
 341         void *bam_strmap_dup(const void*);
 342         void *bam_strmap_init();
 343         void bam_strmap_destroy(void *strmap);
 344
 345         /*!
 346           @abstract Initialize a header structure.
 347           @return   the pointer to the header structure
 348
 349           @discussion This function also modifies the global variable
 350           bam_is_be.
 351          */
 352         bam_header_t *bam_header_init();
 353
 354         /*!
 355           @abstract        Destroy a header structure.
 356           @param  header  pointer to the header
 357          */
 358         void bam_header_destroy(bam_header_t *header);
 359
 360         /*!
 361           @abstract   Read a header structure from BAM.
 362           @param  fp  BAM file handler, opened by bam_open()
 363           @return     pointer to the header structure
 364
 365           @discussion The file position indicator must be placed at the
 366           beginning of the file. Upon success, the position indicator will
 367           be set at the start of the first alignment.
 368          */
 369         bam_header_t *bam_header_read(bamFile fp);
 370
 371         /*!
 372           @abstract      Write a header structure to BAM.
 373           @param  fp     BAM file handler
 374           @param  header pointer to the header structure
 375           @return        always 0 currently
 376          */
 377         int bam_header_write(bamFile fp, const bam_header_t *header);
 378
 379         /*!
 380           @abstract   Read an alignment from BAM.
 381           @param  fp  BAM file handler
 382           @param  b   read alignment; all members are updated.
 383           @return     number of bytes read from the file
 384
 385           @discussion The file position indicator must be
 386           placed right before an alignment. Upon success, this function
 387           will set the position indicator to the start of the next
 388           alignment. This function is not affected by the machine
 389           endianness.
 390          */
 391         int bam_read1(bamFile fp, bam1_t *b);
 392
 393         /*!
 394           @abstract Write an alignment to BAM.
 395           @param  fp       BAM file handler
 396           @param  c        pointer to the bam1_core_t structure
 397           @param  data_len total length of variable size data related to
 398                            the alignment
 399           @param  data     pointer to the concatenated data
 400           @return          number of bytes written to the file
 401
 402           @discussion This function is not affected by the machine
 403           endianness.
 404          */
 405         int bam_write1_core(bamFile fp, const bam1_core_t *c, int data_len, uint8_t *data);
 406
 407         /*!
 408           @abstract   Write an alignment to BAM.
 409           @param  fp  BAM file handler
 410           @param  b   alignment to write
 411           @return     number of bytes written to the file
 412
 413           @abstract It is equivalent to:
 414             bam_write1_core(fp, &b->core, b->data_len, b->data)
 415          */
 416         int bam_write1(bamFile fp, const bam1_t *b);
 417
 418         /*! @function
 419           @abstract  Initiate a pointer to bam1_t struct
 420          */
 421 #define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t)))
 422
 423         /*! @function
 424           @abstract  Free the memory allocated for an alignment.
 425           @param  b  pointer to an alignment
 426          */
 427 #define bam_destroy1(b) do {                                    \
 428                 if (b) { free((b)->data); free(b); }    \
 429         } while (0)
 430
 431         /*!
 432           @abstract       Format a BAM record in the SAM format
 433           @param  header  pointer to the header structure
 434           @param  b       alignment to print
 435           @return         a pointer to the SAM string
 436          */
 437         char *bam_format1(const bam_header_t *header, const bam1_t *b);
 438
 439         char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of);
 440
 441         const char *bam_get_library(bam_header_t *header, const bam1_t *b);
 442
 443         /*! @typedef
 444           @abstract Structure for one alignment covering the pileup position.
 445           @field  b      pointer to the alignment
 446           @field  qpos   position of the read base at the pileup site, 0-based
 447           @field  indel  indel length; 0 for no indel, positive for ins and negative for del
 448           @field  is_del 1 iff the base on the padded read is a deletion
 449           @field  level  the level of the read in the "viewer" mode
 450
 451           @discussion See also bam_plbuf_push() and bam_lplbuf_push(). The
 452           difference between the two functions is that the former does not
 453           set bam_pileup1_t::level, while the later does. Level helps the
 454           implementation of alignment viewers, but calculating this has some
 455           overhead.
 456          */
 457         typedef struct {
 458                 bam1_t *b;
 459                 int32_t qpos;
 460                 int indel, level;
 461                 uint32_t is_del:1, is_head:1, is_tail:1;
 462         } bam_pileup1_t;
 463
 464         struct __bam_plbuf_t;
 465         /*! @abstract pileup buffer */
 466         typedef struct __bam_plbuf_t bam_plbuf_t;
 467
 468         void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask);
 469
 470         /*! @typedef
 471           @abstract    Type of function to be called by bam_plbuf_push().
 472           @param  tid  chromosome ID as is defined in the header
 473           @param  pos  start coordinate of the alignment, 0-based
 474           @param  n    number of elements in pl array
 475           @param  pl   array of alignments
 476           @param  data user provided data
 477           @discussion  See also bam_plbuf_push(), bam_plbuf_init() and bam_pileup1_t.
 478          */
 479         typedef int (*bam_pileup_f)(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data);
 480
 481         /*!
 482           @abstract     Reset a pileup buffer for another pileup process
 483           @param  buf   the pileup buffer to be reset
 484          */
 485         void bam_plbuf_reset(bam_plbuf_t *buf);
 486
 487         /*!
 488           @abstract     Initialize a buffer for pileup.
 489           @param  func  fucntion to be called by bam_pileup_core()
 490           @param  data  user provided data
 491           @return       pointer to the pileup buffer
 492          */
 493         bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data);
 494
 495         /*!
 496           @abstract    Destroy a pileup buffer.
 497           @param  buf  pointer to the pileup buffer
 498          */
 499         void bam_plbuf_destroy(bam_plbuf_t *buf);
 500
 501         /*!
 502           @abstract    Push an alignment to the pileup buffer.
 503           @param  b    alignment to be pushed
 504           @param  buf  pileup buffer
 505           @see         bam_plbuf_init()
 506           @return      always 0 currently
 507
 508           @discussion If all the alignments covering a particular site have
 509           been collected, this function will call the user defined function
 510           as is provided to bam_plbuf_init(). The coordinate of the site and
 511           all the alignments will be transferred to the user defined
 512           function as function parameters.
 513
 514           When all the alignments are pushed to the buffer, this function
 515           needs to be called with b equal to NULL. This will flush the
 516           buffer. A pileup buffer can only be reused when bam_plbuf_reset()
 517           is called.
 518          */
 519         int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf);
 520
 521         int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data);
 522
 523         struct __bam_lplbuf_t;
 524         typedef struct __bam_lplbuf_t bam_lplbuf_t;
 525
 526         void bam_lplbuf_reset(bam_lplbuf_t *buf);
 527
 528         /*! @abstract  bam_plbuf_init() equivalent with level calculated. */
 529         bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data);
 530
 531         /*! @abstract  bam_plbuf_destroy() equivalent with level calculated. */
 532         void bam_lplbuf_destroy(bam_lplbuf_t *tv);
 533
 534         /*! @abstract  bam_plbuf_push() equivalent with level calculated. */
 535         int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf);
 536
 537         struct __bam_index_t;
 538         typedef struct __bam_index_t bam_index_t;
 539
 540         /*!
 541           @abstract   Build index for a BAM file.
 542           @discussion Index file "fn.bai" will be created.
 543           @param  fn  name of the BAM file
 544           @return     always 0 currently
 545          */
 546         int bam_index_build(const char *fn);
 547
 548         /*!
 549           @abstract   Load index from file "fn.bai".
 550           @param  fn  name of the BAM file (NOT the index file)
 551           @return     pointer to the index structure
 552          */
 553         bam_index_t *bam_index_load(const char *fn);
 554
 555         /*!
 556           @abstract    Destroy an index structure.
 557           @param  idx  pointer to the index structure
 558          */
 559         void bam_index_destroy(bam_index_t *idx);
 560
 561         /*! @typedef
 562           @abstract      Type of function to be called by bam_fetch().
 563           @param  b     the alignment
 564           @param  data  user provided data
 565          */
 566         typedef int (*bam_fetch_f)(const bam1_t *b, void *data);
 567
 568         /*!
 569           @abstract Retrieve the alignments that are overlapped with the
 570           specified region.
 571
 572           @discussion A user defined function will be called for each
 573           retrieved alignment ordered by its start position.
 574
 575           @param  fp    BAM file handler
 576           @param  idx   pointer to the alignment index
 577           @param  tid   chromosome ID as is defined in the header
 578           @param  beg   start coordinate, 0-based
 579           @param  end   end coordinate, 0-based
 580           @param  data  user provided data (will be transferred to func)
 581           @param  func  user defined function
 582          */
 583         int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func);
 584
 585         /*!
 586           @abstract       Parse a region in the format: "chr2:100,000-200,000".
 587           @discussion     bam_header_t::hash will be initialized if empty.
 588           @param  header  pointer to the header structure
 589           @param  str     string to be parsed
 590           @param  ref_id  the returned chromosome ID
 591           @param  begin   the returned start coordinate
 592           @param  end     the returned end coordinate
 593           @return         0 on success; -1 on failure
 594          */
 595         int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end);
 596
 597         /*!
 598           @abstract       Retrieve data of a tag
 599           @param  b       pointer to an alignment struct
 600           @param  tag     two-character tag to be retrieved
 601
 602           @return  pointer to the type and data. The first character is the
 603           type that can be 'iIsScCdfAZH'.
 604
 605           @discussion  Use bam_aux2?() series to convert the returned data to
 606           the corresponding type.
 607         */
 608         uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]);
 609
 610         int32_t bam_aux2i(const uint8_t *s);
 611         float bam_aux2f(const uint8_t *s);
 612         double bam_aux2d(const uint8_t *s);
 613         char bam_aux2A(const uint8_t *s);
 614         char *bam_aux2Z(const uint8_t *s);
 615
 616         int bam_aux_del(bam1_t *b, uint8_t *s);
 617         void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);
 618         uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get()
 619
 620         /*!
 621           @abstract Calculate the rightmost coordinate of an alignment on the
 622           reference genome.
 623
 624           @param  c      pointer to the bam1_core_t structure
 625           @param  cigar  the corresponding CIGAR array (from bam1_t::cigar)
 626           @return        the rightmost coordinate, 0-based
 627         */
 628         uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar);
 629
 630         /*!
 631           @abstract      Calculate the length of the query sequence from CIGAR.
 632           @param  c      pointer to the bam1_core_t structure
 633           @param  cigar  the corresponding CIGAR array (from bam1_t::cigar)
 634           @return        length of the query sequence
 635         */
 636         int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar);
 637
 638 #ifdef __cplusplus
 639 }
 640 #endif
 641
 642 /*!
 643   @abstract    Calculate the minimum bin that contains a region [beg,end).
 644   @param  beg  start of the region, 0-based
 645   @param  end  end of the region, 0-based
 646   @return      bin
 647  */
 648 static inline int bam_reg2bin(uint32_t beg, uint32_t end)
 649 {
 650         --end;
 651         if (beg>>14 == end>>14) return 4681 + (beg>>14);
 652         if (beg>>17 == end>>17) return  585 + (beg>>17);
 653         if (beg>>20 == end>>20) return   73 + (beg>>20);
 654         if (beg>>23 == end>>23) return    9 + (beg>>23);
 655         if (beg>>26 == end>>26) return    1 + (beg>>26);
 656         return 0;
 657 }
 658
 659 /*!
 660   @abstract     Copy an alignment
 661   @param  bdst  destination alignment struct
 662   @param  bsrc  source alignment struct
 663   @return       pointer to the destination alignment struct
 664  */
 665 static inline bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
 666 {
 667         uint8_t *data = bdst->data;
 668         int m_data = bdst->m_data;   // backup data and m_data
 669         if (m_data < bsrc->m_data) { // double the capacity
 670                 m_data = bsrc->m_data; kroundup32(m_data);
 671                 data = (uint8_t*)realloc(data, m_data);
 672         }
 673         memcpy(data, bsrc->data, bsrc->data_len); // copy var-len data
 674         *bdst = *bsrc; // copy the rest
 675         // restore the backup
 676         bdst->m_data = m_data;
 677         bdst->data = data;
 678         return bdst;
 679 }
 680
 681 /*!
 682   @abstract     Duplicate an alignment
 683   @param  src   source alignment struct
 684   @return       pointer to the destination alignment struct
 685  */
 686 static inline bam1_t *bam_dup1(const bam1_t *src)
 687 {
 688         bam1_t *b;
 689         b = bam_init1();
 690         *b = *src;
 691         b->m_data = b->data_len;
 692         b->data = (uint8_t*)calloc(b->data_len, 1);
 693         memcpy(b->data, src->data, b->data_len);
 694         return b;
 695 }
 696
 697 #endif