Skip to content

Commit

Permalink
Write output overlaps to a file (not a store); ensure alignment and o…
Browse files Browse the repository at this point in the history 
…verlap outputs are sync'd.
  • Loading branch information
@brianwalenz
brianwalenz committed Jul 15, 2019
1 parent 262edd9 commit 27bda6e
Show file tree
Hide file tree
Showing 5 changed files with 60 additions and 69 deletions.
16 changes: 8 additions & 8 deletions src/overlapAlign/computeAlignments.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -512,15 +512,15 @@ maComputation::computeAlignments(uint32 minOverlapLength,
// Process each overlap.
//

for (uint32 ii=0; ii<_overlapsLen; ii++) {
ovOverlap *ov = _overlaps + ii;
for (uint32 ovlid=0; ovlid<_overlapsLen; ovlid++) {
ovOverlap *ov = _overlaps + ovlid;

if (ov->evalue() == AS_MAX_EVALUE) // Skip overlaps that are flagged as trimmed out
continue; // or otherwise not useful.

if (_verboseAlign > 0) {
fprintf(stderr, "\n");
fprintf(stderr, "computeAlignments()-- ----------------------------------------OVERLAP %u/%u\n", ii, _overlapsLen);
fprintf(stderr, "computeAlignments()-- ----------------------------------------OVERLAP %u/%u\n", ovlid, _overlapsLen);
}

assert(_aID == ov->a_iid);
Expand All @@ -545,7 +545,7 @@ maComputation::computeAlignments(uint32 minOverlapLength,
if ((aend - abgn < _minLength) ||
(bend - bbgn < _minLength)) {
if (_verboseAlign > 0)
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d short.\n", ii, _aID, abgn, aend, _bID, bbgn, bend);
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d short.\n", ovlid, _aID, abgn, aend, _bID, bbgn, bend);
ov->evalue(AS_MAX_EVALUE);
continue;
}
Expand All @@ -560,7 +560,7 @@ maComputation::computeAlignments(uint32 minOverlapLength,
if ((a5 - b5 > _maxEdge) &&
(a3 - b3 > _maxEdge)) {
if (_verboseAlign > 0)
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d contained.\n", ii, _aID, abgn, aend, _bID, bbgn, bend);
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d contained.\n", ovlid, _aID, abgn, aend, _bID, bbgn, bend);
ov->evalue(AS_MAX_EVALUE);
continue;
}
Expand All @@ -574,7 +574,7 @@ maComputation::computeAlignments(uint32 minOverlapLength,
(a3 < b3) &&
(_readData[_aID].trimmedLength - a5 + b5 < thick3)) {
if (_verboseAlign > 0)
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d thick3 %d < %d\n", ii, _aID, abgn, aend, _bID, bbgn, bend, _readData[_aID].trimmedLength - a5 + b5, thick3);
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d thick3 %d < %d\n", ovlid, _aID, abgn, aend, _bID, bbgn, bend, _readData[_aID].trimmedLength - a5 + b5, thick3);
ov->evalue(AS_MAX_EVALUE);
continue;
}
Expand All @@ -588,7 +588,7 @@ maComputation::computeAlignments(uint32 minOverlapLength,
(a3 > b3) &&
(_readData[_aID].trimmedLength - a3 + b3 < thick5)) {
if (_verboseAlign > 0)
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d thick5 %d < %d\n", ii, _aID, abgn, aend, _bID, bbgn, bend, _readData[_aID].trimmedLength - a3 + b3, thick5);
fprintf(stderr, "computeAlignments()-- overlap %4u a=%8u %6d-%-6d b=%8u %6d-%-6d thick5 %d < %d\n", ovlid, _aID, abgn, aend, _bID, bbgn, bend, _readData[_aID].trimmedLength - a3 + b3, thick5);
ov->evalue(AS_MAX_EVALUE);
continue;
}
Expand All @@ -604,7 +604,7 @@ maComputation::computeAlignments(uint32 minOverlapLength,
// Load the b read and compute the alignment.
fetchTrimmedRead(_bID, false, ov->flipped());

computeOverlapAlignment(ov,
computeOverlapAlignment(ovlid,
minOverlapLength,
maxErate,
_overlapSlop,
Expand Down
36 changes: 17 additions & 19 deletions src/overlapAlign/computeOverlapAlignment.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,12 +194,13 @@ computeAlignment(char *aRead, int32 abgn, int32 aend, int32 UNUSED(alen
// which we use to update this end of the overlap.
//
void
maComputation::computeOverlapAlignment(ovOverlap *ovl,
maComputation::computeOverlapAlignment(uint32 ovlid,
uint32 minOverlapLength,
double maxErate,
uint32 overlapSlop,
uint32 maxRepeat) {
ovOverlap ori = *ovl;
ovOverlap *ovl = &_overlaps[ovlid]; // Convenience pointer to the overlap.
ovOverlap ori = _overlaps[ovlid]; // Copy of the original overlap.

int32 alen = _readData[_aID].trimmedLength;
int32 abgn = (int32) ovl->dat.ovl.ahg5;
Expand Down Expand Up @@ -442,17 +443,16 @@ maComputation::computeOverlapAlignment(ovOverlap *ovl,

ovl->erate(editDist / (double)alignLen);

_alignsOvl[_alignsLen] = ovl; // Allocate space for the alignment output.
_alignsA [_alignsLen] = new char [alen + 1];
_alignsB [_alignsLen] = new char [alen + 1];
_alignsA [ovlid] = new char [alen + 1]; // Allocate space for the alignment output.
_alignsB [ovlid] = new char [alen + 1];

for (uint32 ii=0; ii<alen; ii++) {
_alignsA[_alignsLen][ii] = '-';
_alignsB[_alignsLen][ii] = '-';
_alignsA[ovlid][ii] = '-';
_alignsB[ovlid][ii] = '-';
}

_alignsA[_alignsLen][alen] = 0;
_alignsB[_alignsLen][alen] = 0;
_alignsA[ovlid][alen] = 0;
_alignsB[ovlid][alen] = 0;

char *aaln = new char [result.alignmentLength + 1]; // Convert the alignment to a string.
char *baln = new char [result.alignmentLength + 1];
Expand Down Expand Up @@ -483,23 +483,23 @@ maComputation::computeOverlapAlignment(ovOverlap *ovl,
uint32 aa = 0; // Position in sequence A.

for (uint32 ii=0; ii<ovl->dat.ovl.ahg5; ii++, pp++) {
_alignsA[_alignsLen][pp] = _aRead[aa];
_alignsB[_alignsLen][pp] = '-';
_alignsA[ovlid][pp] = _aRead[aa];
_alignsB[ovlid][pp] = '-';
aa++;
}

for (uint32 ii=0; ii<result.alignmentLength; ii++) {
if (aaln[ii] != '-') {
_alignsA[_alignsLen][pp] = aaln[ii];
_alignsB[_alignsLen][pp] = baln[ii];
_alignsA[ovlid][pp] = aaln[ii];
_alignsB[ovlid][pp] = baln[ii];
aa++;
pp++;
}
}

for (uint32 ii=0; ii<ovl->dat.ovl.ahg3; ii++, pp++) {
_alignsA[_alignsLen][pp] = _aRead[aa];
_alignsB[_alignsLen][pp] = '-';
_alignsA[ovlid][pp] = _aRead[aa];
_alignsB[ovlid][pp] = '-';
aa++;
}

Expand All @@ -509,10 +509,8 @@ maComputation::computeOverlapAlignment(ovOverlap *ovl,
delete [] aaln;
delete [] baln;

_alignsA[_alignsLen][alen] = 0;
_alignsB[_alignsLen][alen] = 0;

_alignsLen++;
_alignsA[ovlid][alen] = 0;
_alignsB[ovlid][alen] = 0;
}

// Trash the alignment results.
Expand Down
16 changes: 7 additions & 9 deletions src/overlapAlign/overlapAlign-computation.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,24 +82,25 @@ public:

// Allocate space for the outputs.

_alignsLen = 0;
_alignsMax = _overlapsLen + 1;
_alignsOvl = new ovOverlap * [_overlapsLen];
_alignsA = new char * [_overlapsLen];
_alignsB = new char * [_overlapsLen];

for (uint32 ii=0; ii<_overlapsLen; ii++) {
_alignsA[ii] = NULL;
_alignsB[ii] = NULL;
}
};

~maComputation() {
delete [] _overlaps;
delete [] _aRead;
delete [] _bRead;

for (uint32 ii=0; ii<_alignsLen; ii++) {
for (uint32 ii=0; ii<_overlapsLen; ii++) {
delete [] _alignsA[ii];
delete [] _alignsB[ii];
}

delete [] _alignsOvl;
delete [] _alignsA;
delete [] _alignsB;
};
Expand All @@ -114,7 +115,7 @@ private:
bool trimOverlap_Flipped(ovOverlap *ovl);
bool trimOverlap(ovOverlap *ovl);

void computeOverlapAlignment(ovOverlap *ovl,
void computeOverlapAlignment(uint32 ovlid,
uint32 minOverlapLength,
double maxErate,
uint32 overlapSlop,
Expand Down Expand Up @@ -168,9 +169,6 @@ public:

// Alignment results.

uint32 _alignsLen;
uint32 _alignsMax;
ovOverlap **_alignsOvl;
char **_alignsA;
char **_alignsB;
};
16 changes: 8 additions & 8 deletions src/overlapAlign/overlapAlign-globalData.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,8 @@ public:
ovlStoreName = NULL;
ovlStore = NULL;

outStoreName = NULL;
outStore = NULL;
outFileName = NULL;
outFile = NULL;

clearRangesFileName = NULL;
};
Expand Down Expand Up @@ -113,10 +113,10 @@ public:
ovlStore = new ovStore(ovlStoreName, seqStore);
}

if (outStoreName) {
fprintf(stderr, "Writing overlaps to store '%s'.\n", outStoreName);
if (outFileName) {
fprintf(stderr, "Writing overlaps to file '%s'.\n", outFileName);

outStore = new ovStoreWriter(outStoreName, seqStore);
outFile = new ovFile(seqStore, outFileName, ovFileFullWrite);
}

// Check parameters.
Expand Down Expand Up @@ -153,7 +153,7 @@ public:
delete seqCache;
delete ovlStore;

delete outStore;
delete outFile;
};

// Parameters
Expand Down Expand Up @@ -192,8 +192,8 @@ public:

// Outputs

char *outStoreName;
ovStoreWriter *outStore;
char *outFileName;
ovFile *outFile;

char *clearRangesFileName;
};
Expand Down
45 changes: 20 additions & 25 deletions src/overlapAlign/overlapAlign.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,30 +73,25 @@ overlapWriter(void *G, void *S) {
trGlobalData *g = (trGlobalData *)G;
maComputation *s = (maComputation *)S;

for (uint64 oo=0; oo<s->_overlapsLen; oo++) {
ovOverlap *ovl = s->_overlaps + oo;
// Write the header for our big-pile-o-alignments output.

fprintf(stdout, "\n");
fprintf(stdout, "%6u %6d %6d %s\n", s->_aID, 0, s->_readData[s->_aID].trimmedLength, s->_aRead);

// For each overlap, write the overlap and the alignment string (unless
// it's a garbage overlap).

if (ovl->evalue() == AS_MAX_EVALUE)
for (uint64 oo=0; oo<s->_overlapsLen; oo++) {
if (s->_overlaps[oo].evalue() == AS_MAX_EVALUE)
continue;

g->outStore->writeOverlap(ovl);
}
g->outFile->writeOverlap(&s->_overlaps[oo]);

if (s->_alignsLen > 0) {
fprintf(stdout, "\n");
fprintf(stdout, "%6u %6d %6d %s\n", s->_aID, 0, s->_readData[s->_aID].trimmedLength, s->_aRead);

for (uint32 oo=0; oo<s->_alignsLen; oo++) {
//fprintf(stdout, "\n");
//fprintf(stdout, "%6u %6d %6d\n", s->_alignsOvl[oo]->b_iid, (int32)s->_alignsOvl[oo]->dat.ovl.ahg5, (int32)s->_alignsOvl[oo]->dat.ovl.ahg3);
//fprintf(stdout, "%s\n", s->_alignsA[oo]);
//fprintf(stdout, "%s\n", s->_alignsB[oo]);
fprintf(stdout, "%6u %6d %6d %s\n",
s->_alignsOvl[oo]->b_iid,
(int32)s->_alignsOvl[oo]->dat.ovl.ahg5,
(int32)s->_alignsOvl[oo]->dat.ovl.ahg3,
s->_alignsB[oo]);
}
fprintf(stdout, "%6u %6d %6d %s\n",
s->_overlaps[oo].b_iid,
(int32)s->_overlaps[oo].dat.ovl.ahg5,
(int32)s->_overlaps[oo].dat.ovl.ahg3,
s->_alignsB[oo]);
}

delete s;
Expand Down Expand Up @@ -328,7 +323,7 @@ main(int argc, char **argv) {
}

else if (strcmp(argv[arg], "-align") == 0) {
g->outStoreName = argv[++arg];
g->outFileName = argv[++arg];
}


Expand Down Expand Up @@ -357,7 +352,7 @@ main(int argc, char **argv) {
if (g->seqStoreName == NULL) err.push_back("No sequence store (-S option) supplied.\n");
if (g->ovlStoreName == NULL) err.push_back("No overlap store (-O option) supplied.\n");

if ((g->outStoreName == NULL) &&
if ((g->outFileName == NULL) &&
(g->clearRangesFileName == NULL))
err.push_back("At least one of -trim and -align must be supplied.\n");
#endif
Expand Down Expand Up @@ -420,7 +415,7 @@ main(int argc, char **argv) {
}

else if ((g->clearRangesFileName != NULL) &&
(g->outStoreName == NULL)) {
(g->outFileName == NULL)) {
g->initialize();

fprintf(stderr, "COMPUTE TRIMMING and STOP.\n");
Expand All @@ -430,15 +425,15 @@ main(int argc, char **argv) {
}

else if ((g->clearRangesFileName == NULL) &&
(g->outStoreName != NULL)) {
(g->outFileName != NULL)) {
g->initialize();

fprintf(stderr, "ALIGNING OVERLAPS.\n");
alignOverlaps(g, false);
}

else if ((g->clearRangesFileName != NULL) &&
(g->outStoreName != NULL)) {
(g->outFileName != NULL)) {
g->initialize(sqStore_extend);

fprintf(stderr, "COMPUTE TRIMMING.\n");
Expand Down

0 comments on commit 27bda6e

Please sign in to comment.