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