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 << (int)tagData[index];
 242                     ++index;
 243                     break;
 244
 245                 case('c') :
 246                     out << (int)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             if ( tagData[index] == '\0') break;
 293         }
 294
 295         out << "}";
 296     }
 297
 298     out << "}" << endl;
 299
 300 }
 301
 302 // print BamAlignment in SAM format
 303 void BamTools::PrintSAM(ostream& out, const BamAlignment& a) {
 304
 305     // tab-delimited
 306     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
 307
 308     // write name & alignment flag
 309     out << a.Name << "\t" << a.AlignmentFlag << "\t";
 310
 311     // write reference name
 312     if ( (a.RefID >= 0) && (a.RefID < (int)references.size()) ) out << references[a.RefID].RefName << "\t";
 313     else out << "*\t";
 314
 315     // write position & map quality
 316     out << a.Position+1 << "\t" << a.MapQuality << "\t";
 317
 318     // write CIGAR
 319     const vector<CigarOp>& cigarData = a.CigarData;
 320     if ( cigarData.empty() ) out << "*\t";
 321     else {
 322         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
 323         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
 324         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 325             const CigarOp& op = (*cigarIter);
 326             out << op.Length << op.Type;
 327         }
 328         out << "\t";
 329     }
 330
 331     // write mate reference name, mate position, & insert size
 332     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)references.size()) ) {
 333         if ( a.MateRefID == a.RefID ) out << "=\t";
 334         else out << references[a.MateRefID].RefName << "\t";
 335         out << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
 336     }
 337     else out << "*\t0\t0\t";
 338
 339     // write sequence
 340     if ( a.QueryBases.empty() ) out << "*\t";
 341     else out << a.QueryBases << "\t";
 342
 343     // write qualities
 344     if ( a.Qualities.empty() ) out << "*";
 345     else out << a.Qualities;
 346
 347     // write tag data
 348     const char* tagData = a.TagData.c_str();
 349     const size_t tagDataLength = a.TagData.length();
 350
 351     size_t index = 0;
 352     while ( index < tagDataLength ) {
 353
 354         // write tag name
 355         out << "\t" << a.TagData.substr(index, 2) << ":";
 356         index += 2;
 357
 358         // get data type
 359         char type = a.TagData.at(index);
 360         ++index;
 361
 362         switch (type) {
 363             case('A') :
 364                 out << "A:" << tagData[index];
 365                 ++index;
 366                 break;
 367
 368             case('C') :
 369                 out << "i:" << (int)tagData[index];
 370                 ++index;
 371                 break;
 372
 373             case('c') :
 374                 out << "i:" << (int)tagData[index];
 375                 ++index;
 376                 break;
 377
 378             case('S') :
 379                 out << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
 380                 index += 2;
 381                 break;
 382
 383             case('s') :
 384                 out << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
 385                 index += 2;
 386                 break;
 387
 388             case('I') :
 389                 out << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
 390                 index += 4;
 391                 break;
 392
 393             case('i') :
 394                 out << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
 395                 index += 4;
 396                 break;
 397
 398             case('f') :
 399                 out << "f:" << BgzfData::UnpackFloat(&tagData[index]);
 400                 index += 4;
 401                 break;
 402
 403             case('d') :
 404                 out << "d:" << BgzfData::UnpackDouble(&tagData[index]);
 405                 index += 8;
 406                 break;
 407
 408             case('Z') :
 409             case('H') :
 410                 out << type << ":";
 411                 while (tagData[index]) {
 412                     out << tagData[index];
 413                     ++index;
 414                 }
 415                 ++index;
 416                 break;
 417         }
 418
 419         if ( tagData[index] == '\0') break;
 420     }
 421
 422     out << endl;
 423 }