@@ -48,8 +48,8 @@ public:
48 48 49 49
numOlaps = 0;
50 50 51 -
origLength = UINT32_MAX;
52 -
corrLength = UINT32_MAX; // Makes zeroth read first when sorted by length.
51 +
origLength = 0;
52 +
corrLength = 0;
53 53 54 54
memoryRequired = 0;
55 55
@@ -162,7 +162,6 @@ analyzeLength(tgTig *layout,
162 162
uint32 minOutputCoverage,
163 163
readStatus *status,
164 164
falconConsensus *fc) {
165 -
uint32 readID = layout->tigID();
166 165 167 166
// Estimate the length of the corrected read, using overlap position and depth.
168 167
@@ -197,16 +196,21 @@ analyzeLength(tgTig *layout,
197 196 198 197
// Save our results.
199 198 200 -
status[readID].readID = readID;
199 +
uint32 rid = layout->tigID();
200 + 201 +
status[rid].readID = rid;
201 202 202 -
status[readID].numOlaps = layout->numberOfChildren();
203 +
status[rid].numOlaps = layout->numberOfChildren();
203 204 204 -
status[readID].origLength = layout->length();
205 -
status[readID].corrLength = corLen;
205 +
status[rid].origLength = layout->length();
206 +
status[rid].corrLength = corLen;
206 207 207 -
status[readID].memoryRequired = fc->estimateMemoryUsage(layout->numberOfChildren(),
208 -
basesInOlaps,
209 -
layout->length());
208 +
status[rid].memoryRequired = fc->estimateMemoryUsage(layout->numberOfChildren(),
209 +
basesInOlaps,
210 +
layout->length());
211 + 212 +
//fprintf(stderr, "analyze for readID %u - olaps %u lengths %u %u\n",
213 +
// rid, status[rid].numOlaps, status[rid].origLength, status[rid].corrLength);
210 214
}
211 215 212 216
@@ -511,14 +515,22 @@ main(int argc, char **argv) {
511 515
FILE *stats = AS_UTL_openOutputFile(outName, '.', "stats");
512 516
FILE *log = AS_UTL_openOutputFile(outName, '.', "log");
513 517 518 +
// Initialize. Used to be done in analuzeLength(), but we skip it on
519 +
// empty tigs, and the re-sort below absolutely needs every readStatus
520 +
// element with a valid read id.
521 + 522 +
for (uint32 rr=0; rr<numReads+1; rr++)
523 +
status[rr].readID = rr;
524 + 514 525
// Scan the tigs, computing expected corrected length.
515 526 516 527
for (uint32 ti=1; ti<corStore->numTigs(); ti++) {
517 528
tgTig *layout = corStore->loadTig(ti);
518 529 519 -
analyzeLength(layout, minOutputLength, minOutputCoverage, status, fc);
530 +
if (layout)
531 +
analyzeLength(layout, minOutputLength, minOutputCoverage, status, fc);
520 532 521 -
corStore->unloadTig(layout->tigID());
533 +
corStore->unloadTig(ti);
522 534
}
523 535 524 536
// Sort by expected corrected length, then mark reads for correction until we get the desired
@@ -529,13 +541,13 @@ main(int argc, char **argv) {
529 541
uint64 desiredLength = genomeSize * outCoverage;
530 542
uint64 corrLength = 0;
531 543 532 -
for (uint32 ii=1; ii<numReads+1; ii++) {
533 -
if (status[ii].corrLength == 0)
544 +
for (uint32 rr=1; rr<numReads+1; rr++) {
545 +
if (status[rr].corrLength == 0)
534 546
continue;
535 547 536 -
status[ii].usedForCorrection = (corrLength < desiredLength);
548 +
status[rr].usedForCorrection = (corrLength < desiredLength);
537 549 538 -
corrLength += status[ii].corrLength;
550 +
corrLength += status[rr].corrLength;
539 551
}
540 552 541 553
// Sort by ID.
@@ -547,26 +559,27 @@ main(int argc, char **argv) {
547 559
for (uint32 ti=1; ti<corStore->numTigs(); ti++) {
548 560
tgTig *layout = corStore->loadTig(ti);
549 561 550 -
markEvidence(layout, status);
562 +
if (layout)
563 +
markEvidence(layout, status);
551 564 552 -
corStore->unloadTig(layout->tigID());
565 +
corStore->unloadTig(ti);
553 566
}
554 567 555 568
// And finally, flag any read for correction if it isn't already used as evidence or being corrected.
556 569 557 -
for (uint32 ti=1; ti<numReads+1; ti++) {
558 -
if ((status[ti].usedForCorrection == false) &&
559 -
(status[ti].usedForEvidence == false) &&
560 -
(status[ti].corrLength > minOutputLength))
561 -
status[ti].rescued = true;
570 +
for (uint32 rr=1; rr<numReads+1; rr++) {
571 +
if ((status[rr].usedForCorrection == false) &&
572 +
(status[rr].usedForEvidence == false) &&
573 +
(status[rr].corrLength > minOutputLength))
574 +
status[rr].rescued = true;
562 575
}
563 576 564 577
// Output the list of reads to correct.
565 578 566 -
for (uint32 ti=1; ti<numReads+1; ti++)
567 -
if ((status[ti].usedForCorrection == true) ||
568 -
(status[ti].rescued == true))
569 -
fprintf(roc, "%u\n", status[ti].readID);
579 +
for (uint32 rr=1; rr<numReads+1; rr++)
580 +
if ((status[rr].usedForCorrection == true) ||
581 +
(status[rr].rescued == true))
582 +
fprintf(roc, "%u\n", status[rr].readID);
570 583 571 584
AS_UTL_closeFile(roc);
572 585 RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4