2 #include "BamMultiReader.h"
3 #include "bamtools_pileup.h"
5 using namespace BamTools;
7 struct Pileup::PileupPrivate {
9 // ---------------------
13 BamMultiReader* Reader;
16 bool IsPrintingMapQualities;
22 vector<BamAlignment> CurrentData;
25 // ----------------------
28 PileupPrivate(BamMultiReader* reader, ostream* outStream)
30 , OutStream(outStream)
32 , IsPrintingMapQualities(false)
35 // ----------------------
38 void PrintCurrentData(void);
42 void Pileup::PileupPrivate::PrintCurrentData(void) {
44 // remove any data that ends before CurrentPosition
46 while ( i < CurrentData.size() ) {
47 if ( CurrentData[i].GetEndPosition() < CurrentPosition )
48 CurrentData.erase(CurrentData.begin() + i);
53 // if not data remains, return
54 if ( CurrentData.empty() ) return;
56 // initialize empty strings
58 string baseQuals = "";
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 ) {
67 const BamAlignment& al = (*dataIter);
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) );
74 bases.push_back( toupper(base) );
76 // determine current base quality & store
77 baseQuals.push_back( al.Qualities[CurrentPosition - al.Position] );
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);
87 // print results to OutStream
88 const string& refName = References[CurrentId].RefName;
89 const char refBase = 'N';
91 *OutStream << refName << "\t" << CurrentPosition << "\t" << refBase << "\t" << CurrentData.size() << "\t" << bases << "\t" << baseQuals;
92 if ( IsPrintingMapQualities ) *OutStream << "\t" << mapQuals;
96 bool Pileup::PileupPrivate::Run(void) {
98 // -----------------------------
99 // validate input & output
102 cerr << "Pileup::Run() : Invalid multireader" << endl;
107 cerr << "Pileup::Run() : Invalid output stream" << endl;
111 References = Reader->GetReferenceData();
113 // -----------------------------
114 // process input data
118 if ( !Reader->GetNextAlignment(al) ) {
119 cerr << "Pileup::Run() : Could not read from multireader" << endl;
123 // set initial markers & store first entry
124 CurrentId = al.RefID;
125 CurrentPosition = al.Position;
127 CurrentData.push_back(al);
129 // iterate over remaining data
130 while ( Reader->GetNextAlignment(al) ) {
133 if ( al.RefID == CurrentId ) {
135 // if same position, store and move on
136 if ( al.Position == CurrentPosition )
137 CurrentData.push_back(al);
139 // if less than CurrentPosition - sorting error => ABORT
140 else if ( al.Position < CurrentPosition ) {
141 cerr << "Pileup::Run() : Data not sorted correctly!" << endl;
145 // else print pileup data until 'catching up' to CurrentPosition
147 while ( al.Position > CurrentPosition ) {
151 CurrentData.push_back(al);
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;
161 // else moved forward onto next reference
164 // print any remaining pileup data from previous reference
165 while ( !CurrentData.empty() ) {
170 // store first entry on this new reference, update markers
172 CurrentData.push_back(al);
173 CurrentId = al.RefID;
174 CurrentPosition = al.Position;
178 // ------------------------------------
179 // handle any remaining data entries
181 while ( !CurrentData.empty() ) {
186 // -------------------------
192 // ----------------------------------------------------------
193 // Pileup implementation
195 Pileup::Pileup(BamMultiReader* reader, ostream* outStream) {
196 d = new PileupPrivate(reader, outStream);
199 Pileup::~Pileup(void) {
204 bool Pileup::Run(void) {
208 void Pileup::SetFastaFilename(const string& filename) {
209 d->FastaFilename = filename;
212 void Pileup::SetIsPrintingMapQualities(bool ok) {
213 d->IsPrintingMapQualities = ok;
216 void Pileup::SetRegion(const BamRegion& region) {