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