Skip to content

Commit

Permalink
Replace ovStore::numOverlapsPerFrag() with ovStore::numOverlapsPerRea…
Browse files Browse the repository at this point in the history
…d(). The old

returned the number of overlaps per read in a range, the new returns for all reads.
  • Loading branch information
brianwalenz committed Sep 8, 2017
1 parent e88ad47 commit 1baf7ad
Show file tree
Hide file tree
Showing 5 changed files with 50 additions and 30 deletions.
5 changes: 2 additions & 3 deletions src/bogart/AS_BAT_OverlapCache.C
Original file line number Diff line number Diff line change
Expand Up @@ -231,8 +231,7 @@ OverlapCache::computeOverlapLimit(ovStore *ovlStore, uint64 genomeSize) {

uint32 frstRead = 0;
uint32 lastRead = 0;
uint32 *numPer = ovlStore->numOverlapsPerFrag(frstRead, lastRead);
uint32 totlRead = lastRead - frstRead + 1;
uint32 *numPer = ovlStore->numOverlapsPerRead(RI->numReads());

// Set the minimum number of overlaps per read to twice coverage. Then set the maximum number of
// overlaps per read to a guess of what it will take to fill up memory.
Expand Down Expand Up @@ -268,7 +267,7 @@ OverlapCache::computeOverlapLimit(ovStore *ovlStore, uint64 genomeSize) {
numEqual = 0;
numAbove = 0;

for (uint32 i=0; i<totlRead; i++) {
for (uint32 i=1; i<=RI->numReads(); i++) {
if (numPer[i] < _maxPer) {
numBelow += 1;
olapLoad += numPer[i];
Expand Down
4 changes: 2 additions & 2 deletions src/erateEstimate/erateEstimate.C
Original file line number Diff line number Diff line change
Expand Up @@ -651,12 +651,12 @@ main(int argc, char **argv) {
uint64 *overlapIndex = new uint64 [numIIDs + 1];
uint32 bgn = 0;
uint32 end = 0;
uint32 *overlapLen = ovlStore->numOverlapsPerFrag(bgn, end);
uint32 *overlapLen = ovlStore->numOverlapsPerRead();

overlapIndex[0] = 0;

for (uint32 iid=0; iid<numIIDs; iid++)
overlapIndex[iid+1] = overlapIndex[iid] + overlapLen[iid];
overlapIndex[iid+1] = overlapIndex[iid] + overlapLen[iid+1];

assert(overlapIndex[numIIDs] == numOvls);

Expand Down
51 changes: 31 additions & 20 deletions src/stores/ovStore.C
Original file line number Diff line number Diff line change
Expand Up @@ -487,42 +487,53 @@ ovStore::numOverlapsInRange(void) {



uint32 *
ovStore::numOverlapsPerFrag(uint32 &firstFrag, uint32 &lastFrag) {

if (_firstIIDrequested > _lastIIDrequested)
return(NULL);

firstFrag = _firstIIDrequested;
lastFrag = _lastIIDrequested;
// Fills olapsPerRead with the number of overlaps per read.
// olapsPerRead must be at least 'lastRead+1' in size, not just lastRead-firstRead+1.
// olapsPerRead is not initialized; must be cleared before calling.
//
void
ovStore::numOverlapsPerRead(uint32 *olapsPerRead,
uint32 firstRead,
uint32 lastRead) {

off_t originalPosition = AS_UTL_ftell(_offtFile);

AS_UTL_fseek(_offtFile, (off_t)_firstIIDrequested * sizeof(ovStoreOfft), SEEK_SET);

// Even if we're doing a whole human-size store, this allocation is
// (a) temporary and (b) only 512MB. The only current consumer of
// this code is FragCorrectOVL.c, which doesn't run on the whole
// human, it runs on ~24 pieces, which cuts this down to < 32MB.
AS_UTL_fseek(_offtFile, (off_t)firstRead * sizeof(ovStoreOfft), SEEK_SET);

uint64 len = _lastIIDrequested - _firstIIDrequested + 1;
uint64 len = lastRead - firstRead + 1;

ovStoreOfft *offsets = new ovStoreOfft [len];
uint32 *numolap = new uint32 [len];

uint64 act = AS_UTL_safeRead(_offtFile, offsets, "ovStore::numOverlapsInRange::offsets", sizeof(ovStoreOfft), len);
uint64 act = AS_UTL_safeRead(_offtFile, offsets, "ovStore::numOverlapsPerRead", sizeof(ovStoreOfft), len);

if (len != act)
fprintf(stderr, "AS_OVS_numOverlapsPerFrag()-- short read on offsets! Expected len=" F_U64 " read act=" F_U64 "\n", len, act), exit(1);
fprintf(stderr, "ovStore::numOverlapsPerRead()-- short read on offsets! Expected len=" F_U64 " read act=" F_U64 "\n", len, act), exit(1);

for (uint64 i=0; i<len; i++)
numolap[i] = offsets[i]._numOlaps;
olapsPerRead[firstRead+i] = offsets[i]._numOlaps;

delete [] offsets;

AS_UTL_fseek(_offtFile, originalPosition, SEEK_SET);
}

return(numolap);


uint32 *
ovStore::numOverlapsPerRead(uint32 numReads) {

if ((numReads == 0) && (_gkp != NULL))
numReads = _gkp->gkStore_getNumReads();

assert(numReads > 0);

uint32 *olapsPerRead = new uint32 [numReads+1];

memset(olapsPerRead, 0, (numReads+1) * sizeof(uint32));

numOverlapsPerRead(olapsPerRead, 1, numReads);

return(olapsPerRead);
}


Expand Down
6 changes: 5 additions & 1 deletion src/stores/ovStore.H
Original file line number Diff line number Diff line change
Expand Up @@ -316,7 +316,11 @@ public:
void resetRange(void);

uint64 numOverlapsInRange(void);
uint32 * numOverlapsPerFrag(uint32 &firstFrag, uint32 &lastFrag);

void numOverlapsPerRead(uint32 *olapsPerRead,
uint32 firstRead,
uint32 lastRead);
uint32 *numOverlapsPerRead(uint32 numReads=0);

// Add new evalues for reads between bgnID and endID. No checking of IDs is done, but the number
// of evalues must agree.
Expand Down
14 changes: 10 additions & 4 deletions src/stores/ovStoreDump.C
Original file line number Diff line number Diff line change
Expand Up @@ -254,8 +254,8 @@ dumpStore(ovStore *ovlStore,
uint32 obtDumped = 0;
uint32 merDumped = 0;

uint32 *counts = NULL;
ovStoreHistogram *hist = NULL;
uint32 *counts = NULL;
ovStoreHistogram *hist = NULL;

// Set the range of the reads to dump early so that we can reset it later.

Expand All @@ -271,11 +271,17 @@ dumpStore(ovStore *ovlStore,
}

// If we're dumping counts, and no modifiers, we can just ask the store for the counts
// and set the range to null.
// and set the range to null. Argh, the rest of the code expects counts[] to start at
// bgnID, so we need to rewrite everything.

if ((asCounts) && (dumpType == 0)) {
counts = ovlStore->numOverlapsPerFrag(bgnID, endID);
counts = new uint32 [endID - bgnID + 1];

ovlStore->numOverlapsPerRead(counts, bgnID, endID);
ovlStore->setRange(1, 0);

for (uint32 ii=bgnID; ii<=endID; ii++)
counts[ii - bgnID] = counts[ii];
}

// If we're dumping the erate-vs-length histogram, and no modifiers, grab it from the store and
Expand Down

0 comments on commit 1baf7ad

Please sign in to comment.