src/toolkit/bamtools_stats.cpp

   1 // ***************************************************************************
   2 // bamtools_cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 7 April 2011
   6 // ---------------------------------------------------------------------------
   7 // Prints general alignment statistics for BAM file(s).
   8 // ***************************************************************************
   9
  10 #include "bamtools_stats.h"
  11
  12 #include <api/BamMultiReader.h>
  13 #include <utils/bamtools_options.h>
  14 using namespace BamTools;
  15
  16 #include <cmath>
  17 #include <algorithm>
  18 #include <functional>
  19 #include <iostream>
  20 #include <numeric>
  21 #include <string>
  22 #include <vector>
  23 using namespace std;
  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* m_settings;
  67         unsigned int m_numReads;
  68         unsigned int m_numPaired;
  69         unsigned int m_numProperPair;
  70         unsigned int m_numMapped;
  71         unsigned int m_numBothMatesMapped;
  72         unsigned int m_numForwardStrand;
  73         unsigned int m_numReverseStrand;
  74         unsigned int m_numFirstMate;
  75         unsigned int m_numSecondMate;
  76         unsigned int m_numSingletons;
  77         unsigned int m_numFailedQC;
  78         unsigned int m_numDuplicates;
  79         vector<int> m_insertSizes;
  80 };
  81
  82 StatsTool::StatsToolPrivate::StatsToolPrivate(StatsTool::StatsSettings* settings)
  83     : m_settings(settings)
  84     , m_numReads(0)
  85     , m_numPaired(0)
  86     , m_numProperPair(0)
  87     , m_numMapped(0)
  88     , m_numBothMatesMapped(0)
  89     , m_numForwardStrand(0)
  90     , m_numReverseStrand(0)
  91     , m_numFirstMate(0)
  92     , m_numSecondMate(0)
  93     , m_numSingletons(0)
  94     , m_numFailedQC(0)
  95     , m_numDuplicates(0)
  96 {
  97     m_insertSizes.reserve(100000);
  98 }
  99
 100 // median is of type double because in the case of even number of data elements,
 101 // we need to return the average of middle 2 elements
 102 bool StatsTool::StatsToolPrivate::CalculateMedian(vector<int>& data, double& median) {
 103
 104     // skip if data empty
 105     if ( data.empty() ) return false;
 106
 107     // find middle element
 108     size_t middleIndex = data.size() / 2;
 109     vector<int>::iterator target = data.begin() + middleIndex;
 110     nth_element(data.begin(), target, data.end());
 111
 112     // odd number of elements
 113     if ( (data.size() % 2) != 0) {
 114         median = (double)(*target);
 115         return true;
 116     }
 117
 118     // even number of elements
 119     else {
 120         double rightTarget = (double)(*target);
 121         vector<int>::iterator leftTarget = target - 1;
 122         nth_element(data.begin(), leftTarget, data.end());
 123         median = (double)((rightTarget+*leftTarget)/2.0);
 124         return true;
 125     }
 126 }
 127
 128 // print BAM file alignment stats
 129 void StatsTool::StatsToolPrivate::PrintStats(void) {
 130
 131     cout << endl;
 132     cout << "**********************************************" << endl;
 133     cout << "Stats for BAM file(s): " << endl;
 134     cout << "**********************************************" << endl;
 135     cout << endl;
 136     cout << "Total reads:       " << m_numReads << endl;
 137     cout << "Mapped reads:      " << m_numMapped << "\t(" << ((float)m_numMapped/m_numReads)*100 << "%)" << endl;
 138     cout << "Forward strand:    " << m_numForwardStrand << "\t(" << ((float)m_numForwardStrand/m_numReads)*100 << "%)" << endl;
 139     cout << "Reverse strand:    " << m_numReverseStrand << "\t(" << ((float)m_numReverseStrand/m_numReads)*100 << "%)" << endl;
 140     cout << "Failed QC:         " << m_numFailedQC << "\t(" << ((float)m_numFailedQC/m_numReads)*100 << "%)" << endl;
 141     cout << "Duplicates:        " << m_numDuplicates << "\t(" << ((float)m_numDuplicates/m_numReads)*100 << "%)" << endl;
 142     cout << "Paired-end reads:  " << m_numPaired << "\t(" << ((float)m_numPaired/m_numReads)*100 << "%)" << endl;
 143
 144     if ( m_numPaired != 0 ) {
 145         cout << "'Proper-pairs':    " << m_numProperPair << "\t(" << ((float)m_numProperPair/m_numPaired)*100 << "%)" << endl;
 146         cout << "Both pairs mapped: " << m_numBothMatesMapped << "\t(" << ((float)m_numBothMatesMapped/m_numPaired)*100 << "%)" << endl;
 147         cout << "Read 1:            " << m_numFirstMate << endl;
 148         cout << "Read 2:            " << m_numSecondMate << endl;
 149         cout << "Singletons:        " << m_numSingletons << "\t(" << ((float)m_numSingletons/m_numPaired)*100 << "%)" << endl;
 150     }
 151
 152     if ( m_settings->IsShowingInsertSizeSummary ) {
 153
 154         double avgInsertSize = 0.0;
 155         if ( !m_insertSizes.empty() ) {
 156             avgInsertSize = ( accumulate(m_insertSizes.begin(), m_insertSizes.end(), 0.0) / (double)m_insertSizes.size() );
 157             cout << "Average insert size (absolute value): " << avgInsertSize << endl;
 158         }
 159
 160         double medianInsertSize = 0.0;
 161         if ( CalculateMedian(m_insertSizes, medianInsertSize) )
 162             cout << "Median insert size (absolute value): " << medianInsertSize << endl;
 163     }
 164     cout << endl;
 165 }
 166
 167 // use current input alignment to update BAM file alignment stats
 168 void StatsTool::StatsToolPrivate::ProcessAlignment(const BamAlignment& al) {
 169
 170     // increment total alignment counter
 171     ++m_numReads;
 172
 173     // incrememt counters for pairing-independent flags
 174     if ( al.IsDuplicate() ) ++m_numDuplicates;
 175     if ( al.IsFailedQC()  ) ++m_numFailedQC;
 176     if ( al.IsMapped()    ) ++m_numMapped;
 177
 178     // increment strand counters
 179     if ( al.IsReverseStrand() )
 180         ++m_numReverseStrand;
 181     else
 182         ++m_numForwardStrand;
 183
 184     // if alignment is paired-end
 185     if ( al.IsPaired() ) {
 186
 187         // increment PE counter
 188         ++m_numPaired;
 189
 190         // increment first mate/second mate counters
 191         if ( al.IsFirstMate()  ) ++m_numFirstMate;
 192         if ( al.IsSecondMate() ) ++m_numSecondMate;
 193
 194         // if alignment is mapped, check mate status
 195         if ( al.IsMapped() ) {
 196             // if mate mapped
 197             if ( al.IsMateMapped() )
 198                 ++m_numBothMatesMapped;
 199             // else singleton
 200             else
 201                 ++m_numSingletons;
 202         }
 203
 204         // check for explicit proper pair flag
 205         if ( al.IsProperPair() ) ++m_numProperPair;
 206
 207         // store insert size for first mate
 208         if ( m_settings->IsShowingInsertSizeSummary && al.IsFirstMate() && (al.InsertSize != 0) ) {
 209             int insertSize = abs(al.InsertSize);
 210             m_insertSizes.push_back( insertSize );
 211         }
 212     }
 213 }
 214
 215 bool StatsTool::StatsToolPrivate::Run() {
 216
 217     // set to default input if none provided
 218     if ( !m_settings->HasInput )
 219         m_settings->InputFiles.push_back(Options::StandardIn());
 220
 221     // open the BAM files
 222     BamMultiReader reader;
 223     if ( !reader.Open(m_settings->InputFiles) ) {
 224         cerr << "bamtools stats ERROR: could not open input BAM file(s)... Aborting." << endl;
 225         reader.Close();
 226         return false;
 227     }
 228
 229     // plow through alignments, keeping track of stats
 230     BamAlignment al;
 231     while ( reader.GetNextAlignmentCore(al) )
 232         ProcessAlignment(al);
 233     reader.Close();
 234
 235     // print stats & exit
 236     PrintStats();
 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> ...] [statsOptions]");
 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
 261     delete m_settings;
 262     m_settings = 0;
 263
 264     delete m_impl;
 265     m_impl = 0;
 266 }
 267
 268 int StatsTool::Help(void) {
 269     Options::DisplayHelp();
 270     return 0;
 271 }
 272
 273 int StatsTool::Run(int argc, char* argv[]) {
 274
 275     // parse command line arguments
 276     Options::Parse(argc, argv, 1);
 277
 278     // initialize StatsTool with settings
 279     m_impl = new StatsToolPrivate(m_settings);
 280
 281     // run StatsTool, return success/fail
 282     if ( m_impl->Run() )
 283         return 0;
 284     else
 285         return 1;
 286 }