Estimate the distance between two contigs. * @param numPairs [out] the number of pairs that agree with the * expected distribution * @return the estimated distance */
| 217 | * @return the estimated distance |
| 218 | */ |
| 219 | static int estimateDistance(unsigned len0, unsigned len1, |
| 220 | const Pairs& pairs, const PMF& pmf, |
| 221 | unsigned& numPairs) |
| 222 | { |
| 223 | // The provisional fragment sizes are calculated as if the contigs |
| 224 | // were perfectly adjacent with no overlap or gap. |
| 225 | typedef vector<pair<int, int> > Fragments; |
| 226 | Fragments fragments; |
| 227 | fragments.reserve(pairs.size()); |
| 228 | for (Pairs::const_iterator it = pairs.begin(); |
| 229 | it != pairs.end(); ++it) { |
| 230 | int a0 = it->targetAtQueryStart(); |
| 231 | int a1 = it->mateTargetAtQueryStart(); |
| 232 | if (it->isReverse()) |
| 233 | a0 = len0 - a0; |
| 234 | if (!it->isMateReverse()) |
| 235 | a1 = len1 - a1; |
| 236 | fragments.push_back(opt::rf |
| 237 | ? make_pair(a1, len1 + a0) |
| 238 | : make_pair(a0, len0 + a1)); |
| 239 | } |
| 240 | |
| 241 | // Remove duplicate fragments. |
| 242 | unsigned orig = fragments.size(); |
| 243 | sort(fragments.begin(), fragments.end()); |
| 244 | fragments.erase(unique(fragments.begin(), fragments.end()), |
| 245 | fragments.end()); |
| 246 | numPairs = fragments.size(); |
| 247 | assert((int)orig - (int)numPairs >= 0); |
| 248 | stats.total_frags += orig; |
| 249 | stats.dup_frags += orig - numPairs; |
| 250 | |
| 251 | if (numPairs < opt::npairs) |
| 252 | return INT_MIN; |
| 253 | |
| 254 | vector<int> fragmentSizes; |
| 255 | fragmentSizes.reserve(fragments.size()); |
| 256 | unsigned ma = opt::minAlign; |
| 257 | for (Fragments::const_iterator it = fragments.begin(); |
| 258 | it != fragments.end(); ++it) { |
| 259 | int x = it->second - it->first; |
| 260 | if (!opt::rf && opt::method == MLE |
| 261 | && x <= 2 * int(ma - 1)) { |
| 262 | unsigned align = x / 2; |
| 263 | if (opt::verbose > 0) |
| 264 | #pragma omp critical(cerr) |
| 265 | cerr << PROGRAM ": warning: The observed fragment of " |
| 266 | "size " << x << " bp is shorter than 2*l " |
| 267 | "(l=" << opt::minAlign << ").\n"; |
| 268 | ma = min(ma, align); |
| 269 | } |
| 270 | fragmentSizes.push_back(x); |
| 271 | } |
| 272 | |
| 273 | #pragma omp critical(g_recMA) |
| 274 | g_recMA = min(g_recMA, ma); |
| 275 | switch (opt::method) { |
| 276 | case MLE: |
no test coverage detected