bamtools_sam.cpp

   1 // ***************************************************************************
   2 // bamtools_sam.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 2 June 2010
   7 // ---------------------------------------------------------------------------
   8 // Prints a BAM file in the text-based SAM format.
   9 // ***************************************************************************
  10
  11 #include <cstdlib>
  12 #include <iostream>
  13 #include <sstream>
  14 #include <string>
  15
  16 #include "bamtools_sam.h"
  17 #include "bamtools_options.h"
  18 #include "BamReader.h"
  19 #include "BGZF.h"
  20
  21 using namespace std;
  22 using namespace BamTools;
  23
  24 static RefVector references;
  25
  26 // ---------------------------------------------
  27 // print BamAlignment in SAM format
  28
  29 void PrintSAM(const BamAlignment& a) {
  30
  31     // tab-delimited
  32     // <QNAME> <FLAG> <RNAME> <POS> <MAPQ> <CIGAR> <MRNM> <MPOS> <ISIZE> <SEQ> <QUAL> [ <TAG>:<VTYPE>:<VALUE> [...] ]
  33
  34     ostringstream sb("");
  35
  36     // write name & alignment flag
  37     cout << a.Name << "\t" << a.AlignmentFlag << "\t";
  38
  39     // write reference name
  40     if ( (a.RefID >= 0) && (a.RefID < (int)references.size()) ) cout << references[a.RefID].RefName << "\t";
  41     else cout << "*\t";
  42
  43     // write position & map quality
  44     cout << a.Position+1 << "\t" << a.MapQuality << "\t";
  45
  46     // write CIGAR
  47     const vector<CigarOp>& cigarData = a.CigarData;
  48     if ( cigarData.empty() ) cout << "*\t";
  49     else {
  50         vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
  51         vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();
  52         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
  53             const CigarOp& op = (*cigarIter);
  54             cout << op.Length << op.Type;
  55         }
  56         cout << "\t";
  57     }
  58
  59     // write mate reference name, mate position, & insert size
  60     if ( a.IsPaired() && (a.MateRefID >= 0) && (a.MateRefID < (int)references.size()) ) {
  61         if ( a.MateRefID == a.RefID ) cout << "=\t";
  62         else cout << references[a.MateRefID].RefName << "\t";
  63         cout << a.MatePosition+1 << "\t" << a.InsertSize << "\t";
  64     }
  65     else cout << "*\t0\t0\t";
  66
  67     // write sequence
  68     if ( a.QueryBases.empty() ) cout << "*\t";
  69     else cout << a.QueryBases << "\t";
  70
  71     // write qualities
  72     if ( a.Qualities.empty() ) cout << "*";
  73     else cout << a.Qualities;
  74
  75     // write tag data
  76     const char* tagData = a.TagData.c_str();
  77     const size_t tagDataLength = a.TagData.length();
  78     size_t index = 0;
  79     while ( index < tagDataLength ) {
  80
  81         // write tag name
  82         cout << "\t" << a.TagData.substr(index, 2) << ":";
  83         index += 2;
  84
  85         // get data type
  86         char type = a.TagData.at(index);
  87         ++index;
  88
  89         switch (type) {
  90             case('A') :
  91                 cout << "A:" << tagData[index];
  92                 ++index;
  93                 break;
  94
  95             case('C') :
  96                 cout << "i:" << atoi(&tagData[index]);
  97                 ++index;
  98                 break;
  99
 100             case('c') :
 101                 cout << "i:" << atoi(&tagData[index]);
 102                 ++index;
 103                 break;
 104
 105             case('S') :
 106                 cout << "i:" << BgzfData::UnpackUnsignedShort(&tagData[index]);
 107                 index += 2;
 108                 break;
 109
 110             case('s') :
 111                 cout << "i:" << BgzfData::UnpackSignedShort(&tagData[index]);
 112                 index += 2;
 113                 break;
 114
 115             case('I') :
 116                 cout << "i:" << BgzfData::UnpackUnsignedInt(&tagData[index]);
 117                 index += 4;
 118                 break;
 119
 120             case('i') :
 121                 cout << "i:" << BgzfData::UnpackSignedInt(&tagData[index]);
 122                 index += 4;
 123                 break;
 124
 125             case('f') :
 126                 cout << "f:" << BgzfData::UnpackFloat(&tagData[index]);
 127                 index += 4;
 128                 break;
 129
 130             case('d') :
 131                 cout << "d:" << BgzfData::UnpackDouble(&tagData[index]);
 132                 index += 8;
 133                 break;
 134
 135             case('Z') :
 136             case('H') :
 137                 cout << type << ":";
 138                 while (tagData[index]) {
 139                     cout << tagData[index];
 140                     ++index;
 141                 }
 142                 ++index;
 143                 break;
 144         }
 145     }
 146
 147     // write stream to stdout
 148     cout << sb.str() << endl;
 149 }
 150
 151 // ---------------------------------------------
 152 // SamSettings implementation
 153
 154 struct SamTool::SamSettings {
 155
 156     // flags
 157     bool HasInputBamFilename;
 158     bool HasMaximumOutput;
 159     bool IsOmittingHeader;
 160
 161     // filenames
 162     string InputBamFilename;
 163
 164     // other parameters
 165     int MaximumOutput;
 166
 167     // constructor
 168     SamSettings(void)
 169         : HasInputBamFilename(false)
 170         , HasMaximumOutput(false)
 171         , IsOmittingHeader(false)
 172         , InputBamFilename(Options::StandardIn())
 173     { }
 174 };
 175
 176 // ---------------------------------------------
 177 // SamTool implementation
 178
 179 SamTool::SamTool(void)
 180     : AbstractTool()
 181     , m_settings(new SamSettings)
 182 {
 183     // set program details
 184     Options::SetProgramInfo("bamtools sam", "prints BAM file in SAM text format", "-in <filename>");
 185
 186     // set up options
 187     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 188     Options::AddValueOption("-in", "BAM filename", "the input BAM file", "", m_settings->HasInputBamFilename, m_settings->InputBamFilename, IO_Opts, Options::StandardIn());
 189
 190     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
 191     Options::AddOption("-noheader", "omit the SAM header from output", m_settings->IsOmittingHeader, FilterOpts);
 192     Options::AddValueOption("-num", "N", "maximum number of alignments to output", "", m_settings->HasMaximumOutput, m_settings->MaximumOutput, FilterOpts);
 193 }
 194
 195 SamTool::~SamTool(void) {
 196     delete m_settings;
 197     m_settings = 0;
 198 }
 199
 200 int SamTool::Help(void) {
 201     Options::DisplayHelp();
 202     return 0;
 203 }
 204
 205 int SamTool::Run(int argc, char* argv[]) {
 206
 207     // parse command line arguments
 208     Options::Parse(argc, argv, 1);
 209
 210     // open our BAM reader
 211     BamReader reader;
 212     reader.Open(m_settings->InputBamFilename);
 213
 214     // if header desired, retrieve and print to stdout
 215     if ( !m_settings->IsOmittingHeader ) {
 216         string header = reader.GetHeaderText();
 217         cout << header << endl;
 218     }
 219
 220     // store reference data
 221     references = reader.GetReferenceData();
 222
 223     // print all alignments to stdout in SAM format
 224     if ( !m_settings->HasMaximumOutput ) {
 225         BamAlignment ba;
 226         while( reader.GetNextAlignment(ba) ) {
 227             PrintSAM(ba);
 228         }
 229     }
 230
 231     // print first N alignments to stdout in SAM format
 232     else {
 233         BamAlignment ba;
 234         int alignmentsPrinted = 0;
 235         while ( reader.GetNextAlignment(ba) && (alignmentsPrinted < m_settings->MaximumOutput) ) {
 236             PrintSAM(ba);
 237             ++alignmentsPrinted;
 238         }
 239     }
 240
 241     // clean & exit
 242     reader.Close();
 243     return 0;
 244 }