src/toolkit/bamtools_stats.cpp

   1 // ***************************************************************************
   2 // bamtools_cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 5 October 2010
   7 // ---------------------------------------------------------------------------
   8 // Prints general alignment statistics for BAM file(s).
   9 // ***************************************************************************
  10
  11 #include <cmath>
  12 #include <algorithm>
  13 #include <functional>
  14 #include <iostream>
  15 #include <numeric>
  16 #include <string>
  17 #include <vector>
  18
  19 #include "bamtools_stats.h"
  20 #include "bamtools_options.h"
  21 #include "BamMultiReader.h"
  22 using namespace std;
  23 using namespace BamTools;
  24
  25 // ---------------------------------------------
  26 // StatsSettings implementation
  27
  28 struct StatsTool::StatsSettings {
  29
  30     // flags
  31     bool HasInput;
  32     bool IsShowingInsertSizeSummary;
  33
  34     // filenames
  35     vector<string> InputFiles;
  36
  37     // constructor
  38     StatsSettings(void)
  39         : HasInput(false)
  40         , IsShowingInsertSizeSummary(false)
  41     { }
  42 };
  43
  44 // ---------------------------------------------
  45 // StatsToolPrivate implementation
  46
  47 struct StatsTool::StatsToolPrivate {
  48
  49     // ctor & dtor
  50     public:
  51         StatsToolPrivate(StatsTool::StatsSettings* _settings);
  52         ~StatsToolPrivate(void);
  53
  54     // 'public' interface
  55     public:
  56         bool Run(void);
  57
  58     // internal methods
  59     private:
  60         bool CalculateMedian(vector<int>& data, double& median);
  61         void PrintStats(void);
  62         void ProcessAlignment(const BamAlignment& al);
  63
  64     // data members
  65     private:
  66         StatsTool::StatsSettings* settings;
  67         unsigned int numReads;
  68         unsigned int numPaired;
  69         unsigned int numProperPair;
  70         unsigned int numMapped;
  71         unsigned int numBothMatesMapped;
  72         unsigned int numForwardStrand;
  73         unsigned int numReverseStrand;
  74         unsigned int numFirstMate;
  75         unsigned int numSecondMate;
  76         unsigned int numSingletons;
  77         unsigned int numFailedQC;
  78         unsigned int numDuplicates;
  79         vector<int> insertSizes;
  80 };
  81
  82 StatsTool::StatsToolPrivate::StatsToolPrivate(StatsTool::StatsSettings* _settings)
  83     : settings(_settings)
  84     , numReads(0)
  85     , numPaired(0)
  86     , numProperPair(0)
  87     , numMapped(0)
  88     , numBothMatesMapped(0)
  89     , numForwardStrand(0)
  90     , numReverseStrand(0)
  91     , numFirstMate(0)
  92     , numSecondMate(0)
  93     , numSingletons(0)
  94     , numFailedQC(0)
  95     , numDuplicates(0)
  96 {
  97     insertSizes.reserve(100000);
  98 }
  99
 100 StatsTool::StatsToolPrivate::~StatsToolPrivate(void) { }
 101
 102 // median is of type double because in the case of even number of data elements, we need to return the average of middle 2 elements
 103 bool StatsTool::StatsToolPrivate::CalculateMedian(vector<int>& data, double& median) {
 104
 105     // check that data exists
 106     if ( data.empty() ) return false;
 107
 108     // find middle element
 109     size_t middleIndex = data.size() / 2;
 110     vector<int>::iterator target = data.begin() + middleIndex;
 111     nth_element(data.begin(), target, data.end());
 112
 113     // odd number of elements
 114     if ( (data.size() % 2) != 0) {
 115         median = (double)(*target);
 116         return true;
 117     }
 118
 119     // even number of elements
 120     else {
 121         double rightTarget = (double)(*target);
 122         vector<int>::iterator leftTarget = target - 1;
 123         nth_element(data.begin(), leftTarget, data.end());
 124         median = (double)((rightTarget+*leftTarget)/2.0);
 125         return true;
 126     }
 127 }
 128
 129 // print BAM file alignment stats
 130 void StatsTool::StatsToolPrivate::PrintStats(void) {
 131
 132     cout << endl;
 133     cout << "**********************************************" << endl;
 134     cout << "Stats for BAM file(s): " << endl;
 135     cout << "**********************************************" << endl;
 136     cout << endl;
 137     cout << "Total reads:       " << numReads << endl;
 138     cout << "Mapped reads:      " << numMapped << "\t(" << ((float)numMapped/numReads)*100 << "%)" << endl;
 139     cout << "Forward strand:    " << numForwardStrand << "\t(" << ((float)numForwardStrand/numReads)*100 << "%)" << endl;
 140     cout << "Reverse strand:    " << numReverseStrand << "\t(" << ((float)numReverseStrand/numReads)*100 << "%)" << endl;
 141     cout << "Failed QC:         " << numFailedQC << "\t(" << ((float)numFailedQC/numReads)*100 << "%)" << endl;
 142     cout << "Duplicates:        " << numDuplicates << "\t(" << ((float)numDuplicates/numReads)*100 << "%)" << endl;
 143     cout << "Paired-end reads:  " << numPaired << "\t(" << ((float)numPaired/numReads)*100 << "%)" << endl;
 144
 145     if ( numPaired != 0 ) {
 146         cout << "'Proper-pairs':    " << numProperPair << "\t(" << ((float)numProperPair/numPaired)*100 << "%)" << endl;
 147         cout << "Both pairs mapped: " << numBothMatesMapped << "\t(" << ((float)numBothMatesMapped/numPaired)*100 << "%)" << endl;
 148         cout << "Read 1:            " << numFirstMate << endl;
 149         cout << "Read 2:            " << numSecondMate << endl;
 150         cout << "Singletons:        " << numSingletons << "\t(" << ((float)numSingletons/numPaired)*100 << "%)" << endl;
 151     }
 152
 153     if ( settings->IsShowingInsertSizeSummary ) {
 154
 155         double avgInsertSize = 0.0;
 156         if ( !insertSizes.empty() ) {
 157             avgInsertSize = ( accumulate(insertSizes.begin(), insertSizes.end(), 0.0) / (double)insertSizes.size() );
 158             cout << "Average insert size (absolute value): " << avgInsertSize << endl;
 159         }
 160
 161         double medianInsertSize = 0.0;
 162         if ( CalculateMedian(insertSizes, medianInsertSize) )
 163             cout << "Median insert size (absolute value): " << medianInsertSize << endl;
 164     }
 165     cout << endl;
 166 }
 167
 168 // use current input alignment to update BAM file alignment stats
 169 void StatsTool::StatsToolPrivate::ProcessAlignment(const BamAlignment& al) {
 170
 171     // increment total alignment counter
 172     ++numReads;
 173
 174     // check the paired-independent flags
 175     if ( al.IsDuplicate() ) ++numDuplicates;
 176     if ( al.IsFailedQC()  ) ++numFailedQC;
 177     if ( al.IsMapped()    ) ++numMapped;
 178
 179     // check forward/reverse strand
 180     if ( al.IsReverseStrand() )
 181         ++numReverseStrand;
 182     else
 183         ++numForwardStrand;
 184
 185     // if alignment is paired-end
 186     if ( al.IsPaired() ) {
 187
 188         // increment PE counter
 189         ++numPaired;
 190
 191         // increment first mate/second mate counters
 192         if ( al.IsFirstMate()  ) ++numFirstMate;
 193         if ( al.IsSecondMate() ) ++numSecondMate;
 194
 195         // if alignment is mapped, check mate status
 196         if ( al.IsMapped() ) {
 197             // if mate mapped
 198             if ( al.IsMateMapped() )
 199                 ++numBothMatesMapped;
 200             // else singleton
 201             else
 202                 ++numSingletons;
 203         }
 204
 205         // check for explicit proper pair flag
 206         if ( al.IsProperPair() ) ++numProperPair;
 207
 208         // store insert size for first mate
 209         if ( settings->IsShowingInsertSizeSummary && al.IsFirstMate() && (al.InsertSize != 0) ) {
 210             int insertSize = abs(al.InsertSize);
 211             insertSizes.push_back( insertSize );
 212         }
 213     }
 214 }
 215
 216 bool StatsTool::StatsToolPrivate::Run() {
 217
 218     // opens the BAM files without checking for indexes
 219     BamMultiReader reader;
 220     if ( !reader.Open(settings->InputFiles, false, true) ) {
 221         cerr << "Could not open input BAM file(s)... quitting." << endl;
 222         reader.Close();
 223         return false;
 224     }
 225
 226     // plow through file, keeping track of stats
 227     BamAlignment al;
 228     while ( reader.GetNextAlignmentCore(al) ) {
 229         ProcessAlignment(al);
 230     }
 231
 232     // print stats
 233     PrintStats();
 234
 235     // clean and exit
 236     reader.Close();
 237     return true;
 238 }
 239
 240 // ---------------------------------------------
 241 // StatsTool implementation
 242
 243 StatsTool::StatsTool(void)
 244     : AbstractTool()
 245     , m_settings(new StatsSettings)
 246     , m_impl(0)
 247 {
 248     // set program details
 249     Options::SetProgramInfo("bamtools stats", "prints general alignment statistics", "[-in <filename> -in <filename> ... ]");
 250
 251     // set up options
 252     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 253     Options::AddValueOption("-in", "BAM filename", "the input BAM file", "", m_settings->HasInput,  m_settings->InputFiles,  IO_Opts, Options::StandardIn());
 254
 255     OptionGroup* AdditionalOpts = Options::CreateOptionGroup("Additional Stats");
 256     Options::AddOption("-insert", "summarize insert size data", m_settings->IsShowingInsertSizeSummary, AdditionalOpts);
 257 }
 258
 259 StatsTool::~StatsTool(void) {
 260     delete m_settings;
 261     m_settings = 0;
 262
 263     delete m_impl;
 264     m_impl = 0;
 265 }
 266
 267 int StatsTool::Help(void) {
 268     Options::DisplayHelp();
 269     return 0;
 270 }
 271
 272 int StatsTool::Run(int argc, char* argv[]) {
 273
 274     // parse command line arguments
 275     Options::Parse(argc, argv, 1);
 276
 277     // set to default input if none provided
 278     if ( !m_settings->HasInput )
 279         m_settings->InputFiles.push_back(Options::StandardIn());
 280
 281     // run internal StatsTool implementation, return success/fail
 282     m_impl = new StatsToolPrivate(m_settings);
 283
 284     if ( m_impl->Run() ) return 0;
 285     else return 1;
 286 }