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