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