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.
10 // The 'assembly' aspect of pileup makes this more complicated than the
11 // simpler one-to-one conversion methods for other formats.
12 // ***************************************************************************
15 #include "BamMultiReader.h"
16 #include "bamtools_pileup.h"
18 using namespace BamTools;
20 struct Pileup::PileupPrivate {
22 // ---------------------
26 BamMultiReader* Reader;
29 bool IsPrintingMapQualities;
35 vector<BamAlignment> CurrentData;
38 // ----------------------
41 PileupPrivate(BamMultiReader* reader, ostream* outStream)
43 , OutStream(outStream)
45 , IsPrintingMapQualities(false)
48 // ----------------------
51 void PrintCurrentData(void);
55 void Pileup::PileupPrivate::PrintCurrentData(void) {
57 // remove any data that ends before CurrentPosition
59 while ( i < CurrentData.size() ) {
60 if ( CurrentData[i].GetEndPosition() < CurrentPosition )
61 CurrentData.erase(CurrentData.begin() + i);
66 // if not data remains, return
67 if ( CurrentData.empty() ) return;
69 // initialize empty strings
71 string baseQuals = "";
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 ) {
80 const BamAlignment& al = (*dataIter);
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) );
87 bases.push_back( toupper(base) );
89 // determine current base quality & store
90 baseQuals.push_back( al.Qualities[CurrentPosition - al.Position] );
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);
100 // print results to OutStream
101 const string& refName = References[CurrentId].RefName;
102 const char refBase = 'N';
104 *OutStream << refName << "\t" << CurrentPosition << "\t" << refBase << "\t" << CurrentData.size() << "\t" << bases << "\t" << baseQuals;
105 if ( IsPrintingMapQualities ) *OutStream << "\t" << mapQuals;
109 bool Pileup::PileupPrivate::Run(void) {
111 // -----------------------------
112 // validate input & output
115 cerr << "Pileup::Run() : Invalid multireader" << endl;
120 cerr << "Pileup::Run() : Invalid output stream" << endl;
124 References = Reader->GetReferenceData();
126 // -----------------------------
127 // process input data
131 if ( !Reader->GetNextAlignment(al) ) {
132 cerr << "Pileup::Run() : Could not read from multireader" << endl;
136 // set initial markers & store first entry
137 CurrentId = al.RefID;
138 CurrentPosition = al.Position;
140 CurrentData.push_back(al);
142 // iterate over remaining data
143 while ( Reader->GetNextAlignment(al) ) {
146 if ( al.RefID == CurrentId ) {
148 // if same position, store and move on
149 if ( al.Position == CurrentPosition )
150 CurrentData.push_back(al);
152 // if less than CurrentPosition - sorting error => ABORT
153 else if ( al.Position < CurrentPosition ) {
154 cerr << "Pileup::Run() : Data not sorted correctly!" << endl;
158 // else print pileup data until 'catching up' to CurrentPosition
160 while ( al.Position > CurrentPosition ) {
164 CurrentData.push_back(al);
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;
174 // else moved forward onto next reference
177 // print any remaining pileup data from previous reference
178 while ( !CurrentData.empty() ) {
183 // store first entry on this new reference, update markers
185 CurrentData.push_back(al);
186 CurrentId = al.RefID;
187 CurrentPosition = al.Position;
191 // ------------------------------------
192 // handle any remaining data entries
194 while ( !CurrentData.empty() ) {
199 // -------------------------
205 // ----------------------------------------------------------
206 // Pileup implementation
208 Pileup::Pileup(BamMultiReader* reader, ostream* outStream) {
209 d = new PileupPrivate(reader, outStream);
212 Pileup::~Pileup(void) {
217 bool Pileup::Run(void) {
221 void Pileup::SetFastaFilename(const string& filename) {
222 d->FastaFilename = filename;
225 void Pileup::SetIsPrintingMapQualities(bool ok) {
226 d->IsPrintingMapQualities = ok;
229 void Pileup::SetRegion(const BamRegion& region) {