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