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 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.
40 @copyright Genome Research Ltd.
50 #define BAM_TRUE_OFFSET
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)
62 #define BAM_VIRTUAL_OFFSET16
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)
74 #define BAM_VIRTUAL_OFFSET16
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)
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
96 @discussion Field hash points to null by default. It is a private
102 uint32_t *target_len;
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 */
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
122 #define BAM_CORE_SIZE sizeof(bam1_core_t)
125 * Describing how CIGAR operation/length is packed in a 32-bit integer.
127 #define BAM_CIGAR_SHIFT 4
128 #define BAM_CIGAR_MASK ((1 << BAM_CIGAR_SHIFT) - 1)
133 /*! @abstract CIGAR: match */
135 /*! @abstract CIGAR: insertion to the reference */
137 /*! @abstract CIGAR: deletion from the reference */
139 /*! @abstract CIGAR: skip on the reference (e.g. spliced alignment) */
140 #define BAM_CREF_SKIP 3
141 /*! @abstract CIGAR: clip on the read with clipped sequence present in qseq */
142 #define BAM_CSOFT_CLIP 4
143 /*! @abstract CIGAR: clip on the read with clipped sequence trimmed off */
144 #define BAM_CHARD_CLIP 5
145 /*! @abstract CIGAR: padding */
149 @abstract Structure for core alignment information.
150 @field tid chromosome ID, defined by bam_header_t
151 @field pos 0-based leftmost coordinate
152 @field strand strand; 0 for forward and 1 otherwise
153 @field bin bin calculated by bam_reg2bin()
154 @field qual mapping quality
155 @field l_qname length of the query name
156 @field flag bitwise flag
157 @field n_cigar number of CIGAR operations
158 @field l_qseq length of the query sequence (read)
163 uint32_t bin:16, qual:8, l_qname:8;
164 uint32_t flag:16, n_cigar:16;
172 @abstract Structure for one alignment.
173 @field core core information about the alignment
174 @field l_aux length of auxiliary data
175 @field data_len current length of bam1_t::data
176 @field m_data maximum length of bam1_t::data
177 @field data all variable-length data, concatenated; structure: cigar-qname-seq-qual-aux
178 @field hash hash table for fast retrieval of tag-value pairs; private
182 1. qname is zero tailing and core.l_qname includes the tailing '\0'.
183 2. l_qseq is calculated from the total length of an alignment block
184 on reading or from CIGAR.
188 int l_aux, data_len, m_data;
193 #define bam1_strand(b) (((b)->core.flag&BAM_FREVERSE) != 0)
194 #define bam1_mstrand(b) (((b)->core.flag&BAM_FMREVERSE) != 0)
197 @abstract Get the CIGAR array
198 @param b pointer to an alignment
199 @return pointer to the CIGAR array
201 @discussion In the CIGAR array, each element is a 32-bit integer. The
202 lower 4 bits gives a CIGAR operation and the higher 28 bits keep the
205 #define bam1_cigar(b) ((uint32_t*)((b)->data + (b)->core.l_qname))
208 @abstract Get the name of the query
209 @param b pointer to an alignment
210 @return pointer to the name string, null terminated
212 #define bam1_qname(b) ((char*)((b)->data))
215 @abstract Get query sequence
216 @param b pointer to an alignment
217 @return pointer to sequence
219 @discussion Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G,
220 8 for T and 15 for N. Two bases are packed in one byte with the base
221 at the higher 4 bits having smaller coordinate on the read. It is
222 recommended to use bam1_seqi() macro to get the base.
224 #define bam1_seq(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname)
227 @abstract Get query quality
228 @param b pointer to an alignment
229 @return pointer to quality string
231 #define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + ((b)->core.l_qseq + 1)/2)
234 @abstract Get a base on read
235 @param s Query sequence returned by bam1_seq()
236 @param i The i-th position, 0-based
237 @return 4-bit integer representing the base.
239 #define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
242 @abstract Get query sequence and quality
243 @param b pointer to an alignment
244 @return pointer to the concatenated auxiliary data
246 #define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2)
256 @abstract Round an integer to the next closest power-2 integer.
257 @param x integer to be rounded (in place)
258 @discussion x will be modified.
260 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
264 @abstract Whether the machine is big-endian; modified only in
267 extern int bam_is_be;
269 /*! @abstract Table for converting a nucleotide character to the 4-bit encoding. */
270 extern unsigned char bam_nt16_table[256];
272 /*! @abstract Table for converting a 4-bit encoded nucleotide to a letter. */
273 extern char *bam_nt16_rev_table;
275 extern char bam_nt16_nt4_table[];
281 /*! @abstract TAM file handler */
282 typedef struct __tamFile_t *tamFile;
285 @abstract Open a TAM file, either uncompressed or compressed by gzip/zlib.
286 @param fn TAM file name
287 @return TAM file handler
289 tamFile sam_open(const char *fn);
292 @abstract Close a TAM file handler
293 @param fp TAM file handler
295 void sam_close(tamFile fp);
298 @abstract Read one alignment from a TAM file handler
299 @param fp TAM file handler
300 @param header header information (ordered names of chromosomes)
301 @param b read alignment; all members in b will be updated
302 @return 0 if successful; otherwise negative
304 int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b);
307 @abstract Read header information from a TAB-delimited list file.
308 @param fn_list file name for the list
309 @return a pointer to the header structure
311 @discussion Each line in this file consists of chromosome name and
312 the length of chromosome.
314 bam_header_t *sam_header_read2(const char *fn_list);
316 #define sam_write1(header, b) bam_view1(header, b)
319 @abstract Initialize a header structure.
320 @return the pointer to the header structure
322 @discussion This function also modifies the global variable
325 bam_header_t *bam_header_init();
328 @abstract Destroy a header structure.
329 @param header pointer to the header
331 void bam_header_destroy(bam_header_t *header);
334 @abstract Read a header structure from BAM.
335 @param fp BAM file handler, opened by bam_open()
336 @return pointer to the header structure
338 @discussion The file position indicator must be placed at the
339 beginning of the file. Upon success, the position indicator will
340 be set at the start of the first alignment.
342 bam_header_t *bam_header_read(bamFile fp);
345 @abstract Write a header structure to BAM.
346 @param fp BAM file handler
347 @param header pointer to the header structure
348 @return always 0 currently
350 int bam_header_write(bamFile fp, const bam_header_t *header);
353 @abstract Read an alignment from BAM.
354 @param fp BAM file handler
355 @param b read alignment; all members are updated.
356 @return number of bytes read from the file
358 @discussion The file position indicator must be
359 placed right before an alignment. Upon success, this function
360 will set the position indicator to the start of the next
361 alignment. This function is not affected by the machine
364 int bam_read1(bamFile fp, bam1_t *b);
367 @abstract Write an alignment to BAM.
368 @param fp BAM file handler
369 @param c pointer to the bam1_core_t structure
370 @param data_len total length of variable size data related to
372 @param data pointer to the concatenated data
373 @return number of bytes written to the file
375 @discussion This function is not affected by the machine
378 int bam_write1_core(bamFile fp, const bam1_core_t *c, int data_len, uint8_t *data);
381 @abstract Write an alignment to BAM.
382 @param fp BAM file handler
383 @param b alignment to write
384 @return number of bytes written to the file
386 @abstract It is equivalent to:
387 bam_write1_core(fp, &b->core, b->data_len, b->data)
389 int bam_write1(bamFile fp, const bam1_t *b);
392 @abstract Initiate a pointer to bam1_t struct
394 #define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t)))
397 @abstract Free the memory allocated for an alignment.
398 @param b pointer to an alignment
400 #define bam_destroy1(b) do { \
401 if ((b)->hash) bam_aux_destroy(b); free((b)->data); free(b); \
405 @abstract Print an alignment to the standard output in TAM format.
406 @param header pointer to the header structure
407 @param b alignment to print
409 void bam_view1(const bam_header_t *header, const bam1_t *b);
412 @abstract Merge multiple sorted BAM.
413 @param is_by_qname whether to sort by query name
414 @param out output BAM file name
415 @param n number of files to be merged
416 @param fn names of files to be merged
418 @discussion Padding information may NOT correctly maintained. This
419 function is NOT thread safe.
421 void bam_merge_core(int is_by_qname, const char *out, int n, char * const *fn);
424 @abstract Sort an unsorted BAM file based on the chromosome order
425 and the leftmost position of an alignment
427 @param is_by_qname whether to sort by query name
428 @param fn name of the file to be sorted
429 @param prefix prefix of the output and the temporary files; upon
430 sucessess, prefix.bam will be written.
431 @param max_mem approxiate maximum memory (very inaccurate)
433 @discussion It may create multiple temporary subalignment files
434 and then merge them by calling bam_merge_core(). This function is
437 void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t max_mem);
440 @abstract Structure for one alignment covering the pileup position.
441 @field b pointer to the alignment
442 @field qpos position of the read base at the pileup site, 0-based
443 @field indel indel length; 0 for no indel, positive for ins and negative for del
444 @field is_del 1 iff the base on the padded read is a deletion
445 @field level the level of the read in the "viewer" mode
447 @discussion See also bam_plbuf_push() and bam_lplbuf_push(). The
448 difference between the two functions is that the former does not
449 set bam_pileup1_t::level, while the later does. Level helps the
450 implementation of alignment viewers, but calculating this has some
457 uint32_t is_del:1, is_head:1, is_tail:1;
460 struct __bam_plbuf_t;
461 /*! @abstract pileup buffer */
462 typedef struct __bam_plbuf_t bam_plbuf_t;
465 @abstract Type of function to be called by bam_plbuf_push().
466 @param tid chromosome ID as is defined in the header
467 @param pos start coordinate of the alignment, 0-based
468 @param n number of elements in pl array
469 @param pl array of alignments
470 @param data user provided data
471 @discussion See also bam_plbuf_push(), bam_plbuf_init() and bam_pileup1_t.
473 typedef int (*bam_pileup_f)(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data);
475 void bam_plbuf_reset(bam_plbuf_t *buf);
478 @abstract Initialize a buffer for pileup.
479 @param func fucntion to be called by bam_pileup_core()
480 @param data user provided data
481 @return pointer to the pileup buffer
483 bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data);
486 @abstract Destroy a pileup buffer.
487 @param buf pointer to the pileup buffer
489 void bam_plbuf_destroy(bam_plbuf_t *buf);
492 @abstract Push an alignment to the pileup buffer.
493 @param b alignment to be pushed
494 @param buf pileup buffer
495 @see bam_plbuf_init()
496 @return always 0 currently
498 @discussion If all the alignments covering a particular site have
499 been collected, this function will call the user defined function
500 as is provided to bam_plbuf_init(). The coordinate of the site the
501 all the alignments will be transferred to the user defined
502 function as function parameters.
504 When all the alignments are pushed to the buffer, this function
505 needs to be called with b equal to NULL. This will flush the
506 buffer. A pileup buffer cannot be reused.
508 int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf);
511 @abstract A more convenient interface to bam_plbuf_push()
512 @param fp BAM file handler
513 @param func user defined function
514 @param func_data user provided data
516 @discussion The file position indicator must be placed right
517 before the start of an alignment. See also bam_plbuf_push().
519 int bam_pileup_file(bamFile fp, bam_pileup_f func, void *func_data);
521 struct __bam_lplbuf_t;
522 typedef struct __bam_lplbuf_t bam_lplbuf_t;
524 void bam_lplbuf_reset(bam_lplbuf_t *buf);
526 /*! @abstract bam_plbuf_init() equivalent with level calculated. */
527 bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data);
529 /*! @abstract bam_plbuf_destroy() equivalent with level calculated. */
530 void bam_lplbuf_destroy(bam_lplbuf_t *tv);
532 /*! @abstract bam_plbuf_push() equivalent with level calculated. */
533 int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf);
535 /*! @abstract bam_plbuf_file() equivalent with level calculated. */
536 int bam_lpileup_file(bamFile fp, bam_pileup_f func, void *func_data);
538 struct __bam_index_t;
539 typedef struct __bam_index_t bam_index_t;
542 @abstract Build index for a BAM file.
543 @discussion Index file "fn.bai" will be created.
544 @param fn name of the BAM file
545 @return always 0 currently
547 int bam_index_build(const char *fn);
550 @abstract Load index from file "fn.bai".
551 @param fn name of the BAM file (NOT the index file)
552 @return pointer to the index structure
554 bam_index_t *bam_index_load(const char *fn);
557 @abstract Destroy an index structure.
558 @param idx pointer to the index structure
560 void bam_index_destroy(bam_index_t *idx);
563 @abstract Type of function to be called by bam_fetch().
564 @param b the alignment
565 @param data user provided data
567 typedef int (*bam_fetch_f)(const bam1_t *b, void *data);
570 @abstract Retrieve the alignments that are overlapped with the
573 @discussion A user defined function will be called for each
574 retrieved alignment ordered by its start position.
576 @param fp BAM file handler
577 @param idx pointer to the alignment index
578 @param tid chromosome ID as is defined in the header
579 @param beg start coordinate, 0-based
580 @param end end coordinate, 0-based
581 @param data user provided data (will be transferred to func)
582 @param func user defined function
584 int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func);
587 @abstract Parse a region in the format: "chr2:100,000-200,000".
588 @discussion bam_header_t::hash will be initialized if empty.
589 @param header pointer to the header structure
590 @param str string to be parsed
591 @param ref_id the returned chromosome ID
592 @param begin the returned start coordinate
593 @param end the returned end coordinate
595 void bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end);
597 int32_t bam_aux_geti(bam1_t *b, const char tag[2], int *err);
598 float bam_aux_getf(bam1_t *b, const char tag[2], int *err);
599 char bam_aux_getc(bam1_t *b, const char tag[2], int *err);
600 char *bam_aux_getZH(bam1_t *b, const char tag[2], int *err);
601 void bam_aux_destroy(bam1_t *b);
604 @abstract Calculate the rightmost coordinate of an alignment on the
607 @param c pointer to the bam1_core_t structure
608 @param cigar the corresponding CIGAR array (from bam1_t::cigar)
609 @return the rightmost coordinate, 0-based
611 uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar);
614 @abstract Calculate the length of the query sequence from CIGAR.
615 @param c pointer to the bam1_core_t structure
616 @param cigar the corresponding CIGAR array (from bam1_t::cigar)
617 @return length of the query sequence
619 int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar);
621 int bam_segreg(int32_t pos, const bam1_core_t *c, const uint32_t *cigar, bam_segreg_t *reg);
628 @abstract Calculate the minimum bin that contains a region [beg,end).
629 @param beg start of the region, 0-based
630 @param end end of the region, 0-based
633 static inline int bam_reg2bin(uint32_t beg, uint32_t end)
636 if (beg>>14 == end>>14) return 4681 + (beg>>14);
637 if (beg>>17 == end>>17) return 585 + (beg>>17);
638 if (beg>>20 == end>>20) return 73 + (beg>>20);
639 if (beg>>23 == end>>23) return 9 + (beg>>23);
640 if (beg>>26 == end>>26) return 1 + (beg>>26);
644 static inline void bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
646 uint8_t *data = bdst->data;
647 int m_data = bdst->m_data; // backup data and m_data
648 if (m_data < bsrc->m_data) { // double the capacity
649 m_data = bsrc->m_data; kroundup32(m_data);
650 data = (uint8_t*)realloc(data, m_data);
652 memcpy(data, bsrc->data, bsrc->data_len); // copy var-len data
653 *bdst = *bsrc; // copy the rest
654 // restore the backup
655 bdst->m_data = m_data;