Calculate a k-mer coverage threshold from the given k-mer coverage * histogram. */
| 27 | /** Calculate a k-mer coverage threshold from the given k-mer coverage |
| 28 | * histogram. */ |
| 29 | static inline |
| 30 | float calculateCoverageThreshold(const Histogram& h) |
| 31 | { |
| 32 | float cov = h.firstLocalMinimum(); |
| 33 | if (opt::rank <= 0) { |
| 34 | if (cov == 0) |
| 35 | std::cout << "Unable to determine minimum k-mer coverage\n"; |
| 36 | else |
| 37 | std::cout << "Minimum k-mer coverage is " << cov << std::endl; |
| 38 | } |
| 39 | |
| 40 | for (unsigned iteration = 0; iteration < 100; iteration++) { |
| 41 | Histogram trimmed = h.trimLow((unsigned)roundf(cov)); |
| 42 | if (opt::rank <= 0) |
| 43 | logger(1) << "Coverage: " << cov << "\t" |
| 44 | "Reconstruction: " << trimmed.size() << std::endl; |
| 45 | |
| 46 | unsigned median = trimmed.median(); |
| 47 | float cov1 = sqrt(median); |
| 48 | if (cov1 == cov) { |
| 49 | // The coverage threshold has converged. |
| 50 | if (opt::rank <= 0) |
| 51 | std::cout << "Using a coverage threshold of " |
| 52 | << (unsigned)roundf(cov) << "...\n" |
| 53 | "The median k-mer coverage is " << median << "\n" |
| 54 | "The reconstruction is " << trimmed.size() |
| 55 | << std::endl; |
| 56 | if (!opt::db.empty()) { |
| 57 | addToDb("coverageThreshold", (unsigned)roundf(cov)); |
| 58 | addToDb("medianKcoverage", median); |
| 59 | addToDb("restruction", trimmed.size()); |
| 60 | } |
| 61 | return cov; |
| 62 | } |
| 63 | cov = cov1; |
| 64 | } |
| 65 | if (opt::rank <= 0) |
| 66 | std::cerr << "warning: coverage threshold did not converge" |
| 67 | << std::endl; |
| 68 | return 0; |
| 69 | } |
| 70 | |
| 71 | /** Set the coverage-related parameters e and c from the given k-mer |
| 72 | * coverage histogram. */ |
no test coverage detected