bamtools_convert.cpp

   1 // ***************************************************************************
   2 // bamtools_convert.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 9 July 2010
   7 // ---------------------------------------------------------------------------
   8 // Converts between BAM and a number of other formats
   9 // ***************************************************************************
  10
  11 #include <fstream>
  12 #include <iostream>
  13 #include <sstream>
  14 #include <string>
  15 #include <vector>
  16
  17 #include "bamtools_convert.h"
  18 #include "bamtools_options.h"
  19 #include "bamtools_pileup.h"
  20 #include "bamtools_utilities.h"
  21 #include "BGZF.h"
  22 #include "BamReader.h"
  23 #include "BamMultiReader.h"
  24
  25 using namespace std;
  26 using namespace BamTools;
  27
  28 namespace BamTools {
  29
  30     // format names
  31     static const string FORMAT_BED      = "bed";
  32     static const string FORMAT_BEDGRAPH = "bedgraph";
  33     static const string FORMAT_FASTA    = "fasta";
  34     static const string FORMAT_FASTQ    = "fastq";
  35     static const string FORMAT_JSON     = "json";
  36     static const string FORMAT_SAM      = "sam";
  37     static const string FORMAT_PILEUP   = "pileup";
  38     static const string FORMAT_WIGGLE   = "wig";
  39
  40     // other constants
  41     static const unsigned int FASTA_LINE_MAX = 50;
  42
  43 } // namespace BamTools
  44
  45 struct ConvertTool::ConvertToolPrivate {
  46
  47     // ctor & dtor
  48     public:
  49         ConvertToolPrivate(ConvertTool::ConvertSettings* settings);
  50         ~ConvertToolPrivate(void);
  51
  52     // interface
  53     public:
  54         bool Run(void);
  55
  56     // internal methods
  57     private:
  58         void PrintBed(const BamAlignment& a);
  59         void PrintBedGraph(const BamAlignment& a);
  60         void PrintFasta(const BamAlignment& a);
  61         void PrintFastq(const BamAlignment& a);
  62         void PrintJson(const BamAlignment& a);
  63         void PrintSam(const BamAlignment& a);
  64         void PrintWiggle(const BamAlignment& a);
  65
  66     // data members
  67     private:
  68         ConvertTool::ConvertSettings* m_settings;
  69         RefVector m_references;
  70         ostream m_out;
  71 };
  72
  73 // ---------------------------------------------
  74 // ConvertSettings implementation
  75
  76 struct ConvertTool::ConvertSettings {
  77
  78     // flags
  79     bool HasInputFilenames;
  80     bool HasOutputFilename;
  81     bool HasFormat;
  82     bool HasRegion;
  83
  84     // pileup flags
  85     bool HasFastaFilename;
  86     bool IsPrintingPileupMapQualities;
  87
  88     // options
  89     vector<string> InputFiles;
  90     string OutputFilename;
  91     string Format;
  92     string Region;
  93
  94     // pileup options
  95     string FastaFilename;
  96
  97     // constructor
  98     ConvertSettings(void)
  99         : HasInputFilenames(false)
 100         , HasOutputFilename(false)
 101         , HasFormat(false)
 102         , HasRegion(false)
 103         , HasFastaFilename(false)
 104         , IsPrintingPileupMapQualities(false)
 105         , OutputFilename(Options::StandardOut())
 106     { }
 107 };
 108
 109 // ---------------------------------------------
 110 // ConvertTool implementation
 111
 112 ConvertTool::ConvertTool(void)
 113     : AbstractTool()
 114     , m_settings(new ConvertSettings)
 115     , m_impl(0)
 116 {
 117     // set program details
 118     Options::SetProgramInfo("bamtools convert", "converts between BAM and a number of other formats", "-in <filename> -out <filename> -format <FORMAT>");
 119
 120     // set up options
 121     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 122     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInputFilenames,  m_settings->InputFiles,     IO_Opts, Options::StandardIn());
 123     Options::AddValueOption("-out",    "BAM filename", "the output BAM file",   "", m_settings->HasOutputFilename,  m_settings->OutputFilename, IO_Opts, Options::StandardOut());
 124     Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
 125
 126     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
 127     Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is read automatically if it exists as <filename>.bai. See \'bamtools help index\' for more  details on creating one", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
 128
 129     OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
 130     Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts, "");
 131     Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
 132 }
 133
 134 ConvertTool::~ConvertTool(void) {
 135     delete m_settings;
 136     m_settings = 0;
 137
 138     delete m_impl;
 139     m_impl = 0;
 140 }
 141
 142 int ConvertTool::Help(void) {
 143     Options::DisplayHelp();
 144     return 0;
 145 }
 146
 147 int ConvertTool::Run(int argc, char* argv[]) {
 148
 149     // parse command line arguments
 150     Options::Parse(argc, argv, 1);
 151
 152     // run internal ConvertTool implementation, return success/fail
 153     m_impl = new ConvertToolPrivate(m_settings);
 154
 155     if ( m_impl->Run() )
 156         return 0;
 157     else
 158         return 1;
 159 }
 160
 161 // ---------------------------------------------
 162 // ConvertToolPrivate implementation
 163
 164 ConvertTool::ConvertToolPrivate::ConvertToolPrivate(ConvertTool::ConvertSettings* settings)
 165     : m_settings(settings)
 166     , m_out(cout.rdbuf()) // default output to cout
 167 { }
 168
 169 ConvertTool::ConvertToolPrivate::~ConvertToolPrivate(void) { }
 170
 171 bool ConvertTool::ConvertToolPrivate::Run(void) {
 172
 173     bool convertedOk = true;
 174
 175     // ------------------------------------
 176     // initialize conversion input/output
 177
 178     // set to default input if none provided
 179     if ( !m_settings->HasInputFilenames )
 180         m_settings->InputFiles.push_back(Options::StandardIn());
 181
 182     // open input files
 183     BamMultiReader reader;
 184     reader.Open(m_settings->InputFiles);
 185     m_references = reader.GetReferenceData();
 186
 187     // set region if specified
 188     BamRegion region;
 189     if ( m_settings->HasRegion ) {
 190         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 191             if ( !reader.SetRegion(region) )
 192               cerr << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
 193         }
 194     }
 195
 196     // if output file given
 197     ofstream outFile;
 198     if ( m_settings->HasOutputFilename ) {
 199
 200         // open output file stream
 201         outFile.open(m_settings->OutputFilename.c_str());
 202         if ( !outFile ) {
 203             cerr << "Could not open " << m_settings->OutputFilename << " for output" << endl;
 204             return false;
 205         }
 206
 207         // set m_out to file's streambuf
 208         m_out.rdbuf(outFile.rdbuf());
 209     }
 210
 211     // ------------------------
 212     // pileup is special case
 213     if ( m_settings->Format == FORMAT_PILEUP ) {
 214
 215         // initialize pileup input/output
 216         Pileup pileup(&reader, &m_out);
 217
 218         // ---------------------------
 219         // configure pileup settings
 220
 221         if ( m_settings->HasRegion )
 222             pileup.SetRegion(region);
 223
 224         if ( m_settings->HasFastaFilename )
 225             pileup.SetFastaFilename(m_settings->FastaFilename);
 226
 227         pileup.SetIsPrintingMapQualities( m_settings->IsPrintingPileupMapQualities );
 228
 229         // run pileup
 230         convertedOk = pileup.Run();
 231     }
 232
 233     // -------------------------------------
 234     // else determine 'simpler' format type
 235     else {
 236
 237         bool formatError = false;
 238         void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
 239         if      ( m_settings->Format == FORMAT_BED )      pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
 240         else if ( m_settings->Format == FORMAT_BEDGRAPH ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBedGraph;
 241         else if ( m_settings->Format == FORMAT_FASTA )    pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
 242         else if ( m_settings->Format == FORMAT_FASTQ )    pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
 243         else if ( m_settings->Format == FORMAT_JSON )     pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
 244         else if ( m_settings->Format == FORMAT_SAM )      pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
 245         else if ( m_settings->Format == FORMAT_WIGGLE )   pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintWiggle;
 246         else {
 247             cerr << "Unrecognized format: " << m_settings->Format << endl;
 248             cerr << "Please see help|README (?) for details on supported formats " << endl;
 249             formatError = true;
 250             convertedOk = false;
 251         }
 252
 253         // ------------------------
 254         // do conversion
 255         if ( !formatError ) {
 256             BamAlignment a;
 257             while ( reader.GetNextAlignment(a) ) {
 258                 (this->*pFunction)(a);
 259             }
 260         }
 261     }
 262
 263     // ------------------------
 264     // clean up & exit
 265     reader.Close();
 266     if ( m_settings->HasOutputFilename ) outFile.close();
 267     return convertedOk;
 268 }
 269
 270 // ----------------------------------------------------------
 271 // Conversion/output methods
 272 // ----------------------------------------------------------
 273
 274 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a)      {
 275
 276     // tab-delimited, 0-based half-open
 277     // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
 278     // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
 279
 280     m_out << m_references.at(a.RefID).RefName << "\t"
 281           << a.Position << "\t"
 282           << a.GetEndPosition() + 1 << "\t"
 283           << a.Name << "\t"
 284           << a.MapQuality << "\t"
 285           << (a.IsReverseStrand() ? "-" : "+") << endl;
 286 }
 287
 288 void ConvertTool::ConvertToolPrivate::PrintBedGraph(const BamAlignment& a) {
 289     ;
 290 }
 291
 292 // print BamAlignment in FASTA format
 293 // N.B. - uses QueryBases NOT AlignedBases
 294 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
 295
 296     // >BamAlignment.Name
 297     // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
 298     // ...
 299
 300     // print header
 301     m_out << "> " << a.Name << endl;
 302
 303     // if sequence fits on single line
 304     if ( a.QueryBases.length() <= FASTA_LINE_MAX )
 305         m_out << a.QueryBases << endl;
 306
 307     // else split over multiple lines
 308     else {
 309
 310         size_t position = 0;
 311         size_t seqLength = a.QueryBases.length();
 312
 313         // write subsequences to each line
 314         while ( position < (seqLength - FASTA_LINE_MAX) ) {
 315             m_out << a.QueryBases.substr(position, FASTA_LINE_MAX) << endl;
 316             position += FASTA_LINE_MAX;
 317         }
 318
 319         // write final subsequence
 320         m_out << a.QueryBases.substr(position) << endl;
 321     }
 322 }
 323
 324 // print BamAlignment in FASTQ format
 325 // N.B. - uses QueryBases NOT AlignedBases
 326 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
 327
 328     // @BamAlignment.Name
 329     // BamAlignment.QueryBases
 330     // +
 331     // BamAlignment.Qualities
 332
 333     m_out << "@" << a.Name << endl
 334           << a.QueryBases   << endl
 335           << "+"            << endl
 336           << a.Qualities    << endl;
 337 }
 338
 339 // print BamAlignment in JSON format
 340 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
 341
 342     // write name & alignment flag
 343     m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
 344
 345     // write reference name
 346     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 347         m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
 348
 349     // write position & map quality
 350     m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
 351
 352     // write CIGAR
 353     const vector<CigarOp>& cigarData = a.CigarData;
 354     if ( !cigarData.empty() ) {
 355         m_out << "\"cigar\":[";
 356         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 357         vector<CigarOp>::const_iterator cigarIter = cigarBegin;
 358         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 359         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 360             const CigarOp& op = (*cigarIter);
 361             if (cigarIter != cigarBegin) m_out << ",";
 362             m_out << "\"" << op.Length << op.Type << "\"";
 363         }
 364         m_out << "],";
 365     }
 366
 367     // write mate reference name, mate position, & insert size
 368     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 369         m_out << "\"mate\":{"
 370               << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
 371               << "\"position\":" << a.MatePosition+1
 372               << ",\"insertSize\":" << a.InsertSize << "},";
 373     }
 374
 375     // write sequence
 376     if ( !a.QueryBases.empty() )
 377         m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
 378
 379     // write qualities
 380     if ( !a.Qualities.empty() )
 381         m_out << "\"qualities\":\"" << a.Qualities << "\",";
 382
 383     // write tag data
 384     const char* tagData = a.TagData.c_str();
 385     const size_t tagDataLength = a.TagData.length();
 386     size_t index = 0;
 387     if (index < tagDataLength) {
 388
 389         m_out << "\"tags\":{";
 390
 391         while ( index < tagDataLength ) {
 392
 393             if (index > 0)
 394                 m_out << ",";
 395
 396             // write tag name
 397             m_out << "\"" << a.TagData.substr(index, 2) << "\":";
 398             index += 2;
 399
 400             // get data type
 401             char type = a.TagData.at(index);
 402             ++index;
 403
 404             switch (type) {
 405                 case('A') :
 406                     m_out << "\"" << tagData[index] << "\"";
 407                     ++index;
 408                     break;
 409
 410                 case('C') :
 411                     m_out << (int)tagData[index];
 412                     ++index;
 413                     break;
 414
 415                 case('c') :
 416                     m_out << (int)tagData[index];
 417                     ++index;
 418                     break;
 419
 420                 case('S') :
 421                     m_out << BgzfData::UnpackUnsignedShort(&tagData[index]);
 422                     index += 2;
 423                     break;
 424
 425                 case('s') :
 426                     m_out << BgzfData::UnpackSignedShort(&tagData[index]);
 427                     index += 2;
 428                     break;
 429
 430                 case('I') :
 431                     m_out << BgzfData::UnpackUnsignedInt(&tagData[index]);
 432                     index += 4;
 433                     break;
 434
 435                 case('i') :
 436                     m_out << BgzfData::UnpackSignedInt(&tagData[index]);
 437                     index += 4;
 438                     break;
 439
 440                 case('f') :
 441                     m_out << BgzfData::UnpackFloat(&tagData[index]);
 442                     index += 4;
 443                     break;
 444
 445                 case('d') :
 446                     m_out << BgzfData::UnpackDouble(&tagData[index]);
 447                     index += 8;
 448                     break;
 449
 450                 case('Z') :
 451                 case('H') :
 452                     m_out << "\"";
 453                     while (tagData[index]) {
 454                         m_out << tagData[index];
 455                         ++index;
 456                     }
 457                     m_out << "\"";
 458                     ++index;
 459                     break;
 460             }
 461
 462             if ( tagData[index] == '\0')
 463                 break;
 464         }
 465
 466         m_out << "}";
 467     }
 468
 469     m_out << "}" << endl;
 470
 471 }
 472
 473 // print BamAlignment in SAM format
 474 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
 475
 476     // tab-delimited
 477     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 478
 479     // write name & alignment flag
 480     m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
 481
 482     // write reference name
 483     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 484         m_out << m_references[a.RefID].RefName << "\t";
 485     else
 486         m_out << "*\t";
 487
 488     // write position & map quality
 489     m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
 490
 491     // write CIGAR
 492     const vector<CigarOp>& cigarData = a.CigarData;
 493     if ( cigarData.empty() ) m_out << "*\t";
 494     else {
 495         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 496         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 497         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 498             const CigarOp& op = (*cigarIter);
 499             m_out << op.Length << op.Type;
 500         }
 501         m_out << "\t";
 502     }
 503
 504     // write mate reference name, mate position, & insert size
 505     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 506         if ( a.MateRefID == a.RefID ) m_out << "=\t";
 507         else m_out << m_references[a.MateRefID].RefName << "\t";
 508         m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 509     }
 510     else m_out << "*\t0\t0\t";
 511
 512     // write sequence
 513     if ( a.QueryBases.empty() ) m_out << "*\t";
 514     else m_out << a.QueryBases << "\t";
 515
 516     // write qualities
 517     if ( a.Qualities.empty() ) m_out << "*";
 518     else m_out << a.Qualities;
 519
 520     // write tag data
 521     const char* tagData = a.TagData.c_str();
 522     const size_t tagDataLength = a.TagData.length();
 523
 524     size_t index = 0;
 525     while ( index < tagDataLength ) {
 526
 527         // write tag name
 528         m_out << "\t" << a.TagData.substr(index, 2) << ":";
 529         index += 2;
 530
 531         // get data type
 532         char type = a.TagData.at(index);
 533         ++index;
 534
 535         switch (type) {
 536             case('A') :
 537                 m_out << "A:" << tagData[index];
 538                 ++index;
 539                 break;
 540
 541             case('C') :
 542                 m_out << "i:" << (int)tagData[index];
 543                 ++index;
 544                 break;
 545
 546             case('c') :
 547                 m_out << "i:" << (int)tagData[index];
 548                 ++index;
 549                 break;
 550
 551             case('S') :
 552                 m_out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
 553                 index += 2;
 554                 break;
 555
 556             case('s') :
 557                 m_out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
 558                 index += 2;
 559                 break;
 560
 561             case('I') :
 562                 m_out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
 563                 index += 4;
 564                 break;
 565
 566             case('i') :
 567                 m_out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
 568                 index += 4;
 569                 break;
 570
 571             case('f') :
 572                 m_out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
 573                 index += 4;
 574                 break;
 575
 576             case('d') :
 577                 m_out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
 578                 index += 8;
 579                 break;
 580
 581             case('Z') :
 582             case('H') :
 583                 m_out << type << ":";
 584                 while (tagData[index]) {
 585                     m_out << tagData[index];
 586                     ++index;
 587                 }
 588                 ++index;
 589                 break;
 590         }
 591
 592         if ( tagData[index] == '\0')
 593             break;
 594     }
 595
 596     m_out << endl;
 597 }
 598
 599 void ConvertTool::ConvertToolPrivate::PrintWiggle(const BamAlignment& a) {
 600     ;
 601 }