+731
-535
lines changedFilter options
+731
-535
lines changed Original file line number Diff line number Diff line change
@@ -80,7 +80,8 @@ sub utgcns ($$) {
80 80
print F " -T ../$asm.ctgStore 1 \\\n";
81 81
print F " -P \$jobid \\\n";
82 82
print F " -O ./ctgcns/\$jobid.cns.WORKING \\\n";
83 -
print F " -B ./ctgcns/\$jobid.bams/ \\\n" if (getGlobal("bamOutput") eq "1");
83 +
print F " -norealign \\\n" if ((getGlobal("cnsAlign") eq "0") &&
84 +
(getGlobal("bamOutput") eq "0"));
84 85
print F " -maxcoverage " . getGlobal('cnsMaxCoverage') . " \\\n";
85 86
print F " -e " . getGlobal("cnsErrorRate") . " \\\n";
86 87
print F " -em " . getGlobal("cnsErrorRate") . " \\\n";
@@ -477,39 +478,37 @@ sub consensusLoad ($) {
477 478
$cmd .= "> ./5-consensus/ctgcns.files.ctgStoreLoad.err 2>&1";
478 479 479 480
if (runCommand("unitigging", $cmd)) {
480 -
caExit("failed to load unitig consensus into ctgStore", "$path/ctgcns.files.ctgStoreLoad.err");
481 +
caExit("failed to load unitig consensus into ctgStore", "$path/ctgcns.files.ctgStoreLoad.err");
481 482
}
482 483
unlink "$path/ctgcns.files.ctgStoreLoad.err";
483 484 484 485
stashFile("unitigging/$asm.ctgStore/seqDB.v002.dat");
485 486
stashFile("unitigging/$asm.ctgStore/seqDB.v002.tig");
486 487
}
487 488 488 -
# Consolidate BAM outputs.
489 +
# And extract bam records.
489 490 490 -
if (getGlobal("bamOutput") eq "1") {
491 -
my $samtools = getGlobal("samtools");
491 +
if ((getGlobal("bamOutput") eq "1") &&
492 +
(! fileExists("$asm.contigs.bam"))) {
493 +
my $samtools = getGlobal("samtools");
492 494 493 -
open(F, "find $path/ctgcns -type f -print |");
494 -
open(O, "> $path/ctgcns.bam.list") or caExit("failed to open '$path/ctgcns.bam.list' for writing: $!\n", undef);
495 -
while (<F>) {
496 -
if (m!ctgcns/\d+.bams/\d+.bam!) {
497 -
print O $_;
495 +
$cmd = "$bin/tgTigDisplay \\\n";
496 +
$cmd .= " -S ../$asm.seqStore \\\n";
497 +
$cmd .= " -b \\\n";
498 +
$cmd .= " -o ../$asm.contigs.unsorted.bam \\\n";
499 +
$cmd .= " -L ./5-consensus/ctgcns.files \\\n";
500 +
$cmd .= "> ./5-consensus/ctgcns.files.contigs.bam.err 2>&1";
501 + 502 +
if (runCommand("unitigging", $cmd)) {
503 +
caExit("failed to extract BAM records from ctgcns files", "$path/ctgcns.files.contigs.bam.err");
504 +
}
505 +
unlink "$path/ctgcns.files.contigs.bam.err";
506 + 507 +
if (runCommand(".", "$samtools sort --write-index --threads 2 -o ./$asm.contigs.bam ./$asm.contigs.unsorted.bam > ./ctgcns.bam.err 2>&1")) {
508 +
caExit("failed to sort BAM outputs from ./$asm.contigs.unsorted.bam to ./$asm.contigs.bam", "./ctgcns.bam.err");
498 509
}
499 -
}
500 -
close(O);
501 -
close(F);
502 - 503 -
if (runCommand(".", "$samtools merge -o ./$asm.contigs.unsorted.bam -b $path/ctgcns.bam.list > $path/ctgcns.bam.err 2>&1")) {
504 -
caExit("failed to merge BAM outputs into ./$asm.contigs.unsorted.bam", "$path/ctgcns.bam.err");
505 -
}
506 -
unlink "$path/ctgcns.bam.err";
507 - 508 -
if (runCommand(".", "$samtools sort --write-index --threads 2 -o ./$asm.contigs.bam ./$asm.contigs.unsorted.bam > $path/ctgcns.bam.err 2>&1")) {
509 -
caExit("failed to sort BAM outputs from ./$asm.contigs.unsorted.bam to ./$asm.contigs.bam", "$path/ctgcns.bam.err");
510 -
}
511 -
unlink "./$asm.contigs.unsorted.bam";
512 -
unlink "$path/ctgcns.bam.err";
510 +
unlink "./$asm.contigs.unsorted.bam";
511 +
unlink "./ctgcns.bam.err";
513 512
}
514 513 515 514
# Remove consensus outputs
@@ -1006,7 +1006,8 @@ sub setDefaults () {
1006 1006
setDefault("cnsConsensus", "pbdagcon", "Which consensus algorithm to use; 'pbdagcon' (fast, reliable); 'utgcns' (multialignment output); 'quick' (single read mosaic); default 'pbdagcon'");
1007 1007
setDefault("cnsPartitions", 0, "Attempt to create this many consensus jobs; default '0' = based on the largest tig");
1008 1008 1009 -
setDefault("bamOutput", 1, "Output contigs as BAM");
1009 +
setDefault("cnsAlign", 1, "Realign reads to contigs after consensus");
1010 +
setDefault("bamOutput", 1, "Output contigs as BAM (enables cnsAlign)");
1010 1011
#etDefault("fastaOutput", 1, "Output contigs as FASTA");
1011 1012 1012 1013
##### Correction Options
@@ -19,7 +19,7 @@
19 19
#include "tgStore.H"
20 20 21 21
uint32 MASRmagic = 0x5253414d; // 'MASR', as a big endian integer
22 -
uint32 MASRversion = 1;
22 +
uint32 MASRversion = 2;
23 23 24 24
#define MAX_VERS 1024 // Linked to 10 bits in the header file.
25 25
@@ -125,14 +125,6 @@ tgStore::tgStore(const char *path_,
125 125
assert(0);
126 126
break;
127 127
}
128 - 129 - 130 -
// Fail (again?) if there are no tigs loaded.
131 - 132 -
//if (_tigLen == 0) {
133 -
// fprintf(stderr, "tgStore::tgStore()-- ERROR, didn't find any tigs in the store. Correct version?\n");
134 -
// exit(1);
135 -
//}
136 128
}
137 129 138 130
@@ -643,18 +635,29 @@ tgStore::loadMASR(tgStoreEntry* &R, uint32& L, uint32& M, uint32 V) {
643 635
exit(1);
644 636
}
645 637 646 -
if (MASRversionInFile != MASRversion) {
647 -
fprintf(stderr, "tgStore::loadMASR()-- Failed to open '%s': version number mismatch; file=%d code=%d\n",
648 -
_name, MASRversionInFile, MASRversion);
649 -
exit(1);
650 -
}
651 - 652 -
// Check we're consistent.
653 638
if (L < MASRtotalInFile)
654 639
fprintf(stderr, "tgStore::loadMASR()-- '%s' has more tigs (" F_U32 ") than expected (" F_U32 ").\n",
655 640
_name, MASRtotalInFile, L), exit(1);
656 641 657 -
loadFromFile(R, "MASR", masrLen, F);
642 +
if (MASRversionInFile == 1) {
643 +
tgStoreEntry *Rv1 = new tgStoreEntry [M];
644 +
loadFromFile(Rv1, "MASR", masrLen, F);
645 + 646 +
for (uint32 ii=0; ii<masrLen; ii++)
647 +
R[ii] = Rv1[ii]; // A format conversion!
648 + 649 +
delete [] Rv1;
650 +
}
651 + 652 +
else if (MASRversionInFile == 2) {
653 +
loadFromFile(R, "MASR", masrLen, F);
654 +
}
655 + 656 +
else {
657 +
fprintf(stderr, "tgStore::loadMASR()-- Failed to open '%s': version number mismatch; file=%d code=%d\n",
658 +
_name, MASRversionInFile, MASRversion);
659 +
exit(1);
660 +
}
658 661 659 662
merylutil::closeFile(F, _name);
660 663
}
@@ -98,13 +98,34 @@ public:
98 98
};
99 99 100 100
private:
101 +
struct tgStoreEntryV1 {
102 +
tgTigRecordV1 tigRecord;
103 +
uint64 unusedFlags : 12; // One whole bit for future use.
104 +
uint64 flushNeeded : 1; // If true, this MAR and associated tig are NOT saved to disk.
105 +
uint64 isDeleted : 1; // If true, this MAR has been deleted from the assembly.
106 +
uint64 svID : 10; // 10 -> 1024 versions (HARDCODED in tgStore.C)
107 +
uint64 fileOffset : 40; // 40 -> 1 TB file size; offset in file where MA is stored
108 +
};
101 109
struct tgStoreEntry {
102 -
tgTigRecord tigRecord;
103 -
uint64 unusedFlags : 12; // One whole bit for future use.
104 -
uint64 flushNeeded : 1; // If true, this MAR and associated tig are NOT saved to disk.
105 -
uint64 isDeleted : 1; // If true, this MAR has been deleted from the assembly.
106 -
uint64 svID : 10; // 10 -> 1024 versions (HARDCODED in tgStore.C)
107 -
uint64 fileOffset : 40; // 40 -> 1 TB file size; offset in file where MA is stored
110 +
tgTigRecord tigRecord;
111 +
uint64 unusedFlags : 12; // One whole bit for future use.
112 +
uint64 flushNeeded : 1; // If true, this MAR and associated tig are NOT saved to disk.
113 +
uint64 isDeleted : 1; // If true, this MAR has been deleted from the assembly.
114 +
uint64 svID : 10; // 10 -> 1024 versions (HARDCODED in tgStore.C)
115 +
uint64 fileOffset : 40; // 40 -> 1 TB file size; offset in file where MA is stored
116 + 117 +
tgStoreEntry &
118 +
operator=(const tgStoreEntryV1 &that) {
119 +
tigRecord = that.tigRecord; // Itself an overloaded operator=().
120 + 121 +
unusedFlags = that.unusedFlags;
122 +
flushNeeded = that.flushNeeded;
123 +
isDeleted = that.isDeleted;
124 +
svID = that.svID;
125 +
fileOffset = that.fileOffset;
126 + 127 +
return *this;
128 +
}
108 129
};
109 130 110 131
void writeTigToDisk(tgTig *ma, tgStoreEntry *maRecord);
You can’t perform that action at this time.
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