BamReader.h

   1 // BamReader.h\r
   2 \r
   3 /* The MIT License\r
   4 \r
   5    Copyright (c) 2008 Genome Research Ltd (GRL).\r
   6 \r
   7    Permission is hereby granted, free of charge, to any person obtaining\r
   8    a copy of this software and associated documentation files (the\r
   9    "Software"), to deal in the Software without restriction, including\r
  10    without limitation the rights to use, copy, modify, merge, publish,\r
  11    distribute, sublicense, and/or sell copies of the Software, and to\r
  12    permit persons to whom the Software is furnished to do so, subject to\r
  13    the following conditions:\r
  14 \r
  15    The above copyright notice and this permission notice shall be\r
  16    included in all copies or substantial portions of the Software.\r
  17 \r
  18    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,\r
  19    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF\r
  20    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND\r
  21    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS\r
  22    BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN\r
  23    ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN\r
  24    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\r
  25    SOFTWARE.\r
  26 */\r
  27 \r
  28 /*\r
  29         Implementation of BAM-parsing was translated to C++ directly from Heng Li's SAMtools package \r
  30         (thus the carryover of above MIT license)\r
  31         Contact: Derek Barnett <barnetde@bc.edu>\r
  32 */\r
  33 \r
  34 // Derek Barnett\r
  35 // Marth Lab, Boston College\r
  36 // Last modified: 6 April 2009\r
  37 \r
  38 #ifndef BAMREADER_H\r
  39 #define BAMREADER_H\r
  40 \r
  41 // custom includes\r
  42 #include "BamAlignment.h"\r
  43 #include "STLUtilities.h"\r
  44 \r
  45 // C library includes\r
  46 #include <assert.h>\r
  47 #include <stdint.h>\r
  48 #include <stdio.h>\r
  49 #include <stdlib.h>\r
  50 #include <string.h>\r
  51 \r
  52 // BGZF library includes/defines\r
  53 #include "bgzf.h"\r
  54 typedef BGZF* BamFile;\r
  55 #define bam_open(f_name, mode)      bgzf_open(f_name, mode)\r
  56 #define bam_close(f_ptr)            bgzf_close(f_ptr)\r
  57 #define bam_read(f_ptr, buf, size)  bgzf_read(f_ptr, buf, size)\r
  58 #define bam_write(f_ptr, buf, size) bgzf_write(f_ptr, buf, size)\r
  59 #define bam_tell(f_ptr)             bgzf_tell(f_ptr)\r
  60 #define bam_seek(f_ptr, pos, dir)   bgzf_seek(f_ptr, pos, dir)\r
  61 \r
  62 // size of alignment data block in BAM file (bytes)\r
  63 #define BAM_CORE_SIZE 32\r
  64 \r
  65 // BAM indexing constants\r
  66 #define MAX_BIN           37450 // =(8^6-1)/7+1\r
  67 #define BAM_MIN_CHUNK_GAP 32768 \r
  68 #define BAM_LIDX_SHIFT    14\r
  69 \r
  70 // CIGAR-retrieval mask/shift constants\r
  71 #define BAM_CIGAR_SHIFT 4\r
  72 #define BAM_CIGAR_MASK  ((1 << BAM_CIGAR_SHIFT) - 1)\r
  73 \r
  74 // CIGAR-operation types\r
  75 #define BAM_CMATCH      0\r
  76 #define BAM_CINS        1\r
  77 #define BAM_CDEL        2\r
  78 #define BAM_CREF_SKIP   3\r
  79 #define BAM_CSOFT_CLIP  4\r
  80 #define BAM_CHARD_CLIP  5\r
  81 #define BAM_CPAD        6\r
  82 \r
  83 // --------------------------- //\r
  84 // Bam header info\r
  85 // --------------------------- //\r
  86 \r
  87 // --------------------------- //\r
  88 // BamIndex-related typedefs\r
  89 // --------------------------- //\r
  90 \r
  91 // offset for linear indexing\r
  92 typedef vector<uint64_t> LinearOffsetVector;\r
  93 \r
  94 // chunk boundaries\r
  95 typedef pair<uint64_t, uint64_t> ChunkPair;\r
  96 // list of chunks in a BAM bin\r
  97 typedef vector<ChunkPair> ChunkVector;\r
  98 \r
  99 // BAM bins for a reference sequence\r
 100 // replaces khash - uint32_t is key, ChunkVector is value\r
 101 typedef pair<uint32_t, ChunkVector*> BamBin;\r
 102 typedef vector<BamBin> BinVector;\r
 103 \r
 104 // each reference sequence has a BinVector and LinearOffsetVector\r
 105 typedef pair<BinVector*, LinearOffsetVector*> RefIndex;\r
 106 \r
 107 // full BamIndex defined as: \r
 108 typedef vector<RefIndex*> BamIndex;\r
 109 \r
 110 // ---------------------------------------------------------------------------//\r
 111 \r
 112 class BamReader {\r
 113         \r
 114         public:\r
 115                 // constructors\r
 116                 BamReader(const char* fileName = NULL, const char* indexFilename = NULL);\r
 117 \r
 118         public:\r
 119                 // destructor\r
 120                 ~BamReader(void);\r
 121         \r
 122         // BAM interface methods\r
 123         public:\r
 124 \r
 125                 // ----------------------- //\r
 126                 // File manipulation\r
 127                 // ----------------------- //\r
 128                 \r
 129                 // open BAM file (automatically opens index if provided)\r
 130                 bool Open(void);\r
 131                 \r
 132                 // open BAM index (allows index to be opened separately - i.e. sometime after the BAM file is opened)\r
 133                 bool OpenIndex(void);\r
 134                 \r
 135                 // close BAM file\r
 136                 bool Close(void);\r
 137                 \r
 138                 // get BAM filename\r
 139                 const char* Filename(void) const;\r
 140         \r
 141                 // set BAM filename\r
 142                 void SetFilename(const char*);\r
 143 \r
 144                 // get BAM Index filename\r
 145                 const char* IndexFilename(void) const;\r
 146                 \r
 147                 // set BAM Index filename\r
 148                 void SetIndexFilename(const char*);\r
 149 \r
 150                 // ----------------------- //\r
 151                 // Access BAM header\r
 152                 // ----------------------- //\r
 153                 \r
 154                 // return full header text\r
 155                 const string GetHeaderText(void) const;\r
 156                 \r
 157                 // --------------------------------- //\r
 158                 // Access reference sequence info\r
 159                 // --------------------------------- //\r
 160                 \r
 161                 // return number of reference sequences in BAM file\r
 162                 const int GetReferenceCount(void) const;\r
 163 \r
 164                 // return vector of RefData entries\r
 165                 const RefVector GetReferenceData(void) const;\r
 166 \r
 167                 // get refID from reference name\r
 168                 const int GetRefID(string refName) const;               \r
 169                 \r
 170                 // ----------------------------------------- //\r
 171                 // File position moving\r
 172                 // ----------------------------------------- //\r
 173 \r
 174                 // jumps to 'left' position on refID\r
 175                 // actually jumps before position, so reads that overlap 'left' are included as well\r
 176                 // 'left' defaults to reference begin if not specified\r
 177                 bool Jump(int refID, unsigned int left = 0);\r
 178 \r
 179                 // Jump to beginning of BAM file, clears any region previously set by Jump()\r
 180                 bool Rewind(void);\r
 181                 \r
 182                 // ------------------------------ //\r
 183                 // Access alignments\r
 184                 // ------------------------------ //\r
 185                 \r
 186                 // get next alignment\r
 187                 bool GetNextAlignment(BamAlignment& read);\r
 188 \r
 189                 // allow user to specifiy whether 'AlignedBases' string is calculated when alignment is loaded\r
 190                 void SetCalculateAlignedBases(bool);\r
 191 \r
 192         // internal utility methods\r
 193         private:\r
 194                 int      BinsFromRegion(int, unsigned int, uint16_t[MAX_BIN]);\r
 195                 uint32_t CalculateAlignmentEnd(const unsigned int&, const vector<CigarOp>&);\r
 196                 uint64_t GetOffset(int, unsigned int);\r
 197                 bool     IsOverlap(BamAlignment&);\r
 198                 bool     LoadHeader(void);\r
 199                 bool     LoadIndex(void);\r
 200                 bool     LoadNextAlignment(BamAlignment&);\r
 201 \r
 202         private:        \r
 203                 // main BAM reader components\r
 204                 char*      m_filename;\r
 205                 char*      m_indexFilename;\r
 206                 BamFile    m_file;\r
 207                 BamIndex*  m_index;\r
 208                 RefVector  m_references;\r
 209                 string     m_headerText;\r
 210 \r
 211                 // state flags\r
 212                 bool m_isOpen;                           // BAM file is open for processing\r
 213                 bool m_isIndexLoaded;            // BAM Index data is loaded and available for processing\r
 214                 bool m_isRegionSpecified;        // a region has been specified - specifically, a user has called Jump()\r
 215                 bool m_isCalculateAlignedBases;  // build 'AlignedBases' string when getting an alignment, otherwise skip (default = true)\r
 216 \r
 217                 // region values\r
 218                 int m_currentRefID;\r
 219                 unsigned int m_currentLeft;\r
 220 \r
 221                 // file offset of 1st read in BAM file\r
 222                 int64_t m_alignmentsBeginOffset;\r
 223 \r
 224         private:\r
 225                 // BAM character constants\r
 226                 static const char* DNA_LOOKUP;\r
 227                 static const char* CIGAR_LOOKUP;\r
 228 };\r
 229 \r
 230 #endif /* BAMREADER_H */\r