Resolve ambiguous region using pairwise alignment * (Needleman-Wunsch) ('solutions' contain exactly two paths, from a * source contig to a dest contig) */
| 434 | * source contig to a dest contig) |
| 435 | */ |
| 436 | static ContigPath alignPair(const Graph& g, |
| 437 | const ContigPaths& solutions, ofstream& out) |
| 438 | { |
| 439 | assert(solutions.size() == 2); |
| 440 | assert(solutions[0].size() > 1); |
| 441 | assert(solutions[1].size() > 1); |
| 442 | assert(solutions[0].front() == solutions[1].front()); |
| 443 | assert(solutions[0].back() == solutions[1].back()); |
| 444 | ContigPath fstSol(solutions[0].begin()+1, solutions[0].end()-1); |
| 445 | ContigPath sndSol(solutions[1].begin()+1, solutions[1].end()-1); |
| 446 | |
| 447 | if (fstSol.empty() || sndSol.empty()) { |
| 448 | // This entire sequence may be deleted. |
| 449 | const ContigPath& sol(fstSol.empty() ? sndSol : fstSol); |
| 450 | assert(!sol.empty()); |
| 451 | Sequence consensus(mergePath(g, sol)); |
| 452 | assert(consensus.size() > opt::k - 1); |
| 453 | string::iterator first = consensus.begin() + opt::k - 1; |
| 454 | transform(first, consensus.end(), first, ::tolower); |
| 455 | |
| 456 | unsigned match = opt::k - 1; |
| 457 | float identity = (float)match / consensus.size(); |
| 458 | if (opt::verbose > 2) |
| 459 | cerr << consensus << '\n'; |
| 460 | if (opt::verbose > 1) |
| 461 | cerr << identity |
| 462 | << (identity < opt::identity ? " (too low)\n" : "\n"); |
| 463 | if (identity < opt::identity) |
| 464 | return ContigPath(); |
| 465 | |
| 466 | unsigned coverage = calculatePathProperties(g, sol).coverage; |
| 467 | ContigNode u = outputNewContig(g, |
| 468 | solutions, 1, 1, consensus, coverage, out); |
| 469 | ContigPath path; |
| 470 | path.push_back(solutions.front().front()); |
| 471 | path.push_back(u); |
| 472 | path.push_back(solutions.front().back()); |
| 473 | return path; |
| 474 | } |
| 475 | |
| 476 | Sequence fstPathContig(mergePath(g, fstSol)); |
| 477 | Sequence sndPathContig(mergePath(g, sndSol)); |
| 478 | if (fstPathContig == sndPathContig) { |
| 479 | // These two paths have identical sequence. |
| 480 | if (fstSol.size() == sndSol.size()) { |
| 481 | // A perfect match must be caused by palindrome. |
| 482 | typedef ContigPath::const_iterator It; |
| 483 | pair<It, It> it = mismatch( |
| 484 | fstSol.begin(), fstSol.end(), sndSol.begin()); |
| 485 | assert(it.first != fstSol.end()); |
| 486 | assert(it.second != sndSol.end()); |
| 487 | assert(*it.first |
| 488 | == get(vertex_complement, g, *it.second)); |
| 489 | assert(equal(it.first+1, It(fstSol.end()), it.second+1)); |
| 490 | if (opt::verbose > 1) |
| 491 | cerr << "Palindrome: " |
| 492 | << get(vertex_contig_name, g, *it.first) << '\n'; |
| 493 | return solutions[0]; |
no test coverage detected