BamAux.h

   1 BamConversionMain.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000001656\011307517135\0015027\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0#include <iostream>
   2 #include "BamReader.h"
   3 #include "BamWriter.h"
   4 using namespace BamTools;
   5 using namespace std;
   6
   7 int main(int argc, char* argv[]) {
   8
   9         if(argc != 3) {
  10                 cout << "USAGE: " << argv[0] << " <input BAM file> <output BAM file>" << endl;
  11                 exit(1);
  12         }
  13
  14         // localize our arguments
  15         const char* inputFilename  = argv[1];
  16         const char* outputFilename = argv[2];
  17
  18         // open our BAM reader
  19         BamReader reader;
  20         reader.Open(inputFilename);
  21
  22         // retrieve the SAM header text
  23         string samHeader = reader.GetHeaderText();
  24
  25         // retrieve the reference sequence vector
  26         RefVector referenceSequences = reader.GetReferenceData();
  27
  28         // open the BAM writer
  29         BamWriter writer;
  30         writer.Open(outputFilename, samHeader, referenceSequences);
  31
  32         // copy all of the reads from the input file to the output file
  33         BamAlignment al;
  34         while(reader.GetNextAlignment(al)) writer.SaveAlignment(al);
  35
  36         // close our files
  37         reader.Close();
  38         writer.Close();
  39
  40         return 0;
  41 }
  42 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamDumpMain.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000003172\011307517135\0013602\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
  43 // BamDump.cpp (c) 2009 Derek Barnett\r
  44 // Marth Lab, Department of Biology, Boston College\r
  45 // All rights reserved.\r
  46 // ---------------------------------------------------------------------------\r
  47 // Last modified: 15 July 2009 (DB)\r
  48 // ---------------------------------------------------------------------------\r
  49 // Spits out all alignments in BAM file.\r
  50 //\r
  51 // N.B. - Could result in HUGE text file. This is mostly a debugging tool\r
  52 // for small files.  You have been warned.\r
  53 // ***************************************************************************\r
  54 \r
  55 // Std C/C++ includes\r
  56 #include <cstdlib>\r
  57 #include <iostream>\r
  58 #include <string>\r
  59 using namespace std;\r
  60 \r
  61 // BamTools includes\r
  62 #include "BamReader.h"\r
  63 #include "BamWriter.h"\r
  64 using namespace BamTools;\r
  65 \r
  66 void PrintAlignment(const BamAlignment&);\r
  67 \r
  68 int main(int argc, char* argv[]) {\r
  69 \r
  70         // validate argument count\r
  71         if( argc != 2 ) {\r
  72                 cerr << "USAGE: " << argv[0] << " <input BAM file> " << endl;\r
  73                 exit(1);\r
  74         }\r
  75 \r
  76         string filename = argv[1];\r
  77         cout << "Printing alignments from file: " << filename << endl;\r
  78         \r
  79         BamReader reader;\r
  80         reader.Open(filename);\r
  81         \r
  82         BamAlignment bAlignment;\r
  83         while (reader.GetNextAlignment(bAlignment)) {\r
  84                 PrintAlignment(bAlignment);\r
  85         }\r
  86 \r
  87         reader.Close();\r
  88         return 0;\r
  89 }\r
  90         \r
  91 // Spit out basic BamAlignment data \r
  92 void PrintAlignment(const BamAlignment& alignment) {\r
  93         cout << "---------------------------------" << endl;\r
  94         cout << "Name: "       << alignment.Name << endl;\r
  95         cout << "Aligned to: " << alignment.RefID;\r
  96         cout << ":"            << alignment.Position << endl;\r
  97 }\r
  98 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamReader.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000107374\011307521263\0013300\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
  99 // BamReader.cpp (c) 2009 Derek Barnett, Michael Strömberg\r
 100 // Marth Lab, Department of Biology, Boston College\r
 101 // All rights reserved.\r
 102 // ---------------------------------------------------------------------------\r
 103 // Last modified: 8 December 2009 (DB)\r
 104 // ---------------------------------------------------------------------------\r
 105 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
 106 // Institute.\r
 107 // ---------------------------------------------------------------------------\r
 108 // Provides the basic functionality for reading BAM files\r
 109 // ***************************************************************************\r
 110 \r
 111 // C++ includes\r
 112 #include <algorithm>\r
 113 #include <iterator>\r
 114 #include <string>\r
 115 #include <vector>\r
 116 \r
 117 // BamTools includes\r
 118 #include "BGZF.h"\r
 119 #include "BamReader.h"\r
 120 using namespace BamTools;\r
 121 using namespace std;\r
 122 \r
 123 struct BamReader::BamReaderPrivate {\r
 124 \r
 125     // -------------------------------\r
 126     // data members\r
 127     // -------------------------------\r
 128 \r
 129     // general data\r
 130     BgzfData  mBGZF;\r
 131     string    HeaderText;\r
 132     BamIndex  Index;\r
 133     RefVector References;\r
 134     bool      IsIndexLoaded;\r
 135     int64_t   AlignmentsBeginOffset;\r
 136     string    Filename;\r
 137     string    IndexFilename;\r
 138 \r
 139     // user-specified region values\r
 140     bool IsRegionSpecified;\r
 141     int  CurrentRefID;\r
 142     int  CurrentLeft;\r
 143 \r
 144     // BAM character constants\r
 145     const char* DNA_LOOKUP;\r
 146     const char* CIGAR_LOOKUP;\r
 147 \r
 148     // -------------------------------\r
 149     // constructor & destructor\r
 150     // -------------------------------\r
 151     BamReaderPrivate(void);\r
 152     ~BamReaderPrivate(void);\r
 153 \r
 154     // -------------------------------\r
 155     // "public" interface\r
 156     // -------------------------------\r
 157 \r
 158     // flie operations\r
 159     void Close(void);\r
 160     bool Jump(int refID, int position = 0);\r
 161     void Open(const string& filename, const string& indexFilename = "");\r
 162     bool Rewind(void);\r
 163 \r
 164     // access alignment data\r
 165     bool GetNextAlignment(BamAlignment& bAlignment);\r
 166 \r
 167     // access auxiliary data\r
 168     const string GetHeaderText(void) const;\r
 169     const int GetReferenceCount(void) const;\r
 170     const RefVector GetReferenceData(void) const;\r
 171     const int GetReferenceID(const string& refName) const;\r
 172 \r
 173     // index operations\r
 174     bool CreateIndex(void);\r
 175 \r
 176     // -------------------------------\r
 177     // internal methods\r
 178     // -------------------------------\r
 179 \r
 180     // *** reading alignments and auxiliary data *** //\r
 181 \r
 182     // calculate bins that overlap region ( left to reference end for now )\r
 183     int BinsFromRegion(int refID, int left, uint16_t[MAX_BIN]);\r
 184     // calculates alignment end position based on starting position and provided CIGAR operations\r
 185     int CalculateAlignmentEnd(const int& position, const std::vector<CigarOp>& cigarData);\r
 186     // calculate file offset for first alignment chunk overlapping 'left'\r
 187     int64_t GetOffset(int refID, int left);\r
 188     // checks to see if alignment overlaps current region\r
 189     bool IsOverlap(BamAlignment& bAlignment);\r
 190     // retrieves header text from BAM file\r
 191     void LoadHeaderData(void);\r
 192     // retrieves BAM alignment under file pointer\r
 193     bool LoadNextAlignment(BamAlignment& bAlignment);\r
 194     // builds reference data structure from BAM file\r
 195     void LoadReferenceData(void);\r
 196 \r
 197     // *** index file handling *** //\r
 198 \r
 199     // calculates index for BAM file\r
 200     bool BuildIndex(void);\r
 201     // clear out inernal index data structure\r
 202     void ClearIndex(void);\r
 203     // saves BAM bin entry for index\r
 204     void InsertBinEntry(BamBinMap& binMap, const uint32_t& saveBin, const uint64_t& saveOffset, const uint64_t& lastOffset);\r
 205     // saves linear offset entry for index\r
 206     void InsertLinearOffset(LinearOffsetVector& offsets, const BamAlignment& bAlignment, const uint64_t& lastOffset);\r
 207     // loads index from BAM index file\r
 208     bool LoadIndex(void);\r
 209     // simplifies index by merging 'chunks'\r
 210     void MergeChunks(void);\r
 211     // round-up 32-bit integer to next power-of-2\r
 212     void Roundup32(int& value);\r
 213     // saves index to BAM index file\r
 214     bool WriteIndex(void);\r
 215 };\r
 216 \r
 217 // -----------------------------------------------------\r
 218 // BamReader implementation (wrapper around BRPrivate)\r
 219 // -----------------------------------------------------\r
 220 \r
 221 // constructor\r
 222 BamReader::BamReader(void) {\r
 223     d = new BamReaderPrivate;\r
 224 }\r
 225 \r
 226 // destructor\r
 227 BamReader::~BamReader(void) {\r
 228     delete d;\r
 229     d = 0;\r
 230 }\r
 231 \r
 232 // file operations\r
 233 void BamReader::Close(void) { d->Close(); }\r
 234 bool BamReader::Jump(int refID, int position) { return d->Jump(refID, position); }\r
 235 void BamReader::Open(const string& filename, const string& indexFilename) { d->Open(filename, indexFilename); }\r
 236 bool BamReader::Rewind(void) { return d->Rewind(); }\r
 237 \r
 238 // access alignment data\r
 239 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }\r
 240 \r
 241 // access auxiliary data\r
 242 const string    BamReader::GetHeaderText(void) const { return d->HeaderText; }\r
 243 const int       BamReader::GetReferenceCount(void) const { return d->References.size(); }\r
 244 const RefVector BamReader::GetReferenceData(void) const { return d->References; }\r
 245 const int       BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }\r
 246 \r
 247 // index operations\r
 248 bool BamReader::CreateIndex(void) { return d->CreateIndex(); }\r
 249 \r
 250 // -----------------------------------------------------\r
 251 // BamReaderPrivate implementation\r
 252 // -----------------------------------------------------\r
 253 \r
 254 // constructor\r
 255 BamReader::BamReaderPrivate::BamReaderPrivate(void)\r
 256     : IsIndexLoaded(false)\r
 257     , AlignmentsBeginOffset(0)\r
 258     , IsRegionSpecified(false)\r
 259     , CurrentRefID(0)\r
 260     , CurrentLeft(0)\r
 261     , DNA_LOOKUP("=ACMGRSVTWYHKDBN")\r
 262     , CIGAR_LOOKUP("MIDNSHP")\r
 263 { }\r
 264 \r
 265 // destructor\r
 266 BamReader::BamReaderPrivate::~BamReaderPrivate(void) {\r
 267     Close();\r
 268 }\r
 269 \r
 270 // calculate bins that overlap region ( left to reference end for now )\r
 271 int BamReader::BamReaderPrivate::BinsFromRegion(int refID, int left, uint16_t list[MAX_BIN]) {\r
 272 \r
 273     // get region boundaries\r
 274     uint32_t begin = (unsigned int)left;\r
 275     uint32_t end   = (unsigned int)References.at(refID).RefLength - 1;\r
 276 \r
 277     // initialize list, bin '0' always a valid bin\r
 278     int i = 0;\r
 279     list[i++] = 0;\r
 280 \r
 281     // get rest of bins that contain this region\r
 282     unsigned int k;\r
 283     for (k =    1 + (begin>>26); k <=    1 + (end>>26); ++k) { list[i++] = k; }\r
 284     for (k =    9 + (begin>>23); k <=    9 + (end>>23); ++k) { list[i++] = k; }\r
 285     for (k =   73 + (begin>>20); k <=   73 + (end>>20); ++k) { list[i++] = k; }\r
 286     for (k =  585 + (begin>>17); k <=  585 + (end>>17); ++k) { list[i++] = k; }\r
 287     for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { list[i++] = k; }\r
 288 \r
 289     // return number of bins stored\r
 290     return i;\r
 291 }\r
 292 \r
 293 // populates BAM index data structure from BAM file data\r
 294 bool BamReader::BamReaderPrivate::BuildIndex(void) {\r
 295 \r
 296     // check to be sure file is open\r
 297     if (!mBGZF.IsOpen) { return false; }\r
 298 \r
 299     // move file pointer to beginning of alignments\r
 300     Rewind();\r
 301 \r
 302     // get reference count, reserve index space\r
 303     int numReferences = References.size();\r
 304     for ( int i = 0; i < numReferences; ++i ) {\r
 305         Index.push_back(ReferenceIndex());\r
 306     }\r
 307 \r
 308     // sets default constant for bin, ID, offset, coordinate variables\r
 309     const uint32_t defaultValue = 0xffffffffu;\r
 310 \r
 311     // bin data\r
 312     uint32_t saveBin(defaultValue);\r
 313     uint32_t lastBin(defaultValue);\r
 314 \r
 315     // reference ID data\r
 316     int32_t saveRefID(defaultValue);\r
 317     int32_t lastRefID(defaultValue);\r
 318 \r
 319     // offset data\r
 320     uint64_t saveOffset = mBGZF.Tell();\r
 321     uint64_t lastOffset = saveOffset;\r
 322 \r
 323     // coordinate data\r
 324     int32_t lastCoordinate = defaultValue;\r
 325 \r
 326     BamAlignment bAlignment;\r
 327     while( GetNextAlignment(bAlignment) ) {\r
 328 \r
 329         // change of chromosome, save ID, reset bin\r
 330         if ( lastRefID != bAlignment.RefID ) {\r
 331             lastRefID = bAlignment.RefID;\r
 332             lastBin   = defaultValue;\r
 333         }\r
 334 \r
 335         // if lastCoordinate greater than BAM position - file not sorted properly\r
 336         else if ( lastCoordinate > bAlignment.Position ) {\r
 337             printf("BAM file not properly sorted:\n");\r
 338             printf("Alignment %s : %d > %d on reference (id = %d)", bAlignment.Name.c_str(), lastCoordinate, bAlignment.Position, bAlignment.RefID);\r
 339             exit(1);\r
 340         }\r
 341 \r
 342         // if valid reference && BAM bin spans some minimum cutoff (smaller bin ids span larger regions)\r
 343         if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) {\r
 344 \r
 345             // save linear offset entry (matched to BAM entry refID)\r
 346             ReferenceIndex& refIndex = Index.at(bAlignment.RefID);\r
 347             LinearOffsetVector& offsets = refIndex.Offsets;\r
 348             InsertLinearOffset(offsets, bAlignment, lastOffset);\r
 349         }\r
 350 \r
 351         // if current BamAlignment bin != lastBin, "then possibly write the binning index"\r
 352         if ( bAlignment.Bin != lastBin ) {\r
 353 \r
 354             // if not first time through\r
 355             if ( saveBin != defaultValue ) {\r
 356 \r
 357                 // save Bam bin entry\r
 358                 ReferenceIndex& refIndex = Index.at(saveRefID);\r
 359                 BamBinMap& binMap = refIndex.Bins;\r
 360                 InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);\r
 361             }\r
 362 \r
 363             // update saveOffset\r
 364             saveOffset = lastOffset;\r
 365 \r
 366             // update bin values\r
 367             saveBin = bAlignment.Bin;\r
 368             lastBin = bAlignment.Bin;\r
 369 \r
 370             // update saveRefID\r
 371             saveRefID = bAlignment.RefID;\r
 372 \r
 373             // if invalid RefID, break out (why?)\r
 374             if ( saveRefID < 0 ) { break; }\r
 375         }\r
 376 \r
 377         // make sure that current file pointer is beyond lastOffset\r
 378         if ( mBGZF.Tell() <= (int64_t)lastOffset  ) {\r
 379             printf("Error in BGZF offsets.\n");\r
 380             exit(1);\r
 381         }\r
 382 \r
 383         // update lastOffset\r
 384         lastOffset = mBGZF.Tell();\r
 385 \r
 386         // update lastCoordinate\r
 387         lastCoordinate = bAlignment.Position;\r
 388     }\r
 389 \r
 390     // save any leftover BAM data (as long as refID is valid)\r
 391     if ( saveRefID >= 0 ) {\r
 392         // save Bam bin entry\r
 393         ReferenceIndex& refIndex = Index.at(saveRefID);\r
 394         BamBinMap& binMap = refIndex.Bins;\r
 395         InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);\r
 396     }\r
 397 \r
 398     // simplify index by merging chunks\r
 399     MergeChunks();\r
 400 \r
 401     // iterate over references\r
 402     BamIndex::iterator indexIter = Index.begin();\r
 403     BamIndex::iterator indexEnd  = Index.end();\r
 404     for ( int i = 0; indexIter != indexEnd; ++indexIter, ++i ) {\r
 405 \r
 406         // get reference index data\r
 407         ReferenceIndex& refIndex = (*indexIter);\r
 408         BamBinMap& binMap = refIndex.Bins;\r
 409         LinearOffsetVector& offsets = refIndex.Offsets;\r
 410 \r
 411         // store whether reference has alignments or no\r
 412         References[i].RefHasAlignments = ( binMap.size() > 0 );\r
 413 \r
 414         // sort linear offsets\r
 415         sort(offsets.begin(), offsets.end());\r
 416     }\r
 417 \r
 418 \r
 419     // rewind file pointer to beginning of alignments, return success/fail\r
 420     return Rewind();\r
 421 }\r
 422 \r
 423 // calculates alignment end position based on starting position and provided CIGAR operations\r
 424 int BamReader::BamReaderPrivate::CalculateAlignmentEnd(const int& position, const vector<CigarOp>& cigarData) {\r
 425 \r
 426     // initialize alignment end to starting position\r
 427     int alignEnd = position;\r
 428 \r
 429     // iterate over cigar operations\r
 430     vector<CigarOp>::const_iterator cigarIter = cigarData.begin();\r
 431     vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();\r
 432     for ( ; cigarIter != cigarEnd; ++cigarIter) {\r
 433         char cigarType = (*cigarIter).Type;\r
 434         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {\r
 435             alignEnd += (*cigarIter).Length;\r
 436         }\r
 437     }\r
 438     return alignEnd;\r
 439 }\r
 440 \r
 441 \r
 442 // clear index data structure\r
 443 void BamReader::BamReaderPrivate::ClearIndex(void) {\r
 444     Index.clear(); // sufficient ??\r
 445 }\r
 446 \r
 447 // closes the BAM file\r
 448 void BamReader::BamReaderPrivate::Close(void) {\r
 449     mBGZF.Close();\r
 450     ClearIndex();\r
 451     HeaderText.clear();\r
 452     IsRegionSpecified = false;\r
 453 }\r
 454 \r
 455 // create BAM index from BAM file (keep structure in memory) and write to default index output file\r
 456 bool BamReader::BamReaderPrivate::CreateIndex(void) {\r
 457 \r
 458     // clear out index\r
 459     ClearIndex();\r
 460 \r
 461         // build (& save) index from BAM file\r
 462     bool ok = true;\r
 463     ok &= BuildIndex();\r
 464     ok &= WriteIndex();\r
 465 \r
 466         // return success/fail\r
 467     return ok;\r
 468 }\r
 469 \r
 470 // returns RefID for given RefName (returns References.size() if not found)\r
 471 const int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {\r
 472 \r
 473     // retrieve names from reference data\r
 474     vector<string> refNames;\r
 475     RefVector::const_iterator refIter = References.begin();\r
 476     RefVector::const_iterator refEnd  = References.end();\r
 477     for ( ; refIter != refEnd; ++refIter) {\r
 478         refNames.push_back( (*refIter).RefName );\r
 479     }\r
 480 \r
 481     // return 'index-of' refName ( if not found, returns refNames.size() )\r
 482     return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));\r
 483 }\r
 484 \r
 485 // get next alignment (from specified region, if given)\r
 486 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {\r
 487 \r
 488     // if valid alignment available\r
 489     if ( LoadNextAlignment(bAlignment) ) {\r
 490 \r
 491         // if region not specified, return success\r
 492         if ( !IsRegionSpecified ) { return true; }\r
 493 \r
 494         // load next alignment until region overlap is found\r
 495         while ( !IsOverlap(bAlignment) ) {\r
 496             // if no valid alignment available (likely EOF) return failure\r
 497             if ( !LoadNextAlignment(bAlignment) ) { return false; }\r
 498         }\r
 499 \r
 500         // return success (alignment found that overlaps region)\r
 501         return true;\r
 502     }\r
 503 \r
 504     // no valid alignment\r
 505     else { return false; }\r
 506 }\r
 507 \r
 508 // calculate closest indexed file offset for region specified\r
 509 int64_t BamReader::BamReaderPrivate::GetOffset(int refID, int left) {\r
 510 \r
 511     // calculate which bins overlap this region\r
 512     uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);\r
 513     int numBins = BinsFromRegion(refID, left, bins);\r
 514 \r
 515     // get bins for this reference\r
 516     const ReferenceIndex& refIndex = Index.at(refID);\r
 517     const BamBinMap& binMap        = refIndex.Bins;\r
 518 \r
 519     // get minimum offset to consider\r
 520     const LinearOffsetVector& offsets = refIndex.Offsets;\r
 521     uint64_t minOffset = ( (unsigned int)(left>>BAM_LIDX_SHIFT) >= offsets.size() ) ? 0 : offsets.at(left>>BAM_LIDX_SHIFT);\r
 522 \r
 523     // store offsets to beginning of alignment 'chunks'\r
 524     std::vector<int64_t> chunkStarts;\r
 525 \r
 526     // store all alignment 'chunk' starts for bins in this region\r
 527     for (int i = 0; i < numBins; ++i ) {\r
 528         uint16_t binKey = bins[i];\r
 529 \r
 530         map<uint32_t, ChunkVector>::const_iterator binIter = binMap.find(binKey);\r
 531         if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) {\r
 532 \r
 533             const ChunkVector& chunks = (*binIter).second;\r
 534             std::vector<Chunk>::const_iterator chunksIter = chunks.begin();\r
 535             std::vector<Chunk>::const_iterator chunksEnd  = chunks.end();\r
 536             for ( ; chunksIter != chunksEnd; ++chunksIter) {\r
 537                 const Chunk& chunk = (*chunksIter);\r
 538                 if ( chunk.Stop > minOffset ) {\r
 539                     chunkStarts.push_back( chunk.Start );\r
 540                 }\r
 541             }\r
 542         }\r
 543     }\r
 544 \r
 545     // clean up memory\r
 546     free(bins);\r
 547 \r
 548     // if no alignments found, else return smallest offset for alignment starts\r
 549     if ( chunkStarts.size() == 0 ) { return -1; }\r
 550     else { return *min_element(chunkStarts.begin(), chunkStarts.end()); }\r
 551 }\r
 552 \r
 553 // saves BAM bin entry for index\r
 554 void BamReader::BamReaderPrivate::InsertBinEntry(BamBinMap&      binMap,\r
 555                                                  const uint32_t& saveBin,\r
 556                                                  const uint64_t& saveOffset,\r
 557                                                  const uint64_t& lastOffset)\r
 558 {\r
 559     // look up saveBin\r
 560     BamBinMap::iterator binIter = binMap.find(saveBin);\r
 561 \r
 562     // create new chunk\r
 563     Chunk newChunk(saveOffset, lastOffset);\r
 564 \r
 565     // if entry doesn't exist\r
 566     if ( binIter == binMap.end() ) {\r
 567         ChunkVector newChunks;\r
 568         newChunks.push_back(newChunk);\r
 569         binMap.insert( pair<uint32_t, ChunkVector>(saveBin, newChunks));\r
 570     }\r
 571 \r
 572     // otherwise\r
 573     else {\r
 574         ChunkVector& binChunks = (*binIter).second;\r
 575         binChunks.push_back( newChunk );\r
 576     }\r
 577 }\r
 578 \r
 579 // saves linear offset entry for index\r
 580 void BamReader::BamReaderPrivate::InsertLinearOffset(LinearOffsetVector& offsets,\r
 581                                                      const BamAlignment& bAlignment,\r
 582                                                      const uint64_t&     lastOffset)\r
 583 {\r
 584     // get converted offsets\r
 585     int beginOffset = bAlignment.Position >> BAM_LIDX_SHIFT;\r
 586     int endOffset   = ( CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData) - 1) >> BAM_LIDX_SHIFT;\r
 587 \r
 588     // resize vector if necessary\r
 589     int oldSize = offsets.size();\r
 590     int newSize = endOffset + 1;\r
 591     if ( oldSize < newSize ) {        \r
 592         Roundup32(newSize);\r
 593         offsets.resize(newSize, 0);\r
 594     }\r
 595 \r
 596     // store offset\r
 597     for(int i = beginOffset + 1; i <= endOffset ; ++i) {\r
 598         if ( offsets[i] == 0) {\r
 599             offsets[i] = lastOffset;\r
 600         }\r
 601     }\r
 602 }\r
 603 \r
 604 // returns whether alignment overlaps currently specified region (refID, leftBound)\r
 605 bool BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {\r
 606 \r
 607     // if on different reference sequence, quit\r
 608     if ( bAlignment.RefID != CurrentRefID ) { return false; }\r
 609 \r
 610     // read starts after left boundary\r
 611     if ( bAlignment.Position >= CurrentLeft) { return true; }\r
 612 \r
 613     // return whether alignment end overlaps left boundary\r
 614     return ( CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData) >= CurrentLeft );\r
 615 }\r
 616 \r
 617 // jumps to specified region(refID, leftBound) in BAM file, returns success/fail\r
 618 bool BamReader::BamReaderPrivate::Jump(int refID, int position) {\r
 619 \r
 620     // if data exists for this reference and position is valid    \r
 621     if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {\r
 622 \r
 623                 // set current region\r
 624         CurrentRefID = refID;\r
 625         CurrentLeft  = position;\r
 626         IsRegionSpecified = true;\r
 627 \r
 628                 // calculate offset\r
 629         int64_t offset = GetOffset(CurrentRefID, CurrentLeft);\r
 630 \r
 631                 // if in valid offset, return failure\r
 632         if ( offset == -1 ) { return false; }\r
 633 \r
 634                 // otherwise return success of seek operation\r
 635         else { return mBGZF.Seek(offset); }\r
 636     }\r
 637 \r
 638         // invalid jump request parameters, return failure\r
 639     return false;\r
 640 }\r
 641 \r
 642 // load BAM header data\r
 643 void BamReader::BamReaderPrivate::LoadHeaderData(void) {\r
 644 \r
 645     // check to see if proper BAM header\r
 646     char buffer[4];\r
 647     if (mBGZF.Read(buffer, 4) != 4) {\r
 648         printf("Could not read header type\n");\r
 649         exit(1);\r
 650     }\r
 651 \r
 652     if (strncmp(buffer, "BAM\001", 4)) {\r
 653         printf("wrong header type!\n");\r
 654         exit(1);\r
 655     }\r
 656 \r
 657     // get BAM header text length\r
 658     mBGZF.Read(buffer, 4);\r
 659     const unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);\r
 660 \r
 661     // get BAM header text\r
 662     char* headerText = (char*)calloc(headerTextLength + 1, 1);\r
 663     mBGZF.Read(headerText, headerTextLength);\r
 664     HeaderText = (string)((const char*)headerText);\r
 665 \r
 666     // clean up calloc-ed temp variable\r
 667     free(headerText);\r
 668 }\r
 669 \r
 670 // load existing index data from BAM index file (".bai"), return success/fail\r
 671 bool BamReader::BamReaderPrivate::LoadIndex(void) {\r
 672 \r
 673     // clear out index data\r
 674     ClearIndex();\r
 675 \r
 676     // skip if index file empty\r
 677     if ( IndexFilename.empty() ) { return false; }\r
 678 \r
 679     // open index file, abort on error\r
 680     FILE* indexStream = fopen(IndexFilename.c_str(), "rb");\r
 681     if(!indexStream) {\r
 682         printf("ERROR: Unable to open the BAM index file %s for reading.\n", IndexFilename.c_str() );\r
 683         return false;\r
 684     }\r
 685 \r
 686     // see if index is valid BAM index\r
 687     char magic[4];\r
 688     fread(magic, 1, 4, indexStream);\r
 689     if (strncmp(magic, "BAI\1", 4)) {\r
 690         printf("Problem with index file - invalid format.\n");\r
 691         fclose(indexStream);\r
 692         return false;\r
 693     }\r
 694 \r
 695     // get number of reference sequences\r
 696     uint32_t numRefSeqs;\r
 697     fread(&numRefSeqs, 4, 1, indexStream);\r
 698 \r
 699     // intialize space for BamIndex data structure\r
 700     Index.reserve(numRefSeqs);\r
 701 \r
 702     // iterate over reference sequences\r
 703     for (unsigned int i = 0; i < numRefSeqs; ++i) {\r
 704 \r
 705         // get number of bins for this reference sequence\r
 706         int32_t numBins;\r
 707         fread(&numBins, 4, 1, indexStream);\r
 708 \r
 709         if (numBins > 0) {\r
 710             RefData& refEntry = References[i];\r
 711             refEntry.RefHasAlignments = true;\r
 712         }\r
 713 \r
 714         // intialize BinVector\r
 715         BamBinMap binMap;\r
 716 \r
 717         // iterate over bins for that reference sequence\r
 718         for (int j = 0; j < numBins; ++j) {\r
 719 \r
 720             // get binID\r
 721             uint32_t binID;\r
 722             fread(&binID, 4, 1, indexStream);\r
 723 \r
 724             // get number of regionChunks in this bin\r
 725             uint32_t numChunks;\r
 726             fread(&numChunks, 4, 1, indexStream);\r
 727 \r
 728             // intialize ChunkVector\r
 729             ChunkVector regionChunks;\r
 730             regionChunks.reserve(numChunks);\r
 731 \r
 732             // iterate over regionChunks in this bin\r
 733             for (unsigned int k = 0; k < numChunks; ++k) {\r
 734 \r
 735                 // get chunk boundaries (left, right)\r
 736                 uint64_t left;\r
 737                 uint64_t right;\r
 738                 fread(&left, 8, 1, indexStream);\r
 739                 fread(&right, 8, 1, indexStream);\r
 740 \r
 741                 // save ChunkPair\r
 742                 regionChunks.push_back( Chunk(left, right) );\r
 743             }\r
 744 \r
 745             // sort chunks for this bin\r
 746             sort( regionChunks.begin(), regionChunks.end(), ChunkLessThan );\r
 747 \r
 748             // save binID, chunkVector for this bin\r
 749             binMap.insert( pair<uint32_t, ChunkVector>(binID, regionChunks) );\r
 750         }\r
 751 \r
 752         // load linear index for this reference sequence\r
 753 \r
 754         // get number of linear offsets\r
 755         int32_t numLinearOffsets;\r
 756         fread(&numLinearOffsets, 4, 1, indexStream);\r
 757 \r
 758         // intialize LinearOffsetVector\r
 759         LinearOffsetVector offsets;\r
 760         offsets.reserve(numLinearOffsets);\r
 761 \r
 762         // iterate over linear offsets for this reference sequeence\r
 763         uint64_t linearOffset;\r
 764         for (int j = 0; j < numLinearOffsets; ++j) {\r
 765             // read a linear offset & store\r
 766             fread(&linearOffset, 8, 1, indexStream);\r
 767             offsets.push_back(linearOffset);\r
 768         }\r
 769 \r
 770         // sort linear offsets\r
 771         sort( offsets.begin(), offsets.end() );\r
 772 \r
 773         // store index data for that reference sequence\r
 774         Index.push_back( ReferenceIndex(binMap, offsets) );\r
 775     }\r
 776 \r
 777     // close index file (.bai) and return\r
 778     fclose(indexStream);\r
 779     return true;\r
 780 }\r
 781 \r
 782 // populates BamAlignment with alignment data under file pointer, returns success/fail\r
 783 bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {\r
 784 \r
 785     // read in the 'block length' value, make sure it's not zero\r
 786     char buffer[4];\r
 787     mBGZF.Read(buffer, 4);\r
 788     const unsigned int blockLength = BgzfData::UnpackUnsignedInt(buffer);\r
 789     if ( blockLength == 0 ) { return false; }\r
 790 \r
 791     // keep track of bytes read as method progresses\r
 792     int bytesRead = 4;\r
 793 \r
 794     // read in core alignment data, make sure the right size of data was read\r
 795     char x[BAM_CORE_SIZE];\r
 796     if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }\r
 797     bytesRead += BAM_CORE_SIZE;\r
 798 \r
 799     // set BamAlignment 'core' data and character data lengths\r
 800     unsigned int tempValue;\r
 801     unsigned int queryNameLength;\r
 802     unsigned int numCigarOperations;\r
 803     unsigned int querySequenceLength;\r
 804 \r
 805     bAlignment.RefID    = BgzfData::UnpackSignedInt(&x[0]);\r
 806     bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);\r
 807 \r
 808     tempValue             = BgzfData::UnpackUnsignedInt(&x[8]);\r
 809     bAlignment.Bin        = tempValue >> 16;\r
 810     bAlignment.MapQuality = tempValue >> 8 & 0xff;\r
 811     queryNameLength       = tempValue & 0xff;\r
 812 \r
 813     tempValue                = BgzfData::UnpackUnsignedInt(&x[12]);\r
 814     bAlignment.AlignmentFlag = tempValue >> 16;\r
 815     numCigarOperations       = tempValue & 0xffff;\r
 816 \r
 817     querySequenceLength     = BgzfData::UnpackUnsignedInt(&x[16]);\r
 818     bAlignment.MateRefID    = BgzfData::UnpackSignedInt(&x[20]);\r
 819     bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);\r
 820     bAlignment.InsertSize   = BgzfData::UnpackSignedInt(&x[28]);\r
 821 \r
 822     // calculate lengths/offsets\r
 823     const unsigned int dataLength      = blockLength - BAM_CORE_SIZE;\r
 824     const unsigned int cigarDataOffset = queryNameLength;\r
 825     const unsigned int seqDataOffset   = cigarDataOffset + (numCigarOperations * 4);\r
 826     const unsigned int qualDataOffset  = seqDataOffset + (querySequenceLength+1)/2;\r
 827     const unsigned int tagDataOffset   = qualDataOffset + querySequenceLength;\r
 828     const unsigned int tagDataLen      = dataLength - tagDataOffset;\r
 829 \r
 830     // set up destination buffers for character data\r
 831     char* allCharData   = (char*)calloc(sizeof(char), dataLength);\r
 832     uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);\r
 833     char* seqData       = ((char*)allCharData) + seqDataOffset;\r
 834     char* qualData      = ((char*)allCharData) + qualDataOffset;\r
 835     char* tagData       = ((char*)allCharData) + tagDataOffset;\r
 836 \r
 837     // get character data - make sure proper data size was read\r
 838     if ( mBGZF.Read(allCharData, dataLength) != (signed int)dataLength) { return false; }\r
 839     else {\r
 840 \r
 841         bytesRead += dataLength;\r
 842 \r
 843         // clear out any previous string data\r
 844         bAlignment.Name.clear();\r
 845         bAlignment.QueryBases.clear();\r
 846         bAlignment.Qualities.clear();\r
 847         bAlignment.AlignedBases.clear();\r
 848         bAlignment.CigarData.clear();\r
 849         bAlignment.TagData.clear();\r
 850 \r
 851         // save name\r
 852         bAlignment.Name = (string)((const char*)(allCharData));\r
 853 \r
 854         // save query sequence\r
 855         for (unsigned int i = 0; i < querySequenceLength; ++i) {\r
 856             char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];\r
 857             bAlignment.QueryBases.append( 1, singleBase );\r
 858         }\r
 859 \r
 860         // save sequence length\r
 861         bAlignment.Length = bAlignment.QueryBases.length();\r
 862 \r
 863         // save qualities, convert from numeric QV to FASTQ character\r
 864         for (unsigned int i = 0; i < querySequenceLength; ++i) {\r
 865             char singleQuality = (char)(qualData[i]+33);\r
 866             bAlignment.Qualities.append( 1, singleQuality );\r
 867         }\r
 868 \r
 869         // save CIGAR-related data;\r
 870         int k = 0;\r
 871         for (unsigned int i = 0; i < numCigarOperations; ++i) {\r
 872 \r
 873             // build CigarOp struct\r
 874             CigarOp op;\r
 875             op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);\r
 876             op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];\r
 877 \r
 878             // save CigarOp\r
 879             bAlignment.CigarData.push_back(op);\r
 880 \r
 881             // build AlignedBases string\r
 882             switch (op.Type) {\r
 883 \r
 884                 case ('M') :\r
 885                 case ('I') : bAlignment.AlignedBases.append( bAlignment.QueryBases.substr(k, op.Length) ); // for 'M', 'I' - write bases\r
 886                 case ('S') : k += op.Length;                                                               // for 'S' - skip over query bases\r
 887                              break;\r
 888 \r
 889                 case ('D') : bAlignment.AlignedBases.append( op.Length, '-' );  // for 'D' - write gap character\r
 890                              break;\r
 891 \r
 892                 case ('P') : bAlignment.AlignedBases.append( op.Length, '*' );  // for 'P' - write padding character;\r
 893                              break;\r
 894 \r
 895                 case ('N') : bAlignment.AlignedBases.append( op.Length, 'N' );  // for 'N' - write N's, skip bases in query sequence\r
 896                              k += op.Length;\r
 897                              break;\r
 898 \r
 899                 case ('H') : break;                                             // for 'H' - do nothing, move to next op\r
 900 \r
 901                 default    : printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here\r
 902                              exit(1);\r
 903             }\r
 904         }\r
 905 \r
 906         // read in the tag data\r
 907         bAlignment.TagData.resize(tagDataLen);\r
 908         memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLen);\r
 909     }\r
 910 \r
 911     free(allCharData);\r
 912     return true;\r
 913 }\r
 914 \r
 915 // loads reference data from BAM file\r
 916 void BamReader::BamReaderPrivate::LoadReferenceData(void) {\r
 917 \r
 918     // get number of reference sequences\r
 919     char buffer[4];\r
 920     mBGZF.Read(buffer, 4);\r
 921     const unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);\r
 922     if (numberRefSeqs == 0) { return; }\r
 923     References.reserve((int)numberRefSeqs);\r
 924 \r
 925     // iterate over all references in header\r
 926     for (unsigned int i = 0; i != numberRefSeqs; ++i) {\r
 927 \r
 928         // get length of reference name\r
 929         mBGZF.Read(buffer, 4);\r
 930         const unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);\r
 931         char* refName = (char*)calloc(refNameLength, 1);\r
 932 \r
 933         // get reference name and reference sequence length\r
 934         mBGZF.Read(refName, refNameLength);\r
 935         mBGZF.Read(buffer, 4);\r
 936         const int refLength = BgzfData::UnpackSignedInt(buffer);\r
 937 \r
 938         // store data for reference\r
 939         RefData aReference;\r
 940         aReference.RefName   = (string)((const char*)refName);\r
 941         aReference.RefLength = refLength;\r
 942         References.push_back(aReference);\r
 943 \r
 944         // clean up calloc-ed temp variable\r
 945         free(refName);\r
 946     }\r
 947 }\r
 948 \r
 949 // merges 'alignment chunks' in BAM bin (used for index building)\r
 950 void BamReader::BamReaderPrivate::MergeChunks(void) {\r
 951 \r
 952     // iterate over reference enties\r
 953     BamIndex::iterator indexIter = Index.begin();\r
 954     BamIndex::iterator indexEnd  = Index.end();\r
 955     for ( ; indexIter != indexEnd; ++indexIter ) {\r
 956 \r
 957         // get BAM bin map for this reference\r
 958         ReferenceIndex& refIndex = (*indexIter);\r
 959         BamBinMap& bamBinMap = refIndex.Bins;\r
 960 \r
 961         // iterate over BAM bins\r
 962         BamBinMap::iterator binIter = bamBinMap.begin();\r
 963         BamBinMap::iterator binEnd  = bamBinMap.end();\r
 964         for ( ; binIter != binEnd; ++binIter ) {\r
 965 \r
 966             // get chunk vector for this bin\r
 967             ChunkVector& binChunks = (*binIter).second;\r
 968             if ( binChunks.size() == 0 ) { continue; }\r
 969 \r
 970             ChunkVector mergedChunks;\r
 971             mergedChunks.push_back( binChunks[0] );\r
 972 \r
 973             // iterate over chunks\r
 974             int i = 0;\r
 975             ChunkVector::iterator chunkIter = binChunks.begin();\r
 976             ChunkVector::iterator chunkEnd  = binChunks.end();\r
 977             for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {\r
 978 \r
 979                 // get 'currentChunk' based on numeric index\r
 980                 Chunk& currentChunk = mergedChunks[i];\r
 981 \r
 982                 // get iteratorChunk based on vector iterator\r
 983                 Chunk& iteratorChunk = (*chunkIter);\r
 984 \r
 985                 // if currentChunk.Stop(shifted) == iterator Chunk.Start(shifted)\r
 986                 if ( currentChunk.Stop>>16 == iteratorChunk.Start>>16 ) {\r
 987 \r
 988                     // set currentChunk.Stop to iteratorChunk.Stop\r
 989                     currentChunk.Stop = iteratorChunk.Stop;\r
 990                 }\r
 991 \r
 992                 // otherwise\r
 993                 else {\r
 994                     // set currentChunk + 1 to iteratorChunk\r
 995                     mergedChunks.push_back(iteratorChunk);\r
 996                     ++i;\r
 997                 }\r
 998             }\r
 999 \r
1000             // saved merged chunk vector\r
1001             (*binIter).second = mergedChunks;\r
1002         }\r
1003     }\r
1004 }\r
1005 \r
1006 // opens BAM file (and index)\r
1007 void BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename) {\r
1008 \r
1009     Filename = filename;\r
1010     IndexFilename = indexFilename;\r
1011 \r
1012     // open the BGZF file for reading, retrieve header text & reference data\r
1013     mBGZF.Open(filename, "rb");\r
1014     LoadHeaderData();\r
1015     LoadReferenceData();\r
1016 \r
1017     // store file offset of first alignment\r
1018     AlignmentsBeginOffset = mBGZF.Tell();\r
1019 \r
1020     // open index file & load index data (if exists)\r
1021     if ( !IndexFilename.empty() ) {\r
1022         LoadIndex();\r
1023     }\r
1024 }\r
1025 \r
1026 // returns BAM file pointer to beginning of alignment data\r
1027 bool BamReader::BamReaderPrivate::Rewind(void) {\r
1028 \r
1029     // find first reference that has alignments in the BAM file\r
1030     int refID = 0;\r
1031     int refCount = References.size();\r
1032     for ( ; refID < refCount; ++refID ) {\r
1033         if ( References.at(refID).RefHasAlignments ) { break; }\r
1034     }\r
1035 \r
1036     // store default bounds for first alignment\r
1037     CurrentRefID = refID;\r
1038     CurrentLeft = 0;\r
1039     IsRegionSpecified = false;\r
1040 \r
1041     // return success/failure of seek\r
1042     return mBGZF.Seek(AlignmentsBeginOffset);\r
1043 }\r
1044 \r
1045 // rounds value up to next power-of-2 (used in index building)\r
1046 void BamReader::BamReaderPrivate::Roundup32(int& value) {    \r
1047     --value;\r
1048     value |= value >> 1;\r
1049     value |= value >> 2;\r
1050     value |= value >> 4;\r
1051     value |= value >> 8;\r
1052     value |= value >> 16;\r
1053     ++value;\r
1054 }\r
1055 \r
1056 // saves index data to BAM index file (".bai"), returns success/fail\r
1057 bool BamReader::BamReaderPrivate::WriteIndex(void) {\r
1058 \r
1059     IndexFilename = Filename + ".bai";\r
1060     FILE* indexStream = fopen(IndexFilename.c_str(), "wb");\r
1061     if ( indexStream == 0 ) {\r
1062         printf("ERROR: Could not open file to save index\n");\r
1063         return false;\r
1064     }\r
1065 \r
1066     // write BAM index header\r
1067     fwrite("BAI\1", 1, 4, indexStream);\r
1068 \r
1069     // write number of reference sequences\r
1070     int32_t numReferenceSeqs = Index.size();\r
1071     fwrite(&numReferenceSeqs, 4, 1, indexStream);\r
1072 \r
1073     // iterate over reference sequences\r
1074     BamIndex::const_iterator indexIter = Index.begin();\r
1075     BamIndex::const_iterator indexEnd  = Index.end();\r
1076     for ( ; indexIter != indexEnd; ++ indexIter ) {\r
1077 \r
1078         // get reference index data\r
1079         const ReferenceIndex& refIndex = (*indexIter);\r
1080         const BamBinMap& binMap = refIndex.Bins;\r
1081         const LinearOffsetVector& offsets = refIndex.Offsets;\r
1082 \r
1083         // write number of bins\r
1084         int32_t binCount = binMap.size();\r
1085         fwrite(&binCount, 4, 1, indexStream);\r
1086 \r
1087         // iterate over bins\r
1088         BamBinMap::const_iterator binIter = binMap.begin();\r
1089         BamBinMap::const_iterator binEnd  = binMap.end();\r
1090         for ( ; binIter != binEnd; ++binIter ) {\r
1091 \r
1092             // get bin data (key and chunk vector)\r
1093             const uint32_t& binKey = (*binIter).first;\r
1094             const ChunkVector& binChunks = (*binIter).second;\r
1095 \r
1096             // save BAM bin key\r
1097             fwrite(&binKey, 4, 1, indexStream);\r
1098 \r
1099             // save chunk count\r
1100             int32_t chunkCount = binChunks.size();\r
1101             fwrite(&chunkCount, 4, 1, indexStream);\r
1102 \r
1103             // iterate over chunks\r
1104             ChunkVector::const_iterator chunkIter = binChunks.begin();\r
1105             ChunkVector::const_iterator chunkEnd  = binChunks.end();\r
1106             for ( ; chunkIter != chunkEnd; ++chunkIter ) {\r
1107 \r
1108                 // get current chunk data\r
1109                 const Chunk& chunk    = (*chunkIter);\r
1110                 const uint64_t& start = chunk.Start;\r
1111                 const uint64_t& stop  = chunk.Stop;\r
1112 \r
1113                 // save chunk offsets\r
1114                 fwrite(&start, 8, 1, indexStream);\r
1115                 fwrite(&stop,  8, 1, indexStream);\r
1116             }\r
1117         }\r
1118 \r
1119         // write linear offsets size\r
1120         int32_t offsetSize = offsets.size();\r
1121         fwrite(&offsetSize, 4, 1, indexStream);\r
1122 \r
1123         // iterate over linear offsets\r
1124         LinearOffsetVector::const_iterator offsetIter = offsets.begin();\r
1125         LinearOffsetVector::const_iterator offsetEnd  = offsets.end();\r
1126         for ( ; offsetIter != offsetEnd; ++offsetIter ) {\r
1127 \r
1128             // write linear offset value\r
1129             const uint64_t& linearOffset = (*offsetIter);\r
1130             fwrite(&linearOffset, 8, 1, indexStream);\r
1131         }\r
1132     }\r
1133 \r
1134     // flush buffer, close file, and return success\r
1135     fflush(indexStream);\r
1136     fclose(indexStream);\r
1137     return true;\r
1138 }\r
1139 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamReader.h\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000005312\011307520655\0012736\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
1140 // BamReader.h (c) 2009 Derek Barnett, Michael Strömberg\r
1141 // Marth Lab, Department of Biology, Boston College\r
1142 // All rights reserved.\r
1143 // ---------------------------------------------------------------------------\r
1144 // Last modified: 8 December 2009 (DB)\r
1145 // ---------------------------------------------------------------------------\r
1146 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
1147 // Institute.\r
1148 // ---------------------------------------------------------------------------\r
1149 // Provides the basic functionality for reading BAM files\r
1150 // ***************************************************************************\r
1151 \r
1152 #ifndef BAMREADER_H\r
1153 #define BAMREADER_H\r
1154 \r
1155 // C++ includes\r
1156 #include <string>\r
1157 \r
1158 // BamTools includes\r
1159 #include "BamAux.h"\r
1160 \r
1161 namespace BamTools {\r
1162 \r
1163 class BamReader {\r
1164 \r
1165     // constructor / destructor\r
1166     public:\r
1167         BamReader(void);\r
1168         ~BamReader(void);\r
1169 \r
1170     // public interface\r
1171     public:\r
1172 \r
1173         // ----------------------\r
1174         // BAM file operations\r
1175         // ----------------------\r
1176 \r
1177         // close BAM file\r
1178         void Close(void);\r
1179         // performs random-access jump to reference, position\r
1180         bool Jump(int refID, int position = 0);\r
1181         // opens BAM file (and optional BAM index file, if provided)\r
1182         void Open(const std::string& filename, const std::string& indexFilename = "");\r
1183         // returns file pointer to beginning of alignments\r
1184         bool Rewind(void);\r
1185 \r
1186         // ----------------------\r
1187         // access alignment data\r
1188         // ----------------------\r
1189 \r
1190         // retrieves next available alignment (returns success/fail)\r
1191         bool GetNextAlignment(BamAlignment& bAlignment);\r
1192 \r
1193         // ----------------------\r
1194         // access auxiliary data\r
1195         // ----------------------\r
1196 \r
1197         // returns SAM header text\r
1198         const std::string GetHeaderText(void) const;\r
1199         // returns number of reference sequences\r
1200         const int GetReferenceCount(void) const;\r
1201         // returns vector of reference objects\r
1202         const BamTools::RefVector GetReferenceData(void) const;\r
1203         // returns reference id (used for BamReader::Jump()) for the given reference name\r
1204         const int GetReferenceID(const std::string& refName) const;\r
1205 \r
1206         // ----------------------\r
1207         // BAM index operations\r
1208         // ----------------------\r
1209 \r
1210         // creates index for BAM file, saves to file (default = bamFilename + ".bai")\r
1211         bool CreateIndex(void);\r
1212 \r
1213     // private implementation\r
1214     private:\r
1215         struct BamReaderPrivate;\r
1216         BamReaderPrivate* d;\r
1217 };\r
1218 \r
1219 } // namespace BamTools\r
1220 \r
1221 #endif // BAMREADER_H\r
1222 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamTrimMain.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000005351\011307520655\0013612\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
1223 // BamTrimMain.cpp (c) 2009 Derek Barnett\r
1224 // Marth Lab, Department of Biology, Boston College\r
1225 // All rights reserved.\r
1226 // ---------------------------------------------------------------------------\r
1227 // Last modified: 15 July 2009 (DB)\r
1228 // ---------------------------------------------------------------------------\r
1229 // Basic example of reading/writing BAM files. Pulls alignments overlapping \r
1230 // the range specified by user from one BAM file and writes to a new BAM file.\r
1231 // ***************************************************************************\r
1232 \r
1233 // Std C/C++ includes\r
1234 #include <cstdlib>\r
1235 #include <iostream>\r
1236 #include <string>\r
1237 using namespace std;\r
1238 \r
1239 // BamTools includes\r
1240 #include "BamReader.h"\r
1241 #include "BamWriter.h"\r
1242 using namespace BamTools;\r
1243 \r
1244 int main(int argc, char* argv[]) {\r
1245 \r
1246         // validate argument count\r
1247         if( argc != 7 ) {\r
1248                 cerr << "USAGE: " << argv[0] << " <input BAM file> <input BAM index file> <output BAM file> <reference name> <leftBound> <rightBound> " << endl;\r
1249                 exit(1);\r
1250         }\r
1251 \r
1252         // store arguments\r
1253         string inBamFilename  = argv[1];\r
1254         string indexFilename  = argv[2];\r
1255         string outBamFilename = argv[3];\r
1256         string referenceName  = argv[4];\r
1257         string leftBound_str  = argv[5];\r
1258         string rightBound_str = argv[6];\r
1259 \r
1260         // open our BAM reader\r
1261         BamReader reader;\r
1262         reader.Open(inBamFilename, indexFilename);\r
1263 \r
1264         // get header & reference information\r
1265         string header = reader.GetHeaderText();\r
1266         RefVector references = reader.GetReferenceData();\r
1267         \r
1268         // open our BAM writer\r
1269         BamWriter writer;\r
1270         writer.Open(outBamFilename, header, references);\r
1271 \r
1272         // get reference ID from name\r
1273         int refID = 0;\r
1274         RefVector::const_iterator refIter = references.begin();\r
1275         RefVector::const_iterator refEnd  = references.end();\r
1276         for ( ; refIter != refEnd; ++refIter ) {\r
1277                 if ( (*refIter).RefName == referenceName ) { break; }\r
1278                 ++refID;\r
1279         }\r
1280         \r
1281         // validate ID\r
1282         if ( refIter == refEnd ) {\r
1283                 cerr << "Reference: " << referenceName << " not found." << endl;\r
1284                 exit(1);\r
1285         }\r
1286 \r
1287         // convert boundary arguments to numeric values\r
1288         int leftBound  = (int) atoi( leftBound_str.c_str()  );\r
1289         int rightBound = (int) atoi( rightBound_str.c_str() );\r
1290         \r
1291         // attempt jump to range of interest\r
1292         if ( reader.Jump(refID, leftBound) ) {\r
1293                 \r
1294                 // while data exists and alignment begin before right bound\r
1295                 BamAlignment bAlignment;\r
1296                 while ( reader.GetNextAlignment(bAlignment) && (bAlignment.Position <= rightBound) ) {\r
1297                         // save alignment to archive\r
1298                         writer.SaveAlignment(bAlignment);\r
1299                 }\r
1300         } \r
1301         \r
1302         // if jump failed\r
1303         else {\r
1304                 cerr << "Could not jump to ref:pos " << referenceName << ":" << leftBound << endl;\r
1305                 exit(1);\r
1306         }\r
1307 \r
1308         // clean up and exit\r
1309         reader.Close();\r
1310         writer.Close(); \r
1311         return 0;\r
1312 }\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamWriter.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000022463\011307516613\0013350\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
1313 // BamWriter.cpp (c) 2009 Michael Strömberg, Derek Barnett\r
1314 // Marth Lab, Department of Biology, Boston College\r
1315 // All rights reserved.\r
1316 // ---------------------------------------------------------------------------\r
1317 // Last modified: 8 December 2009 (DB)\r
1318 // ---------------------------------------------------------------------------\r
1319 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
1320 // Institute.\r
1321 // ---------------------------------------------------------------------------\r
1322 // Provides the basic functionality for producing BAM files\r
1323 // ***************************************************************************\r
1324 \r
1325 // BGZF includes\r
1326 #include "BGZF.h"\r
1327 #include "BamWriter.h"\r
1328 using namespace BamTools;\r
1329 using namespace std;\r
1330 \r
1331 struct BamWriter::BamWriterPrivate {\r
1332 \r
1333     // data members\r
1334     BgzfData mBGZF;\r
1335 \r
1336     // constructor / destructor\r
1337     BamWriterPrivate(void) { }\r
1338     ~BamWriterPrivate(void) {\r
1339         mBGZF.Close();\r
1340     }\r
1341 \r
1342     // "public" interface\r
1343     void Close(void);\r
1344     void Open(const std::string& filename, const std::string& samHeader, const BamTools::RefVector& referenceSequences);\r
1345     void SaveAlignment(const BamTools::BamAlignment& al);\r
1346 \r
1347     // internal methods\r
1348     void CreatePackedCigar(const std::vector<CigarOp>& cigarOperations, std::string& packedCigar);\r
1349     void EncodeQuerySequence(const std::string& query, std::string& encodedQuery);\r
1350 };\r
1351 \r
1352 // -----------------------------------------------------\r
1353 // BamWriter implementation\r
1354 // -----------------------------------------------------\r
1355 \r
1356 // constructor\r
1357 BamWriter::BamWriter(void) {\r
1358     d = new BamWriterPrivate;\r
1359 }\r
1360 \r
1361 // destructor\r
1362 BamWriter::~BamWriter(void) {\r
1363     delete d;\r
1364     d = 0;\r
1365 }\r
1366 \r
1367 // closes the alignment archive\r
1368 void BamWriter::Close(void) {\r
1369     d->Close();\r
1370 }\r
1371 \r
1372 // opens the alignment archive\r
1373 void BamWriter::Open(const string& filename, const string& samHeader, const RefVector& referenceSequences) {\r
1374     d->Open(filename, samHeader, referenceSequences);\r
1375 }\r
1376 \r
1377 // saves the alignment to the alignment archive\r
1378 void BamWriter::SaveAlignment(const BamAlignment& al) {\r
1379     d->SaveAlignment(al);\r
1380 }\r
1381 \r
1382 // -----------------------------------------------------\r
1383 // BamWriterPrivate implementation\r
1384 // -----------------------------------------------------\r
1385 \r
1386 // closes the alignment archive\r
1387 void BamWriter::BamWriterPrivate::Close(void) {\r
1388     mBGZF.Close();\r
1389 }\r
1390 \r
1391 // creates a cigar string from the supplied alignment\r
1392 void BamWriter::BamWriterPrivate::CreatePackedCigar(const vector<CigarOp>& cigarOperations, string& packedCigar) {\r
1393 \r
1394     // initialize\r
1395     const unsigned int numCigarOperations = cigarOperations.size();\r
1396     packedCigar.resize(numCigarOperations * BT_SIZEOF_INT);\r
1397 \r
1398     // pack the cigar data into the string\r
1399     unsigned int* pPackedCigar = (unsigned int*)packedCigar.data();\r
1400 \r
1401     unsigned int cigarOp;\r
1402     vector<CigarOp>::const_iterator coIter;\r
1403     for(coIter = cigarOperations.begin(); coIter != cigarOperations.end(); coIter++) {\r
1404 \r
1405         switch(coIter->Type) {\r
1406             case 'M':\r
1407                   cigarOp = BAM_CMATCH;\r
1408                   break;\r
1409             case 'I':\r
1410                   cigarOp = BAM_CINS;\r
1411                   break;\r
1412             case 'D':\r
1413                   cigarOp = BAM_CDEL;\r
1414                   break;\r
1415             case 'N':\r
1416                   cigarOp = BAM_CREF_SKIP;\r
1417                   break;\r
1418             case 'S':\r
1419                   cigarOp = BAM_CSOFT_CLIP;\r
1420                   break;\r
1421             case 'H':\r
1422                   cigarOp = BAM_CHARD_CLIP;\r
1423                   break;\r
1424             case 'P':\r
1425                   cigarOp = BAM_CPAD;\r
1426                   break;\r
1427             default:\r
1428                   printf("ERROR: Unknown cigar operation found: %c\n", coIter->Type);\r
1429                   exit(1);\r
1430         }\r
1431 \r
1432         *pPackedCigar = coIter->Length << BAM_CIGAR_SHIFT | cigarOp;\r
1433         pPackedCigar++;\r
1434     }\r
1435 }\r
1436 \r
1437 // encodes the supplied query sequence into 4-bit notation\r
1438 void BamWriter::BamWriterPrivate::EncodeQuerySequence(const string& query, string& encodedQuery) {\r
1439 \r
1440     // prepare the encoded query string\r
1441     const unsigned int queryLen = query.size();\r
1442     const unsigned int encodedQueryLen = (unsigned int)((queryLen / 2.0) + 0.5);\r
1443     encodedQuery.resize(encodedQueryLen);\r
1444     char* pEncodedQuery = (char*)encodedQuery.data();\r
1445     const char* pQuery = (const char*)query.data();\r
1446 \r
1447     unsigned char nucleotideCode;\r
1448     bool useHighWord = true;\r
1449 \r
1450     while(*pQuery) {\r
1451 \r
1452         switch(*pQuery) {\r
1453             case '=':\r
1454                     nucleotideCode = 0;\r
1455                     break;\r
1456             case 'A':\r
1457                     nucleotideCode = 1;\r
1458                     break;\r
1459             case 'C':\r
1460                     nucleotideCode = 2;\r
1461                     break;\r
1462             case 'G':\r
1463                     nucleotideCode = 4;\r
1464                     break;\r
1465             case 'T':\r
1466                     nucleotideCode = 8;\r
1467                     break;\r
1468             case 'N':\r
1469                     nucleotideCode = 15;\r
1470                     break;\r
1471             default:\r
1472                     printf("ERROR: Only the following bases are supported in the BAM format: {=, A, C, G, T, N}. Found [%c]\n", *pQuery);\r
1473                     exit(1);\r
1474         }\r
1475 \r
1476         // pack the nucleotide code\r
1477         if(useHighWord) {\r
1478             *pEncodedQuery = nucleotideCode << 4;\r
1479             useHighWord = false;\r
1480         } else {\r
1481             *pEncodedQuery |= nucleotideCode;\r
1482             pEncodedQuery++;\r
1483             useHighWord = true;\r
1484         }\r
1485 \r
1486         // increment the query position\r
1487         pQuery++;\r
1488     }\r
1489 }\r
1490 \r
1491 // opens the alignment archive\r
1492 void BamWriter::BamWriterPrivate::Open(const string& filename, const string& samHeader, const RefVector& referenceSequences) {\r
1493 \r
1494     // open the BGZF file for writing\r
1495     mBGZF.Open(filename, "wb");\r
1496 \r
1497     // ================\r
1498     // write the header\r
1499     // ================\r
1500 \r
1501     // write the BAM signature\r
1502     const unsigned char SIGNATURE_LENGTH = 4;\r
1503     const char* BAM_SIGNATURE = "BAM\1";\r
1504     mBGZF.Write(BAM_SIGNATURE, SIGNATURE_LENGTH);\r
1505 \r
1506     // write the SAM header text length\r
1507     const unsigned int samHeaderLen = samHeader.size();\r
1508     mBGZF.Write((char*)&samHeaderLen, BT_SIZEOF_INT);\r
1509 \r
1510     // write the SAM header text\r
1511     if(samHeaderLen > 0) {\r
1512         mBGZF.Write(samHeader.data(), samHeaderLen);\r
1513     }\r
1514 \r
1515     // write the number of reference sequences\r
1516     const unsigned int numReferenceSequences = referenceSequences.size();\r
1517     mBGZF.Write((char*)&numReferenceSequences, BT_SIZEOF_INT);\r
1518 \r
1519     // =============================\r
1520     // write the sequence dictionary\r
1521     // =============================\r
1522 \r
1523     RefVector::const_iterator rsIter;\r
1524     for(rsIter = referenceSequences.begin(); rsIter != referenceSequences.end(); rsIter++) {\r
1525 \r
1526         // write the reference sequence name length\r
1527         const unsigned int referenceSequenceNameLen = rsIter->RefName.size() + 1;\r
1528         mBGZF.Write((char*)&referenceSequenceNameLen, BT_SIZEOF_INT);\r
1529 \r
1530         // write the reference sequence name\r
1531         mBGZF.Write(rsIter->RefName.c_str(), referenceSequenceNameLen);\r
1532 \r
1533         // write the reference sequence length\r
1534         mBGZF.Write((char*)&rsIter->RefLength, BT_SIZEOF_INT);\r
1535     }\r
1536 }\r
1537 \r
1538 // saves the alignment to the alignment archive\r
1539 void BamWriter::BamWriterPrivate::SaveAlignment(const BamAlignment& al) {\r
1540 \r
1541     // initialize\r
1542     const unsigned int nameLen            = al.Name.size() + 1;\r
1543     const unsigned int queryLen           = al.QueryBases.size();\r
1544     const unsigned int numCigarOperations = al.CigarData.size();\r
1545 \r
1546     // create our packed cigar string\r
1547     string packedCigar;\r
1548     CreatePackedCigar(al.CigarData, packedCigar);\r
1549     const unsigned int packedCigarLen = packedCigar.size();\r
1550 \r
1551     // encode the query\r
1552     string encodedQuery;\r
1553     EncodeQuerySequence(al.QueryBases, encodedQuery);\r
1554     const unsigned int encodedQueryLen = encodedQuery.size();\r
1555 \r
1556     // store the tag data length\r
1557     const unsigned int tagDataLength = al.TagData.size() + 1;\r
1558 \r
1559     // assign the BAM core data\r
1560     unsigned int buffer[8];\r
1561     buffer[0] = al.RefID;\r
1562     buffer[1] = al.Position;\r
1563     buffer[2] = (al.Bin << 16) | (al.MapQuality << 8) | nameLen;\r
1564     buffer[3] = (al.AlignmentFlag << 16) | numCigarOperations;\r
1565     buffer[4] = queryLen;\r
1566     buffer[5] = al.MateRefID;\r
1567     buffer[6] = al.MatePosition;\r
1568     buffer[7] = al.InsertSize;\r
1569 \r
1570     // write the block size\r
1571     const unsigned int dataBlockSize = nameLen + packedCigarLen + encodedQueryLen + queryLen + tagDataLength;\r
1572     const unsigned int blockSize = BAM_CORE_SIZE + dataBlockSize;\r
1573     mBGZF.Write((char*)&blockSize, BT_SIZEOF_INT);\r
1574 \r
1575     // write the BAM core\r
1576     mBGZF.Write((char*)&buffer, BAM_CORE_SIZE);\r
1577 \r
1578     // write the query name\r
1579     mBGZF.Write(al.Name.c_str(), nameLen);\r
1580 \r
1581     // write the packed cigar\r
1582     mBGZF.Write(packedCigar.data(), packedCigarLen);\r
1583 \r
1584     // write the encoded query sequence\r
1585     mBGZF.Write(encodedQuery.data(), encodedQueryLen);\r
1586 \r
1587     // write the base qualities\r
1588     string baseQualities = al.Qualities;\r
1589     char* pBaseQualities = (char*)al.Qualities.data();\r
1590     for(unsigned int i = 0; i < queryLen; i++) { pBaseQualities[i] -= 33; }\r
1591     mBGZF.Write(pBaseQualities, queryLen);\r
1592 \r
1593     // write the read group tag\r
1594     mBGZF.Write(al.TagData.data(), tagDataLength);\r
1595 }\r
1596 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BamWriter.h\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000003034\011307516613\0013006\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
1597 // BamWriter.h (c) 2009 Michael Strömberg, Derek Barnett\r
1598 // Marth Lab, Department of Biology, Boston College\r
1599 // All rights reserved.\r
1600 // ---------------------------------------------------------------------------\r
1601 // Last modified: 8 December 2009 (DB)\r
1602 // ---------------------------------------------------------------------------\r
1603 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
1604 // Institute.\r
1605 // ---------------------------------------------------------------------------\r
1606 // Provides the basic functionality for producing BAM files\r
1607 // ***************************************************************************\r
1608 \r
1609 #ifndef BAMWRITER_H\r
1610 #define BAMWRITER_H\r
1611 \r
1612 // C++ includes\r
1613 #include <string>\r
1614 \r
1615 // BamTools includes\r
1616 #include "BamAux.h"\r
1617 \r
1618 namespace BamTools {\r
1619 \r
1620 class BamWriter {\r
1621 \r
1622     // constructor/destructor\r
1623     public:\r
1624         BamWriter(void);\r
1625         ~BamWriter(void);\r
1626 \r
1627     // public interface\r
1628     public:\r
1629         // closes the alignment archive\r
1630         void Close(void);\r
1631         // opens the alignment archive\r
1632         void Open(const std::string& filename, const std::string& samHeader, const BamTools::RefVector& referenceSequences);\r
1633         // saves the alignment to the alignment archive\r
1634         void SaveAlignment(const BamTools::BamAlignment& al);\r
1635 \r
1636     // private implementation\r
1637     private:\r
1638         struct BamWriterPrivate;\r
1639         BamWriterPrivate* d;\r
1640 };\r
1641 \r
1642 } // namespace BamTools\r
1643 \r
1644 #endif // BAMWRITER_H\r
1645 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BGZF.cpp\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000024215\011307516613\0012201\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
1646 // BGZF.cpp (c) 2009 Derek Barnett, Michael Strömberg\r
1647 // Marth Lab, Department of Biology, Boston College\r
1648 // All rights reserved.\r
1649 // ---------------------------------------------------------------------------\r
1650 // Last modified: 8 December 2009 (DB)\r
1651 // ---------------------------------------------------------------------------\r
1652 // BGZF routines were adapted from the bgzf.c code developed at the Broad\r
1653 // Institute.\r
1654 // ---------------------------------------------------------------------------\r
1655 // Provides the basic functionality for reading & writing BGZF files\r
1656 // ***************************************************************************\r
1657 \r
1658 #include <algorithm>\r
1659 #include "BGZF.h"\r
1660 using namespace BamTools;\r
1661 using std::string;\r
1662 using std::min;\r
1663 \r
1664 BgzfData::BgzfData(void)\r
1665     : UncompressedBlockSize(DEFAULT_BLOCK_SIZE)\r
1666     , CompressedBlockSize(MAX_BLOCK_SIZE)\r
1667     , BlockLength(0)\r
1668     , BlockOffset(0)\r
1669     , BlockAddress(0)\r
1670     , IsOpen(false)\r
1671     , IsWriteOnly(false)\r
1672     , Stream(NULL)\r
1673     , UncompressedBlock(NULL)\r
1674     , CompressedBlock(NULL)\r
1675 {\r
1676     try {\r
1677         CompressedBlock   = new char[CompressedBlockSize];\r
1678         UncompressedBlock = new char[UncompressedBlockSize];\r
1679     } catch( std::bad_alloc& ba ) {\r
1680         printf("ERROR: Unable to allocate memory for our BGZF object.\n");\r
1681         exit(1);\r
1682     }\r
1683 }\r
1684 \r
1685 // destructor\r
1686 BgzfData::~BgzfData(void) {\r
1687     if(CompressedBlock)   delete [] CompressedBlock;\r
1688     if(UncompressedBlock) delete [] UncompressedBlock;\r
1689 }\r
1690 \r
1691 // closes BGZF file\r
1692 void BgzfData::Close(void) {\r
1693 \r
1694     if (!IsOpen) { return; }\r
1695     IsOpen = false;\r
1696 \r
1697     // flush the BGZF block\r
1698     if ( IsWriteOnly ) { FlushBlock(); }\r
1699 \r
1700     // flush and close\r
1701     fflush(Stream);\r
1702     fclose(Stream);\r
1703 }\r
1704 \r
1705 // compresses the current block\r
1706 int BgzfData::DeflateBlock(void) {\r
1707 \r
1708     // initialize the gzip header\r
1709     char* buffer = CompressedBlock;\r
1710     unsigned int bufferSize = CompressedBlockSize;\r
1711 \r
1712     memset(buffer, 0, 18);\r
1713     buffer[0]  = GZIP_ID1;\r
1714     buffer[1]  = (char)GZIP_ID2;\r
1715     buffer[2]  = CM_DEFLATE;\r
1716     buffer[3]  = FLG_FEXTRA;\r
1717     buffer[9]  = (char)OS_UNKNOWN;\r
1718     buffer[10] = BGZF_XLEN;\r
1719     buffer[12] = BGZF_ID1;\r
1720     buffer[13] = BGZF_ID2;\r
1721     buffer[14] = BGZF_LEN;\r
1722 \r
1723     // loop to retry for blocks that do not compress enough\r
1724     int inputLength = BlockOffset;\r
1725     int compressedLength = 0;\r
1726 \r
1727     while(true) {\r
1728 \r
1729         z_stream zs;\r
1730         zs.zalloc    = NULL;\r
1731         zs.zfree     = NULL;\r
1732         zs.next_in   = (Bytef*)UncompressedBlock;\r
1733         zs.avail_in  = inputLength;\r
1734         zs.next_out  = (Bytef*)&buffer[BLOCK_HEADER_LENGTH];\r
1735         zs.avail_out = bufferSize - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;\r
1736 \r
1737         // initialize the zlib compression algorithm\r
1738         if(deflateInit2(&zs, Z_DEFAULT_COMPRESSION, Z_DEFLATED, GZIP_WINDOW_BITS, Z_DEFAULT_MEM_LEVEL, Z_DEFAULT_STRATEGY) != Z_OK) {\r
1739             printf("ERROR: zlib deflate initialization failed.\n");\r
1740             exit(1);\r
1741         }\r
1742 \r
1743         // compress the data\r
1744         int status = deflate(&zs, Z_FINISH);\r
1745         if(status != Z_STREAM_END) {\r
1746 \r
1747             deflateEnd(&zs);\r
1748 \r
1749             // reduce the input length and try again\r
1750             if(status == Z_OK) {\r
1751                 inputLength -= 1024;\r
1752                 if(inputLength < 0) {\r
1753                     printf("ERROR: input reduction failed.\n");\r
1754                     exit(1);\r
1755                 }\r
1756                 continue;\r
1757             }\r
1758 \r
1759             printf("ERROR: zlib deflate failed.\n");\r
1760             exit(1);\r
1761         }\r
1762 \r
1763         // finalize the compression routine\r
1764         if(deflateEnd(&zs) != Z_OK) {\r
1765             printf("ERROR: deflate end failed.\n");\r
1766             exit(1);\r
1767         }\r
1768 \r
1769         compressedLength = zs.total_out;\r
1770         compressedLength += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;\r
1771 \r
1772         if(compressedLength > MAX_BLOCK_SIZE) {\r
1773             printf("ERROR: deflate overflow.\n");\r
1774             exit(1);\r
1775         }\r
1776 \r
1777         break;\r
1778     }\r
1779 \r
1780     // store the compressed length\r
1781     BgzfData::PackUnsignedShort(&buffer[16], (unsigned short)(compressedLength - 1));\r
1782 \r
1783     // store the CRC32 checksum\r
1784     unsigned int crc = crc32(0, NULL, 0);\r
1785     crc = crc32(crc, (Bytef*)UncompressedBlock, inputLength);\r
1786     BgzfData::PackUnsignedInt(&buffer[compressedLength - 8], crc);\r
1787     BgzfData::PackUnsignedInt(&buffer[compressedLength - 4], inputLength);\r
1788 \r
1789     // ensure that we have less than a block of data left\r
1790     int remaining = BlockOffset - inputLength;\r
1791     if(remaining > 0) {\r
1792         if(remaining > inputLength) {\r
1793             printf("ERROR: remainder too large.\n");\r
1794             exit(1);\r
1795         }\r
1796         memcpy(UncompressedBlock, UncompressedBlock + inputLength, remaining);\r
1797     }\r
1798 \r
1799     BlockOffset = remaining;\r
1800     return compressedLength;\r
1801 }\r
1802 \r
1803 // flushes the data in the BGZF block\r
1804 void BgzfData::FlushBlock(void) {\r
1805 \r
1806     // flush all of the remaining blocks\r
1807     while(BlockOffset > 0) {\r
1808 \r
1809         // compress the data block\r
1810         int blockLength = DeflateBlock();\r
1811 \r
1812         // flush the data to our output stream\r
1813         int numBytesWritten = fwrite(CompressedBlock, 1, blockLength, Stream);\r
1814 \r
1815         if(numBytesWritten != blockLength) {\r
1816             printf("ERROR: Expected to write %u bytes during flushing, but wrote %u bytes.\n", blockLength, numBytesWritten);\r
1817             exit(1);\r
1818         }\r
1819 \r
1820         BlockAddress += blockLength;\r
1821     }\r
1822 }\r
1823 \r
1824 // de-compresses the current block\r
1825 int BgzfData::InflateBlock(const int& blockLength) {\r
1826 \r
1827     // Inflate the block in m_BGZF.CompressedBlock into m_BGZF.UncompressedBlock\r
1828     z_stream zs;\r
1829     zs.zalloc    = NULL;\r
1830     zs.zfree     = NULL;\r
1831     zs.next_in   = (Bytef*)CompressedBlock + 18;\r
1832     zs.avail_in  = blockLength - 16;\r
1833     zs.next_out  = (Bytef*)UncompressedBlock;\r
1834     zs.avail_out = UncompressedBlockSize;\r
1835 \r
1836     int status = inflateInit2(&zs, GZIP_WINDOW_BITS);\r
1837     if (status != Z_OK) {\r
1838         printf("inflateInit failed\n");\r
1839         exit(1);\r
1840     }\r
1841 \r
1842     status = inflate(&zs, Z_FINISH);\r
1843     if (status != Z_STREAM_END) {\r
1844         inflateEnd(&zs);\r
1845         printf("inflate failed\n");\r
1846         exit(1);\r
1847     }\r
1848 \r
1849     status = inflateEnd(&zs);\r
1850     if (status != Z_OK) {\r
1851         printf("inflateEnd failed\n");\r
1852         exit(1);\r
1853     }\r
1854 \r
1855     return zs.total_out;\r
1856 }\r
1857 \r
1858 void BgzfData::Open(const string& filename, const char* mode) {\r
1859 \r
1860     if ( strcmp(mode, "rb") == 0 ) {\r
1861         IsWriteOnly = false;\r
1862     } else if ( strcmp(mode, "wb") == 0) {\r
1863         IsWriteOnly = true;\r
1864     } else {\r
1865         printf("ERROR: Unknown file mode: %s\n", mode);\r
1866         exit(1);\r
1867     }\r
1868 \r
1869     Stream = fopen(filename.c_str(), mode);\r
1870     if(!Stream) {\r
1871         printf("ERROR: Unable to open the BAM file %s\n", filename.c_str() );\r
1872         exit(1);\r
1873     }\r
1874     IsOpen = true;\r
1875 }\r
1876 \r
1877 int BgzfData::Read(char* data, const unsigned int dataLength) {\r
1878 \r
1879    if (dataLength == 0) { return 0; }\r
1880 \r
1881    char* output = data;\r
1882    unsigned int numBytesRead = 0;\r
1883    while (numBytesRead < dataLength) {\r
1884 \r
1885        int bytesAvailable = BlockLength - BlockOffset;\r
1886        if (bytesAvailable <= 0) {\r
1887            if ( ReadBlock() != 0 ) { return -1; }\r
1888            bytesAvailable = BlockLength - BlockOffset;\r
1889            if ( bytesAvailable <= 0 ) { break; }\r
1890        }\r
1891 \r
1892        char* buffer   = UncompressedBlock;\r
1893        int copyLength = min( (int)(dataLength-numBytesRead), bytesAvailable );\r
1894        memcpy(output, buffer + BlockOffset, copyLength);\r
1895 \r
1896        BlockOffset  += copyLength;\r
1897        output       += copyLength;\r
1898        numBytesRead += copyLength;\r
1899    }\r
1900 \r
1901    if ( BlockOffset == BlockLength ) {\r
1902        BlockAddress = ftell(Stream);\r
1903        BlockOffset  = 0;\r
1904        BlockLength  = 0;\r
1905    }\r
1906 \r
1907    return numBytesRead;\r
1908 }\r
1909 \r
1910 int BgzfData::ReadBlock(void) {\r
1911 \r
1912     char    header[BLOCK_HEADER_LENGTH];\r
1913     int64_t blockAddress = ftell(Stream);\r
1914 \r
1915     int count = fread(header, 1, sizeof(header), Stream);\r
1916     if (count == 0) {\r
1917         BlockLength = 0;\r
1918         return 0;\r
1919     }\r
1920 \r
1921     if (count != sizeof(header)) {\r
1922         printf("read block failed - count != sizeof(header)\n");\r
1923         return -1;\r
1924     }\r
1925 \r
1926     if (!BgzfData::CheckBlockHeader(header)) {\r
1927         printf("read block failed - CheckBlockHeader() returned false\n");\r
1928         return -1;\r
1929     }\r
1930 \r
1931     int blockLength = BgzfData::UnpackUnsignedShort(&header[16]) + 1;\r
1932     char* compressedBlock = CompressedBlock;\r
1933     memcpy(compressedBlock, header, BLOCK_HEADER_LENGTH);\r
1934     int remaining = blockLength - BLOCK_HEADER_LENGTH;\r
1935 \r
1936     count = fread(&compressedBlock[BLOCK_HEADER_LENGTH], 1, remaining, Stream);\r
1937     if (count != remaining) {\r
1938         printf("read block failed - count != remaining\n");\r
1939         return -1;\r
1940     }\r
1941 \r
1942     count = InflateBlock(blockLength);\r
1943     if (count < 0) { return -1; }\r
1944 \r
1945     if ( BlockLength != 0 ) {\r
1946         BlockOffset = 0;\r
1947     }\r
1948 \r
1949     BlockAddress = blockAddress;\r
1950     BlockLength  = count;\r
1951     return 0;\r
1952 }\r
1953 \r
1954 bool BgzfData::Seek(int64_t position) {\r
1955 \r
1956     int     blockOffset  = (position & 0xFFFF);\r
1957     int64_t blockAddress = (position >> 16) & 0xFFFFFFFFFFFFLL;\r
1958 \r
1959     if (fseek(Stream, blockAddress, SEEK_SET) != 0) {\r
1960         printf("ERROR: Unable to seek in BAM file\n");\r
1961         exit(1);\r
1962     }\r
1963 \r
1964     BlockLength  = 0;\r
1965     BlockAddress = blockAddress;\r
1966     BlockOffset  = blockOffset;\r
1967     return true;\r
1968 }\r
1969 \r
1970 int64_t BgzfData::Tell(void) {\r
1971     return ( (BlockAddress << 16) | (BlockOffset & 0xFFFF) );\r
1972 }\r
1973 \r
1974 // writes the supplied data into the BGZF buffer\r
1975 unsigned int BgzfData::Write(const char* data, const unsigned int dataLen) {\r
1976 \r
1977     // initialize\r
1978     unsigned int numBytesWritten = 0;\r
1979     const char* input = data;\r
1980     unsigned int blockLength = UncompressedBlockSize;\r
1981 \r
1982     // copy the data to the buffer\r
1983     while(numBytesWritten < dataLen) {\r
1984         unsigned int copyLength = min(blockLength - BlockOffset, dataLen - numBytesWritten);\r
1985         char* buffer = UncompressedBlock;\r
1986         memcpy(buffer + BlockOffset, input, copyLength);\r
1987 \r
1988         BlockOffset     += copyLength;\r
1989         input           += copyLength;\r
1990         numBytesWritten += copyLength;\r
1991 \r
1992         if(BlockOffset == blockLength) {\r
1993             FlushBlock();\r
1994         }\r
1995     }\r
1996 \r
1997     return numBytesWritten;\r
1998 }\r
1999 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0BGZF.h\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000013625\011307516613\0011651\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0// ***************************************************************************\r
2000 // BGZF.h (c) 2009 Derek Barnett, Michael Strömberg\r
2001 // Marth Lab, Department of Biology, Boston College\r
2002 // All rights reserved.\r
2003 // ---------------------------------------------------------------------------\r
2004 // Last modified: 8 December 2009 (DB)\r
2005 // ---------------------------------------------------------------------------\r
2006 // BGZF routines were adapted from the bgzf.c code developed at the Broad\r
2007 // Institute.\r
2008 // ---------------------------------------------------------------------------\r
2009 // Provides the basic functionality for reading & writing BGZF files\r
2010 // ***************************************************************************\r
2011 \r
2012 #ifndef BGZF_H\r
2013 #define BGZF_H\r
2014 \r
2015 // 'C' includes\r
2016 #include <cstdio>\r
2017 #include <cstdlib>\r
2018 #include <cstring>\r
2019 \r
2020 // C++ includes\r
2021 #include <string>\r
2022 \r
2023 // zlib includes\r
2024 #include "zlib.h"\r
2025 \r
2026 // Platform-specific type definitions\r
2027 #ifdef _MSC_VER\r
2028         typedef char                 int8_t;\r
2029         typedef unsigned char       uint8_t;\r
2030         typedef short               int16_t;\r
2031         typedef unsigned short     uint16_t;\r
2032         typedef int                 int32_t;\r
2033         typedef unsigned int       uint32_t;\r
2034         typedef long long           int64_t;\r
2035         typedef unsigned long long uint64_t;\r
2036 #else\r
2037         #include <stdint.h>\r
2038 #endif\r
2039 \r
2040 namespace BamTools {\r
2041 \r
2042 // zlib constants\r
2043 const int GZIP_ID1   = 31;\r
2044 const int GZIP_ID2   = 139;\r
2045 const int CM_DEFLATE = 8;\r
2046 const int FLG_FEXTRA = 4;\r
2047 const int OS_UNKNOWN = 255;\r
2048 const int BGZF_XLEN  = 6;\r
2049 const int BGZF_ID1   = 66;\r
2050 const int BGZF_ID2   = 67;\r
2051 const int BGZF_LEN   = 2;\r
2052 const int GZIP_WINDOW_BITS    = -15;\r
2053 const int Z_DEFAULT_MEM_LEVEL = 8;\r
2054 \r
2055 // BZGF constants\r
2056 const int BLOCK_HEADER_LENGTH = 18;\r
2057 const int BLOCK_FOOTER_LENGTH = 8;\r
2058 const int MAX_BLOCK_SIZE      = 65536;\r
2059 const int DEFAULT_BLOCK_SIZE  = 65536;\r
2060 \r
2061 struct BgzfData {\r
2062 \r
2063     // data members\r
2064     unsigned int UncompressedBlockSize;\r
2065     unsigned int CompressedBlockSize;\r
2066     unsigned int BlockLength;\r
2067     unsigned int BlockOffset;\r
2068     uint64_t BlockAddress;\r
2069     bool     IsOpen;\r
2070     bool     IsWriteOnly;\r
2071     FILE*    Stream;\r
2072     char*    UncompressedBlock;\r
2073     char*    CompressedBlock;\r
2074 \r
2075     // constructor & destructor\r
2076     BgzfData(void);\r
2077     ~BgzfData(void);\r
2078 \r
2079     // closes BGZF file\r
2080     void Close(void);\r
2081     // compresses the current block\r
2082     int DeflateBlock(void);\r
2083     // flushes the data in the BGZF block\r
2084     void FlushBlock(void);\r
2085     // de-compresses the current block\r
2086     int InflateBlock(const int& blockLength);\r
2087     // opens the BGZF file for reading (mode is either "rb" for reading, or "wb" for writing\r
2088     void Open(const std::string& filename, const char* mode);\r
2089     // reads BGZF data into a byte buffer\r
2090     int Read(char* data, const unsigned int dataLength);\r
2091     // reads BGZF block\r
2092     int ReadBlock(void);\r
2093     // seek to position in BAM file\r
2094     bool Seek(int64_t position);\r
2095     // get file position in BAM file\r
2096     int64_t Tell(void);\r
2097     // writes the supplied data into the BGZF buffer\r
2098     unsigned int Write(const char* data, const unsigned int dataLen);\r
2099 \r
2100     // checks BGZF block header\r
2101     static inline bool CheckBlockHeader(char* header);\r
2102     // packs an unsigned integer into the specified buffer\r
2103     static inline void PackUnsignedInt(char* buffer, unsigned int value);\r
2104     // packs an unsigned short into the specified buffer\r
2105     static inline void PackUnsignedShort(char* buffer, unsigned short value);\r
2106     // unpacks a buffer into a signed int\r
2107     static inline signed int UnpackSignedInt(char* buffer);\r
2108     // unpacks a buffer into a unsigned int\r
2109     static inline unsigned int UnpackUnsignedInt(char* buffer);\r
2110     // unpacks a buffer into a unsigned short\r
2111     static inline unsigned short UnpackUnsignedShort(char* buffer);\r
2112 };\r
2113 \r
2114 // -------------------------------------------------------------\r
2115 \r
2116 inline\r
2117 bool BgzfData::CheckBlockHeader(char* header) {\r
2118 \r
2119     return (header[0] == GZIP_ID1 &&\r
2120             header[1] == (char)GZIP_ID2 &&\r
2121             header[2] == Z_DEFLATED &&\r
2122             (header[3] & FLG_FEXTRA) != 0 &&\r
2123             BgzfData::UnpackUnsignedShort(&header[10]) == BGZF_XLEN &&\r
2124             header[12] == BGZF_ID1 &&\r
2125             header[13] == BGZF_ID2 &&\r
2126             BgzfData::UnpackUnsignedShort(&header[14]) == BGZF_LEN );\r
2127 }\r
2128 \r
2129 // packs an unsigned integer into the specified buffer\r
2130 inline\r
2131 void BgzfData::PackUnsignedInt(char* buffer, unsigned int value) {\r
2132     buffer[0] = (char)value;\r
2133     buffer[1] = (char)(value >> 8);\r
2134     buffer[2] = (char)(value >> 16);\r
2135     buffer[3] = (char)(value >> 24);\r
2136 }\r
2137 \r
2138 // packs an unsigned short into the specified buffer\r
2139 inline\r
2140 void BgzfData::PackUnsignedShort(char* buffer, unsigned short value) {\r
2141     buffer[0] = (char)value;\r
2142     buffer[1] = (char)(value >> 8);\r
2143 }\r
2144 \r
2145 // unpacks a buffer into a signed int\r
2146 inline\r
2147 signed int BgzfData::UnpackSignedInt(char* buffer) {\r
2148     union { signed int value; unsigned char valueBuffer[sizeof(signed int)]; } un;\r
2149     un.value = 0;\r
2150     un.valueBuffer[0] = buffer[0];\r
2151     un.valueBuffer[1] = buffer[1];\r
2152     un.valueBuffer[2] = buffer[2];\r
2153     un.valueBuffer[3] = buffer[3];\r
2154     return un.value;\r
2155 }\r
2156 \r
2157 // unpacks a buffer into an unsigned int\r
2158 inline\r
2159 unsigned int BgzfData::UnpackUnsignedInt(char* buffer) {\r
2160     union { unsigned int value; unsigned char valueBuffer[sizeof(unsigned int)]; } un;\r
2161     un.value = 0;\r
2162     un.valueBuffer[0] = buffer[0];\r
2163     un.valueBuffer[1] = buffer[1];\r
2164     un.valueBuffer[2] = buffer[2];\r
2165     un.valueBuffer[3] = buffer[3];\r
2166     return un.value;\r
2167 }\r
2168 \r
2169 // unpacks a buffer into an unsigned short\r
2170 inline\r
2171 unsigned short BgzfData::UnpackUnsignedShort(char* buffer) {\r
2172     union { unsigned short value; unsigned char valueBuffer[sizeof(unsigned short)];} un;\r
2173     un.value = 0;\r
2174     un.valueBuffer[0] = buffer[0];\r
2175     un.valueBuffer[1] = buffer[1];\r
2176     return un.value;\r
2177 }\r
2178 \r
2179 } // namespace BamTools\r
2180 \r
2181 #endif // BGZF_H\r
2182 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0Makefile\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000001044\011307517214\0012376\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0CXX=         g++\r
2183 CXXFLAGS=       -Wall -O3\r
2184 PROG=           BamConversion BamDump BamTrim\r
2185 LIBS=           -lz\r
2186 \r
2187 all: $(PROG)\r
2188 \r
2189 BamConversion: BGZF.o BamReader.o BamWriter.o BamConversionMain.o\r
2190         $(CXX) $(CXXFLAGS) -o $@ BGZF.o BamReader.o BamWriter.o BamConversionMain.o $(LIBS)\r
2191 \r
2192 BamDump: BGZF.o BamReader.o BamDumpMain.o\r
2193         $(CXX) $(CXXFLAGS) -o $@ BGZF.o BamReader.o BamDumpMain.o $(LIBS)\r
2194 \r
2195 BamTrim: BGZF.o BamReader.o BamWriter.o BamTrimMain.o\r
2196         $(CXX) $(CXXFLAGS) -o $@ BGZF.o BamReader.o BamWriter.o BamTrimMain.o $(LIBS)\r
2197 \r
2198 clean:\r
2199         rm -fr gmon.out *.o *.a a.out *~\r
2200 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0README\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\00000664\00001074\00001074\000000003607\011307522302\0011617\0 0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0ustar  \0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0barnett\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0------------------------------------------------------------
2201 README : BAMTOOLS
2202 ------------------------------------------------------------
2203
2204 BamTools: a C++ API for reading/writing BAM files.
2205
2206 I.   Introduction
2207 II.  Usage
2208 III. Contact
2209
2210 ------------------------------------------------------------
2211
2212 I. Introduction:
2213
2214 The API consists of 2 main modules - BamReader and BamWriter. As you would expect,
2215 BamReader provides read-access to BAM files, while BamWriter does the writing of BAM
2216 files. BamReader provides an interface for random-access (jumping) in a BAM file,
2217 as well as generating BAM index files.
2218
2219 An additional file, BamAux.h, is included as well.
2220 This file contains the common data structures and typedefs used throught the API.
2221
2222 BGZF.h & BGZF.cpp contain our implementation of the Broad Institute's
2223 BGZF compression format.
2224
2225 BamConversion, BamDump, and BamTrim are 3 'toy' examples on how the API can be used.
2226
2227 ------------------------------------------------------------
2228
2229 II. Usage :
2230
2231 To use this API, you simply need to do 3 things:
2232
2233     1 - Drop the BamTools files somewhere the compiler can find them.
2234         (i.e. in your source tree, or somewhere else in your include path)
2235
2236     2 - Import BamTools API with the following lines of code
2237         #include "BamReader.h"     // as needed
2238         #include "BamWriter.h"     // as needed
2239         using namespace BamTools;
2240
2241     3 - Compile with '-lz' ('l' as in Lima) to access ZLIB compression library
2242             (For VS users, I can provide you zlib headers - just contact me).
2243
2244 See any included programs and Makefile for more specific compiling/usage examples.
2245 See comments in the header files for more detailed API documentation.
2246
2247 ------------------------------------------------------------
2248
2249 III. Contact :
2250
2251 Feel free to contact me with any questions, comments, suggestions, bug reports, etc.
2252   - Derek Barnett
2253
2254 http://sourceforge.net/projects/bamtools
2255 \0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0