1 /*
2  * aligner_metrics.h
3  */
5 #ifndef ALIGNER_METRICS_H_
6 #define ALIGNER_METRICS_H_
8 #include <math.h>
9 #include <iostream>
10 #include <seqan/sequence.h>
11 #include "alphabet.h"
12 #include "timer.h"
14 using namespace std;
16 /**
17  * Borrowed from http://www.johndcook.com/standard_deviation.html,
18  * which in turn is borrowed from Knuth.
19  */
20 class RunningStat {
21 public:
22         RunningStat() : m_n(0), m_tot(0.0) { }
24         void clear() {
25                 m_n = 0;
26                 m_tot = 0.0;
27         }
29         void push(float x) {
30                 m_n++;
31                 m_tot += x;
32                 // See Knuth TAOCP vol 2, 3rd edition, page 232
33                 if (m_n == 1) {
34                         m_oldM = m_newM = x;
35                         m_oldS = 0.0;
36                 } else {
37                         m_newM = m_oldM + (x - m_oldM)/m_n;
38                         m_newS = m_oldS + (x - m_oldM)*(x - m_newM);
39                         // set up for next iteration
40                         m_oldM = m_newM;
41                         m_oldS = m_newS;
42                 }
43         }
45         int num() const {
46                 return m_n;
47         }
49         double tot() const {
50                 return m_tot;
51         }
53         double mean() const {
54                 return (m_n > 0) ? m_newM : 0.0;
55         }
57         double variance() const {
58                 return ( (m_n > 1) ? m_newS/(m_n - 1) : 0.0 );
59         }
61         double stddev() const {
62                 return sqrt(variance());
63         }
65 private:
66         int m_n;
67         double m_tot;
68         double m_oldM, m_newM, m_oldS, m_newS;
69 };
71 /**
72  * Encapsulates a set of metrics that we would like an aligner to keep
73  * track of, so that we can possibly use it to diagnose performance
74  * issues.
75  */
76 class AlignerMetrics {
78 public:
80         AlignerMetrics() :
81                 curBacktracks_(0),
82                 curBwtOps_(0),
83                 first_(true),
84                 curIsLowEntropy_(false),
85                 curIsHomoPoly_(false),
87                 curNumNs_(0),
116                 timer_(cout, "", false)
117                 { }
119         void printSummary() {
120                 if(!first_) {
122                 }
123                 cout << "AlignerMetrics:" << endl;
124                 cout << "  # Reads:             " << reads_ << endl;
125                 float hopct = (reads_ > 0) ? (((float)homoReads_)/((float)reads_)) : (0.0);
126                 hopct *= 100.0;
127                 cout << "  % homo-polymeric:    " << (hopct) << endl;
128                 float lopct = (reads_ > 0) ? ((float)lowEntReads_/(float)(reads_)) : (0.0);
129                 lopct *= 100.0;
130                 cout << "  % low-entropy:       " << (lopct) << endl;
131                 float unpct = (reads_ > 0) ? ((float)unalignedReads_/(float)(reads_)) : (0.0);
132                 unpct *= 100.0;
133                 cout << "  % unaligned:         " << (unpct) << endl;
134                 float npct = (reads_ > 0) ? ((float)threeOrMoreNReads_/(float)(reads_)) : (0.0);
135                 npct *= 100.0;
136                 cout << "  % with 3 or more Ns: " << (npct) << endl;
137                 cout << endl;
138                 cout << "  Total BWT ops:    avg: " << bwtOpsPerRead_.mean() << ", stddev: " << bwtOpsPerRead_.stddev() << endl;
139                 cout << "  Total Backtracks: avg: " << backtracksPerRead_.mean() << ", stddev: " << backtracksPerRead_.stddev() << endl;
140                 time_t elapsed = timer_.elapsed();
141                 cout << "  BWT ops per second:    " << (bwtOpsPerRead_.tot()/elapsed) << endl;
142                 cout << "  Backtracks per second: " << (backtracksPerRead_.tot()/elapsed) << endl;
143                 cout << endl;
144                 cout << "  Homo-poly:" << endl;
145                 cout << "    BWT ops:    avg: " << bwtOpsPerHomoRead_.mean() << ", stddev: " << bwtOpsPerHomoRead_.stddev() << endl;
146                 cout << "    Backtracks: avg: " << backtracksPerHomoRead_.mean() << ", stddev: " << backtracksPerHomoRead_.stddev() << endl;
147                 cout << "  Low-entropy:" << endl;
148                 cout << "    BWT ops:    avg: " << bwtOpsPerLoEntRead_.mean() << ", stddev: " << bwtOpsPerLoEntRead_.stddev() << endl;
149                 cout << "    Backtracks: avg: " << backtracksPerLoEntRead_.mean() << ", stddev: " << backtracksPerLoEntRead_.stddev() << endl;
150                 cout << "  High-entropy:" << endl;
151                 cout << "    BWT ops:    avg: " << bwtOpsPerHiEntRead_.mean() << ", stddev: " << bwtOpsPerHiEntRead_.stddev() << endl;
152                 cout << "    Backtracks: avg: " << backtracksPerHiEntRead_.mean() << ", stddev: " << backtracksPerHiEntRead_.stddev() << endl;
153                 cout << endl;
154                 cout << "  Unaligned:" << endl;
155                 cout << "    BWT ops:    avg: " << bwtOpsPerUnalignedRead_.mean() << ", stddev: " << bwtOpsPerUnalignedRead_.stddev() << endl;
156                 cout << "    Backtracks: avg: " << backtracksPerUnalignedRead_.mean() << ", stddev: " << backtracksPerUnalignedRead_.stddev() << endl;
157                 cout << "  Aligned:" << endl;
158                 cout << "    BWT ops:    avg: " << bwtOpsPerAlignedRead_.mean() << ", stddev: " << bwtOpsPerAlignedRead_.stddev() << endl;
159                 cout << "    Backtracks: avg: " << backtracksPerAlignedRead_.mean() << ", stddev: " << backtracksPerAlignedRead_.stddev() << endl;
160                 cout << endl;
161                 cout << "  0 Ns:" << endl;
162                 cout << "    BWT ops:    avg: " << bwtOpsPer0nRead_.mean() << ", stddev: " << bwtOpsPer0nRead_.stddev() << endl;
163                 cout << "    Backtracks: avg: " << backtracksPer0nRead_.mean() << ", stddev: " << backtracksPer0nRead_.stddev() << endl;
164                 cout << "  1 N:" << endl;
165                 cout << "    BWT ops:    avg: " << bwtOpsPer1nRead_.mean() << ", stddev: " << bwtOpsPer1nRead_.stddev() << endl;
166                 cout << "    Backtracks: avg: " << backtracksPer1nRead_.mean() << ", stddev: " << backtracksPer1nRead_.stddev() << endl;
167                 cout << "  2 Ns:" << endl;
168                 cout << "    BWT ops:    avg: " << bwtOpsPer2nRead_.mean() << ", stddev: " << bwtOpsPer2nRead_.stddev() << endl;
169                 cout << "    Backtracks: avg: " << backtracksPer2nRead_.mean() << ", stddev: " << backtracksPer2nRead_.stddev() << endl;
170                 cout << "  >2 Ns:" << endl;
171                 cout << "    BWT ops:    avg: " << bwtOpsPer3orMoreNRead_.mean() << ", stddev: " << bwtOpsPer3orMoreNRead_.stddev() << endl;
172                 cout << "    Backtracks: avg: " << backtracksPer3orMoreNRead_.mean() << ", stddev: " << backtracksPer3orMoreNRead_.stddev() << endl;
173                 cout << endl;
174         }
176         /**
177          *
178          */
180                 if(!first_) {
182                 }
183                 first_ = false;
184                 float ent = entropyDna5(read);
185                 curIsLowEntropy_ = (ent < 0.75f);
186                 curIsHomoPoly_ = (ent < 0.001f);
187                 curHadRanges_ = false;
188                 curBwtOps_ = 0;
189                 curBacktracks_ = 0;
190                 // Count Ns
191                 curNumNs_ = 0;
192                 const size_t len = seqan::length(read);
193                 for(size_t i = 0; i < len; i++) {
194                         if((int)read[i] == 4) curNumNs_++;
195                 }
196         }
198         /**
199          *
200          */
201         void setReadHasRange() {
202                 curHadRanges_ = true;
203         }
205         /**
206          * Commit the running statistics for this read to
207          */
208         void finishRead() {
211                 else if(curIsLowEntropy_) lowEntReads_++;
217                 // Drill down by entropy
218                 if(curIsHomoPoly_) {
221                 } else if(curIsLowEntropy_) {
224                 } else {
227                 }
228                 // Drill down by whether it aligned
232                 } else {
235                 }
236                 if(curNumNs_ == 0) {
240                 } else if(curNumNs_ == 1) {
244                 } else if(curNumNs_ == 2) {
248                 } else {
252                 }
253         }
255         // Running-total of the number of backtracks and BWT ops for the
256         // current read
257         uint32_t curBacktracks_;
258         uint32_t curBwtOps_;
260 protected:
262         bool first_;
264         // true iff the current read is low entropy
265         bool curIsLowEntropy_;
266         // true if current read is all 1 char (or very close)
267         bool curIsHomoPoly_;
268         // true iff the current read has had one or more ranges reported
270         // number of Ns in current read
271         int curNumNs_;
273         // # reads
275         // # homo-poly reads
277         // # low-entropy reads
279         // # high-entropy reads
281         // # reads with alignments
283         // # reads without alignments
285         // # reads with 3 or more Ns
287         // # reads with < 3 Ns
290         // Distribution of BWT operations per read
294         // Distribution of BWT operations per homo-poly read
298         // Distribution of BWT operations per low-entropy read
302         // Distribution of BWT operations per high-entropy read
306         // Distribution of BWT operations per read that "aligned" (for
307         // which a range was arrived at - range may not have necessarily
308         // lead to an alignment)
312         // Distribution of BWT operations per read that didn't align
316         // Distribution of BWT operations/backtracks per read with no Ns