3 Copyright (c) 2008 Genome Research Ltd (GRL).
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:
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
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
26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
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.
40 @copyright Genome Research Ltd.
49 #define BAM_VIRTUAL_OFFSET16
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)
61 #define BAM_TRUE_OFFSET
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 */
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 hash hash table for fast name lookup
77 @field rg2lib hash table for @RG-ID -> LB lookup
78 @field l_text length of the plain text in the header
79 @field text plain text
81 @discussion Field hash points to null by default. It is a private
93 /*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
95 /*! @abstract the read is mapped in a proper pair */
96 #define BAM_FPROPER_PAIR 2
97 /*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
99 /*! @abstract the mate is unmapped */
100 #define BAM_FMUNMAP 8
101 /*! @abstract the read is mapped to the reverse strand */
102 #define BAM_FREVERSE 16
103 /*! @abstract the mate is mapped to the reverse strand */
104 #define BAM_FMREVERSE 32
105 /*! @abstract this is read1 */
106 #define BAM_FREAD1 64
107 /*! @abstract this is read2 */
108 #define BAM_FREAD2 128
109 /*! @abstract not primary alignment */
110 #define BAM_FSECONDARY 256
111 /*! @abstract QC failure */
112 #define BAM_FQCFAIL 512
113 /*! @abstract optical or PCR duplicate */
114 #define BAM_FDUP 1024
116 /*! @abstract defautl mask for pileup */
117 #define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP)
119 #define BAM_CORE_SIZE sizeof(bam1_core_t)
122 * Describing how CIGAR operation/length is packed in a 32-bit integer.
124 #define BAM_CIGAR_SHIFT 4
125 #define BAM_CIGAR_MASK ((1 << BAM_CIGAR_SHIFT) - 1)
130 /*! @abstract CIGAR: match */
132 /*! @abstract CIGAR: insertion to the reference */
134 /*! @abstract CIGAR: deletion from the reference */
136 /*! @abstract CIGAR: skip on the reference (e.g. spliced alignment) */
137 #define BAM_CREF_SKIP 3
138 /*! @abstract CIGAR: clip on the read with clipped sequence present in qseq */
139 #define BAM_CSOFT_CLIP 4
140 /*! @abstract CIGAR: clip on the read with clipped sequence trimmed off */
141 #define BAM_CHARD_CLIP 5
142 /*! @abstract CIGAR: padding */
146 @abstract Structure for core alignment information.
147 @field tid chromosome ID, defined by bam_header_t
148 @field pos 0-based leftmost coordinate
149 @field strand strand; 0 for forward and 1 otherwise
150 @field bin bin calculated by bam_reg2bin()
151 @field qual mapping quality
152 @field l_qname length of the query name
153 @field flag bitwise flag
154 @field n_cigar number of CIGAR operations
155 @field l_qseq length of the query sequence (read)
160 uint32_t bin:16, qual:8, l_qname:8;
161 uint32_t flag:16, n_cigar:16;
169 @abstract Structure for one alignment.
170 @field core core information about the alignment
171 @field l_aux length of auxiliary data
172 @field data_len current length of bam1_t::data
173 @field m_data maximum length of bam1_t::data
174 @field data all variable-length data, concatenated; structure: cigar-qname-seq-qual-aux
178 1. qname is zero tailing and core.l_qname includes the tailing '\0'.
179 2. l_qseq is calculated from the total length of an alignment block
180 on reading or from CIGAR.
184 int l_aux, data_len, m_data;
188 #define bam1_strand(b) (((b)->core.flag&BAM_FREVERSE) != 0)
189 #define bam1_mstrand(b) (((b)->core.flag&BAM_FMREVERSE) != 0)
192 @abstract Get the CIGAR array
193 @param b pointer to an alignment
194 @return pointer to the CIGAR array
196 @discussion In the CIGAR array, each element is a 32-bit integer. The
197 lower 4 bits gives a CIGAR operation and the higher 28 bits keep the
200 #define bam1_cigar(b) ((uint32_t*)((b)->data + (b)->core.l_qname))
203 @abstract Get the name of the query
204 @param b pointer to an alignment
205 @return pointer to the name string, null terminated
207 #define bam1_qname(b) ((char*)((b)->data))
210 @abstract Get query sequence
211 @param b pointer to an alignment
212 @return pointer to sequence
214 @discussion Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G,
215 8 for T and 15 for N. Two bases are packed in one byte with the base
216 at the higher 4 bits having smaller coordinate on the read. It is
217 recommended to use bam1_seqi() macro to get the base.
219 #define bam1_seq(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname)
222 @abstract Get query quality
223 @param b pointer to an alignment
224 @return pointer to quality string
226 #define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + ((b)->core.l_qseq + 1)/2)
229 @abstract Get a base on read
230 @param s Query sequence returned by bam1_seq()
231 @param i The i-th position, 0-based
232 @return 4-bit integer representing the base.
234 #define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
237 @abstract Get query sequence and quality
238 @param b pointer to an alignment
239 @return pointer to the concatenated auxiliary data
241 #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)
245 @abstract Round an integer to the next closest power-2 integer.
246 @param x integer to be rounded (in place)
247 @discussion x will be modified.
249 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
253 @abstract Whether the machine is big-endian; modified only in
256 extern int bam_is_be;
258 /*! @abstract Table for converting a nucleotide character to the 4-bit encoding. */
259 extern unsigned char bam_nt16_table[256];
261 /*! @abstract Table for converting a 4-bit encoded nucleotide to a letter. */
262 extern char *bam_nt16_rev_table;
264 extern char bam_nt16_nt4_table[];
270 /*! @abstract TAM file handler */
271 typedef struct __tamFile_t *tamFile;
274 @abstract Open a SAM file for reading, either uncompressed or compressed by gzip/zlib.
275 @param fn SAM file name
276 @return SAM file handler
278 tamFile sam_open(const char *fn);
281 @abstract Close a SAM file handler
282 @param fp SAM file handler
284 void sam_close(tamFile fp);
287 @abstract Read one alignment from a SAM file handler
288 @param fp SAM file handler
289 @param header header information (ordered names of chromosomes)
290 @param b read alignment; all members in b will be updated
291 @return 0 if successful; otherwise negative
293 int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b);
296 @abstract Read header information from a TAB-delimited list file.
297 @param fn_list file name for the list
298 @return a pointer to the header structure
300 @discussion Each line in this file consists of chromosome name and
301 the length of chromosome.
303 bam_header_t *sam_header_read2(const char *fn_list);
306 @abstract Read header from a SAM file (if present)
307 @param fp SAM file handler
308 @return pointer to header struct; 0 if no @SQ lines available
310 bam_header_t *sam_header_read(tamFile fp);
313 @abstract Parse @SQ lines a update a header struct
314 @param h pointer to the header struct to be updated
315 @return number of target sequences
317 @discussion bam_header_t::{n_targets,target_len,target_name} will
318 be destroyed in the first place.
320 int sam_header_parse(bam_header_t *h);
323 @abstract Parse @RG lines a update a header struct
324 @param h pointer to the header struct to be updated
325 @return number of @RG lines
327 @discussion bam_header_t::rg2lib will be destroyed in the first
330 int sam_header_parse_rg(bam_header_t *h);
332 #define sam_write1(header, b) bam_view1(header, b)
334 int bam_strmap_put(void *strmap, const char *rg, const char *lib);
335 const char *bam_strmap_get(const void *strmap, const char *rg);
336 void *bam_strmap_dup(const void*);
337 void *bam_strmap_init();
338 void bam_strmap_destroy(void *strmap);
341 @abstract Initialize a header structure.
342 @return the pointer to the header structure
344 @discussion This function also modifies the global variable
347 bam_header_t *bam_header_init();
350 @abstract Destroy a header structure.
351 @param header pointer to the header
353 void bam_header_destroy(bam_header_t *header);
356 @abstract Read a header structure from BAM.
357 @param fp BAM file handler, opened by bam_open()
358 @return pointer to the header structure
360 @discussion The file position indicator must be placed at the
361 beginning of the file. Upon success, the position indicator will
362 be set at the start of the first alignment.
364 bam_header_t *bam_header_read(bamFile fp);
367 @abstract Write a header structure to BAM.
368 @param fp BAM file handler
369 @param header pointer to the header structure
370 @return always 0 currently
372 int bam_header_write(bamFile fp, const bam_header_t *header);
375 @abstract Read an alignment from BAM.
376 @param fp BAM file handler
377 @param b read alignment; all members are updated.
378 @return number of bytes read from the file
380 @discussion The file position indicator must be
381 placed right before an alignment. Upon success, this function
382 will set the position indicator to the start of the next
383 alignment. This function is not affected by the machine
386 int bam_read1(bamFile fp, bam1_t *b);
389 @abstract Write an alignment to BAM.
390 @param fp BAM file handler
391 @param c pointer to the bam1_core_t structure
392 @param data_len total length of variable size data related to
394 @param data pointer to the concatenated data
395 @return number of bytes written to the file
397 @discussion This function is not affected by the machine
400 int bam_write1_core(bamFile fp, const bam1_core_t *c, int data_len, uint8_t *data);
403 @abstract Write an alignment to BAM.
404 @param fp BAM file handler
405 @param b alignment to write
406 @return number of bytes written to the file
408 @abstract It is equivalent to:
409 bam_write1_core(fp, &b->core, b->data_len, b->data)
411 int bam_write1(bamFile fp, const bam1_t *b);
414 @abstract Initiate a pointer to bam1_t struct
416 #define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t)))
419 @abstract Free the memory allocated for an alignment.
420 @param b pointer to an alignment
422 #define bam_destroy1(b) do { \
423 free((b)->data); free(b); \
427 @abstract Format a BAM record in the SAM format
428 @param header pointer to the header structure
429 @param b alignment to print
430 @return a pointer to the SAM string
432 char *bam_format1(const bam_header_t *header, const bam1_t *b);
434 char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int is_hex);
437 @abstract Structure for one alignment covering the pileup position.
438 @field b pointer to the alignment
439 @field qpos position of the read base at the pileup site, 0-based
440 @field indel indel length; 0 for no indel, positive for ins and negative for del
441 @field is_del 1 iff the base on the padded read is a deletion
442 @field level the level of the read in the "viewer" mode
444 @discussion See also bam_plbuf_push() and bam_lplbuf_push(). The
445 difference between the two functions is that the former does not
446 set bam_pileup1_t::level, while the later does. Level helps the
447 implementation of alignment viewers, but calculating this has some
454 uint32_t is_del:1, is_head:1, is_tail:1;
457 struct __bam_plbuf_t;
458 /*! @abstract pileup buffer */
459 typedef struct __bam_plbuf_t bam_plbuf_t;
461 void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask);
464 @abstract Type of function to be called by bam_plbuf_push().
465 @param tid chromosome ID as is defined in the header
466 @param pos start coordinate of the alignment, 0-based
467 @param n number of elements in pl array
468 @param pl array of alignments
469 @param data user provided data
470 @discussion See also bam_plbuf_push(), bam_plbuf_init() and bam_pileup1_t.
472 typedef int (*bam_pileup_f)(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data);
475 @abstract Reset a pileup buffer for another pileup process
476 @param buf the pileup buffer to be reset
478 void bam_plbuf_reset(bam_plbuf_t *buf);
481 @abstract Initialize a buffer for pileup.
482 @param func fucntion to be called by bam_pileup_core()
483 @param data user provided data
484 @return pointer to the pileup buffer
486 bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data);
489 @abstract Destroy a pileup buffer.
490 @param buf pointer to the pileup buffer
492 void bam_plbuf_destroy(bam_plbuf_t *buf);
495 @abstract Push an alignment to the pileup buffer.
496 @param b alignment to be pushed
497 @param buf pileup buffer
498 @see bam_plbuf_init()
499 @return always 0 currently
501 @discussion If all the alignments covering a particular site have
502 been collected, this function will call the user defined function
503 as is provided to bam_plbuf_init(). The coordinate of the site and
504 all the alignments will be transferred to the user defined
505 function as function parameters.
507 When all the alignments are pushed to the buffer, this function
508 needs to be called with b equal to NULL. This will flush the
509 buffer. A pileup buffer can only be reused when bam_plbuf_reset()
512 int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf);
514 int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data);
516 struct __bam_lplbuf_t;
517 typedef struct __bam_lplbuf_t bam_lplbuf_t;
519 void bam_lplbuf_reset(bam_lplbuf_t *buf);
521 /*! @abstract bam_plbuf_init() equivalent with level calculated. */
522 bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data);
524 /*! @abstract bam_plbuf_destroy() equivalent with level calculated. */
525 void bam_lplbuf_destroy(bam_lplbuf_t *tv);
527 /*! @abstract bam_plbuf_push() equivalent with level calculated. */
528 int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf);
530 struct __bam_index_t;
531 typedef struct __bam_index_t bam_index_t;
534 @abstract Build index for a BAM file.
535 @discussion Index file "fn.bai" will be created.
536 @param fn name of the BAM file
537 @return always 0 currently
539 int bam_index_build(const char *fn);
542 @abstract Load index from file "fn.bai".
543 @param fn name of the BAM file (NOT the index file)
544 @return pointer to the index structure
546 bam_index_t *bam_index_load(const char *fn);
549 @abstract Destroy an index structure.
550 @param idx pointer to the index structure
552 void bam_index_destroy(bam_index_t *idx);
555 @abstract Type of function to be called by bam_fetch().
556 @param b the alignment
557 @param data user provided data
559 typedef int (*bam_fetch_f)(const bam1_t *b, void *data);
562 @abstract Retrieve the alignments that are overlapped with the
565 @discussion A user defined function will be called for each
566 retrieved alignment ordered by its start position.
568 @param fp BAM file handler
569 @param idx pointer to the alignment index
570 @param tid chromosome ID as is defined in the header
571 @param beg start coordinate, 0-based
572 @param end end coordinate, 0-based
573 @param data user provided data (will be transferred to func)
574 @param func user defined function
576 int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func);
579 @abstract Parse a region in the format: "chr2:100,000-200,000".
580 @discussion bam_header_t::hash will be initialized if empty.
581 @param header pointer to the header structure
582 @param str string to be parsed
583 @param ref_id the returned chromosome ID
584 @param begin the returned start coordinate
585 @param end the returned end coordinate
586 @return 0 on success; -1 on failure
588 int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end);
591 @abstract Retrieve data of a tag
592 @param b pointer to an alignment struct
593 @param tag two-character tag to be retrieved
595 @return pointer to the type and data. The first character is the
596 type that can be 'iIsScCdfAZH'.
598 @discussion Use bam_aux2?() series to convert the returned data to
599 the corresponding type.
601 uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]);
603 int32_t bam_aux2i(const uint8_t *s);
604 float bam_aux2f(const uint8_t *s);
605 double bam_aux2d(const uint8_t *s);
606 char bam_aux2A(const uint8_t *s);
607 char *bam_aux2Z(const uint8_t *s);
609 int bam_aux_del(bam1_t *b, uint8_t *s);
610 void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);
611 uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get()
614 @abstract Calculate the rightmost coordinate of an alignment on the
617 @param c pointer to the bam1_core_t structure
618 @param cigar the corresponding CIGAR array (from bam1_t::cigar)
619 @return the rightmost coordinate, 0-based
621 uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar);
624 @abstract Calculate the length of the query sequence from CIGAR.
625 @param c pointer to the bam1_core_t structure
626 @param cigar the corresponding CIGAR array (from bam1_t::cigar)
627 @return length of the query sequence
629 int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar);
637 int bam_segreg(int32_t pos, const bam1_core_t *c, const uint32_t *cigar, bam_segreg_t *reg);
644 @abstract Calculate the minimum bin that contains a region [beg,end).
645 @param beg start of the region, 0-based
646 @param end end of the region, 0-based
649 static inline int bam_reg2bin(uint32_t beg, uint32_t end)
652 if (beg>>14 == end>>14) return 4681 + (beg>>14);
653 if (beg>>17 == end>>17) return 585 + (beg>>17);
654 if (beg>>20 == end>>20) return 73 + (beg>>20);
655 if (beg>>23 == end>>23) return 9 + (beg>>23);
656 if (beg>>26 == end>>26) return 1 + (beg>>26);
661 @abstract Copy an alignment
662 @param bdst destination alignment struct
663 @param bsrc source alignment struct
664 @return pointer to the destination alignment struct
666 static inline bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
668 uint8_t *data = bdst->data;
669 int m_data = bdst->m_data; // backup data and m_data
670 if (m_data < bsrc->m_data) { // double the capacity
671 m_data = bsrc->m_data; kroundup32(m_data);
672 data = (uint8_t*)realloc(data, m_data);
674 memcpy(data, bsrc->data, bsrc->data_len); // copy var-len data
675 *bdst = *bsrc; // copy the rest
676 // restore the backup
677 bdst->m_data = m_data;
683 @abstract Duplicate an alignment
684 @param src source alignment struct
685 @return pointer to the destination alignment struct
687 static inline bam1_t *bam_dup1(const bam1_t *src)
692 b->m_data = b->data_len;
693 b->data = (uint8_t*)calloc(b->data_len, 1);
694 memcpy(b->data, src->data, b->data_len);