text: /* @include ajseqbam **********************************************************
  **
  ** AJAX BAM format sequence processing functions
  **
  ** These functions control all aspects of AJAX BAM file processing
  **
  ** @author Copyright (C) 2010 Peter Rice ported from samtools
  ** @version $Revision: 1.21 $
  ** @modified 2010-2011 Peter Rice
  ** @modified $Date: 2012/07/02 17:24:52 $ by $Author: rice $
  ** @@
  **
  ** This library is free software; you can redistribute it and/or
  ** modify it under the terms of the GNU Lesser General Public
  ** License as published by the Free Software Foundation; either
  ** version 2.1 of the License, or (at your option) any later version.
  **
  ** This library is distributed in the hope that it will be useful,
  ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  ** Lesser General Public License for more details.
  **
  ** You should have received a copy of the GNU Lesser General Public
  ** License along with this library; if not, write to the Free Software
  ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  ** MA  02110-1301,  USA.
  **
  ******************************************************************************/
  
  /* The MIT License
  **
  **   Copyright (c) 2008 Genome Research Ltd (GRL).
  **
  **   Permission is hereby granted, free of charge, to any person obtaining
  **   a copy of this software and associated documentation files (the
  **   "Software"), to deal in the Software without restriction, including
  **   without limitation the rights to use, copy, modify, merge, publish,
  **   distribute, sublicense, and/or sell copies of the Software, and to
  **   permit persons to whom the Software is furnished to do so, subject to
  **   the following conditions:
  **
  **   The above copyright notice and this permission notice shall be
  **   included in all copies or substantial portions of the Software.
  **
  **   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  **   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  **   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  **   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  **   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  **   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  **   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  **   SOFTWARE.
  */
  
  /* Contact: Heng Li <lh3@sanger.ac.uk> */
  
  /*
  ** much modified for EMBOSS by Peter Rice pmr@ebi.ac.uk May 2010
  ** lists changed to AjPList
  ** hashes changed to AjPTable
  ** strings changed to AjPStr
  ** fixed-length datatypes int32_t etc. changed to EMBOSS types
  */
  
  #ifndef AJSEQBAM_H
  #define AJSEQBAM_H
  
  /*
  **  BAM library provides I/O and various operations on manipulating files
  **  in the BAM (Binary Alignment/Mapping) or SAM (Sequence Alignment/Map)
  **  format. It now supports importing from or exporting to TAM, sorting,
  **  merging, generating pileup, and quickly retrieval of reads overlapped
  **  with a specified region.
  **
  **  copyright Genome Research Ltd.
  */
  
  /* ========================================================================= */
  /* ============================= include files ============================= */
  /* ========================================================================= */
  
  #include "ajdefine.h"
  #include "ajlist.h"
  #include "ajtable.h"
  
  #include <stdlib.h>
  #include <string.h>
  #include <stdio.h>
  
  #include "zlib.h"
  
  AJ_BEGIN_DECLS
  
  
  
  
  /* ========================================================================= */
  /* =============================== constants =============================== */
  /* ========================================================================= */
  
  
  
  
  #ifdef WIN32
  #define inline __inline
  #endif /* WIN32 */
  
  /* The read is paired in sequencing, no matter whether it is mapped in a pair */
  #define BAM_FPAIRED        1
  
  /* The read is mapped in a proper pair */
  #define BAM_FPROPER_PAIR   2
  
  /* The read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
  #define BAM_FUNMAP         4
  
  /* The mate is unmapped */
  #define BAM_FMUNMAP        8
  
  /* The read is mapped to the reverse strand */
  #define BAM_FREVERSE      16
  
  /* The mate is mapped to the reverse strand */
  #define BAM_FMREVERSE     32
  
  /* This is read1 */
  #define BAM_FREAD1        64
  
  /* This is read2 */
  #define BAM_FREAD2       128
  
  /* Not primary alignment */
  #define BAM_FSECONDARY   256
  
  /* QC failure */
  #define BAM_FQCFAIL      512
  
  /* Optical or PCR duplicate */
  #define BAM_FDUP        1024
  
  
  #define BAM_OFDEC          0
  #define BAM_OFHEX          1
  #define BAM_OFSTR          2
  
  
  /* Default mask for pileup */
  #define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP)
  
  #define BAM_CORE_SIZE   sizeof(AjOSeqBamCore)
  
  /*
  ** Describing how CIGAR operation/length is packed in a 32-bit integer.
  */
  
  #define BAM_CIGAR_SHIFT 4
  #define BAM_CIGAR_MASK  ((1 << BAM_CIGAR_SHIFT) - 1)
  
  /*
  **  CIGAR operations.
  */
  
  /* CIGAR: M match or mismatch */
  #define BAM_CMATCH      0
  
  /* CIGAR: I insertion to the reference */
  #define BAM_CINS        1
  
  /* CIGAR: D deletion from the reference */
  #define BAM_CDEL        2
  
  /* CIGAR: N skip on the reference (e.g. spliced alignment) */
  #define BAM_CREF_SKIP   3
  
  /* CIGAR: S clip on the read with clipped sequence present in qseq */
  #define BAM_CSOFT_CLIP  4
  
  /* CIGAR: H clip on the read with clipped sequence trimmed off */
  #define BAM_CHARD_CLIP  5
  
  /* CIGAR: P padding */
  #define BAM_CPAD        6
  
  /* CIGAR: = match */
  #define BAM_CEQUAL        7
  
  /* CIGAR: X mismatch */
  #define BAM_CDIFF        8
  
  
  extern const char* cigarcode;
  extern const char* bam_nt16_rev_table;
  
  
  
  
  /* ========================================================================= */
  /* ============================== public data ============================== */
  /* ========================================================================= */
  
  
  
  
  /* @data AjPSeqBamBgzf ********************************************************
  **
  ** BGZF file handling object
  **
  ** @alias AjOSeqBamBgzf
  ** @alias AjSSeqBamBgzf
  **
  ** @attr file [FILE*] File object
  ** @attr cache [AjPTable] Block cache
  ** @attr uncompressed_block [void*] Uncompressed block data
  ** @attr compressed_block [void*] Compressed block data
  ** @attr error [const char*] Error description
  ** @attr block_address [ajlong] Block offset
  ** @attr file_descriptor [int] File descriptor
  ** @attr cache_size [int] Cache size
  ** @attr uncompressed_block_size [int] Uncompressed block size
  ** @attr compressed_block_size [int] Compressed block size
  ** @attr block_length [int] Block length
  ** @attr block_offset [int] Block offset
  ** @attr open_mode [char] Open_mode 'r' or 'w'
  ** @attr owned_file [char] Boolean
  ** @attr is_uncompressed [char] Boolean
  ** @attr Padding [char[5]] Padding
  **
  ******************************************************************************/
  
  typedef struct AjSSeqBamBgzf
  {
      FILE* file;
      AjPTable cache;
      void* uncompressed_block;
      void* compressed_block;
      const char* error;
      ajlong block_address;
      int file_descriptor;
      int cache_size;
      int uncompressed_block_size;
      int compressed_block_size;
      int block_length;
      int block_offset;
      char open_mode;
      char owned_file;
      char is_uncompressed;
      char Padding[5];
  } AjOSeqBamBgzf;
  
  #define AjPSeqBamBgzf AjOSeqBamBgzf*
  
  
  
  
  #define BAM_VIRTUAL_OFFSET16
  
  
  
  
  /* #abstract BAM file handler */
  
  /* @data AjPSeqBamHeader ******************************************************
  **
  ** BAM alignment file header data
  **
  ** @attr  target_name [char**] names of the reference sequences
  ** @attr  target_len  [ajuint*] lengths of the reference sequences
  ** @attr  dict        [AjPList] header dictionary
  ** @attr  hash        [AjPTable] hash table for fast name lookup
  ** @attr  rg2lib      [AjPTable] hash table for @RG-ID -> LB lookup
  ** @attr  text        [char*] plain text
  ** @attr  n_targets   [ajint] number of reference sequences
  ** @attr  l_text      [ajint] length of the plain text in the header
  **
  ** @@
  ** discussion Field hash points to null by default. It is a private
  **  member.
  ******************************************************************************/
  
  typedef struct AjSSeqBamheader
  {
      char **target_name;
      ajuint *target_len;
      AjPList dict;
      AjPTable hash;
      AjPTable rg2lib;
      char *text;
      ajint n_targets;
      ajint l_text;
  } AjOSeqBamHeader;
  
  #define AjPSeqBamHeader AjOSeqBamHeader*
  
  
  
  
  /* @data AjPSeqBamCore ********************************************************
  **
  ** Structure for core alignment information.
  **
  ** @attr  tid     [ajint]  read ID, defined by AjPSeqBamheader
  ** @attr  pos     [ajint]  0-based leftmost coordinate
  ** @attr  bin     [ajushort]  bin calculated by ajSeqBamReg2bin()
  ** @attr  qual    [unsigned char]  mapping quality
  ** @attr  l_qname [unsigned char]  length of the query name
  ** @attr  flag    [ajushort]  bitwise flag
  ** @attr  n_cigar [ajushort]  number of CIGAR operations
  ** @attr  l_qseq  [ajint]  length of the query sequence (read)
  ** @attr  mtid    [ajint]  paired read (mate) ID
  ** @attr  mpos    [ajint]  paired read (mate) position
  ** @attr  isize   [ajint]  insert size for paired reads
  ******************************************************************************/
  
  typedef struct AjSBamSeqCore
  {
      ajint tid;
      ajint pos;
      ajushort bin;
      unsigned char qual;
      unsigned char l_qname;
      ajushort flag;
      ajushort n_cigar;
      ajint l_qseq;
      ajint mtid;
      ajint mpos;
      ajint isize;
  } AjOSeqBamCore;
  
  #define AjPSeqBamCore AjOSeqBamCore*
  
  
  
  
  /* @data AjPSeqBam ************************************************************
  **
  ** Structure for one alignment.
  **
  ** @alias AjSSeqBam
  ** @alias AjOSeqBam
  **
  ** @attr  core      [AjOSeqBamCore]  core information about the alignment
  ** @attr  data      [unsigned char*] all variable-length data, concatenated;
  **                             structure: cigar-qname-seq-qual-aux
  ** @attr  l_aux      [int]  length of auxiliary data
  ** @attr  data_len   [int]  current length of data
  ** @attr  m_data     [int]  maximum reserved size of data
  ** @attr  Padding    [int]  Padding to alignment boundary
  **
  ** @@
  ** discussion Notes:
  **
  **   1. qname is zero tailed and core.l_qname includes the tailing '\0'.
  **
  **   2. l_qseq is calculated from the total length of an alignment block
  **      on reading or from CIGAR.
  ******************************************************************************/
  
  typedef struct AjSSeqBam
  {
      AjOSeqBamCore core;
      unsigned char *data;
      int l_aux;
      int data_len;
      int m_data;
      int Padding;
  } AjOSeqBam;
  
  #define AjPSeqBam AjOSeqBam*
  
  
  
  /* ========================================================================= */
  /* =========================== public functions ============================ */
  /* ========================================================================= */
  
  
  
  
  #define MAJSEQBAMSTRAND(b) (((b)->core.flag&BAM_FREVERSE) != 0)
  #define MAJSEQBAMMSTRAND(b) (((b)->core.flag&BAM_FMREVERSE) != 0)
  
  
  
  
  /*
  **  Get the CIGAR array
  **  param  b  pointer to an alignment
  **  return    pointer to the CIGAR array
  **
  **  In the CIGAR array, each element is a 32-bit integer. The
  **  lower 4 bits gives a CIGAR operation and the higher 28 bits keep the
  **  length of a CIGAR.
  */
  #define MAJSEQBAMCIGAR(b) ((ajuint*)((b)->data + (b)->core.l_qname))
  
  
  /*
  **  Get the name of the query
  **  param  b  pointer to an alignment
  **  return    pointer to the name string, null terminated
  */
  #define MAJSEQBAMQNAME(b) ((char*)((b)->data))
  
  
  /*
  **  Get query sequence
  **  param  b  pointer to an alignment
  **  return    pointer to sequence
  **
  **  Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G,
  **  8 for T and 15 for N. Two bases are packed in one byte with the base
  **  at the higher 4 bits having smaller coordinate on the read. It is
  **  recommended to use bam1_seqi() macro to get the base.
  */
  #define MAJSEQBAMSEQ(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname)
  
  
  /*
  **  Get query quality
  **  param  b  pointer to an alignment
  **  return    pointer to quality string
  */
  #define MAJSEQBAMQUAL(b) ((b)->data + (b)->core.n_cigar*4 +             \
                            (b)->core.l_qname + ((b)->core.l_qseq + 1)/2)
  
  /*
  **  Get a base on read
  **  param  s  Query sequence returned by bam1_seq()
  **  param  i  The i-th position, 0-based
  **  return    4-bit integer representing the base.
  */
  #define MAJSEQBAMSEQI(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
  
  /*
  **  Get pointer to the list of auxiliary data
  **  param  b  pointer to an alignment
  **  return    pointer to the concatenated auxiliary data
  */
  #define MAJSEQBAMAUX(b) ((b)->data + (b)->core.n_cigar*4 + \
  			 (b)->core.l_qname + \
  		         (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2)
  
  #ifndef kroundup32
  /*
  **  Round an integer to the next closest power-2 integer.
  **  param  x  integer to be rounded (in place)
  **  x will be modified.
  */
  #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4,    \
                         (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
  #endif /* !kroundup32 */
  
  
  
  
  /*
  ** Prototype definitions
  */
  
  void ajSeqBamAuxAppend(AjPSeqBam b, const char tag[2],
  		       char type, int len, const unsigned char* data);
  
  AjPSeqBamBgzf ajSeqBamBgzfNew(FILE* file, const char* mode);
  
  AjPSeqBamBgzf ajSeqBamBgzfOpenfd(int fd, const char *mode);
  AjPSeqBamBgzf ajSeqBamBgzfOpenC(const char* path, const char *mode);
  
  ajuint ajSeqBamCalend(const AjPSeqBamCore c, const ajuint *cigar);
  
  int ajSeqBamBgzfClose(AjPSeqBamBgzf fp);
  int ajSeqBamBgzfEof(AjPSeqBamBgzf fp);
  int ajSeqBamBgzfFlush(AjPSeqBamBgzf fp);
  int ajSeqBamBgzfRead(AjPSeqBamBgzf fp, void* data, int length);
  ajlong ajSeqBamBgzfSeek(AjPSeqBamBgzf fp, ajlong pos, int where);
  AjBool ajSeqBamBgzfSetInfile(AjPSeqBamBgzf gzfile, AjPFile outf);
  AjBool ajSeqBamBgzfSetOutfile(AjPSeqBamBgzf gzfile, AjPFile outf);
  int ajSeqBamBgzfWrite(AjPSeqBamBgzf fp, const void* data, int length);
  
  void ajSeqBamDel(AjPSeqBam *Pbam);
  
  AjPSeqBamHeader ajSeqBamHeaderNew(void);
  AjPSeqBamHeader ajSeqBamHeaderNewN(ajint n);
  AjPSeqBamHeader ajSeqBamHeaderNewTextC(const char* txt);
  AjPSeqBamHeader ajSeqBamHeaderRead(AjPSeqBamBgzf gzfile);
  void ajSeqBamHeaderDel(AjPSeqBamHeader *Pheader);
  int ajSeqBamHeaderWrite(AjPSeqBamBgzf fp, const AjPSeqBamHeader header);
  AjPTable ajSeqBamHeaderGetRefseqTags(const AjPSeqBamHeader header);
  AjPTable ajSeqBamHeaderGetReadgroupTags(const AjPSeqBamHeader header);
  const char* ajSeqBamHeaderGetSortorder(const AjPSeqBamHeader header);
  void ajSeqBamHeaderSetTextC(AjPSeqBamHeader header, const char* txt);
  
  int ajSeqBamRead(AjPSeqBamBgzf fp, AjPSeqBam b);
  
  int ajSeqBamReg2bin(ajuint beg, ajuint end);
  
  int ajSeqBamValidate(const AjPSeqBamHeader header, const AjPSeqBam b);
  int ajSeqBamWrite(AjPSeqBamBgzf fp, const AjPSeqBam b);
  
  const char *ajSeqBamGetLibrary(AjPSeqBamHeader header, const AjPSeqBam b);
  
  AjPTable ajSeqBamHeaderLineParse(const char *headerLine);
  
  /*
  ** End of prototype definitions
  */
  
  
  
  
  AJ_END_DECLS
  
  #endif /* !AJSEQBAM_H */
  