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: 11 June 2011
   7 // ---------------------------------------------------------------------------
   8 // Converts between BAM and a number of other formats
   9 // ***************************************************************************
  10
  11 #include "bamtools_convert.h"
  12
  13 #include <api/BamConstants.h>
  14 #include <api/BamMultiReader.h>
  15 #include <utils/bamtools_fasta.h>
  16 #include <utils/bamtools_options.h>
  17 #include <utils/bamtools_pileup_engine.h>
  18 #include <utils/bamtools_utilities.h>
  19 using namespace BamTools;
  20
  21 #include <fstream>
  22 #include <iostream>
  23 #include <sstream>
  24 #include <string>
  25 #include <vector>
  26 using namespace std;
  27
  28 namespace BamTools {
  29
  30 // ---------------------------------------------
  31 // ConvertTool constants
  32
  33 // supported conversion format command-line names
  34 static const string FORMAT_BED    = "bed";
  35 static const string FORMAT_FASTA  = "fasta";
  36 static const string FORMAT_FASTQ  = "fastq";
  37 static const string FORMAT_JSON   = "json";
  38 static const string FORMAT_SAM    = "sam";
  39 static const string FORMAT_PILEUP = "pileup";
  40 static const string FORMAT_YAML   = "yaml";
  41
  42 // other constants
  43 static const unsigned int FASTA_LINE_MAX = 50;
  44
  45 // ---------------------------------------------
  46 // ConvertPileupFormatVisitor declaration
  47
  48 class ConvertPileupFormatVisitor : public PileupVisitor {
  49
  50     // ctor & dtor
  51     public:
  52         ConvertPileupFormatVisitor(const RefVector& references,
  53                                    const string& fastaFilename,
  54                                    const bool isPrintingMapQualities,
  55                                    ostream* out);
  56         ~ConvertPileupFormatVisitor(void);
  57
  58     // PileupVisitor interface implementation
  59     public:
  60         void Visit(const PileupPosition& pileupData);
  61
  62     // data members
  63     private:
  64         Fasta     m_fasta;
  65         bool      m_hasFasta;
  66         bool      m_isPrintingMapQualities;
  67         ostream*  m_out;
  68         RefVector m_references;
  69 };
  70
  71 } // namespace BamTools
  72
  73 // ---------------------------------------------
  74 // ConvertSettings implementation
  75
  76 struct ConvertTool::ConvertSettings {
  77
  78     // flag
  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         , FastaFilename("")
 109     { }
 110 };
 111
 112 // ---------------------------------------------
 113 // ConvertToolPrivate implementation
 114
 115 struct ConvertTool::ConvertToolPrivate {
 116
 117     // ctor & dtor
 118     public:
 119         ConvertToolPrivate(ConvertTool::ConvertSettings* settings)
 120             : m_settings(settings)
 121             , m_out(cout.rdbuf())
 122         { }
 123
 124         ~ConvertToolPrivate(void) { }
 125
 126     // interface
 127     public:
 128         bool Run(void);
 129
 130     // internal methods
 131     private:
 132         void PrintBed(const BamAlignment& a);
 133         void PrintFasta(const BamAlignment& a);
 134         void PrintFastq(const BamAlignment& a);
 135         void PrintJson(const BamAlignment& a);
 136         void PrintSam(const BamAlignment& a);
 137         void PrintYaml(const BamAlignment& a);
 138
 139         // special case - uses the PileupEngine
 140         bool RunPileupConversion(BamMultiReader* reader);
 141
 142     // data members
 143     private:
 144         ConvertTool::ConvertSettings* m_settings;
 145         RefVector m_references;
 146         ostream m_out;
 147 };
 148
 149 bool ConvertTool::ConvertToolPrivate::Run(void) {
 150
 151     // ------------------------------------
 152     // initialize conversion input/output
 153
 154     // set to default input if none provided
 155     if ( !m_settings->HasInput )
 156         m_settings->InputFiles.push_back(Options::StandardIn());
 157
 158     // open input files
 159     BamMultiReader reader;
 160     if ( !reader.Open(m_settings->InputFiles) ) {
 161         cerr << "bamtools convert ERROR: could not open input BAM file(s)... Aborting." << endl;
 162         return false;
 163     }
 164
 165     // if input is not stdin & a region is provided, look for index files
 166     if ( m_settings->HasInput && m_settings->HasRegion ) {
 167         if ( !reader.LocateIndexes() ) {
 168             cerr << "bamtools convert ERROR: could not locate index file(s)... Aborting." << endl;
 169             return false;
 170         }
 171     }
 172
 173     // retrieve reference data
 174     m_references = reader.GetReferenceData();
 175
 176     // set region if specified
 177     BamRegion region;
 178     if ( m_settings->HasRegion ) {
 179         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 180
 181             if ( reader.HasIndexes() ) {
 182                 if ( !reader.SetRegion(region) ) {
 183                     cerr << "bamtools convert ERROR: set region failed. Check that REGION describes a valid range" << endl;
 184                     reader.Close();
 185                     return false;
 186                 }
 187             }
 188
 189         } else {
 190             cerr << "bamtools convert ERROR: could not parse REGION: " << m_settings->Region << endl;
 191             cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 192                  << endl;
 193             reader.Close();
 194             return false;
 195         }
 196     }
 197
 198     // if output file given
 199     ofstream outFile;
 200     if ( m_settings->HasOutput ) {
 201
 202         // open output file stream
 203         outFile.open(m_settings->OutputFilename.c_str());
 204         if ( !outFile ) {
 205             cerr << "bamtools convert ERROR: could not open " << m_settings->OutputFilename
 206                  << " for output" << endl;
 207             return false;
 208         }
 209
 210         // set m_out to file's streambuf
 211         m_out.rdbuf(outFile.rdbuf());
 212     }
 213
 214     // -------------------------------------
 215     // do conversion based on format
 216
 217      bool convertedOk = true;
 218
 219     // pileup is special case
 220     // conversion not done per alignment, like the other formats
 221     if ( m_settings->Format == FORMAT_PILEUP )
 222         convertedOk = RunPileupConversion(&reader);
 223
 224     // all other formats
 225     else {
 226
 227         bool formatError = false;
 228
 229         // set function pointer to proper conversion method
 230         void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
 231         if      ( m_settings->Format == FORMAT_BED )   pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
 232         else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
 233         else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
 234         else if ( m_settings->Format == FORMAT_JSON )  pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
 235         else if ( m_settings->Format == FORMAT_SAM )   pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
 236         else if ( m_settings->Format == FORMAT_YAML )  pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
 237         else {
 238             cerr << "bamtools convert ERROR: unrecognized format: " << m_settings->Format << endl;
 239             cerr << "Please see documentation for list of supported formats " << endl;
 240             formatError = true;
 241             convertedOk = false;
 242         }
 243
 244         // if format selected ok
 245         if ( !formatError ) {
 246
 247             // if SAM format & not omitting header, print SAM header first
 248             if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
 249                 m_out << reader.GetHeaderText();
 250
 251             // iterate through file, doing conversion
 252             BamAlignment a;
 253             while ( reader.GetNextAlignment(a) )
 254                 (this->*pFunction)(a);
 255
 256             // set flag for successful conversion
 257             convertedOk = true;
 258         }
 259     }
 260
 261     // ------------------------
 262     // clean up & exit
 263     reader.Close();
 264     if ( m_settings->HasOutput )
 265         outFile.close();
 266     return convertedOk;
 267 }
 268
 269 // ----------------------------------------------------------
 270 // Conversion/output methods
 271 // ----------------------------------------------------------
 272
 273 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a) {
 274
 275     // tab-delimited, 0-based half-open
 276     // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
 277     // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
 278
 279     m_out << m_references.at(a.RefID).RefName << "\t"
 280           << a.Position << "\t"
 281           << a.GetEndPosition() + 1 << "\t"
 282           << a.Name << "\t"
 283           << a.MapQuality << "\t"
 284           << (a.IsReverseStrand() ? "-" : "+") << endl;
 285 }
 286
 287 // print BamAlignment in FASTA format
 288 // N.B. - uses QueryBases NOT AlignedBases
 289 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
 290
 291     // >BamAlignment.Name
 292     // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
 293     // ...
 294     //
 295     // N.B. - QueryBases are reverse-complemented if aligned to reverse strand
 296
 297     // print header
 298     m_out << ">" << a.Name << endl;
 299
 300     // handle reverse strand alignment - bases
 301     string sequence = a.QueryBases;
 302     if ( a.IsReverseStrand() )
 303         Utilities::ReverseComplement(sequence);
 304
 305     // if sequence fits on single line
 306     if ( sequence.length() <= FASTA_LINE_MAX )
 307         m_out << sequence << endl;
 308
 309     // else split over multiple lines
 310     else {
 311
 312         size_t position = 0;
 313         size_t seqLength = sequence.length(); // handle reverse strand alignment - bases & qualitiesth();
 314
 315         // write subsequences to each line
 316         while ( position < (seqLength - FASTA_LINE_MAX) ) {
 317             m_out << sequence.substr(position, FASTA_LINE_MAX) << endl;
 318             position += FASTA_LINE_MAX;
 319         }
 320
 321         // write final subsequence
 322         m_out << sequence.substr(position) << endl;
 323     }
 324 }
 325
 326 // print BamAlignment in FASTQ format
 327 // N.B. - uses QueryBases NOT AlignedBases
 328 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
 329
 330     // @BamAlignment.Name
 331     // BamAlignment.QueryBases
 332     // +
 333     // BamAlignment.Qualities
 334     //
 335     // N.B. - QueryBases are reverse-complemented (& Qualities reversed) if aligned to reverse strand .
 336     //        Name is appended "/1" or "/2" if paired-end, to reflect which mate this entry is.
 337
 338     // handle paired-end alignments
 339     string name = a.Name;
 340     if ( a.IsPaired() )
 341         name.append( (a.IsFirstMate() ? "/1" : "/2") );
 342
 343     // handle reverse strand alignment - bases & qualities
 344     string qualities = a.Qualities;
 345     string sequence  = a.QueryBases;
 346     if ( a.IsReverseStrand() ) {
 347         Utilities::Reverse(qualities);
 348         Utilities::ReverseComplement(sequence);
 349     }
 350
 351     // write to output stream
 352     m_out << "@" << name << endl
 353           << sequence    << endl
 354           << "+"         << endl
 355           << qualities   << endl;
 356 }
 357
 358 // print BamAlignment in JSON format
 359 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
 360
 361     // write name & alignment flag
 362     m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
 363
 364     // write reference name
 365     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 366         m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
 367
 368     // write position & map quality
 369     m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
 370
 371     // write CIGAR
 372     const vector<CigarOp>& cigarData = a.CigarData;
 373     if ( !cigarData.empty() ) {
 374         m_out << "\"cigar\":[";
 375         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 376         vector<CigarOp>::const_iterator cigarIter  = cigarBegin;
 377         vector<CigarOp>::const_iterator cigarEnd   = cigarData.end();
 378         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 379             const CigarOp& op = (*cigarIter);
 380             if (cigarIter != cigarBegin)
 381                 m_out << ",";
 382             m_out << "\"" << op.Length << op.Type << "\"";
 383         }
 384         m_out << "],";
 385     }
 386
 387     // write mate reference name, mate position, & insert size
 388     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 389         m_out << "\"mate\":{"
 390               << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
 391               << "\"position\":" << a.MatePosition+1
 392               << ",\"insertSize\":" << a.InsertSize << "},";
 393     }
 394
 395     // write sequence
 396     if ( !a.QueryBases.empty() )
 397         m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
 398
 399     // write qualities
 400     if ( !a.Qualities.empty() ) {
 401         string::const_iterator s = a.Qualities.begin();
 402         m_out << "\"qualities\":[" << static_cast<short>(*s) - 33;
 403         ++s;
 404         for ( ; s != a.Qualities.end(); ++s )
 405             m_out << "," << static_cast<short>(*s) - 33;
 406         m_out << "],";
 407     }
 408
 409     // write alignment's source BAM file
 410     m_out << "\"filename\":" << a.Filename << ",";
 411
 412     // write tag data
 413     const char* tagData = a.TagData.c_str();
 414     const size_t tagDataLength = a.TagData.length();
 415     size_t index = 0;
 416     if ( index < tagDataLength ) {
 417
 418         m_out << "\"tags\":{";
 419
 420         while ( index < tagDataLength ) {
 421
 422             if ( index > 0 )
 423                 m_out << ",";
 424
 425             // write tag name
 426             m_out << "\"" << a.TagData.substr(index, 2) << "\":";
 427             index += 2;
 428
 429             // get data type
 430             char type = a.TagData.at(index);
 431             ++index;
 432             switch ( type ) {
 433                 case (Constants::BAM_TAG_TYPE_ASCII) :
 434                     m_out << "\"" << tagData[index] << "\"";
 435                     ++index;
 436                     break;
 437
 438                 case (Constants::BAM_TAG_TYPE_INT8)  :
 439                 case (Constants::BAM_TAG_TYPE_UINT8) :
 440                     m_out << (int)tagData[index];
 441                     ++index;
 442                     break;
 443
 444                 case (Constants::BAM_TAG_TYPE_INT16) :
 445                     m_out << BamTools::UnpackSignedShort(&tagData[index]);
 446                     index += sizeof(int16_t);
 447                     break;
 448
 449                 case (Constants::BAM_TAG_TYPE_UINT16) :
 450                     m_out << BamTools::UnpackUnsignedShort(&tagData[index]);
 451                     index += sizeof(uint16_t);
 452                     break;
 453
 454                 case (Constants::BAM_TAG_TYPE_INT32) :
 455                     m_out << BamTools::UnpackSignedInt(&tagData[index]);
 456                     index += sizeof(int32_t);
 457                     break;
 458
 459                 case (Constants::BAM_TAG_TYPE_UINT32) :
 460                     m_out << BamTools::UnpackUnsignedInt(&tagData[index]);
 461                     index += sizeof(uint32_t);
 462                     break;
 463
 464                 case (Constants::BAM_TAG_TYPE_FLOAT) :
 465                     m_out << BamTools::UnpackFloat(&tagData[index]);
 466                     index += sizeof(float);
 467                     break;
 468
 469                 case (Constants::BAM_TAG_TYPE_HEX)    :
 470                 case (Constants::BAM_TAG_TYPE_STRING) :
 471                     m_out << "\"";
 472                     while (tagData[index]) {
 473                         if (tagData[index] == '\"')
 474                             m_out << "\\\""; // escape for json
 475                         else
 476                             m_out << tagData[index];
 477                         ++index;
 478                     }
 479                     m_out << "\"";
 480                     ++index;
 481                     break;
 482             }
 483
 484             if ( tagData[index] == '\0')
 485                 break;
 486         }
 487
 488         m_out << "}";
 489     }
 490
 491     m_out << "}" << endl;
 492 }
 493
 494 // print BamAlignment in SAM format
 495 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
 496
 497     // tab-delimited
 498     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 499
 500     // write name & alignment flag
 501     m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
 502
 503     // write reference name
 504     if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
 505         m_out << m_references[a.RefID].RefName << "\t";
 506     else
 507         m_out << "*\t";
 508
 509     // write position & map quality
 510     m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
 511
 512     // write CIGAR
 513     const vector<CigarOp>& cigarData = a.CigarData;
 514     if ( cigarData.empty() ) m_out << "*\t";
 515     else {
 516         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 517         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 518         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 519             const CigarOp& op = (*cigarIter);
 520             m_out << op.Length << op.Type;
 521         }
 522         m_out << "\t";
 523     }
 524
 525     // write mate reference name, mate position, & insert size
 526     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
 527         if ( a.MateRefID == a.RefID )
 528             m_out << "=\t";
 529         else
 530             m_out << m_references[a.MateRefID].RefName << "\t";
 531         m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 532     }
 533     else
 534         m_out << "*\t0\t0\t";
 535
 536     // write sequence
 537     if ( a.QueryBases.empty() )
 538         m_out << "*\t";
 539     else
 540         m_out << a.QueryBases << "\t";
 541
 542     // write qualities
 543     if ( a.Qualities.empty() )
 544         m_out << "*";
 545     else
 546         m_out << a.Qualities;
 547
 548     // write tag data
 549     const char* tagData = a.TagData.c_str();
 550     const size_t tagDataLength = a.TagData.length();
 551
 552     size_t index = 0;
 553     while ( index < tagDataLength ) {
 554
 555         // write tag name
 556         string tagName = a.TagData.substr(index, 2);
 557         m_out << "\t" << tagName << ":";
 558         index += 2;
 559
 560         // get data type
 561         char type = a.TagData.at(index);
 562         ++index;
 563         switch ( type ) {
 564             case (Constants::BAM_TAG_TYPE_ASCII) :
 565                 m_out << "A:" << tagData[index];
 566                 ++index;
 567                 break;
 568
 569             case (Constants::BAM_TAG_TYPE_INT8)  :
 570             case (Constants::BAM_TAG_TYPE_UINT8) :
 571                 m_out << "i:" << (int)tagData[index];
 572                 ++index;
 573                 break;
 574
 575             case (Constants::BAM_TAG_TYPE_INT16) :
 576                 m_out << "i:" << BamTools::UnpackSignedShort(&tagData[index]);
 577                 index += sizeof(int16_t);
 578                 break;
 579
 580             case (Constants::BAM_TAG_TYPE_UINT16) :
 581                 m_out << "i:" << BamTools::UnpackUnsignedShort(&tagData[index]);
 582                 index += sizeof(uint16_t);
 583                 break;
 584
 585             case (Constants::BAM_TAG_TYPE_INT32) :
 586                 m_out << "i:" << BamTools::UnpackSignedInt(&tagData[index]);
 587                 index += sizeof(int32_t);
 588                 break;
 589
 590             case (Constants::BAM_TAG_TYPE_UINT32) :
 591                 m_out << "i:" << BamTools::UnpackUnsignedInt(&tagData[index]);
 592                 index += sizeof(uint32_t);
 593                 break;
 594
 595             case (Constants::BAM_TAG_TYPE_FLOAT) :
 596                 m_out << "f:" << BamTools::UnpackFloat(&tagData[index]);
 597                 index += sizeof(float);
 598                 break;
 599
 600             case (Constants::BAM_TAG_TYPE_HEX)    :
 601             case (Constants::BAM_TAG_TYPE_STRING) :
 602                 m_out << type << ":";
 603                 while (tagData[index]) {
 604                     m_out << tagData[index];
 605                     ++index;
 606                 }
 607                 ++index;
 608                 break;
 609         }
 610
 611         if ( tagData[index] == '\0')
 612             break;
 613     }
 614
 615     m_out << endl;
 616 }
 617
 618 // Print BamAlignment in YAML format
 619 void ConvertTool::ConvertToolPrivate::PrintYaml(const BamAlignment& a) {
 620
 621     // write alignment name
 622     m_out << "---" << endl;
 623     m_out << a.Name << ":" << endl;
 624
 625     // write alignment data
 626     m_out << "   " << "AlndBases: "     << a.AlignedBases << endl;
 627     m_out << "   " << "Qualities: "     << a.Qualities << endl;
 628     m_out << "   " << "Name: "          << a.Name << endl;
 629     m_out << "   " << "Length: "        << a.Length << endl;
 630     m_out << "   " << "TagData: "       << a.TagData << endl;
 631     m_out << "   " << "RefID: "         << a.RefID << endl;
 632     m_out << "   " << "RefName: "       << m_references[a.RefID].RefName << endl;
 633     m_out << "   " << "Position: "      << a.Position << endl;
 634     m_out << "   " << "Bin: "           << a.Bin << endl;
 635     m_out << "   " << "MapQuality: "    << a.MapQuality << endl;
 636     m_out << "   " << "AlignmentFlag: " << a.AlignmentFlag << endl;
 637     m_out << "   " << "MateRefID: "     << a.MateRefID << endl;
 638     m_out << "   " << "MatePosition: "  << a.MatePosition << endl;
 639     m_out << "   " << "InsertSize: "    << a.InsertSize << endl;
 640     m_out << "   " << "Filename: "      << a.Filename << endl;
 641
 642     // write Cigar data
 643     const vector<CigarOp>& cigarData = a.CigarData;
 644     if ( !cigarData.empty() ) {
 645         m_out << "   " <<  "Cigar: ";
 646         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 647         vector<CigarOp>::const_iterator cigarIter  = cigarBegin;
 648         vector<CigarOp>::const_iterator cigarEnd   = cigarData.end();
 649         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 650             const CigarOp& op = (*cigarIter);
 651             m_out << op.Length << op.Type;
 652         }
 653         m_out << endl;
 654     }
 655 }
 656
 657 bool ConvertTool::ConvertToolPrivate::RunPileupConversion(BamMultiReader* reader) {
 658
 659     // check for valid BamMultiReader
 660     if ( reader == 0 ) return false;
 661
 662     // set up our pileup format 'visitor'
 663     ConvertPileupFormatVisitor* v = new ConvertPileupFormatVisitor(m_references,
 664                                                                    m_settings->FastaFilename,
 665                                                                    m_settings->IsPrintingPileupMapQualities,
 666                                                                    &m_out);
 667
 668     // set up PileupEngine
 669     PileupEngine pileup;
 670     pileup.AddVisitor(v);
 671
 672     // iterate through data
 673     BamAlignment al;
 674     while ( reader->GetNextAlignment(al) )
 675         pileup.AddAlignment(al);
 676     pileup.Flush();
 677
 678     // clean up
 679     delete v;
 680     v = 0;
 681
 682     // return success
 683     return true;
 684 }
 685
 686 // ---------------------------------------------
 687 // ConvertTool implementation
 688
 689 ConvertTool::ConvertTool(void)
 690     : AbstractTool()
 691     , m_settings(new ConvertSettings)
 692     , m_impl(0)
 693 {
 694     // set program details
 695     Options::SetProgramInfo("bamtools convert", "converts BAM to a number of other formats", "-format <FORMAT> [-in <filename> -in <filename> ...] [-out <filename>] [-region <REGION>] [format-specific options]");
 696
 697     // set up options
 698     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 699     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInput,   m_settings->InputFiles,     IO_Opts, Options::StandardIn());
 700     Options::AddValueOption("-out",    "BAM filename", "the output BAM file",   "", m_settings->HasOutput,  m_settings->OutputFilename, IO_Opts, Options::StandardOut());
 701     Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
 702     Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is used automatically if it exists. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
 703
 704     OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
 705     Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts);
 706     Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
 707
 708     OptionGroup* SamOpts = Options::CreateOptionGroup("SAM Options");
 709     Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingSamHeader, SamOpts);
 710 }
 711
 712 ConvertTool::~ConvertTool(void) {
 713
 714     delete m_settings;
 715     m_settings = 0;
 716
 717     delete m_impl;
 718     m_impl = 0;
 719 }
 720
 721 int ConvertTool::Help(void) {
 722     Options::DisplayHelp();
 723     return 0;
 724 }
 725
 726 int ConvertTool::Run(int argc, char* argv[]) {
 727
 728     // parse command line arguments
 729     Options::Parse(argc, argv, 1);
 730
 731     // initialize ConvertTool with settings
 732     m_impl = new ConvertToolPrivate(m_settings);
 733
 734     // run ConvertTool, return success/fail
 735     if ( m_impl->Run() )
 736         return 0;
 737     else
 738         return 1;
 739 }
 740
 741 // ---------------------------------------------
 742 // ConvertPileupFormatVisitor implementation
 743
 744 ConvertPileupFormatVisitor::ConvertPileupFormatVisitor(const RefVector& references,
 745                                                        const string& fastaFilename,
 746                                                        const bool isPrintingMapQualities,
 747                                                        ostream* out)
 748     : PileupVisitor()
 749     , m_hasFasta(false)
 750     , m_isPrintingMapQualities(isPrintingMapQualities)
 751     , m_out(out)
 752     , m_references(references)
 753 {
 754     // set up Fasta reader if file is provided
 755     if ( !fastaFilename.empty() ) {
 756
 757         // check for FASTA index
 758         string indexFilename = "";
 759         if ( Utilities::FileExists(fastaFilename + ".fai") )
 760             indexFilename = fastaFilename + ".fai";
 761
 762         // open FASTA file
 763         if ( m_fasta.Open(fastaFilename, indexFilename) )
 764             m_hasFasta = true;
 765     }
 766 }
 767
 768 ConvertPileupFormatVisitor::~ConvertPileupFormatVisitor(void) {
 769     // be sure to close Fasta reader
 770     if ( m_hasFasta ) {
 771         m_fasta.Close();
 772         m_hasFasta = false;
 773     }
 774 }
 775
 776 void ConvertPileupFormatVisitor::Visit(const PileupPosition& pileupData ) {
 777
 778     // skip if no alignments at this position
 779     if ( pileupData.PileupAlignments.empty() ) return;
 780
 781     // retrieve reference name
 782     const string& referenceName = m_references[pileupData.RefId].RefName;
 783     const int& position   = pileupData.Position;
 784
 785     // retrieve reference base from FASTA file, if one provided; otherwise default to 'N'
 786     char referenceBase('N');
 787     if ( m_hasFasta && (pileupData.Position < m_references[pileupData.RefId].RefLength) ) {
 788         if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position, referenceBase ) ) {
 789             cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
 790             return;
 791         }
 792     }
 793
 794     // get count of alleles at this position
 795     const int numberAlleles = pileupData.PileupAlignments.size();
 796
 797     // -----------------------------------------------------------
 798     // build strings based on alleles at this positionInAlignment
 799
 800     stringstream bases("");
 801     stringstream baseQualities("");
 802     stringstream mapQualities("");
 803
 804     // iterate over alignments at this pileup position
 805     vector<PileupAlignment>::const_iterator pileupIter = pileupData.PileupAlignments.begin();
 806     vector<PileupAlignment>::const_iterator pileupEnd  = pileupData.PileupAlignments.end();
 807     for ( ; pileupIter != pileupEnd; ++pileupIter ) {
 808         const PileupAlignment pa = (*pileupIter);
 809         const BamAlignment& ba = pa.Alignment;
 810
 811         // if beginning of read segment
 812         if ( pa.IsSegmentBegin )
 813             bases << '^' << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
 814
 815         // if current base is not a DELETION
 816         if ( !pa.IsCurrentDeletion ) {
 817
 818             // get base at current position
 819             char base = ba.QueryBases.at(pa.PositionInAlignment);
 820
 821             // if base matches reference
 822             if ( base == '=' ||
 823                  toupper(base) == toupper(referenceBase) ||
 824                  tolower(base) == tolower(referenceBase) )
 825             {
 826                 base = ( ba.IsReverseStrand() ? ',' : '.' );
 827             }
 828
 829             // mismatches reference
 830             else base = ( ba.IsReverseStrand() ? tolower(base) : toupper(base) );
 831
 832             // store base
 833             bases << base;
 834
 835             // if next position contains insertion
 836             if ( pa.IsNextInsertion ) {
 837                 bases << '+' << pa.InsertionLength;
 838                 for (int i = 1; i <= pa.InsertionLength; ++i) {
 839                     char insertedBase = (char)ba.QueryBases.at(pa.PositionInAlignment + i);
 840                     bases << (ba.IsReverseStrand() ? (char)tolower(insertedBase) : (char)toupper(insertedBase) );
 841                 }
 842             }
 843
 844             // if next position contains DELETION
 845             else if ( pa.IsNextDeletion ) {
 846                 bases << '-' << pa.DeletionLength;
 847                 for (int i = 1; i <= pa.DeletionLength; ++i) {
 848                     char deletedBase('N');
 849                     if ( m_hasFasta && (pileupData.Position+i < m_references[pileupData.RefId].RefLength) ) {
 850                         if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position+i, deletedBase ) ) {
 851                             cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
 852                             return;
 853                         }
 854                     }
 855                     bases << (ba.IsReverseStrand() ? (char)tolower(deletedBase) : (char)toupper(deletedBase) );
 856                 }
 857             }
 858         }
 859
 860         // otherwise, DELETION
 861         else bases << '*';
 862
 863         // if end of read segment
 864         if ( pa.IsSegmentEnd )
 865             bases << '$';
 866
 867         // store current base quality
 868         baseQualities << ba.Qualities.at(pa.PositionInAlignment);
 869
 870         // save alignment map quality if desired
 871         if ( m_isPrintingMapQualities )
 872             mapQualities << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
 873     }
 874
 875     // ----------------------
 876     // print results
 877
 878     // tab-delimited
 879     // <refName> <1-based pos> <refBase> <numberAlleles> <bases> <qualities> [mapQuals]
 880
 881     const string TAB = "\t";
 882     *m_out << referenceName       << TAB
 883            << position + 1        << TAB
 884            << referenceBase       << TAB
 885            << numberAlleles       << TAB
 886            << bases.str()         << TAB
 887            << baseQualities.str() << TAB
 888            << mapQualities.str()  << endl;
 889 }