Forms contigs based on the consensus of each base and outputs them * to the file specified by the -o option. */
| 366 | /** Forms contigs based on the consensus of each base and outputs them |
| 367 | * to the file specified by the -o option. */ |
| 368 | static void consensus(const string& outPath, const string& pileupPath) |
| 369 | { |
| 370 | ofstream outFile(outPath.c_str()); |
| 371 | assert_good(outFile, outPath); |
| 372 | |
| 373 | ofstream pileupFile; |
| 374 | ostream& pileupOut |
| 375 | = pileupPath.empty() || pileupPath == "-" ? cout |
| 376 | : (pileupFile.open(pileupPath.c_str()), pileupFile); |
| 377 | assert_good(pileupOut, pileupPath); |
| 378 | |
| 379 | unsigned numIgnored = 0; |
| 380 | for (ContigMap::const_iterator it = g_contigs.begin(); |
| 381 | it != g_contigs.end(); ++it) { |
| 382 | const ContigCount& contig = it->second; |
| 383 | unsigned seqLength = it->second.counts.size(); |
| 384 | |
| 385 | Sequence outSeq(seqLength, 'N'); |
| 386 | unsigned sumBest = 0; |
| 387 | unsigned sumSecond = 0; |
| 388 | for (unsigned x = 0; x < seqLength; x++) { |
| 389 | char c = selectBase( |
| 390 | it->second.counts[x], sumBest, sumSecond); |
| 391 | outSeq[x] = islower(contig.seq[x]) ? tolower(c) : c; |
| 392 | } |
| 393 | |
| 394 | if (outSeq.find_first_of("ACGT") != string::npos) { |
| 395 | // Check that the average percent agreement was enough to |
| 396 | // write the contig to file. |
| 397 | float percentAgreement |
| 398 | = sumBest / (float)(sumBest + sumSecond); |
| 399 | if (isnan(percentAgreement) || percentAgreement < .9) { |
| 400 | numIgnored++; |
| 401 | if (opt::csToNt) { |
| 402 | if (opt::verbose > 0) |
| 403 | cerr << "warning: Contig " << it->first |
| 404 | << " has less than 90% agreement " |
| 405 | "and will not be converted.\n"; |
| 406 | } else |
| 407 | continue; |
| 408 | } else { |
| 409 | if (opt::csToNt) |
| 410 | fixUnknown(outSeq, contig.seq); |
| 411 | ostringstream comment; |
| 412 | comment << outSeq.length() << ' ' << contig.coverage; |
| 413 | if (!contig.comment.empty()) |
| 414 | comment << ' ' << contig.comment; |
| 415 | outFile << FastaRecord( |
| 416 | it->first, comment.str(), outSeq); |
| 417 | assert(outFile.good()); |
| 418 | } |
| 419 | |
| 420 | if (opt::verbose > 1) { |
| 421 | // ID pos reference genotype A C G T |
| 422 | if (opt::csToNt) |
| 423 | for (unsigned i = 0; i < seqLength - 1; i++) |
| 424 | cout << it->first << '\t' << 1+i |
| 425 | << '\t' << contig.seq[i] |
no test coverage detected