bamtools_convert.cpp

   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 June 2010
   7 // ---------------------------------------------------------------------------
   8 // Converts between BAM and a number of other formats
   9 // ***************************************************************************
  10
  11 #include <iostream>
  12 #include <sstream>
  13 #include <string>
  14 #include <vector>
  15
  16 #include "bamtools_convert.h"
  17 //#include "bamtools_format.h"
  18 #include "bamtools_options.h"
  19 #include "bamtools_utilities.h"
  20 #include "BGZF.h"
  21 #include "BamReader.h"
  22 #include "BamMultiReader.h"
  23
  24 using namespace std;
  25 using namespace BamTools;
  26
  27 static RefVector references;
  28
  29 namespace BamTools {
  30
  31     static const string FORMAT_FASTA = "fasta";
  32     static const string FORMAT_FASTQ = "fastq";
  33     static const string FORMAT_JSON  = "json";
  34     static const string FORMAT_SAM   = "sam";
  35
  36     void PrintFASTA(ostream& out, const BamAlignment& a);
  37     void PrintFASTQ(ostream& out, const BamAlignment& a);
  38     void PrintJSON(ostream& out, const BamAlignment& a);
  39     void PrintSAM(ostream& out, const BamAlignment& a);
  40
  41 } // namespace BamTools
  42
  43 // ---------------------------------------------
  44 // ConvertSettings implementation
  45
  46 struct ConvertTool::ConvertSettings {
  47
  48     // flags
  49     bool HasInputBamFilenames;
  50     bool HasOutputBamFilename;
  51     bool HasFormat;
  52     bool HasRegion;
  53
  54     // filenames
  55     //string InputFilename;
  56     vector<string> InputFiles;
  57     string OutputFilename;
  58     string Format;
  59
  60     // region
  61     string Region;
  62
  63     // constructor
  64     ConvertSettings(void)
  65         : HasInputBamFilenames(false)
  66         , HasOutputBamFilename(false)
  67         //, InputFilename(Options::StandardIn())
  68         , OutputFilename(Options::StandardOut())
  69     { }
  70 };
  71
  72 // ---------------------------------------------
  73 // ConvertTool implementation
  74
  75 ConvertTool::ConvertTool(void)
  76     : AbstractTool()
  77     , m_settings(new ConvertSettings)
  78 {
  79     // set program details
  80     Options::SetProgramInfo("bamtools convert", "converts between BAM and a number of other formats", "-in <filename> -out <filename> -format <FORMAT>");
  81
  82     // set up options
  83     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  84     //Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInputBamFilename,  m_settings->InputFilename,  IO_Opts, Options::StandardIn());
  85     Options::AddValueOption("-in",  "BAM filename", "the input BAM file(s)", "", m_settings->HasInputBamFilenames,  m_settings->InputFiles, IO_Opts, Options::StandardIn());
  86     Options::AddValueOption("-out",    "BAM filename", "the output BAM file",   "", m_settings->HasOutputBamFilename, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
  87     Options::AddValueOption("-format", "FORMAT", "the output file format - see README for recognized formats", "", m_settings->HasFormat, m_settings->Format, IO_Opts);
  88     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  89     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);
  90 }
  91
  92 ConvertTool::~ConvertTool(void) {
  93     delete m_settings;
  94     m_settings = 0;
  95 }
  96
  97 int ConvertTool::Help(void) {
  98     Options::DisplayHelp();
  99     return 0;
 100 }
 101
 102 int ConvertTool::Run(int argc, char* argv[]) {
 103
 104     bool convertedOk = true;
 105
 106     // parse command line arguments
 107     Options::Parse(argc, argv, 1);
 108
 109     // open files
 110     BamMultiReader reader;
 111     reader.Open(m_settings->InputFiles);
 112     references = reader.GetReferenceData();
 113
 114     BamRegion region;
 115     if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 116         if ( !reader.SetRegion(region) ) {
 117            cerr << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
 118         }
 119     }
 120
 121     // ----------------------------------------
 122     // do conversion,depending on desired output format
 123
 124     // FASTA
 125     if ( m_settings->Format == FORMAT_FASTA ) {
 126         //cout << "Converting to FASTA" << endl;
 127     }
 128
 129     // FASTQ
 130     else if ( m_settings->Format == FORMAT_FASTQ) {
 131         //cout << "Converting to FASTQ" << endl;
 132     }
 133
 134     // JSON
 135     else if ( m_settings->Format == FORMAT_JSON ) {
 136         //cout << "Converting to JSON" << endl;
 137         BamAlignment alignment;
 138         while ( reader.GetNextAlignment(alignment) ) {
 139             PrintJSON(cout, alignment);
 140         }
 141
 142     }
 143
 144     // SAM
 145     else if ( m_settings->Format == FORMAT_SAM ) {
 146         BamAlignment alignment;
 147         while ( reader.GetNextAlignment(alignment) ) {
 148             PrintSAM(cout, alignment);
 149         }
 150     }
 151
 152     // uncrecognized format
 153     else {
 154         cerr << "Unrecognized format: " << m_settings->Format << endl;
 155         cerr << "Please see help|README (?) for details on supported formats " << endl;
 156         convertedOk = false;
 157     }
 158
 159     // ------------------------
 160     // clean up & exit
 161     reader.Close();
 162     return (int)convertedOk;
 163 }
 164
 165 // ----------------------------------------------------------
 166 // Conversion/output methods
 167 // ----------------------------------------------------------
 168
 169 // print BamAlignment in FASTA format
 170 void BamTools::PrintFASTA(ostream& out, const BamAlignment& a) {
 171
 172 }
 173
 174 // print BamAlignment in FASTQ format
 175 void BamTools::PrintFASTQ(ostream& out, const BamAlignment& a) {
 176
 177 }
 178
 179 // print BamAlignment in JSON format
 180 void BamTools::PrintJSON(ostream& out, const BamAlignment& a) {
 181
 182     // tab-delimited
 183     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 184
 185     // write name & alignment flag
 186     out << "{\"name\":\"" << a.Name << "\",\"alignmentFlag\":\""
 187         << a.AlignmentFlag << "\",";
 188
 189     // write reference name
 190     if ( (a.RefID >= 0) && (a.RefID < (int)references.size()) ) out << "\"reference\":\"" << references[a.RefID].RefName << "\",";
 191     //else out << "*\t";
 192
 193     // write position & map quality
 194     out << "\"position\":" << a.Position+1 << ",\"mapQuality\":" << a.MapQuality << ",";
 195
 196     // write CIGAR
 197     const vector<CigarOp>& cigarData = a.CigarData;
 198     if ( !cigarData.empty() ) {
 199         out << "\"cigar\":[";
 200         vector<CigarOp>::const_iterator cigarBegin = cigarData.begin();
 201         vector<CigarOp>::const_iterator cigarIter = cigarBegin;
 202         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 203         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 204             const CigarOp& op = (*cigarIter);
 205             if (cigarIter != cigarBegin) out << ",";
 206             out << "\"" << op.Length << op.Type << "\"";
 207         }
 208         out << "],";
 209     }
 210
 211     // write mate reference name, mate position, & insert size
 212     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)references.size()) ) {
 213         out << "\"mate\":{"
 214             << "\"reference\":\"" << references[a.MateRefID].RefName << "\","
 215             << "\"position\":" << a.MatePosition+1
 216             << ",\"insertSize\":" << a.InsertSize << "},";
 217     }
 218
 219     // write sequence
 220     if ( !a.QueryBases.empty() ) {
 221         out << "\"queryBases\":\"" << a.QueryBases << "\",";
 222     }
 223
 224     // write qualities
 225     if ( !a.Qualities.empty() ) {
 226         out << "\"qualities\":\"" << a.Qualities << "\",";
 227     }
 228
 229     // write tag data
 230     const char* tagData = a.TagData.c_str();
 231     const size_t tagDataLength = a.TagData.length();
 232     size_t index = 0;
 233     if (index < tagDataLength) {
 234
 235         out << "\"tags\":{";
 236
 237         while ( index < tagDataLength ) {
 238
 239             if (index > 0)
 240                 out << ",";
 241
 242             // write tag name
 243             out << "\"" << a.TagData.substr(index, 2) << "\":";
 244             index += 2;
 245
 246             // get data type
 247             char type = a.TagData.at(index);
 248             ++index;
 249
 250             switch (type) {
 251                 case('A') :
 252                     out << "\"" << tagData[index] << "\"";
 253                     ++index;
 254                     break;
 255
 256                 case('C') :
 257                     out << (int)tagData[index];
 258                     ++index;
 259                     break;
 260
 261                 case('c') :
 262                     out << (int)tagData[index];
 263                     ++index;
 264                     break;
 265
 266                 case('S') :
 267                     out << BgzfData::UnpackUnsignedShort(&tagData[index]);
 268                     index += 2;
 269                     break;
 270
 271                 case('s') :
 272                     out << BgzfData::UnpackSignedShort(&tagData[index]);
 273                     index += 2;
 274                     break;
 275
 276                 case('I') :
 277                     out << BgzfData::UnpackUnsignedInt(&tagData[index]);
 278                     index += 4;
 279                     break;
 280
 281                 case('i') :
 282                     out << BgzfData::UnpackSignedInt(&tagData[index]);
 283                     index += 4;
 284                     break;
 285
 286                 case('f') :
 287                     out << BgzfData::UnpackFloat(&tagData[index]);
 288                     index += 4;
 289                     break;
 290
 291                 case('d') :
 292                     out << BgzfData::UnpackDouble(&tagData[index]);
 293                     index += 8;
 294                     break;
 295
 296                 case('Z') :
 297                 case('H') :
 298                     out << "\"";
 299                     while (tagData[index]) {
 300                         out << tagData[index];
 301                         ++index;
 302                     }
 303                     out << "\"";
 304                     ++index;
 305                     break;
 306             }
 307
 308             if ( tagData[index] == '\0') break;
 309         }
 310
 311         out << "}";
 312     }
 313
 314     out << "}" << endl;
 315
 316 }
 317
 318 // print BamAlignment in SAM format
 319 void BamTools::PrintSAM(ostream& out, const BamAlignment& a) {
 320
 321     // tab-delimited
 322     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 323
 324     // write name & alignment flag
 325     out << a.Name << "\t" << a.AlignmentFlag << "\t";
 326
 327     // write reference name
 328     if ( (a.RefID >= 0) && (a.RefID < (int)references.size()) ) out << references[a.RefID].RefName << "\t";
 329     else out << "*\t";
 330
 331     // write position & map quality
 332     out << a.Position+1 << "\t" << a.MapQuality << "\t";
 333
 334     // write CIGAR
 335     const vector<CigarOp>& cigarData = a.CigarData;
 336     if ( cigarData.empty() ) out << "*\t";
 337     else {
 338         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 339         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 340         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 341             const CigarOp& op = (*cigarIter);
 342             out << op.Length << op.Type;
 343         }
 344         out << "\t";
 345     }
 346
 347     // write mate reference name, mate position, & insert size
 348     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)references.size()) ) {
 349         if ( a.MateRefID == a.RefID ) out << "=\t";
 350         else out << references[a.MateRefID].RefName << "\t";
 351         out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 352     }
 353     else out << "*\t0\t0\t";
 354
 355     // write sequence
 356     if ( a.QueryBases.empty() ) out << "*\t";
 357     else out << a.QueryBases << "\t";
 358
 359     // write qualities
 360     if ( a.Qualities.empty() ) out << "*";
 361     else out << a.Qualities;
 362
 363     // write tag data
 364     const char* tagData = a.TagData.c_str();
 365     const size_t tagDataLength = a.TagData.length();
 366
 367     size_t index = 0;
 368     while ( index < tagDataLength ) {
 369
 370         // write tag name
 371         out << "\t" << a.TagData.substr(index, 2) << ":";
 372         index += 2;
 373
 374         // get data type
 375         char type = a.TagData.at(index);
 376         ++index;
 377
 378         switch (type) {
 379             case('A') :
 380                 out << "A:" << tagData[index];
 381                 ++index;
 382                 break;
 383
 384             case('C') :
 385                 out << "i:" << (int)tagData[index];
 386                 ++index;
 387                 break;
 388
 389             case('c') :
 390                 out << "i:" << (int)tagData[index];
 391                 ++index;
 392                 break;
 393
 394             case('S') :
 395                 out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
 396                 index += 2;
 397                 break;
 398
 399             case('s') :
 400                 out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
 401                 index += 2;
 402                 break;
 403
 404             case('I') :
 405                 out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
 406                 index += 4;
 407                 break;
 408
 409             case('i') :
 410                 out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
 411                 index += 4;
 412                 break;
 413
 414             case('f') :
 415                 out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
 416                 index += 4;
 417                 break;
 418
 419             case('d') :
 420                 out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
 421                 index += 8;
 422                 break;
 423
 424             case('Z') :
 425             case('H') :
 426                 out << type << ":";
 427                 while (tagData[index]) {
 428                     out << tagData[index];
 429                     ++index;
 430                 }
 431                 ++index;
 432                 break;
 433         }
 434
 435         if ( tagData[index] == '\0') break;
 436     }
 437
 438     out << endl;
 439 }