@@ -228,6 +228,8 @@ main_clean_hairpins(int argc, const char **argv) {
228228 }
229229 }
230230
231+ vector<double > hist (20 , 0.0 );
232+
231233 if (!check_first || hairpin_fraction > max_hairpin_rate) {
232234
233235 // Input: paired-end reads with end1 and end2
@@ -265,6 +267,9 @@ main_clean_hairpins(int argc, const char **argv) {
265267 const double min_len = min (end_one.seq .length (), end_two.seq .length ());
266268 const double percent_match = sim/min_len;
267269
270+ const int hist_bin = hist.size ()*((sim - 0.001 )/(min_len + 0.001 ));
271+ hist[hist_bin]++;
272+
268273 if (percent_match > cutoff) {
269274
270275 sum_percent_match_bad += percent_match;
@@ -287,6 +292,13 @@ main_clean_hairpins(int argc, const char **argv) {
287292 << sum_percent_match_good/(n_reads - n_bad_reads) << endl
288293 << " mean_percent_match_hairpin: "
289294 << sum_percent_match_bad/n_bad_reads << endl;
295+
296+ const auto total = accumulate (cbegin (hist), cend (hist), 0.0 );
297+ transform (cbegin (hist), cend (hist), begin (hist),
298+ [&total](const double t) {return t/total;});
299+ for (auto i = 0 ; i < std::size (hist); ++i)
300+ cout << i << ' \t ' << std::setprecision (3 ) << hist[i] << endl;
301+ cout << endl;
290302 }
291303 }
292304 catch (const runtime_error &e) {
0 commit comments