src/api/internal/BamWriter_p.cpp

   1 // ***************************************************************************
   2 // BamWriter_p.cpp (c) 2010 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 6 October 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides the basic functionality for producing BAM files
   8 // ***************************************************************************
   9
  10 #include <api/BamAlignment.h>
  11 #include <api/BamConstants.h>
  12 #include <api/internal/BamException_p.h>
  13 #include <api/internal/BamWriter_p.h>
  14 using namespace BamTools;
  15 using namespace BamTools::Internal;
  16
  17 #include <cstdlib>
  18 #include <cstring>
  19 using namespace std;
  20
  21 // ctor
  22 BamWriterPrivate::BamWriterPrivate(void)
  23     : m_isBigEndian( BamTools::SystemIsBigEndian() )
  24 { }
  25
  26 // dtor
  27 BamWriterPrivate::~BamWriterPrivate(void) {
  28     Close();
  29 }
  30
  31 // calculates minimum bin for a BAM alignment interval
  32 uint32_t BamWriterPrivate::CalculateMinimumBin(const int begin, int end) const {
  33     --end;
  34     if ( (begin >> 14) == (end >> 14) ) return 4681 + (begin >> 14);
  35     if ( (begin >> 17) == (end >> 17) ) return  585 + (begin >> 17);
  36     if ( (begin >> 20) == (end >> 20) ) return   73 + (begin >> 20);
  37     if ( (begin >> 23) == (end >> 23) ) return    9 + (begin >> 23);
  38     if ( (begin >> 26) == (end >> 26) ) return    1 + (begin >> 26);
  39     return 0;
  40 }
  41
  42 // closes the alignment archive
  43 void BamWriterPrivate::Close(void) {
  44
  45     // skip if file not open
  46     if ( !IsOpen() ) return;
  47
  48     // close output stream
  49     try {
  50         m_stream.Close();
  51     } catch ( BamException& e ) {
  52         m_errorString = e.what();
  53     }
  54 }
  55
  56 // creates a cigar string from the supplied alignment
  57 void BamWriterPrivate::CreatePackedCigar(const vector<CigarOp>& cigarOperations, string& packedCigar) {
  58
  59     // initialize
  60     const size_t numCigarOperations = cigarOperations.size();
  61     packedCigar.resize(numCigarOperations * Constants::BAM_SIZEOF_INT);
  62
  63     // pack the cigar data into the string
  64     unsigned int* pPackedCigar = (unsigned int*)packedCigar.data();
  65
  66     // iterate over cigar operations
  67     vector<CigarOp>::const_iterator coIter = cigarOperations.begin();
  68     vector<CigarOp>::const_iterator coEnd  = cigarOperations.end();
  69     for ( ; coIter != coEnd; ++coIter ) {
  70
  71         // store op in packedCigar
  72         uint8_t cigarOp;
  73         switch ( coIter->Type ) {
  74             case (Constants::BAM_CIGAR_MATCH_CHAR)    : cigarOp = Constants::BAM_CIGAR_MATCH;    break;
  75             case (Constants::BAM_CIGAR_INS_CHAR)      : cigarOp = Constants::BAM_CIGAR_INS;      break;
  76             case (Constants::BAM_CIGAR_DEL_CHAR)      : cigarOp = Constants::BAM_CIGAR_DEL;      break;
  77             case (Constants::BAM_CIGAR_REFSKIP_CHAR)  : cigarOp = Constants::BAM_CIGAR_REFSKIP;  break;
  78             case (Constants::BAM_CIGAR_SOFTCLIP_CHAR) : cigarOp = Constants::BAM_CIGAR_SOFTCLIP; break;
  79             case (Constants::BAM_CIGAR_HARDCLIP_CHAR) : cigarOp = Constants::BAM_CIGAR_HARDCLIP; break;
  80             case (Constants::BAM_CIGAR_PAD_CHAR)      : cigarOp = Constants::BAM_CIGAR_PAD;      break;
  81             case (Constants::BAM_CIGAR_SEQMATCH_CHAR) : cigarOp = Constants::BAM_CIGAR_SEQMATCH; break;
  82             case (Constants::BAM_CIGAR_MISMATCH_CHAR) : cigarOp = Constants::BAM_CIGAR_MISMATCH; break;
  83             default:
  84                 const string message = string("invalid CIGAR operation type") + coIter->Type;
  85                 throw BamException("BamWriter::CreatePackedCigar", message);
  86         }
  87
  88         *pPackedCigar = coIter->Length << Constants::BAM_CIGAR_SHIFT | cigarOp;
  89         pPackedCigar++;
  90     }
  91 }
  92
  93 // encodes the supplied query sequence into 4-bit notation
  94 void BamWriterPrivate::EncodeQuerySequence(const string& query, string& encodedQuery) {
  95
  96     // prepare the encoded query string
  97     const size_t queryLength = query.size();
  98     const size_t encodedQueryLength = static_cast<size_t>((queryLength+1)/2);
  99     encodedQuery.resize(encodedQueryLength);
 100     char* pEncodedQuery = (char*)encodedQuery.data();
 101     const char* pQuery = (const char*)query.data();
 102
 103     // walk through original query sequence, encoding its bases
 104     unsigned char nucleotideCode;
 105     bool useHighWord = true;
 106     while ( *pQuery ) {
 107         switch ( *pQuery ) {
 108             case (Constants::BAM_DNA_EQUAL) : nucleotideCode = Constants::BAM_BASECODE_EQUAL; break;
 109             case (Constants::BAM_DNA_A)     : nucleotideCode = Constants::BAM_BASECODE_A;     break;
 110             case (Constants::BAM_DNA_C)     : nucleotideCode = Constants::BAM_BASECODE_C;     break;
 111             case (Constants::BAM_DNA_M)     : nucleotideCode = Constants::BAM_BASECODE_M;     break;
 112             case (Constants::BAM_DNA_G)     : nucleotideCode = Constants::BAM_BASECODE_G;     break;
 113             case (Constants::BAM_DNA_R)     : nucleotideCode = Constants::BAM_BASECODE_R;     break;
 114             case (Constants::BAM_DNA_S)     : nucleotideCode = Constants::BAM_BASECODE_S;     break;
 115             case (Constants::BAM_DNA_V)     : nucleotideCode = Constants::BAM_BASECODE_V;     break;
 116             case (Constants::BAM_DNA_T)     : nucleotideCode = Constants::BAM_BASECODE_T;     break;
 117             case (Constants::BAM_DNA_W)     : nucleotideCode = Constants::BAM_BASECODE_W;     break;
 118             case (Constants::BAM_DNA_Y)     : nucleotideCode = Constants::BAM_BASECODE_Y;     break;
 119             case (Constants::BAM_DNA_H)     : nucleotideCode = Constants::BAM_BASECODE_H;     break;
 120             case (Constants::BAM_DNA_K)     : nucleotideCode = Constants::BAM_BASECODE_K;     break;
 121             case (Constants::BAM_DNA_D)     : nucleotideCode = Constants::BAM_BASECODE_D;     break;
 122             case (Constants::BAM_DNA_B)     : nucleotideCode = Constants::BAM_BASECODE_B;     break;
 123             case (Constants::BAM_DNA_N)     : nucleotideCode = Constants::BAM_BASECODE_N;     break;
 124             default:
 125                 const string message = string("invalid base: ") + *pQuery;
 126                 throw BamException("BamWriter::EncodeQuerySequence", message);
 127         }
 128
 129         // pack the nucleotide code
 130         if ( useHighWord ) {
 131             *pEncodedQuery = nucleotideCode << 4;
 132             useHighWord = false;
 133         } else {
 134             *pEncodedQuery |= nucleotideCode;
 135             ++pEncodedQuery;
 136             useHighWord = true;
 137         }
 138
 139         // increment the query position
 140         ++pQuery;
 141     }
 142 }
 143
 144 // returns a description of the last error that occurred
 145 std::string BamWriterPrivate::GetErrorString(void) const {
 146     return m_errorString;
 147 }
 148
 149 // returns whether BAM file is open for writing or not
 150 bool BamWriterPrivate::IsOpen(void) const {
 151     return m_stream.IsOpen;
 152 }
 153
 154 // opens the alignment archive
 155 bool BamWriterPrivate::Open(const string& filename,
 156                             const string& samHeaderText,
 157                             const RefVector& referenceSequences)
 158 {
 159     try {
 160
 161         // open the BGZF file for writing, return failure if error
 162         m_stream.Open(filename, "wb");
 163
 164         // write BAM file 'metadata' components
 165         WriteMagicNumber();
 166         WriteSamHeaderText(samHeaderText);
 167         WriteReferences(referenceSequences);
 168
 169         // return success
 170         return true;
 171
 172     } catch ( BamException& e ) {
 173         m_errorString = e.what();
 174         return false;
 175     }
 176 }
 177
 178 // saves the alignment to the alignment archive
 179 bool BamWriterPrivate::SaveAlignment(const BamAlignment& al) {
 180
 181     try {
 182
 183         // if BamAlignment contains only the core data and a raw char data buffer
 184         // (as a result of BamReader::GetNextAlignmentCore())
 185         if ( al.SupportData.HasCoreOnly )
 186             WriteCoreAlignment(al);
 187
 188         // otherwise, BamAlignment should contain character in the standard fields: Name, QueryBases, etc
 189         // (resulting from BamReader::GetNextAlignment() *OR* being generated directly by client code)
 190         else WriteAlignment(al);
 191
 192         // if we get here, everything OK
 193         return true;
 194
 195     } catch ( BamException& e ) {
 196         m_errorString = e.what();
 197         return false;
 198     }
 199 }
 200
 201 void BamWriterPrivate::SetWriteCompressed(bool ok) {
 202     // modifying compression is not allowed if BAM file is open
 203     if ( !IsOpen() )
 204         m_stream.SetWriteCompressed(ok);
 205 }
 206
 207 void BamWriterPrivate::WriteAlignment(const BamAlignment& al) {
 208
 209     // calculate char lengths
 210     const unsigned int nameLength         = al.Name.size() + 1;
 211     const unsigned int numCigarOperations = al.CigarData.size();
 212     const unsigned int queryLength        = al.QueryBases.size();
 213     const unsigned int tagDataLength      = al.TagData.size();
 214
 215     // no way to tell if BamAlignment.Bin is already defined (no default, invalid value)
 216     // force calculation of Bin before storing
 217     const int endPosition = al.GetEndPosition();
 218     const uint32_t alignmentBin = CalculateMinimumBin(al.Position, endPosition);
 219
 220     // create our packed cigar string
 221     string packedCigar;
 222     CreatePackedCigar(al.CigarData, packedCigar);
 223     const unsigned int packedCigarLength = packedCigar.size();
 224
 225     // encode the query
 226     string encodedQuery;
 227     EncodeQuerySequence(al.QueryBases, encodedQuery);
 228     const unsigned int encodedQueryLength = encodedQuery.size();
 229
 230     // write the block size
 231     const unsigned int dataBlockSize = nameLength +
 232                                        packedCigarLength +
 233                                        encodedQueryLength +
 234                                        queryLength +
 235                                        tagDataLength;
 236     unsigned int blockSize = Constants::BAM_CORE_SIZE + dataBlockSize;
 237     if ( m_isBigEndian ) BamTools::SwapEndian_32(blockSize);
 238     m_stream.Write((char*)&blockSize, Constants::BAM_SIZEOF_INT);
 239
 240     // assign the BAM core data
 241     uint32_t buffer[Constants::BAM_CORE_BUFFER_SIZE];
 242     buffer[0] = al.RefID;
 243     buffer[1] = al.Position;
 244     buffer[2] = (alignmentBin << 16) | (al.MapQuality << 8) | nameLength;
 245     buffer[3] = (al.AlignmentFlag << 16) | numCigarOperations;
 246     buffer[4] = queryLength;
 247     buffer[5] = al.MateRefID;
 248     buffer[6] = al.MatePosition;
 249     buffer[7] = al.InsertSize;
 250
 251     // swap BAM core endian-ness, if necessary
 252     if ( m_isBigEndian ) {
 253         for ( int i = 0; i < 8; ++i )
 254             BamTools::SwapEndian_32(buffer[i]);
 255     }
 256
 257     // write the BAM core
 258     m_stream.Write((char*)&buffer, Constants::BAM_CORE_SIZE);
 259
 260     // write the query name
 261     m_stream.Write(al.Name.c_str(), nameLength);
 262
 263     // write the packed cigar
 264     if ( m_isBigEndian ) {
 265         char* cigarData = new char[packedCigarLength]();
 266         memcpy(cigarData, packedCigar.data(), packedCigarLength);
 267         if ( m_isBigEndian ) {
 268             for ( size_t i = 0; i < packedCigarLength; ++i )
 269                 BamTools::SwapEndian_32p(&cigarData[i]);
 270         }
 271         m_stream.Write(cigarData, packedCigarLength);
 272         delete[] cigarData; // TODO: cleanup on Write exception thrown?
 273     }
 274     else
 275         m_stream.Write(packedCigar.data(), packedCigarLength);
 276
 277     // write the encoded query sequence
 278     m_stream.Write(encodedQuery.data(), encodedQueryLength);
 279
 280     // write the base qualities
 281     char* pBaseQualities = (char*)al.Qualities.data();
 282     for ( size_t i = 0; i < queryLength; ++i )
 283         pBaseQualities[i] -= 33; // FASTQ conversion
 284     m_stream.Write(pBaseQualities, queryLength);
 285
 286     // write the read group tag
 287     if ( m_isBigEndian ) {
 288
 289         char* tagData = new char[tagDataLength]();
 290         memcpy(tagData, al.TagData.data(), tagDataLength);
 291
 292         size_t i = 0;
 293         while ( i < tagDataLength ) {
 294
 295             i += Constants::BAM_TAG_TAGSIZE;  // skip tag chars (e.g. "RG", "NM", etc.)
 296             const char type = tagData[i];     // get tag type at position i
 297             ++i;
 298
 299             switch ( type ) {
 300
 301                 case(Constants::BAM_TAG_TYPE_ASCII) :
 302                 case(Constants::BAM_TAG_TYPE_INT8)  :
 303                 case(Constants::BAM_TAG_TYPE_UINT8) :
 304                     ++i;
 305                     break;
 306
 307                 case(Constants::BAM_TAG_TYPE_INT16)  :
 308                 case(Constants::BAM_TAG_TYPE_UINT16) :
 309                     BamTools::SwapEndian_16p(&tagData[i]);
 310                     i += sizeof(uint16_t);
 311                     break;
 312
 313                 case(Constants::BAM_TAG_TYPE_FLOAT)  :
 314                 case(Constants::BAM_TAG_TYPE_INT32)  :
 315                 case(Constants::BAM_TAG_TYPE_UINT32) :
 316                     BamTools::SwapEndian_32p(&tagData[i]);
 317                     i += sizeof(uint32_t);
 318                     break;
 319
 320                 case(Constants::BAM_TAG_TYPE_HEX) :
 321                 case(Constants::BAM_TAG_TYPE_STRING) :
 322                     // no endian swapping necessary for hex-string/string data
 323                     while ( tagData[i] )
 324                         ++i;
 325                     // increment one more for null terminator
 326                     ++i;
 327                     break;
 328
 329                 case(Constants::BAM_TAG_TYPE_ARRAY) :
 330
 331                 {
 332                     // read array type
 333                     const char arrayType = tagData[i];
 334                     ++i;
 335
 336                     // swap endian-ness of number of elements in place, then retrieve for loop
 337                     BamTools::SwapEndian_32p(&tagData[i]);
 338                     int32_t numElements;
 339                     memcpy(&numElements, &tagData[i], sizeof(uint32_t));
 340                     i += sizeof(uint32_t);
 341
 342                     // swap endian-ness of array elements
 343                     for ( int j = 0; j < numElements; ++j ) {
 344                         switch (arrayType) {
 345                             case (Constants::BAM_TAG_TYPE_INT8)  :
 346                             case (Constants::BAM_TAG_TYPE_UINT8) :
 347                                 // no endian-swapping necessary
 348                                 ++i;
 349                                 break;
 350                             case (Constants::BAM_TAG_TYPE_INT16)  :
 351                             case (Constants::BAM_TAG_TYPE_UINT16) :
 352                                 BamTools::SwapEndian_16p(&tagData[i]);
 353                                 i += sizeof(uint16_t);
 354                                 break;
 355                             case (Constants::BAM_TAG_TYPE_FLOAT)  :
 356                             case (Constants::BAM_TAG_TYPE_INT32)  :
 357                             case (Constants::BAM_TAG_TYPE_UINT32) :
 358                                 BamTools::SwapEndian_32p(&tagData[i]);
 359                                 i += sizeof(uint32_t);
 360                                 break;
 361                             default:
 362                                 delete[] tagData;
 363                                 const string message = string("invalid binary array type: ") + arrayType;
 364                                 throw BamException("BamWriter::SaveAlignment", message);
 365                         }
 366                     }
 367
 368                     break;
 369                 }
 370
 371                 default :
 372                     delete[] tagData;
 373                     const string message = string("invalid tag type: ") + type;
 374                     throw BamException("BamWriter::SaveAlignment", message);
 375             }
 376         }
 377
 378         m_stream.Write(tagData, tagDataLength);
 379         delete[] tagData; // TODO: cleanup on Write exception thrown?
 380     }
 381     else
 382         m_stream.Write(al.TagData.data(), tagDataLength);
 383 }
 384
 385 void BamWriterPrivate::WriteCoreAlignment(const BamAlignment& al) {
 386
 387     // write the block size
 388     unsigned int blockSize = al.SupportData.BlockLength;
 389     if ( m_isBigEndian ) BamTools::SwapEndian_32(blockSize);
 390     m_stream.Write((char*)&blockSize, Constants::BAM_SIZEOF_INT);
 391
 392     // re-calculate bin (in case BamAlignment's position has been previously modified)
 393     const uint32_t alignmentBin = CalculateMinimumBin(al.Position, al.GetEndPosition());
 394
 395     // assign the BAM core data
 396     uint32_t buffer[Constants::BAM_CORE_BUFFER_SIZE];
 397     buffer[0] = al.RefID;
 398     buffer[1] = al.Position;
 399     buffer[2] = (alignmentBin << 16) | (al.MapQuality << 8) | al.SupportData.QueryNameLength;
 400     buffer[3] = (al.AlignmentFlag << 16) | al.SupportData.NumCigarOperations;
 401     buffer[4] = al.SupportData.QuerySequenceLength;
 402     buffer[5] = al.MateRefID;
 403     buffer[6] = al.MatePosition;
 404     buffer[7] = al.InsertSize;
 405
 406     // swap BAM core endian-ness, if necessary
 407     if ( m_isBigEndian ) {
 408         for ( int i = 0; i < 8; ++i )
 409             BamTools::SwapEndian_32(buffer[i]);
 410     }
 411
 412     // write the BAM core
 413     m_stream.Write((char*)&buffer, Constants::BAM_CORE_SIZE);
 414
 415     // write the raw char data
 416     m_stream.Write((char*)al.SupportData.AllCharData.data(),
 417                    al.SupportData.BlockLength-Constants::BAM_CORE_SIZE);
 418 }
 419
 420 void BamWriterPrivate::WriteMagicNumber(void) {
 421     // write BAM file 'magic number'
 422     m_stream.Write(Constants::BAM_HEADER_MAGIC, Constants::BAM_HEADER_MAGIC_LENGTH);
 423 }
 424
 425 void BamWriterPrivate::WriteReferences(const BamTools::RefVector& referenceSequences) {
 426
 427     // write the number of reference sequences
 428     uint32_t numReferenceSequences = referenceSequences.size();
 429     if ( m_isBigEndian ) BamTools::SwapEndian_32(numReferenceSequences);
 430     m_stream.Write((char*)&numReferenceSequences, Constants::BAM_SIZEOF_INT);
 431
 432     // foreach reference sequence
 433     RefVector::const_iterator rsIter = referenceSequences.begin();
 434     RefVector::const_iterator rsEnd  = referenceSequences.end();
 435     for ( ; rsIter != rsEnd; ++rsIter ) {
 436
 437         // write the reference sequence name length
 438         uint32_t referenceSequenceNameLen = rsIter->RefName.size() + 1;
 439         if ( m_isBigEndian ) BamTools::SwapEndian_32(referenceSequenceNameLen);
 440         m_stream.Write((char*)&referenceSequenceNameLen, Constants::BAM_SIZEOF_INT);
 441
 442         // write the reference sequence name
 443         m_stream.Write(rsIter->RefName.c_str(), referenceSequenceNameLen);
 444
 445         // write the reference sequence length
 446         int32_t referenceLength = rsIter->RefLength;
 447         if ( m_isBigEndian ) BamTools::SwapEndian_32(referenceLength);
 448         m_stream.Write((char*)&referenceLength, Constants::BAM_SIZEOF_INT);
 449     }
 450 }
 451
 452 void BamWriterPrivate::WriteSamHeaderText(const std::string& samHeaderText) {
 453
 454     // write the SAM header  text length
 455     uint32_t samHeaderLen = samHeaderText.size();
 456     if ( m_isBigEndian ) BamTools::SwapEndian_32(samHeaderLen);
 457     m_stream.Write((char*)&samHeaderLen, Constants::BAM_SIZEOF_INT);
 458
 459     // write the SAM header text
 460     if ( samHeaderLen > 0 )
 461         m_stream.Write(samHeaderText.data(), samHeaderLen);
 462 }