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: 7 April 2011
7 // ---------------------------------------------------------------------------
8 // Converts between BAM and a number of other formats
9 // ***************************************************************************
11 #include "bamtools_convert.h"
13 #include <api/BamConstants.h>
14 #include <api/BamMultiReader.h>
15 #include <utils/bamtools_fasta.h>
16 #include <utils/bamtools_options.h>
17 #include <utils/bamtools_pileup_engine.h>
18 #include <utils/bamtools_utilities.h>
19 using namespace BamTools;
30 // ---------------------------------------------
31 // ConvertTool constants
33 // supported conversion format command-line names
34 static const string FORMAT_BED = "bed";
35 static const string FORMAT_FASTA = "fasta";
36 static const string FORMAT_FASTQ = "fastq";
37 static const string FORMAT_JSON = "json";
38 static const string FORMAT_SAM = "sam";
39 static const string FORMAT_PILEUP = "pileup";
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 : m_settings(settings)
121 , m_out(cout.rdbuf())
124 ~ConvertToolPrivate(void) { }
132 void PrintBed(const BamAlignment& a);
133 void PrintFasta(const BamAlignment& a);
134 void PrintFastq(const BamAlignment& a);
135 void PrintJson(const BamAlignment& a);
136 void PrintSam(const BamAlignment& a);
137 void PrintYaml(const BamAlignment& a);
139 // special case - uses the PileupEngine
140 bool RunPileupConversion(BamMultiReader* reader);
144 ConvertTool::ConvertSettings* m_settings;
145 RefVector m_references;
149 bool ConvertTool::ConvertToolPrivate::Run(void) {
151 // ------------------------------------
152 // initialize conversion input/output
154 // set to default input if none provided
155 if ( !m_settings->HasInput )
156 m_settings->InputFiles.push_back(Options::StandardIn());
159 BamMultiReader reader;
160 if ( !reader.Open(m_settings->InputFiles) ) {
161 cerr << "bamtools convert ERROR: could not open input BAM file(s)... Aborting." << endl;
165 // if input is not stdin & a region is provided, look for index files
166 if ( m_settings->HasInput && m_settings->HasRegion ) {
167 if ( !reader.LocateIndexes() ) {
168 cerr << "bamtools convert ERROR: could not locate index file(s)... Aborting." << endl;
173 // retrieve reference data
174 m_references = reader.GetReferenceData();
176 // set region if specified
178 if ( m_settings->HasRegion ) {
179 if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
181 if ( reader.HasIndexes() ) {
182 if ( !reader.SetRegion(region) ) {
183 cerr << "bamtools convert ERROR: set region failed. Check that REGION describes a valid range" << endl;
190 cerr << "bamtools convert ERROR: could not parse REGION: " << m_settings->Region << endl;
191 cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
198 // if output file given
200 if ( m_settings->HasOutput ) {
202 // open output file stream
203 outFile.open(m_settings->OutputFilename.c_str());
205 cerr << "bamtools convert ERROR: could not open " << m_settings->OutputFilename
206 << " for output" << endl;
210 // set m_out to file's streambuf
211 m_out.rdbuf(outFile.rdbuf());
214 // -------------------------------------
215 // do conversion based on format
217 bool convertedOk = true;
219 // pileup is special case
220 // conversion not done per alignment, like the other formats
221 if ( m_settings->Format == FORMAT_PILEUP )
222 convertedOk = RunPileupConversion(&reader);
227 bool formatError = false;
229 // set function pointer to proper conversion method
230 void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
231 if ( m_settings->Format == FORMAT_BED ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
232 else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
233 else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
234 else if ( m_settings->Format == FORMAT_JSON ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
235 else if ( m_settings->Format == FORMAT_SAM ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
236 else if ( m_settings->Format == FORMAT_YAML ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
238 cerr << "bamtools convert ERROR: unrecognized format: " << m_settings->Format << endl;
239 cerr << "Please see documentation for list of supported formats " << endl;
244 // if format selected ok
245 if ( !formatError ) {
247 // if SAM format & not omitting header, print SAM header first
248 if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
249 m_out << reader.GetHeaderText();
251 // iterate through file, doing conversion
253 while ( reader.GetNextAlignment(a) )
254 (this->*pFunction)(a);
256 // set flag for successful conversion
261 // ------------------------
264 if ( m_settings->HasOutput )
269 // ----------------------------------------------------------
270 // Conversion/output methods
271 // ----------------------------------------------------------
273 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a) {
275 // tab-delimited, 0-based half-open
276 // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
277 // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
279 m_out << m_references.at(a.RefID).RefName << "\t"
280 << a.Position << "\t"
281 << a.GetEndPosition() + 1 << "\t"
283 << a.MapQuality << "\t"
284 << (a.IsReverseStrand() ? "-" : "+") << endl;
287 // print BamAlignment in FASTA format
288 // N.B. - uses QueryBases NOT AlignedBases
289 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
291 // >BamAlignment.Name
292 // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
295 // N.B. - QueryBases are reverse-complemented if aligned to reverse strand
298 m_out << "> " << a.Name << endl;
300 // handle reverse strand alignment - bases
301 string sequence = a.QueryBases;
302 if ( a.IsReverseStrand() )
303 Utilities::ReverseComplement(sequence);
305 // if sequence fits on single line
306 if ( sequence.length() <= FASTA_LINE_MAX )
307 m_out << sequence << endl;
309 // else split over multiple lines
313 size_t seqLength = sequence.length(); // handle reverse strand alignment - bases & qualitiesth();
315 // write subsequences to each line
316 while ( position < (seqLength - FASTA_LINE_MAX) ) {
317 m_out << sequence.substr(position, FASTA_LINE_MAX) << endl;
318 position += FASTA_LINE_MAX;
321 // write final subsequence
322 m_out << sequence.substr(position) << endl;
326 // print BamAlignment in FASTQ format
327 // N.B. - uses QueryBases NOT AlignedBases
328 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
330 // @BamAlignment.Name
331 // BamAlignment.QueryBases
333 // BamAlignment.Qualities
335 // N.B. - QueryBases are reverse-complemented (& Qualities reversed) if aligned to reverse strand .
336 // Name is appended "/1" or "/2" if paired-end, to reflect which mate this entry is.
338 // handle paired-end alignments
339 string name = a.Name;
341 name.append( (a.IsFirstMate() ? "/1" : "/2") );
343 // handle reverse strand alignment - bases & qualities
344 string qualities = a.Qualities;
345 string sequence = a.QueryBases;
346 if ( a.IsReverseStrand() ) {
347 Utilities::Reverse(qualities);
348 Utilities::ReverseComplement(sequence);
351 // write to output stream
352 m_out << "@" << name << endl
355 << qualities << endl;
358 // print BamAlignment in JSON format
359 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
361 // write name & alignment flag
362 m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
364 // write reference name
365 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
366 m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
368 // write position & map quality
369 m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
372 const vector<CigarOp>& cigarData = a.CigarData;
373 if ( !cigarData.empty() ) {
374 m_out << "\"cigar\":[";
375 vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
376 vector<CigarOp>::const_iterator cigarIter = cigarBegin;
377 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
378 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
379 const CigarOp& op = (*cigarIter);
380 if (cigarIter != cigarBegin)
382 m_out << "\"" << op.Length << op.Type << "\"";
387 // write mate reference name, mate position, & insert size
388 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
389 m_out << "\"mate\":{"
390 << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
391 << "\"position\":" << a.MatePosition+1
392 << ",\"insertSize\":" << a.InsertSize << "},";
396 if ( !a.QueryBases.empty() )
397 m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
400 if ( !a.Qualities.empty() ) {
401 string::const_iterator s = a.Qualities.begin();
402 m_out << "\"qualities\":[" << static_cast<short>(*s) - 33;
404 for ( ; s != a.Qualities.end(); ++s )
405 m_out << "," << static_cast<short>(*s) - 33;
409 // write alignment's source BAM file
410 m_out << "\"filename\":" << a.Filename << ",";
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);
433 case (Constants::BAM_TAG_TYPE_ASCII) :
434 m_out << "\"" << tagData[index] << "\"";
438 case (Constants::BAM_TAG_TYPE_INT8) :
439 case (Constants::BAM_TAG_TYPE_UINT8) :
440 m_out << (int)tagData[index];
444 case (Constants::BAM_TAG_TYPE_INT16) :
445 m_out << BamTools::UnpackSignedShort(&tagData[index]);
446 index += sizeof(int16_t);
449 case (Constants::BAM_TAG_TYPE_UINT16) :
450 m_out << BamTools::UnpackUnsignedShort(&tagData[index]);
451 index += sizeof(uint16_t);
454 case (Constants::BAM_TAG_TYPE_INT32) :
455 m_out << BamTools::UnpackSignedInt(&tagData[index]);
456 index += sizeof(int32_t);
459 case (Constants::BAM_TAG_TYPE_UINT32) :
460 m_out << BamTools::UnpackUnsignedInt(&tagData[index]);
461 index += sizeof(uint32_t);
464 case (Constants::BAM_TAG_TYPE_FLOAT) :
465 m_out << BamTools::UnpackFloat(&tagData[index]);
466 index += sizeof(float);
469 case (Constants::BAM_TAG_TYPE_HEX) :
470 case (Constants::BAM_TAG_TYPE_STRING) :
472 while (tagData[index]) {
473 if (tagData[index] == '\"')
474 m_out << "\\\""; // escape for json
476 m_out << tagData[index];
484 if ( tagData[index] == '\0')
491 m_out << "}" << endl;
494 // print BamAlignment in SAM format
495 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
498 // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
500 // write name & alignment flag
501 m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
503 // write reference name
504 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
505 m_out << m_references[a.RefID].RefName << "\t";
509 // write position & map quality
510 m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
513 const vector<CigarOp>& cigarData = a.CigarData;
514 if ( cigarData.empty() ) m_out << "*\t";
516 vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
517 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
518 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
519 const CigarOp& op = (*cigarIter);
520 m_out << op.Length << op.Type;
525 // write mate reference name, mate position, & insert size
526 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
527 if ( a.MateRefID == a.RefID )
530 m_out << m_references[a.MateRefID].RefName << "\t";
531 m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
534 m_out << "*\t0\t0\t";
537 if ( a.QueryBases.empty() )
540 m_out << a.QueryBases << "\t";
543 if ( a.Qualities.empty() )
546 m_out << a.Qualities;
549 const char* tagData = a.TagData.c_str();
550 const size_t tagDataLength = a.TagData.length();
553 while ( index < tagDataLength ) {
556 string tagName = a.TagData.substr(index, 2);
557 m_out << "\t" << tagName << ":";
561 char type = a.TagData.at(index);
564 case (Constants::BAM_TAG_TYPE_ASCII) :
565 m_out << "A:" << tagData[index];
569 case (Constants::BAM_TAG_TYPE_INT8) :
570 case (Constants::BAM_TAG_TYPE_UINT8) :
571 m_out << "i:" << (int)tagData[index];
575 case (Constants::BAM_TAG_TYPE_INT16) :
576 m_out << "i:" << BamTools::UnpackSignedShort(&tagData[index]);
577 index += sizeof(int16_t);
580 case (Constants::BAM_TAG_TYPE_UINT16) :
581 m_out << "i:" << BamTools::UnpackUnsignedShort(&tagData[index]);
582 index += sizeof(uint16_t);
585 case (Constants::BAM_TAG_TYPE_INT32) :
586 m_out << "i:" << BamTools::UnpackSignedInt(&tagData[index]);
587 index += sizeof(int32_t);
590 case (Constants::BAM_TAG_TYPE_UINT32) :
591 m_out << "i:" << BamTools::UnpackUnsignedInt(&tagData[index]);
592 index += sizeof(uint32_t);
595 case (Constants::BAM_TAG_TYPE_FLOAT) :
596 m_out << "f:" << BamTools::UnpackFloat(&tagData[index]);
597 index += sizeof(float);
600 case (Constants::BAM_TAG_TYPE_HEX) :
601 case (Constants::BAM_TAG_TYPE_STRING) :
602 m_out << type << ":";
603 while (tagData[index]) {
604 m_out << tagData[index];
611 if ( tagData[index] == '\0')
618 // Print BamAlignment in YAML format
619 void ConvertTool::ConvertToolPrivate::PrintYaml(const BamAlignment& a) {
621 // write alignment name
622 m_out << "---" << endl;
623 m_out << a.Name << ":" << endl;
625 // write alignment data
626 m_out << " " << "AlndBases: " << a.AlignedBases << endl;
627 m_out << " " << "Qualities: " << a.Qualities << endl;
628 m_out << " " << "Name: " << a.Name << endl;
629 m_out << " " << "Length: " << a.Length << endl;
630 m_out << " " << "TagData: " << a.TagData << endl;
631 m_out << " " << "RefID: " << a.RefID << endl;
632 m_out << " " << "RefName: " << m_references[a.RefID].RefName << endl;
633 m_out << " " << "Position: " << a.Position << endl;
634 m_out << " " << "Bin: " << a.Bin << endl;
635 m_out << " " << "MapQuality: " << a.MapQuality << endl;
636 m_out << " " << "AlignmentFlag: " << a.AlignmentFlag << endl;
637 m_out << " " << "MateRefID: " << a.MateRefID << endl;
638 m_out << " " << "MatePosition: " << a.MatePosition << endl;
639 m_out << " " << "InsertSize: " << a.InsertSize << endl;
640 m_out << " " << "Filename: " << a.Filename << endl;
643 const vector<CigarOp>& cigarData = a.CigarData;
644 if ( !cigarData.empty() ) {
645 m_out << " " << "Cigar: ";
646 vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
647 vector<CigarOp>::const_iterator cigarIter = cigarBegin;
648 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
649 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
650 const CigarOp& op = (*cigarIter);
651 m_out << op.Length << op.Type;
657 bool ConvertTool::ConvertToolPrivate::RunPileupConversion(BamMultiReader* reader) {
659 // check for valid BamMultiReader
660 if ( reader == 0 ) return false;
662 // set up our pileup format 'visitor'
663 ConvertPileupFormatVisitor* v = new ConvertPileupFormatVisitor(m_references,
664 m_settings->FastaFilename,
665 m_settings->IsPrintingPileupMapQualities,
668 // set up PileupEngine
670 pileup.AddVisitor(v);
672 // iterate through data
674 while ( reader->GetNextAlignment(al) )
675 pileup.AddAlignment(al);
686 // ---------------------------------------------
687 // ConvertTool implementation
689 ConvertTool::ConvertTool(void)
691 , m_settings(new ConvertSettings)
694 // set program details
695 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]");
698 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
699 Options::AddValueOption("-in", "BAM filename", "the input BAM file(s)", "", m_settings->HasInput, m_settings->InputFiles, IO_Opts, Options::StandardIn());
700 Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
701 Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
702 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);
704 OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
705 Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts);
706 Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
708 OptionGroup* SamOpts = Options::CreateOptionGroup("SAM Options");
709 Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingSamHeader, SamOpts);
712 ConvertTool::~ConvertTool(void) {
721 int ConvertTool::Help(void) {
722 Options::DisplayHelp();
726 int ConvertTool::Run(int argc, char* argv[]) {
728 // parse command line arguments
729 Options::Parse(argc, argv, 1);
731 // initialize ConvertTool with settings
732 m_impl = new ConvertToolPrivate(m_settings);
734 // run ConvertTool, return success/fail
741 // ---------------------------------------------
742 // ConvertPileupFormatVisitor implementation
744 ConvertPileupFormatVisitor::ConvertPileupFormatVisitor(const RefVector& references,
745 const string& fastaFilename,
746 const bool isPrintingMapQualities,
750 , m_isPrintingMapQualities(isPrintingMapQualities)
752 , m_references(references)
754 // set up Fasta reader if file is provided
755 if ( !fastaFilename.empty() ) {
757 // check for FASTA index
758 string indexFilename = "";
759 if ( Utilities::FileExists(fastaFilename + ".fai") )
760 indexFilename = fastaFilename + ".fai";
763 if ( m_fasta.Open(fastaFilename, indexFilename) )
768 ConvertPileupFormatVisitor::~ConvertPileupFormatVisitor(void) {
769 // be sure to close Fasta reader
776 void ConvertPileupFormatVisitor::Visit(const PileupPosition& pileupData ) {
778 // skip if no alignments at this position
779 if ( pileupData.PileupAlignments.empty() ) return;
781 // retrieve reference name
782 const string& referenceName = m_references[pileupData.RefId].RefName;
783 const int& position = pileupData.Position;
785 // retrieve reference base from FASTA file, if one provided; otherwise default to 'N'
786 char referenceBase('N');
787 if ( m_hasFasta && (pileupData.Position < m_references[pileupData.RefId].RefLength) ) {
788 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position, referenceBase ) ) {
789 cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
794 // get count of alleles at this position
795 const int numberAlleles = pileupData.PileupAlignments.size();
797 // -----------------------------------------------------------
798 // build strings based on alleles at this positionInAlignment
800 stringstream bases("");
801 stringstream baseQualities("");
802 stringstream mapQualities("");
804 // iterate over alignments at this pileup position
805 vector<PileupAlignment>::const_iterator pileupIter = pileupData.PileupAlignments.begin();
806 vector<PileupAlignment>::const_iterator pileupEnd = pileupData.PileupAlignments.end();
807 for ( ; pileupIter != pileupEnd; ++pileupIter ) {
808 const PileupAlignment pa = (*pileupIter);
809 const BamAlignment& ba = pa.Alignment;
811 // if beginning of read segment
812 if ( pa.IsSegmentBegin )
813 bases << '^' << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
815 // if current base is not a DELETION
816 if ( !pa.IsCurrentDeletion ) {
818 // get base at current position
819 char base = ba.QueryBases.at(pa.PositionInAlignment);
821 // if base matches reference
823 toupper(base) == toupper(referenceBase) ||
824 tolower(base) == tolower(referenceBase) )
826 base = ( ba.IsReverseStrand() ? ',' : '.' );
829 // mismatches reference
830 else base = ( ba.IsReverseStrand() ? tolower(base) : toupper(base) );
835 // if next position contains insertion
836 if ( pa.IsNextInsertion ) {
837 bases << '+' << pa.InsertionLength;
838 for (int i = 1; i <= pa.InsertionLength; ++i) {
839 char insertedBase = (char)ba.QueryBases.at(pa.PositionInAlignment + i);
840 bases << (ba.IsReverseStrand() ? (char)tolower(insertedBase) : (char)toupper(insertedBase) );
844 // if next position contains DELETION
845 else if ( pa.IsNextDeletion ) {
846 bases << '-' << pa.DeletionLength;
847 for (int i = 1; i <= pa.DeletionLength; ++i) {
848 char deletedBase('N');
849 if ( m_hasFasta && (pileupData.Position+i < m_references[pileupData.RefId].RefLength) ) {
850 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position+i, deletedBase ) ) {
851 cerr << "bamtools convert ERROR: pileup conversion - could not read reference base from FASTA file" << endl;
855 bases << (ba.IsReverseStrand() ? (char)tolower(deletedBase) : (char)toupper(deletedBase) );
860 // otherwise, DELETION
863 // if end of read segment
864 if ( pa.IsSegmentEnd )
867 // store current base quality
868 baseQualities << ba.Qualities.at(pa.PositionInAlignment);
870 // save alignment map quality if desired
871 if ( m_isPrintingMapQualities )
872 mapQualities << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
875 // ----------------------
879 // <refName> <1-based pos> <refBase> <numberAlleles> <bases> <qualities> [mapQuals]
881 const string TAB = "\t";
882 *m_out << referenceName << TAB
883 << position + 1 << TAB
884 << referenceBase << TAB
885 << numberAlleles << TAB
886 << bases.str() << TAB
887 << baseQualities.str() << TAB
888 << mapQualities.str() << endl;