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: 4 October 2010
7 // ---------------------------------------------------------------------------
8 // Converts between BAM and a number of other formats
9 // ***************************************************************************
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"
22 #include "BamMultiReader.h"
24 using namespace BamTools;
28 // ---------------------------------------------
29 // ConvertTool constants
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 static const string FORMAT_YAML = "yaml";
43 static const unsigned int FASTA_LINE_MAX = 50;
45 // ---------------------------------------------
46 // ConvertPileupFormatVisitor declaration
48 class ConvertPileupFormatVisitor : public PileupVisitor {
52 ConvertPileupFormatVisitor(const RefVector& references,
53 const string& fastaFilename,
54 const bool isPrintingMapQualities,
56 ~ConvertPileupFormatVisitor(void);
58 // PileupVisitor interface implementation
60 void Visit(const PileupPosition& pileupData);
66 bool m_isPrintingMapQualities;
68 RefVector m_references;
71 } // namespace BamTools
73 // ---------------------------------------------
74 // ConvertSettings implementation
76 struct ConvertTool::ConvertSettings {
85 bool HasFastaFilename;
86 bool IsOmittingSamHeader;
87 bool IsPrintingPileupMapQualities;
90 vector<string> InputFiles;
91 string OutputFilename;
104 , HasFastaFilename(false)
105 , IsOmittingSamHeader(false)
106 , IsPrintingPileupMapQualities(false)
107 , OutputFilename(Options::StandardOut())
112 // ---------------------------------------------
113 // ConvertToolPrivate implementation
115 struct ConvertTool::ConvertToolPrivate {
119 ConvertToolPrivate(ConvertTool::ConvertSettings* settings);
120 ~ConvertToolPrivate(void);
128 void PrintBed(const BamAlignment& a);
129 void PrintBedGraph(const BamAlignment& a);
130 void PrintFasta(const BamAlignment& a);
131 void PrintFastq(const BamAlignment& a);
132 void PrintJson(const BamAlignment& a);
133 void PrintSam(const BamAlignment& a);
134 void PrintWiggle(const BamAlignment& a);
135 void PrintYaml(const BamAlignment& a);
137 // special case - uses the PileupEngine
138 bool RunPileupConversion(BamMultiReader* reader);
142 ConvertTool::ConvertSettings* m_settings;
143 RefVector m_references;
147 ConvertTool::ConvertToolPrivate::ConvertToolPrivate(ConvertTool::ConvertSettings* settings)
148 : m_settings(settings)
149 , m_out(cout.rdbuf()) // default output to cout
152 ConvertTool::ConvertToolPrivate::~ConvertToolPrivate(void) { }
154 bool ConvertTool::ConvertToolPrivate::Run(void) {
156 // ------------------------------------
157 // initialize conversion input/output
159 // set to default input if none provided
160 if ( !m_settings->HasInput )
161 m_settings->InputFiles.push_back(Options::StandardIn());
164 BamMultiReader reader;
165 if ( !m_settings->HasInput ) { // don't attempt to open index for stdin
166 if ( !reader.Open(m_settings->InputFiles, false) ) {
167 cerr << "Could not open input files" << endl;
171 if ( !reader.Open(m_settings->InputFiles, true) ) {
172 if ( !reader.Open(m_settings->InputFiles, false) ) {
173 cerr << "Could not open input files" << endl;
176 cerr << "Opened reader without index file, jumping is disabled." << endl;
180 m_references = reader.GetReferenceData();
182 // set region if specified
184 if ( m_settings->HasRegion ) {
185 if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
186 if ( !reader.SetRegion(region) ) {
187 cerr << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
191 cerr << "Could not parse REGION: " << m_settings->Region << endl;
196 // if output file given
198 if ( m_settings->HasOutput ) {
200 // open output file stream
201 outFile.open(m_settings->OutputFilename.c_str());
203 cerr << "Could not open " << m_settings->OutputFilename << " for output" << endl;
207 // set m_out to file's streambuf
208 m_out.rdbuf(outFile.rdbuf());
211 // -------------------------------------
212 // do conversion based on format
214 bool convertedOk = true;
216 // pileup is special case
217 // conversion not done per alignment, like the other formats
218 if ( m_settings->Format == FORMAT_PILEUP )
219 convertedOk = RunPileupConversion(&reader);
224 bool formatError = false;
226 // set function pointer to proper conversion method
227 void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
228 if ( m_settings->Format == FORMAT_BED ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
229 else if ( m_settings->Format == FORMAT_BEDGRAPH ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBedGraph;
230 else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
231 else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
232 else if ( m_settings->Format == FORMAT_JSON ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
233 else if ( m_settings->Format == FORMAT_SAM ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
234 else if ( m_settings->Format == FORMAT_WIGGLE ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintWiggle;
235 else if ( m_settings->Format == FORMAT_YAML ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
237 cerr << "Unrecognized format: " << m_settings->Format << endl;
238 cerr << "Please see help|README (?) for details on supported formats " << endl;
243 // if format selected ok
244 if ( !formatError ) {
246 // if SAM format & not omitting header, print SAM header first
247 if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
248 m_out << reader.GetHeaderText();
250 // iterate through file, doing conversion
252 while ( reader.GetNextAlignment(a) )
253 (this->*pFunction)(a);
255 // set flag for successful conversion
260 // ------------------------
264 if ( m_settings->HasOutput ) outFile.close();
268 // ----------------------------------------------------------
269 // Conversion/output methods
270 // ----------------------------------------------------------
272 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a) {
274 // tab-delimited, 0-based half-open
275 // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
276 // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
278 m_out << m_references.at(a.RefID).RefName << "\t"
279 << a.Position << "\t"
280 << a.GetEndPosition() + 1 << "\t"
282 << a.MapQuality << "\t"
283 << (a.IsReverseStrand() ? "-" : "+") << endl;
286 void ConvertTool::ConvertToolPrivate::PrintBedGraph(const BamAlignment& a) {
290 // print BamAlignment in FASTA format
291 // N.B. - uses QueryBases NOT AlignedBases
292 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
294 // >BamAlignment.Name
295 // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
298 // N.B. - QueryBases are reverse-complemented if aligned to reverse strand
301 m_out << "> " << a.Name << endl;
303 // handle reverse strand alignment - bases
304 string sequence = a.QueryBases;
305 if ( a.IsReverseStrand() )
306 Utilities::ReverseComplement(sequence);
308 // if sequence fits on single line
309 if ( sequence.length() <= FASTA_LINE_MAX )
310 m_out << sequence << endl;
312 // else split over multiple lines
316 size_t seqLength = sequence.length(); // handle reverse strand alignment - bases & qualitiesth();
318 // write subsequences to each line
319 while ( position < (seqLength - FASTA_LINE_MAX) ) {
320 m_out << sequence.substr(position, FASTA_LINE_MAX) << endl;
321 position += FASTA_LINE_MAX;
324 // write final subsequence
325 m_out << sequence.substr(position) << endl;
329 // print BamAlignment in FASTQ format
330 // N.B. - uses QueryBases NOT AlignedBases
331 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
333 // @BamAlignment.Name
334 // BamAlignment.QueryBases
336 // BamAlignment.Qualities
338 // N.B. - QueryBases are reverse-complemented (& Qualities reversed) if aligned to reverse strand .
339 // Name is appended "/1" or "/2" if paired-end, to reflect which mate this entry is.
341 // handle paired-end alignments
342 string name = a.Name;
344 name.append( (a.IsFirstMate() ? "/1" : "/2") );
346 // handle reverse strand alignment - bases & qualities
347 string qualities = a.Qualities;
348 string sequence = a.QueryBases;
349 if ( a.IsReverseStrand() ) {
350 Utilities::Reverse(qualities);
351 Utilities::ReverseComplement(sequence);
354 // write to output stream
355 m_out << "@" << name << endl
358 << qualities << endl;
361 // print BamAlignment in JSON format
362 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
364 // write name & alignment flag
365 m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
367 // write reference name
368 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
369 m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
371 // write position & map quality
372 m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
375 const vector<CigarOp>& cigarData = a.CigarData;
376 if ( !cigarData.empty() ) {
377 m_out << "\"cigar\":[";
378 vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
379 vector<CigarOp>::const_iterator cigarIter = cigarBegin;
380 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
381 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
382 const CigarOp& op = (*cigarIter);
383 if (cigarIter != cigarBegin) m_out << ",";
384 m_out << "\"" << op.Length << op.Type << "\"";
389 // write mate reference name, mate position, & insert size
390 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
391 m_out << "\"mate\":{"
392 << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
393 << "\"position\":" << a.MatePosition+1
394 << ",\"insertSize\":" << a.InsertSize << "},";
398 if ( !a.QueryBases.empty() )
399 m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
402 if ( !a.Qualities.empty() ) {
403 string::const_iterator s = a.Qualities.begin();
404 m_out << "\"qualities\":[" << static_cast<short>(*s) - 33;
406 for (; s != a.Qualities.end(); ++s) {
407 m_out << "," << static_cast<short>(*s) - 33;
413 const char* tagData = a.TagData.c_str();
414 const size_t tagDataLength = a.TagData.length();
416 if (index < tagDataLength) {
418 m_out << "\"tags\":{";
420 while ( index < tagDataLength ) {
426 m_out << "\"" << a.TagData.substr(index, 2) << "\":";
430 char type = a.TagData.at(index);
435 m_out << "\"" << tagData[index] << "\"";
440 m_out << (int)tagData[index];
445 m_out << (int)tagData[index];
450 m_out << BgzfData::UnpackUnsignedShort(&tagData[index]);
455 m_out << BgzfData::UnpackSignedShort(&tagData[index]);
460 m_out << BgzfData::UnpackUnsignedInt(&tagData[index]);
465 m_out << BgzfData::UnpackSignedInt(&tagData[index]);
470 m_out << BgzfData::UnpackFloat(&tagData[index]);
475 m_out << BgzfData::UnpackDouble(&tagData[index]);
482 while (tagData[index]) {
483 if (tagData[index] == '\"')
484 m_out << "\\\""; // escape for json
486 m_out << tagData[index];
494 if ( tagData[index] == '\0')
501 m_out << "}" << endl;
504 // print BamAlignment in SAM format
505 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
508 // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
510 // write name & alignment flag
511 m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
513 // write reference name
514 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
515 m_out << m_references[a.RefID].RefName << "\t";
519 // write position & map quality
520 m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
523 const vector<CigarOp>& cigarData = a.CigarData;
524 if ( cigarData.empty() ) m_out << "*\t";
526 vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
527 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
528 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
529 const CigarOp& op = (*cigarIter);
530 m_out << op.Length << op.Type;
535 // write mate reference name, mate position, & insert size
536 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
537 if ( a.MateRefID == a.RefID ) m_out << "=\t";
538 else m_out << m_references[a.MateRefID].RefName << "\t";
539 m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
541 else m_out << "*\t0\t0\t";
544 if ( a.QueryBases.empty() ) m_out << "*\t";
545 else m_out << a.QueryBases << "\t";
548 if ( a.Qualities.empty() ) m_out << "*";
549 else m_out << a.Qualities;
552 const char* tagData = a.TagData.c_str();
553 const size_t tagDataLength = a.TagData.length();
556 while ( index < tagDataLength ) {
559 string tagName = a.TagData.substr(index, 2);
560 m_out << "\t" << tagName << ":";
564 char type = a.TagData.at(index);
568 m_out << "A:" << tagData[index];
573 m_out << "i:" << (int)tagData[index];
578 m_out << "i:" << (int)tagData[index];
583 m_out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
588 m_out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
593 m_out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
598 m_out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
603 m_out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
608 m_out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
614 m_out << type << ":";
615 while (tagData[index]) {
616 m_out << tagData[index];
623 if ( tagData[index] == '\0')
630 void ConvertTool::ConvertToolPrivate::PrintWiggle(const BamAlignment& a) {
634 // Print BamAlignment in YAML format
635 void ConvertTool::ConvertToolPrivate::PrintYaml(const BamAlignment& a) {
637 // write alignment name
638 m_out << "---" << endl;
639 m_out << a.Name << ":" << endl;
641 // write alignment data
642 m_out << " " << "AlndBases: " << a.AlignedBases << endl;
643 m_out << " " << "Qualities: " << a.Qualities << endl;
644 m_out << " " << "Name: " << a.Name << endl;
645 m_out << " " << "Length: " << a.Length << endl;
646 m_out << " " << "TagData: " << a.TagData << endl;
647 m_out << " " << "RefID: " << a.RefID << endl;
648 m_out << " " << "RefName: " << m_references[a.RefID].RefName << endl;
649 m_out << " " << "Position: " << a.Position << endl;
650 m_out << " " << "Bin: " << a.Bin << endl;
651 m_out << " " << "MapQuality: " << a.MapQuality << endl;
652 m_out << " " << "AlignmentFlag: " << a.AlignmentFlag << endl;
653 m_out << " " << "MateRefID: " << a.MateRefID << endl;
654 m_out << " " << "MatePosition: " << a.MatePosition << endl;
655 m_out << " " << "InsertSize: " << a.InsertSize << endl;
658 const vector<CigarOp>& cigarData = a.CigarData;
659 if ( !cigarData.empty() ) {
660 m_out << " " << "Cigar: ";
661 vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
662 vector<CigarOp>::const_iterator cigarIter = cigarBegin;
663 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
664 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
665 const CigarOp& op = (*cigarIter);
666 m_out << op.Length << op.Type;
672 bool ConvertTool::ConvertToolPrivate::RunPileupConversion(BamMultiReader* reader) {
674 // check for valid BamMultiReader
675 if ( reader == 0 ) return false;
677 // set up our pileup format 'visitor'
678 ConvertPileupFormatVisitor* v = new ConvertPileupFormatVisitor(m_references,
679 m_settings->FastaFilename,
680 m_settings->IsPrintingPileupMapQualities,
683 // set up PileupEngine
685 pileup.AddVisitor(v);
687 // iterate through data
689 while ( reader->GetNextAlignment(al) )
690 pileup.AddAlignment(al);
701 // ---------------------------------------------
702 // ConvertTool implementation
704 ConvertTool::ConvertTool(void)
706 , m_settings(new ConvertSettings)
709 // set program details
710 Options::SetProgramInfo("bamtools convert", "converts BAM to a number of other formats", "-format <FORMAT> [-in <filename> -in <filename> ...] [-out <filename>] [-region <REGION>] [format-specific options]");
713 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
714 Options::AddValueOption("-in", "BAM filename", "the input BAM file(s)", "", m_settings->HasInput, m_settings->InputFiles, IO_Opts, Options::StandardIn());
715 Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
716 Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
717 Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is used automatically if it exists. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
719 OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
720 Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts);
721 Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
723 OptionGroup* SamOpts = Options::CreateOptionGroup("SAM Options");
724 Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingSamHeader, SamOpts);
727 ConvertTool::~ConvertTool(void) {
735 int ConvertTool::Help(void) {
736 Options::DisplayHelp();
740 int ConvertTool::Run(int argc, char* argv[]) {
742 // parse command line arguments
743 Options::Parse(argc, argv, 1);
745 // run internal ConvertTool implementation, return success/fail
746 m_impl = new ConvertToolPrivate(m_settings);
754 // ---------------------------------------------
755 // ConvertPileupFormatVisitor implementation
757 ConvertPileupFormatVisitor::ConvertPileupFormatVisitor(const RefVector& references,
758 const string& fastaFilename,
759 const bool isPrintingMapQualities,
763 , m_isPrintingMapQualities(isPrintingMapQualities)
765 , m_references(references)
767 // set up Fasta reader if file is provided
768 if ( !fastaFilename.empty() ) {
770 // check for FASTA index
771 string indexFilename = "";
772 if ( Utilities::FileExists(fastaFilename + ".fai") )
773 indexFilename = fastaFilename + ".fai";
776 if ( m_fasta.Open(fastaFilename, indexFilename) )
781 ConvertPileupFormatVisitor::~ConvertPileupFormatVisitor(void) {
782 // be sure to close Fasta reader
789 void ConvertPileupFormatVisitor::Visit(const PileupPosition& pileupData ) {
791 // skip if no alignments at this position
792 if ( pileupData.PileupAlignments.empty() ) return;
794 // retrieve reference name
795 const string& referenceName = m_references[pileupData.RefId].RefName;
796 const int& position = pileupData.Position;
798 // retrieve reference base from FASTA file, if one provided; otherwise default to 'N'
799 char referenceBase('N');
800 if ( m_hasFasta && (pileupData.Position < m_references[pileupData.RefId].RefLength) ) {
801 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position, referenceBase ) ) {
802 cerr << "Pileup error : Could not read reference base from FASTA file" << endl;
807 // get count of alleles at this position
808 const int numberAlleles = pileupData.PileupAlignments.size();
810 // -----------------------------------------------------------
811 // build strings based on alleles at this positionInAlignment
813 stringstream bases("");
814 stringstream baseQualities("");
815 stringstream mapQualities("");
817 // iterate over alignments at this pileup position
818 vector<PileupAlignment>::const_iterator pileupIter = pileupData.PileupAlignments.begin();
819 vector<PileupAlignment>::const_iterator pileupEnd = pileupData.PileupAlignments.end();
820 for ( ; pileupIter != pileupEnd; ++pileupIter ) {
821 const PileupAlignment pa = (*pileupIter);
822 const BamAlignment& ba = pa.Alignment;
824 // if beginning of read segment
825 if ( pa.IsSegmentBegin )
826 bases << '^' << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
828 // if current base is not a DELETION
829 if ( !pa.IsCurrentDeletion ) {
831 // get base at current position
832 char base = ba.QueryBases.at(pa.PositionInAlignment);
834 // if base matches reference
836 toupper(base) == toupper(referenceBase) ||
837 tolower(base) == tolower(referenceBase) )
839 base = (ba.IsReverseStrand() ? ',' : '.' );
842 // mismatches reference
843 else base = (ba.IsReverseStrand() ? tolower(base) : toupper(base) );
848 // if next position contains insertion
849 if ( pa.IsNextInsertion ) {
850 bases << '+' << pa.InsertionLength;
851 for (int i = 1; i <= pa.InsertionLength; ++i) {
852 char insertedBase = (char)ba.QueryBases.at(pa.PositionInAlignment + i);
853 bases << (ba.IsReverseStrand() ? (char)tolower(insertedBase) : (char)toupper(insertedBase) );
857 // if next position contains DELETION
858 else if ( pa.IsNextDeletion ) {
859 bases << '-' << pa.DeletionLength;
860 for (int i = 1; i <= pa.DeletionLength; ++i) {
861 char deletedBase('N');
862 if ( m_hasFasta && (pileupData.Position+i < m_references[pileupData.RefId].RefLength) ) {
863 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position+i, deletedBase ) ) {
864 cerr << "Pileup error : Could not read reference base from FASTA file" << endl;
868 bases << (ba.IsReverseStrand() ? (char)tolower(deletedBase) : (char)toupper(deletedBase) );
873 // otherwise, DELETION
876 // if end of read segment
877 if ( pa.IsSegmentEnd ) bases << '$';
879 // store current base quality
880 baseQualities << ba.Qualities.at(pa.PositionInAlignment);
882 // save alignment map quality if desired
883 if ( m_isPrintingMapQualities )
884 mapQualities << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
887 // ----------------------
891 // <refName> <1-based pos> <refBase> <numberAlleles> <bases> <qualities> [mapQuals]
893 const string TAB = "\t";
894 *m_out << referenceName << TAB
895 << position + 1 << TAB
896 << referenceBase << TAB
897 << numberAlleles << TAB
898 << bases.str() << TAB
899 << baseQualities.str() << TAB
900 << mapQualities.str() << endl;