bamtools_pileup.cpp

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