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: 16 September 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";
42 static const unsigned int FASTA_LINE_MAX = 50;
44 // ---------------------------------------------
45 // ConvertPileupFormatVisitor declaration
47 class ConvertPileupFormatVisitor : public PileupVisitor {
51 ConvertPileupFormatVisitor(const RefVector& references,
52 const string& fastaFilename,
53 const bool isPrintingMapQualities,
55 ~ConvertPileupFormatVisitor(void);
57 // PileupVisitor interface implementation
59 // void Visit(const int& refId, const int& position, const vector<BamAlignment>& alignments);
61 void Visit(const PileupPosition& pileupData);
67 bool m_isPrintingMapQualities;
69 RefVector m_references;
72 } // namespace BamTools
74 // ---------------------------------------------
75 // ConvertSettings implementation
77 struct ConvertTool::ConvertSettings {
86 bool HasFastaFilename;
87 bool IsOmittingSamHeader;
88 bool IsPrintingPileupMapQualities;
91 vector<string> InputFiles;
92 string OutputFilename;
100 ConvertSettings(void)
105 , HasFastaFilename(false)
106 , IsOmittingSamHeader(false)
107 , IsPrintingPileupMapQualities(false)
108 , OutputFilename(Options::StandardOut())
113 // ---------------------------------------------
114 // ConvertToolPrivate implementation
116 struct ConvertTool::ConvertToolPrivate {
120 ConvertToolPrivate(ConvertTool::ConvertSettings* settings);
121 ~ConvertToolPrivate(void);
129 void PrintBed(const BamAlignment& a);
130 void PrintBedGraph(const BamAlignment& a);
131 void PrintFasta(const BamAlignment& a);
132 void PrintFastq(const BamAlignment& a);
133 void PrintJson(const BamAlignment& a);
134 void PrintSam(const BamAlignment& a);
135 void PrintWiggle(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 reader.Open(m_settings->InputFiles);
166 m_references = reader.GetReferenceData();
168 // set region if specified
170 if ( m_settings->HasRegion ) {
171 if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
172 if ( !reader.SetRegion(region) ) {
173 cerr << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
177 cerr << "Could not parse REGION: " << m_settings->Region << endl;
182 // if output file given
184 if ( m_settings->HasOutput ) {
186 // open output file stream
187 outFile.open(m_settings->OutputFilename.c_str());
189 cerr << "Could not open " << m_settings->OutputFilename << " for output" << endl;
193 // set m_out to file's streambuf
194 m_out.rdbuf(outFile.rdbuf());
197 // -------------------------------------
198 // do conversion based on format
200 bool convertedOk = true;
202 // pileup is special case
203 // conversion not done per alignment, like the other formats
204 if ( m_settings->Format == FORMAT_PILEUP )
205 convertedOk = RunPileupConversion(&reader);
210 bool formatError = false;
212 // set function pointer to proper conversion method
213 void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
214 if ( m_settings->Format == FORMAT_BED ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
215 else if ( m_settings->Format == FORMAT_BEDGRAPH ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBedGraph;
216 else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
217 else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
218 else if ( m_settings->Format == FORMAT_JSON ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
219 else if ( m_settings->Format == FORMAT_SAM ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
220 else if ( m_settings->Format == FORMAT_WIGGLE ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintWiggle;
222 cerr << "Unrecognized format: " << m_settings->Format << endl;
223 cerr << "Please see help|README (?) for details on supported formats " << endl;
228 // if format selected ok
229 if ( !formatError ) {
231 // if SAM format & not omitting header, print SAM header first
232 if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
233 m_out << reader.GetHeaderText();
235 // iterate through file, doing conversion
237 while ( reader.GetNextAlignment(a) )
238 (this->*pFunction)(a);
240 // set flag for successful conversion
245 // ------------------------
249 if ( m_settings->HasOutput ) outFile.close();
253 // ----------------------------------------------------------
254 // Conversion/output methods
255 // ----------------------------------------------------------
257 void ConvertTool::ConvertToolPrivate::PrintBed(const BamAlignment& a) {
259 // tab-delimited, 0-based half-open
260 // (e.g. a 50-base read aligned to pos 10 could have BED coordinates (10, 60) instead of BAM coordinates (10, 59) )
261 // <chromName> <chromStart> <chromEnd> <readName> <score> <strand>
263 m_out << m_references.at(a.RefID).RefName << "\t"
264 << a.Position << "\t"
265 << a.GetEndPosition() + 1 << "\t"
267 << a.MapQuality << "\t"
268 << (a.IsReverseStrand() ? "-" : "+") << endl;
271 void ConvertTool::ConvertToolPrivate::PrintBedGraph(const BamAlignment& a) {
275 // print BamAlignment in FASTA format
276 // N.B. - uses QueryBases NOT AlignedBases
277 void ConvertTool::ConvertToolPrivate::PrintFasta(const BamAlignment& a) {
279 // >BamAlignment.Name
280 // BamAlignment.QueryBases (up to FASTA_LINE_MAX bases per line)
284 m_out << "> " << a.Name << endl;
286 // if sequence fits on single line
287 if ( a.QueryBases.length() <= FASTA_LINE_MAX )
288 m_out << a.QueryBases << endl;
290 // else split over multiple lines
294 size_t seqLength = a.QueryBases.length();
296 // write subsequences to each line
297 while ( position < (seqLength - FASTA_LINE_MAX) ) {
298 m_out << a.QueryBases.substr(position, FASTA_LINE_MAX) << endl;
299 position += FASTA_LINE_MAX;
302 // write final subsequence
303 m_out << a.QueryBases.substr(position) << endl;
307 // print BamAlignment in FASTQ format
308 // N.B. - uses QueryBases NOT AlignedBases
309 void ConvertTool::ConvertToolPrivate::PrintFastq(const BamAlignment& a) {
311 // @BamAlignment.Name
312 // BamAlignment.QueryBases
314 // BamAlignment.Qualities
316 m_out << "@" << a.Name << endl
317 << a.QueryBases << endl
319 << a.Qualities << endl;
322 // print BamAlignment in JSON format
323 void ConvertTool::ConvertToolPrivate::PrintJson(const BamAlignment& a) {
325 // write name & alignment flag
326 m_out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\"" << a.AlignmentFlag << "\",";
328 // write reference name
329 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
330 m_out << "\"reference\":\"" << m_references[a.RefID].RefName << "\",";
332 // write position & map quality
333 m_out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
336 const vector<CigarOp>& cigarData = a.CigarData;
337 if ( !cigarData.empty() ) {
338 m_out << "\"cigar\":[";
339 vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
340 vector<CigarOp>::const_iterator cigarIter = cigarBegin;
341 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
342 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
343 const CigarOp& op = (*cigarIter);
344 if (cigarIter != cigarBegin) m_out << ",";
345 m_out << "\"" << op.Length << op.Type << "\"";
350 // write mate reference name, mate position, & insert size
351 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
352 m_out << "\"mate\":{"
353 << "\"reference\":\"" << m_references[a.MateRefID].RefName << "\","
354 << "\"position\":" << a.MatePosition+1
355 << ",\"insertSize\":" << a.InsertSize << "},";
359 if ( !a.QueryBases.empty() )
360 m_out << "\"queryBases\":\"" << a.QueryBases << "\",";
363 if ( !a.Qualities.empty() ) {
364 string::const_iterator s = a.Qualities.begin();
365 m_out << "\"qualities\":[" << static_cast<short>(*s) - 33;
367 for (; s != a.Qualities.end(); ++s) {
368 m_out << "," << static_cast<short>(*s) - 33;
374 const char* tagData = a.TagData.c_str();
375 const size_t tagDataLength = a.TagData.length();
377 if (index < tagDataLength) {
379 m_out << "\"tags\":{";
381 while ( index < tagDataLength ) {
387 m_out << "\"" << a.TagData.substr(index, 2) << "\":";
391 char type = a.TagData.at(index);
396 m_out << "\"" << tagData[index] << "\"";
401 m_out << (int)tagData[index];
406 m_out << (int)tagData[index];
411 m_out << BgzfData::UnpackUnsignedShort(&tagData[index]);
416 m_out << BgzfData::UnpackSignedShort(&tagData[index]);
421 m_out << BgzfData::UnpackUnsignedInt(&tagData[index]);
426 m_out << BgzfData::UnpackSignedInt(&tagData[index]);
431 m_out << BgzfData::UnpackFloat(&tagData[index]);
436 m_out << BgzfData::UnpackDouble(&tagData[index]);
443 while (tagData[index]) {
444 m_out << tagData[index];
452 if ( tagData[index] == '\0')
459 m_out << "}" << endl;
462 // print BamAlignment in SAM format
463 void ConvertTool::ConvertToolPrivate::PrintSam(const BamAlignment& a) {
466 // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
468 // write name & alignment flag
469 m_out << a.Name << "\t" << a.AlignmentFlag << "\t";
471 // write reference name
472 if ( (a.RefID >= 0) && (a.RefID < (int)m_references.size()) )
473 m_out << m_references[a.RefID].RefName << "\t";
477 // write position & map quality
478 m_out << a.Position+1 << "\t" << a.MapQuality << "\t";
481 const vector<CigarOp>& cigarData = a.CigarData;
482 if ( cigarData.empty() ) m_out << "*\t";
484 vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
485 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
486 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
487 const CigarOp& op = (*cigarIter);
488 m_out << op.Length << op.Type;
493 // write mate reference name, mate position, & insert size
494 if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)m_references.size()) ) {
495 if ( a.MateRefID == a.RefID ) m_out << "=\t";
496 else m_out << m_references[a.MateRefID].RefName << "\t";
497 m_out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
499 else m_out << "*\t0\t0\t";
502 if ( a.QueryBases.empty() ) m_out << "*\t";
503 else m_out << a.QueryBases << "\t";
506 if ( a.Qualities.empty() ) m_out << "*";
507 else m_out << a.Qualities;
510 const char* tagData = a.TagData.c_str();
511 const size_t tagDataLength = a.TagData.length();
514 while ( index < tagDataLength ) {
517 string tagName = a.TagData.substr(index, 2);
518 m_out << "\t" << tagName << ":";
522 char type = a.TagData.at(index);
526 m_out << "A:" << tagData[index];
531 m_out << "i:" << (int)tagData[index];
536 m_out << "i:" << (int)tagData[index];
541 m_out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
546 m_out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
551 m_out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
556 m_out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
561 m_out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
566 m_out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
572 m_out << type << ":";
573 while (tagData[index]) {
574 m_out << tagData[index];
581 if ( tagData[index] == '\0')
588 void ConvertTool::ConvertToolPrivate::PrintWiggle(const BamAlignment& a) {
592 bool ConvertTool::ConvertToolPrivate::RunPileupConversion(BamMultiReader* reader) {
594 // check for valid BamMultiReader
595 if ( reader == 0 ) return false;
597 // set up our pileup format 'visitor'
598 ConvertPileupFormatVisitor* v = new ConvertPileupFormatVisitor(m_references,
599 m_settings->FastaFilename,
600 m_settings->IsPrintingPileupMapQualities,
603 // set up PileupEngine
605 pileup.AddVisitor(v);
607 // iterate through data
609 while ( reader->GetNextAlignment(al) )
610 pileup.AddAlignment(al);
621 // ---------------------------------------------
622 // ConvertTool implementation
624 ConvertTool::ConvertTool(void)
626 , m_settings(new ConvertSettings)
629 // set program details
630 Options::SetProgramInfo("bamtools convert", "converts BAM to a number of other formats", "-format <FORMAT> [-in <filename> -in <filename> ...] [-out <filename>] [other options]");
633 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
634 Options::AddValueOption("-in", "BAM filename", "the input BAM file(s)", "", m_settings->HasInput, m_settings->InputFiles, IO_Opts, Options::StandardIn());
635 Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
636 Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
638 OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
639 Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is read automatically if it exists as <filename>.bai. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
641 OptionGroup* PileupOpts = Options::CreateOptionGroup("Pileup Options");
642 Options::AddValueOption("-fasta", "FASTA filename", "FASTA reference file", "", m_settings->HasFastaFilename, m_settings->FastaFilename, PileupOpts, "");
643 Options::AddOption("-mapqual", "print the mapping qualities", m_settings->IsPrintingPileupMapQualities, PileupOpts);
645 OptionGroup* SamOpts = Options::CreateOptionGroup("SAM Options");
646 Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingSamHeader, SamOpts);
649 ConvertTool::~ConvertTool(void) {
657 int ConvertTool::Help(void) {
658 Options::DisplayHelp();
662 int ConvertTool::Run(int argc, char* argv[]) {
664 // parse command line arguments
665 Options::Parse(argc, argv, 1);
667 // run internal ConvertTool implementation, return success/fail
668 m_impl = new ConvertToolPrivate(m_settings);
676 // ---------------------------------------------
677 // ConvertPileupFormatVisitor implementation
679 ConvertPileupFormatVisitor::ConvertPileupFormatVisitor(const RefVector& references,
680 const string& fastaFilename,
681 const bool isPrintingMapQualities,
685 , m_isPrintingMapQualities(isPrintingMapQualities)
687 , m_references(references)
689 // set up Fasta reader if file is provided
690 if ( !fastaFilename.empty() ) {
692 // check for FASTA index
693 string indexFilename = "";
694 if ( Utilities::FileExists(fastaFilename + ".fai") )
695 indexFilename = fastaFilename + ".fai";
698 if ( m_fasta.Open(fastaFilename, indexFilename) )
703 ConvertPileupFormatVisitor::~ConvertPileupFormatVisitor(void) {
704 // be sure to close Fasta reader
711 void ConvertPileupFormatVisitor::Visit(const PileupPosition& pileupData ) {
713 // skip if no alignments at this position
714 if ( pileupData.PileupAlignments.empty() ) return;
716 // retrieve reference name
717 const string& referenceName = m_references[pileupData.RefId].RefName;
718 const int& position = pileupData.Position;
720 // retrieve reference base from FASTA file, if one provided; otherwise default to 'N'
721 char referenceBase('N');
722 if ( m_hasFasta && (pileupData.Position < m_references[pileupData.RefId].RefLength) ) {
723 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position, referenceBase ) ) {
724 cerr << "Pileup error : Could not read reference base from FASTA file" << endl;
729 // get count of alleles at this position
730 const int numberAlleles = pileupData.PileupAlignments.size();
732 // -----------------------------------------------------------
733 // build strings based on alleles at this positionInAlignment
735 stringstream bases("");
736 stringstream baseQualities("");
737 stringstream mapQualities("");
739 // iterate over alignments at this pileup position
740 vector<PileupAlignment>::const_iterator pileupIter = pileupData.PileupAlignments.begin();
741 vector<PileupAlignment>::const_iterator pileupEnd = pileupData.PileupAlignments.end();
742 for ( ; pileupIter != pileupEnd; ++pileupIter ) {
743 const PileupAlignment pa = (*pileupIter);
744 const BamAlignment& ba = pa.Alignment;
746 // if beginning of read segment
747 if ( pa.IsSegmentBegin )
748 bases << '^' << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
750 // if current base is not a DELETION
751 if ( !pa.IsCurrentDeletion ) {
753 // get base at current position
754 char base = ba.QueryBases.at(pa.PositionInAlignment);
756 // if base matches reference
758 toupper(base) == toupper(referenceBase) ||
759 tolower(base) == tolower(referenceBase) )
761 base = (ba.IsReverseStrand() ? ',' : '.' );
764 // mismatches reference
765 else base = (ba.IsReverseStrand() ? tolower(base) : toupper(base) );
770 // if next position contains insertion
771 if ( pa.IsNextInsertion ) {
772 bases << '+' << pa.InsertionLength;
773 for (int i = 1; i <= pa.InsertionLength; ++i) {
774 char insertedBase = (char)ba.QueryBases.at(pa.PositionInAlignment + i);
775 bases << (ba.IsReverseStrand() ? (char)tolower(insertedBase) : (char)toupper(insertedBase) );
779 // if next position contains DELETION
780 else if ( pa.IsNextDeletion ) {
781 bases << '-' << pa.DeletionLength;
782 for (int i = 1; i <= pa.DeletionLength; ++i) {
783 char deletedBase('N');
784 if ( m_hasFasta && (pileupData.Position+i < m_references[pileupData.RefId].RefLength) ) {
785 if ( !m_fasta.GetBase(pileupData.RefId, pileupData.Position+i, deletedBase ) ) {
786 cerr << "Pileup error : Could not read reference base from FASTA file" << endl;
790 bases << (ba.IsReverseStrand() ? (char)tolower(deletedBase) : (char)toupper(deletedBase) );
795 // otherwise, DELETION
798 // if end of read segment
799 if ( pa.IsSegmentEnd ) bases << '$';
801 // store current base quality
802 baseQualities << ba.Qualities.at(pa.PositionInAlignment);
804 // save alignment map quality if desired
805 if ( m_isPrintingMapQualities )
806 mapQualities << ( ((int)ba.MapQuality > 93) ? (char)126 : (char)((int)ba.MapQuality+33) );
809 // ----------------------
813 // <refName> <1-based pos> <refBase> <numberAlleles> <bases> <qualities> [mapQuals]
815 const string TAB = "\t";
816 *m_out << referenceName << TAB
817 << position + 1 << TAB
818 << referenceBase << TAB
819 << numberAlleles << TAB
820 << bases.str() << TAB
821 << baseQualities.str() << TAB
822 << mapQualities.str() << endl;