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