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 HasInputFilenames;
  80     bool HasOutputFilename;
  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         : HasInputFilenames(false)
 101         , HasOutputFilename(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 between BAM and a number of other formats", "-in <filename> -out <filename> -format <FORMAT>");
 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->HasInputFilenames,  m_settings->InputFiles,     IO_Opts, Options::StandardIn());
 125     Options::AddValueOption("-out",    "BAM filename", "the output BAM file",   "", m_settings->HasOutputFilename,  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->HasInputFilenames )
 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->HasOutputFilename ) {
 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->HasOutputFilename ) 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         m_out << "\"qualities\":\"" << a.Qualities << "\",";
 393
 394     // write tag data
 395     const char* tagData = a.TagData.c_str();
 396     const size_t tagDataLength = a.TagData.length();
 397     size_t index = 0;
 398     if (index < tagDataLength) {
 399
 400         m_out << "\"tags\":{";
 401
 402         while ( index < tagDataLength ) {
 403
 404             if (index > 0)
 405                 m_out << ",";
 406
 407             // write tag name
 408             m_out << "\"" << a.TagData.substr(index, 2) << "\":";
 409             index += 2;
 410
 411             // get data type
 412             char type = a.TagData.at(index);
 413             ++index;
 414
 415             switch (type) {
 416                 case('A') :
 417                     m_out << "\"" << tagData[index] << "\"";
 418                     ++index;
 419                     break;
 420
 421                 case('C') :
 422                     m_out << (int)tagData[index];
 423                     ++index;
 424                     break;
 425
 426                 case('c') :
 427                     m_out << (int)tagData[index];
 428                     ++index;
 429                     break;
 430
 431                 case('S') :
 432                     m_out << BgzfData::UnpackUnsignedShort(&tagData[index]);
 433                     index += 2;
 434                     break;
 435
 436                 case('s') :
 437                     m_out << BgzfData::UnpackSignedShort(&tagData[index]);
 438                     index += 2;
 439                     break;
 440
 441                 case('I') :
 442                     m_out << BgzfData::UnpackUnsignedInt(&tagData[index]);
 443                     index += 4;
 444                     break;
 445
 446                 case('i') :
 447                     m_out << BgzfData::UnpackSignedInt(&tagData[index]);
 448                     index += 4;
 449                     break;
 450
 451                 case('f') :
 452                     m_out << BgzfData::UnpackFloat(&tagData[index]);
 453                     index += 4;
 454                     break;
 455
 456                 case('d') :
 457                     m_out << BgzfData::UnpackDouble(&tagData[index]);
 458                     index += 8;
 459                     break;
 460
 461                 case('Z') :
 462                 case('H') :
 463                     m_out << "\"";
 464                     while (tagData[index]) {
 465                         m_out << tagData[index];
 466                         ++index;
 467                     }
 468                     m_out << "\"";
 469                     ++index;
 470                     break;
 471             }
 472
 473             if ( tagData[index] == '\0')
 474                 break;
 475         }
 476
 477         m_out << "}";
 478     }
 479
 480     m_out << "}" << endl;
 481
 482 }
 483
 484 // print BamAlignment in SAM format
 485 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
 486
 487     // tab-delimited
 488     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 489
 490     // write name & alignment flag
 491     m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
 492
 493     // write reference name
 494     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 495         m_out << m_references[a.RefID].RefName << "\t";
 496     else
 497         m_out << "*\t";
 498
 499     // write position & map quality
 500     m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
 501
 502     // write CIGAR
 503     const vector<CigarOp>& cigarData = a.CigarData;
 504     if ( cigarData.empty() ) m_out << "*\t";
 505     else {
 506         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 507         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 508         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 509             const CigarOp& op = (*cigarIter);
 510             m_out << op.Length << op.Type;
 511         }
 512         m_out << "\t";
 513     }
 514
 515     // write mate reference name, mate position, & insert size
 516     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 517         if ( a.MateRefID == a.RefID ) m_out << "=\t";
 518         else m_out << m_references[a.MateRefID].RefName << "\t";
 519         m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 520     }
 521     else m_out << "*\t0\t0\t";
 522
 523     // write sequence
 524     if ( a.QueryBases.empty() ) m_out << "*\t";
 525     else m_out << a.QueryBases << "\t";
 526
 527     // write qualities
 528     if ( a.Qualities.empty() ) m_out << "*";
 529     else m_out << a.Qualities;
 530
 531     // write tag data
 532     const char* tagData = a.TagData.c_str();
 533     const size_t tagDataLength = a.TagData.length();
 534
 535     size_t index = 0;
 536     while ( index < tagDataLength ) {
 537
 538         // write tag name
 539         string tagName = a.TagData.substr(index, 2);
 540         m_out << "\t" << tagName << ":";
 541         index += 2;
 542
 543         // get data type
 544         char type = a.TagData.at(index);
 545         ++index;
 546         switch (type) {
 547             case('A') :
 548                 m_out << "A:" << tagData[index];
 549                 ++index;
 550                 break;
 551
 552             case('C') :
 553                 m_out << "i:" << (int)tagData[index];
 554                 ++index;
 555                 break;
 556
 557             case('c') :
 558                 m_out << "i:" << (int)tagData[index];
 559                 ++index;
 560                 break;
 561
 562             case('S') :
 563                 m_out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
 564                 index += 2;
 565                 break;
 566
 567             case('s') :
 568                 m_out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
 569                 index += 2;
 570                 break;
 571
 572             case('I') :
 573                 m_out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
 574                 index += 4;
 575                 break;
 576
 577             case('i') :
 578                 m_out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
 579                 index += 4;
 580                 break;
 581
 582             case('f') :
 583                 m_out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
 584                 index += 4;
 585                 break;
 586
 587             case('d') :
 588                 m_out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
 589                 index += 8;
 590                 break;
 591
 592             case('Z') :
 593             case('H') :
 594                 m_out << type << ":";
 595                 while (tagData[index]) {
 596                     m_out << tagData[index];
 597                     ++index;
 598                 }
 599                 ++index;
 600                 break;
 601         }
 602
 603         if ( tagData[index] == '\0')
 604             break;
 605     }
 606
 607     m_out << endl;
 608 }
 609
 610 void ConvertTool::ConvertToolPrivate::PrintWiggle(const BamAlignment& a) {
 611     ;
 612 }