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