src/utils/bamtools_pileup.cpp

   1 // ***************************************************************************
   2 // bamtools_pileup.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 13 July 2010
   7 // ---------------------------------------------------------------------------
   8 // Provides pileup conversion functionality.
   9 //
  10 // The 'assembly' aspect of pileup makes this more complicated than the
  11 // simpler one-to-one conversion methods for other formats.
  12 // ***************************************************************************
  13
  14 #include <vector>
  15 #include "BamMultiReader.h"
  16 #include "bamtools_pileup.h"
  17 using namespace std;
  18 using namespace BamTools;
  19
  20 struct Pileup::PileupPrivate {
  21
  22     // ---------------------
  23     // data members
  24
  25     // IO & settings
  26     BamMultiReader* Reader;
  27     ostream* OutStream;
  28     string FastaFilename;
  29     bool IsPrintingMapQualities;
  30     BamRegion Region;
  31
  32     // parsing data
  33     int CurrentId;
  34     int CurrentPosition;
  35     vector<BamAlignment> CurrentData;
  36     RefVector References;
  37
  38     // ----------------------
  39     // ctor
  40
  41     PileupPrivate(BamMultiReader* reader, ostream* outStream)
  42         : Reader(reader)
  43         , OutStream(outStream)
  44         , FastaFilename("")
  45         , IsPrintingMapQualities(false)
  46     { }
  47
  48     // ----------------------
  49     // internal methods
  50
  51     void PrintCurrentData(void);
  52     bool Run(void);
  53 };
  54
  55 void Pileup::PileupPrivate::PrintCurrentData(void) {
  56
  57     // remove any data that ends before CurrentPosition
  58     size_t i = 0;
  59     while ( i < CurrentData.size() ) {
  60         if ( CurrentData[i].GetEndPosition() < CurrentPosition )
  61             CurrentData.erase(CurrentData.begin() + i);
  62         else
  63             ++i;
  64     }
  65
  66     // if not data remains, return
  67     if ( CurrentData.empty() ) return;
  68
  69     // initialize empty strings
  70     string bases     = "";
  71     string baseQuals = "";
  72     string mapQuals  = "";
  73
  74     // iterate over alignments
  75     vector<BamAlignment>::const_iterator dataIter = CurrentData.begin();
  76     vector<BamAlignment>::const_iterator dataEnd  = CurrentData.end();
  77     for ( ; dataIter != dataEnd; ++dataIter ) {
  78
  79         // retrieve alignment
  80         const BamAlignment& al = (*dataIter);
  81
  82         // determine current base character & store
  83         const char base = al.AlignedBases[CurrentPosition -al.Position];
  84         if ( al.IsReverseStrand() )
  85             bases.push_back( tolower(base) );
  86         else
  87             bases.push_back( toupper(base) );
  88
  89         // determine current base quality & store
  90         baseQuals.push_back( al.Qualities[CurrentPosition - al.Position] );
  91
  92         // if using mapQuals, determine current mapQual & store
  93         if ( IsPrintingMapQualities ) {
  94             int mapQuality = (int)(al.MapQuality + 33);
  95             if ( mapQuality > 126 ) mapQuality = 126;
  96             mapQuals.push_back((char)mapQuality);
  97         }
  98     }
  99
 100     // print results to OutStream
 101     const string& refName = References[CurrentId].RefName;
 102     const char refBase = 'N';
 103
 104     *OutStream << refName << "\t" << CurrentPosition << "\t" << refBase << "\t" << CurrentData.size() << "\t" << bases << "\t" << baseQuals;
 105     if ( IsPrintingMapQualities ) *OutStream << "\t" << mapQuals;
 106     *OutStream << endl;
 107 }
 108
 109 bool Pileup::PileupPrivate::Run(void) {
 110
 111     // -----------------------------
 112     // validate input & output
 113
 114     if ( !Reader ) {
 115         cerr << "Pileup::Run() : Invalid multireader" << endl;
 116         return false;
 117     }
 118
 119     if ( !OutStream) {
 120         cerr << "Pileup::Run() : Invalid output stream" << endl;
 121         return false;
 122     }
 123
 124     References = Reader->GetReferenceData();
 125
 126     // -----------------------------
 127     // process input data
 128
 129     // get first entry
 130     BamAlignment al;
 131     if ( !Reader->GetNextAlignment(al) ) {
 132         cerr << "Pileup::Run() : Could not read from multireader" << endl;
 133         return false;
 134     }
 135
 136     // set initial markers & store first entry
 137     CurrentId = al.RefID;
 138     CurrentPosition = al.Position;
 139     CurrentData.clear();
 140     CurrentData.push_back(al);
 141
 142     // iterate over remaining data
 143     while ( Reader->GetNextAlignment(al) ) {
 144
 145         // if same reference
 146         if ( al.RefID == CurrentId ) {
 147
 148             // if same position, store and move on
 149             if ( al.Position == CurrentPosition )
 150                 CurrentData.push_back(al);
 151
 152             // if less than CurrentPosition - sorting error => ABORT
 153             else if ( al.Position < CurrentPosition ) {
 154                 cerr << "Pileup::Run() : Data not sorted correctly!" << endl;
 155                 return false;
 156             }
 157
 158             // else print pileup data until 'catching up' to CurrentPosition
 159             else {
 160                 while ( al.Position > CurrentPosition ) {
 161                     PrintCurrentData();
 162                     ++CurrentPosition;
 163                 }
 164                 CurrentData.push_back(al);
 165             }
 166         }
 167
 168         // if reference ID less than CurrentID - sorting error => ABORT
 169         else if ( al.RefID < CurrentId ) {
 170             cerr << "Pileup::Run() : Data not sorted correctly!" << endl;
 171             return false;
 172         }
 173
 174         // else moved forward onto next reference
 175         else {
 176
 177             // print any remaining pileup data from previous reference
 178             while ( !CurrentData.empty() ) {
 179                 PrintCurrentData();
 180                 ++CurrentPosition;
 181             }
 182
 183             // store first entry on this new reference, update markers
 184             CurrentData.clear();
 185             CurrentData.push_back(al);
 186             CurrentId = al.RefID;
 187             CurrentPosition = al.Position;
 188         }
 189     }
 190
 191     // ------------------------------------
 192     // handle  any remaining data entries
 193
 194     while ( !CurrentData.empty() ) {
 195         PrintCurrentData();
 196         ++CurrentPosition;
 197     }
 198
 199     // -------------------------
 200     // return success
 201
 202     return true;
 203 }
 204
 205 // ----------------------------------------------------------
 206 // Pileup implementation
 207
 208 Pileup::Pileup(BamMultiReader* reader, ostream* outStream) {
 209     d = new PileupPrivate(reader, outStream);
 210 }
 211
 212 Pileup::~Pileup(void) {
 213     delete d;
 214     d = 0;
 215 }
 216
 217 bool Pileup::Run(void) {
 218     return d->Run();
 219 }
 220
 221 void Pileup::SetFastaFilename(const string& filename) {
 222     d->FastaFilename = filename;
 223 }
 224
 225 void Pileup::SetIsPrintingMapQualities(bool ok) {
 226     d->IsPrintingMapQualities = ok;
 227 }
 228
 229 void Pileup::SetRegion(const BamRegion& region) {
 230     d->Region = region;
 231 }