src/api/internal/BamWriter_p.cpp

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