BamAux.h

   1 // ***************************************************************************\r
   2 // BamAux.h (c) 2009 Derek Barnett, Michael Strömberg\r
   3 // Marth Lab, Department of Biology, Boston College\r
   4 // All rights reserved.\r
   5 // ---------------------------------------------------------------------------\r
   6 // Last modified: 11 December 2009 (DB)\r
   7 // ---------------------------------------------------------------------------\r
   8 // Provides the basic constants, data structures, etc. for using BAM files\r
   9 // ***************************************************************************\r
  10 \r
  11 #ifndef BAMAUX_H\r
  12 #define BAMAUX_H\r
  13 \r
  14 // C inclues\r
  15 #include <cstdlib>\r
  16 #include <cstring>\r
  17 \r
  18 // C++ includes\r
  19 #include <exception>\r
  20 #include <map>\r
  21 #include <string>\r
  22 #include <utility>\r
  23 #include <vector>\r
  24 \r
  25 // Platform-specific type definitions\r
  26 #ifndef BAMTOOLS_TYPES\r
  27 #define BAMTOOLS_TYPES\r
  28         #ifdef _MSC_VER\r
  29                         typedef char                 int8_t;\r
  30                         typedef unsigned char       uint8_t;\r
  31                         typedef short               int16_t;\r
  32                         typedef unsigned short     uint16_t;\r
  33                         typedef int                 int32_t;\r
  34                         typedef unsigned int       uint32_t;\r
  35                         typedef long long           int64_t;\r
  36                         typedef unsigned long long uint64_t;\r
  37         #else\r
  38                         #include <stdint.h>\r
  39         #endif\r
  40 #endif // BAMTOOLS_TYPES\r
  41 \r
  42 namespace BamTools {\r
  43 \r
  44 // BAM constants\r
  45 const int BAM_CORE_SIZE   = 32;\r
  46 const int BAM_CMATCH      = 0;\r
  47 const int BAM_CINS        = 1;\r
  48 const int BAM_CDEL        = 2;\r
  49 const int BAM_CREF_SKIP   = 3;\r
  50 const int BAM_CSOFT_CLIP  = 4;\r
  51 const int BAM_CHARD_CLIP  = 5;\r
  52 const int BAM_CPAD        = 6;\r
  53 const int BAM_CIGAR_SHIFT = 4;\r
  54 const int BAM_CIGAR_MASK  = ((1 << BAM_CIGAR_SHIFT) - 1);\r
  55 \r
  56 // BAM index constants\r
  57 const int MAX_BIN           = 37450;    // =(8^6-1)/7+1\r
  58 const int BAM_MIN_CHUNK_GAP = 32768;\r
  59 const int BAM_LIDX_SHIFT    = 14;\r
  60 \r
  61 // Explicit variable sizes\r
  62 const int BT_SIZEOF_INT = 4;\r
  63 \r
  64 struct CigarOp;\r
  65 \r
  66 struct BamAlignment {\r
  67 \r
  68     // Queries against alignment flag\r
  69     public:\r
  70         // Returns true if this read is a PCR duplicate (determined by external app)\r
  71         bool IsDuplicate(void) const { return ( (AlignmentFlag & DUPLICATE) != 0 ); }\r
  72         // Returns true if this read failed quality control (determined by external app)\r
  73         bool IsFailedQC(void) const { return ( (AlignmentFlag & QC_FAILED) != 0 ); }\r
  74         // Returns true if alignment is first mate on read\r
  75         bool IsFirstMate(void) const { return ( (AlignmentFlag & READ_1) != 0 ); }\r
  76         // Returns true if alignment is mapped\r
  77         bool IsMapped(void) const { return ( (AlignmentFlag & UNMAPPED) == 0 ); }\r
  78         // Returns true if alignment's mate is mapped\r
  79         bool IsMateMapped(void) const { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }\r
  80         // Returns true if alignment's mate mapped to reverse strand\r
  81         bool IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }\r
  82         // Returns true if alignment part of paired-end read\r
  83         bool IsPaired(void) const { return ( (AlignmentFlag & PAIRED) != 0 ); }\r
  84         // Returns true if this position is primary alignment (determined by external app)\r
  85         bool IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY) == 0 ); }\r
  86         // Returns true if alignment is part of read that satisfied paired-end resolution (determined by external app)\r
  87         bool IsProperPair(void) const { return ( (AlignmentFlag & PROPER_PAIR) != 0 ); }\r
  88         // Returns true if alignment mapped to reverse strand\r
  89         bool IsReverseStrand(void) const { return ( (AlignmentFlag & REVERSE) != 0 ); }\r
  90         // Returns true if alignment is second mate on read\r
  91         bool IsSecondMate(void) const { return ( (AlignmentFlag & READ_2) != 0 ); }\r
  92 \r
  93     public:\r
  94 \r
  95         // get "RG" tag data\r
  96         bool GetReadGroup(std::string& readGroup) const {\r
  97 \r
  98             if ( TagData.empty() ) { return false; }\r
  99 \r
 100             // localize the tag data\r
 101             char* pTagData = (char*)TagData.data();\r
 102             const unsigned int tagDataLen = TagData.size();\r
 103             unsigned int numBytesParsed = 0;\r
 104 \r
 105             bool foundReadGroupTag = false;\r
 106             while( numBytesParsed < tagDataLen ) {\r
 107 \r
 108                 const char* pTagType = pTagData;\r
 109                 const char* pTagStorageType = pTagData + 2;\r
 110                 pTagData       += 3;\r
 111                 numBytesParsed += 3;\r
 112 \r
 113                 // check the current tag\r
 114                 if ( std::strncmp(pTagType, "RG", 2) == 0 ) {\r
 115                     foundReadGroupTag = true;\r
 116                     break;\r
 117                 }\r
 118 \r
 119                 // get the storage class and find the next tag\r
 120                 SkipToNextTag( *pTagStorageType, pTagData, numBytesParsed );\r
 121             }\r
 122 \r
 123             // return if the read group tag was not present\r
 124             if ( !foundReadGroupTag ) { return false; }\r
 125 \r
 126             // assign the read group\r
 127             const unsigned int readGroupLen = std::strlen(pTagData);\r
 128             readGroup.resize(readGroupLen);\r
 129             std::memcpy( (char*)readGroup.data(), pTagData, readGroupLen );\r
 130             return true;\r
 131         }\r
 132 \r
 133         // get "NM" tag data - contributed by Aaron Quinlan\r
 134         bool GetEditDistance(uint8_t& editDistance) const {\r
 135 \r
 136             if ( TagData.empty() ) { return false; }\r
 137 \r
 138             // localize the tag data\r
 139             char* pTagData = (char*)TagData.data();\r
 140             const unsigned int tagDataLen = TagData.size();\r
 141             unsigned int numBytesParsed = 0;\r
 142 \r
 143             bool foundEditDistanceTag = false;\r
 144             while( numBytesParsed < tagDataLen ) {\r
 145 \r
 146                 const char* pTagType = pTagData;\r
 147                 const char* pTagStorageType = pTagData + 2;\r
 148                 pTagData       += 3;\r
 149                 numBytesParsed += 3;\r
 150 \r
 151                 // check the current tag\r
 152                 if ( strncmp(pTagType, "NM", 2) == 0 ) {\r
 153                     foundEditDistanceTag = true;\r
 154                     break;\r
 155                 }\r
 156 \r
 157                 // get the storage class and find the next tag\r
 158                 SkipToNextTag( *pTagStorageType, pTagData, numBytesParsed );\r
 159             }\r
 160             // return if the edit distance tag was not present\r
 161             if ( !foundEditDistanceTag ) { return false; }\r
 162 \r
 163             // assign the editDistance value\r
 164             memcpy(&editDistance, pTagData, 1);\r
 165             return true;\r
 166         }\r
 167 \r
 168     private:\r
 169         static void SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {\r
 170             switch(storageType) {\r
 171 \r
 172                 case 'A':\r
 173                 case 'c':\r
 174                 case 'C':\r
 175                         ++numBytesParsed;\r
 176                         ++pTagData;\r
 177                         break;\r
 178 \r
 179                 case 's':\r
 180                 case 'S':\r
 181                 case 'f':\r
 182                         numBytesParsed += 2;\r
 183                         pTagData       += 2;\r
 184                         break;\r
 185 \r
 186                 case 'i':\r
 187                 case 'I':\r
 188                         numBytesParsed += 4;\r
 189                         pTagData       += 4;\r
 190                         break;\r
 191 \r
 192                 case 'Z':\r
 193                 case 'H':\r
 194                         while(*pTagData) {\r
 195                             ++numBytesParsed;\r
 196                             ++pTagData;\r
 197                         }\r
 198                         break;\r
 199 \r
 200                 default:\r
 201                         printf("ERROR: Unknown tag storage class encountered: [%c]\n", *pTagData);\r
 202                         exit(1);\r
 203             }\r
 204         }\r
 205 \r
 206     // Data members\r
 207     public:\r
 208         std::string  Name;              // Read name\r
 209         int32_t      Length;            // Query length\r
 210         std::string  QueryBases;        // 'Original' sequence (as reported from sequencing machine)\r
 211         std::string  AlignedBases;      // 'Aligned' sequence (includes any indels, padding, clipping)\r
 212         std::string  Qualities;         // FASTQ qualities (ASCII characters, not numeric values)\r
 213         std::string  TagData;           // Tag data (accessor methods will pull the requested information out)\r
 214         int32_t      RefID;             // ID number for reference sequence\r
 215         int32_t      Position;          // Position (0-based) where alignment starts\r
 216         uint16_t     Bin;               // Bin in BAM file where this alignment resides\r
 217         uint16_t     MapQuality;        // Mapping quality score\r
 218         uint32_t     AlignmentFlag;     // Alignment bit-flag - see Is<something>() methods for available queries\r
 219         std::vector<CigarOp> CigarData; // CIGAR operations for this alignment\r
 220         int32_t      MateRefID;         // ID number for reference sequence where alignment's mate was aligned\r
 221         int32_t      MatePosition;      // Position (0-based) where alignment's mate starts\r
 222         int32_t      InsertSize;        // Mate-pair insert size\r
 223 \r
 224     // Alignment flag query constants\r
 225     private:\r
 226         enum { PAIRED        = 1,\r
 227                PROPER_PAIR   = 2,\r
 228                UNMAPPED      = 4,\r
 229                MATE_UNMAPPED = 8,\r
 230                REVERSE       = 16,\r
 231                MATE_REVERSE  = 32,\r
 232                READ_1        = 64,\r
 233                READ_2        = 128,\r
 234                SECONDARY     = 256,\r
 235                QC_FAILED     = 512,\r
 236                DUPLICATE     = 1024\r
 237              };\r
 238 };\r
 239 \r
 240 // ----------------------------------------------------------------\r
 241 // Auxiliary data structs & typedefs\r
 242 \r
 243 struct CigarOp {\r
 244     char     Type;   // Operation type (MIDNSHP)\r
 245     uint32_t Length; // Operation length (number of bases)\r
 246 };\r
 247 \r
 248 struct RefData {\r
 249     // data members\r
 250     std::string RefName;          // Name of reference sequence\r
 251     int         RefLength;        // Length of reference sequence\r
 252     bool        RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence\r
 253     // constructor\r
 254     RefData(void)\r
 255         : RefLength(0)\r
 256         , RefHasAlignments(false)\r
 257     { }\r
 258 };\r
 259 \r
 260 typedef std::vector<RefData> RefVector;\r
 261 typedef std::vector<BamAlignment> BamAlignmentVector;\r
 262 \r
 263 // ----------------------------------------------------------------\r
 264 // Indexing structs & typedefs\r
 265 \r
 266 struct Chunk {\r
 267     // data members\r
 268     uint64_t Start;\r
 269     uint64_t Stop;\r
 270     // constructor\r
 271     Chunk(const uint64_t& start = 0, const uint64_t& stop = 0)\r
 272         : Start(start)\r
 273         , Stop(stop)\r
 274     { }\r
 275 };\r
 276 \r
 277 inline\r
 278 bool ChunkLessThan(const Chunk& lhs, const Chunk& rhs) {\r
 279     return lhs.Start < rhs.Start;\r
 280 }\r
 281 \r
 282 typedef std::vector<Chunk> ChunkVector;\r
 283 typedef std::map<uint32_t, ChunkVector> BamBinMap;\r
 284 typedef std::vector<uint64_t> LinearOffsetVector;\r
 285 \r
 286 struct ReferenceIndex {\r
 287     // data members\r
 288     BamBinMap Bins;\r
 289     LinearOffsetVector Offsets;\r
 290     // constructor\r
 291     ReferenceIndex(const BamBinMap& binMap = BamBinMap(),\r
 292                    const LinearOffsetVector& offsets = LinearOffsetVector())\r
 293         : Bins(binMap)\r
 294         , Offsets(offsets)\r
 295     { }\r
 296 };\r
 297 \r
 298 typedef std::vector<ReferenceIndex> BamIndex;\r
 299 \r
 300 } // namespace BamTools\r
 301 \r
 302 #endif // BAMAUX_H\r