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