diff --git a/esl_gencode.c b/esl_gencode.c index d0c6205..e869e03 100644 --- a/esl_gencode.c +++ b/esl_gencode.c @@ -726,7 +726,7 @@ esl_gencode_WorkstateDestroy(ESL_GENCODE_WORKSTATE *wrk) if(wrk->orf_block != NULL) { - esl_sq_DestroyBlock(wrk->orf_block); + esl_gencode_OrfBlockDestroy(wrk->orf_block); wrk->orf_block = NULL; } @@ -779,11 +779,135 @@ esl_gencode_WorkstateCreate(ESL_GETOPTS *go, ESL_GENCODE *gcode) * 6. Functions for processing ORFs *****************************************************************/ +/* Residue storage of an ORF block comes in chunks of at least this many bytes. + * Kept small: when a block is created and freed for every query, 1 MB chunks + * raised peak memory about 30%, and 512 kB chunks sometimes did; 64 kB to + * 256 kB did not. */ +#define eslORF_CHUNK (1 << 16) + +/* Function: esl_gencode_OrfBlockCreate() + * Synopsis: Create an empty block of ORFs. + * + * Purpose: Create a block with room for ORFs; it grows as needed. + * + * Returns: the new block; on allocation failure. + */ +ESL_ORF_BLOCK * +esl_gencode_OrfBlockCreate(int count) +{ + ESL_ORF_BLOCK *block = NULL; + int status; + + if (count < 1) count = 1; + ESL_ALLOC(block, sizeof(ESL_ORF_BLOCK)); + block->list = NULL; + block->chunk = NULL; + block->chunkalloc = NULL; + block->count = 0; + block->listSize = count; + block->nchunks = 0; + block->cur = 0; + block->used = 0; + + ESL_ALLOC(block->list, sizeof(ESL_ORF) * count); + ESL_ALLOC(block->chunk, sizeof(ESL_DSQ *)); + ESL_ALLOC(block->chunkalloc, sizeof(int64_t)); + block->chunk[0] = NULL; + ESL_ALLOC(block->chunk[0], sizeof(ESL_DSQ) * eslORF_CHUNK); + block->chunkalloc[0] = eslORF_CHUNK; + block->nchunks = 1; + return block; + + ERROR: + esl_gencode_OrfBlockDestroy(block); + return NULL; +} + +/* Function: esl_gencode_OrfBlockReuse() + * Synopsis: Empty a block of ORFs, keeping its memory. + */ +void +esl_gencode_OrfBlockReuse(ESL_ORF_BLOCK *block) +{ + block->count = 0; + block->cur = 0; + block->used = 0; +} + +/* Function: esl_gencode_OrfBlockDestroy() + * Synopsis: Free a block of ORFs. + */ +void +esl_gencode_OrfBlockDestroy(ESL_ORF_BLOCK *block) +{ + int c; + + if (block == NULL) return; + if (block->chunk != NULL) + for (c = 0; c < block->nchunks; c++) free(block->chunk[c]); + free(block->chunk); + free(block->chunkalloc); + free(block->list); + free(block); +} + +/* orf_block_add() + * Add the finished ORF to . + */ +static int +orf_block_add(ESL_ORF_BLOCK *block, const ESL_SQ *psq) +{ + ESL_ORF *orf; + int64_t nbytes = psq->n + 2; /* residues and the two sentinels */ + int status; + + if (block->count == block->listSize) + { + ESL_REALLOC(block->list, sizeof(ESL_ORF) * block->listSize * 2); + block->listSize *= 2; + } + + if (block->used + nbytes > block->chunkalloc[block->cur]) + { /* move to the next chunk. Chunks after hold no ORFs, so one of them can be created or resized */ + int next = (block->used > 0 ? block->cur + 1 : block->cur); + if (next == block->nchunks) + { + ESL_REALLOC(block->chunk, sizeof(ESL_DSQ *) * (block->nchunks + 1)); + ESL_REALLOC(block->chunkalloc, sizeof(int64_t) * (block->nchunks + 1)); + block->chunk[next] = NULL; + block->chunkalloc[next] = 0; + block->nchunks++; + } + if (block->chunkalloc[next] < nbytes) + { + int64_t newalloc = ESL_MAX(nbytes, eslORF_CHUNK); + ESL_REALLOC(block->chunk[next], sizeof(ESL_DSQ) * newalloc); + block->chunkalloc[next] = newalloc; + } + block->cur = next; + block->used = 0; + } + + orf = block->list + block->count; + orf->dsq = block->chunk[block->cur] + block->used; + memcpy(orf->dsq, psq->dsq, sizeof(ESL_DSQ) * nbytes); + orf->n = psq->n; + orf->start = psq->start; + orf->end = psq->end; + block->used += nbytes; + block->count++; + return eslOK; + + ERROR: + return status; +} + int esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq) { int status = eslOK; + int named = FALSE; /* TRUE if was given a name and description for output */ ESL_SQ *psq = wrk->psq[wrk->frame]; psq->end = (wrk->is_revcomp ? wrk->apos+1 : wrk->apos-1); if (wrk->in_orf[wrk->frame] && psq->n >= wrk->minlen) @@ -796,6 +920,7 @@ esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq) if (wrk->orf_block == NULL) { /* esl-translate: name and description go to FASTA output */ esl_sq_FormatName(psq, "orf%d", wrk->orfcount); esl_sq_FormatDesc(psq, "source=%s coords=%" PRId64 "..%" PRId64 " length=%" PRId64 " frame=%d desc=%s", psq->source, psq->start, psq->end, psq->n, wrk->frame + 1 + (wrk->is_revcomp ? 3 : 0), sq->desc); + named = TRUE; } /* if we do not have a block to write ORFs to then write ORFs to file */ if (wrk->orf_block == NULL) @@ -804,20 +929,26 @@ esl_gencode_ProcessOrf(ESL_GENCODE_WORKSTATE *wrk, ESL_SQ *sq) } else { - if (wrk->orf_block->count == wrk->orf_block->listSize) - { - status = esl_sq_BlockGrowTo(wrk->orf_block, wrk->orf_block->listSize + 128, TRUE, psq->abc); - if (status != eslOK) ESL_XEXCEPTION(eslEMEM, "Cannot increase size of ORF sequence block"); - } - - esl_sq_Copy(psq, &(wrk->orf_block->list[wrk->orf_block->count])); - - wrk->orf_block->count++; + if ((status = orf_block_add(wrk->orf_block, psq)) != eslOK) ESL_XEXCEPTION(eslEMEM, "Cannot add to the ORF block"); } } - esl_sq_Reuse(psq); - esl_sq_SetSource(psq, sq->name); + /* Ready for the next ORF. Unless it was named for output above, only + * its residues and coordinates have changed since it was last reset, and + * its source is still the one esl_gencode_ProcessStart() set. + */ + if (named) + { + esl_sq_Reuse(psq); + esl_sq_SetSource(psq, sq->name); + } + else + { + psq->dsq[1] = eslDSQ_SENTINEL; + psq->n = 0; + psq->start = 0; + psq->end = 0; + } wrk->in_orf[wrk->frame] = FALSE; ERROR: diff --git a/esl_gencode.h b/esl_gencode.h index b27051d..cb08be0 100644 --- a/esl_gencode.h +++ b/esl_gencode.h @@ -22,6 +22,35 @@ typedef struct { const ESL_ALPHABET *aa_abc; // A reference to amino alphabet that caller is maintaining } ESL_GENCODE; +/* A translated ORF in an : its residues, and where it lies on + * the nucleotide sequence it came from. + */ +typedef struct { + ESL_DSQ *dsq; // residues 1..n, with sentinel bytes at 0 and n+1; points into the block's residue storage + int64_t n; // length in residues + int64_t start; // nucleotide coord of the ORF's first base + int64_t end; // nucleotide coord of its last base; < start on the reverse strand +} ESL_ORF; + +/* The ORFs translated from one or more nucleotide sequences. Their residues + * are kept end to end in a few large chunks, so an ORF costs its residues + * and one . + */ +typedef struct { + ESL_ORF *list; // ORFs 0..count-1 + int count; // number of ORFs + int listSize; // number of ORFs is allocated for + ESL_DSQ **chunk; // residue storage, chunks 0..nchunks-1; a chunk doesn't move while it holds ORFs + int64_t *chunkalloc; // allocated size of each chunk, in bytes + int nchunks; // number of chunks + int cur; // chunk being filled + int64_t used; // bytes used in chunk +} ESL_ORF_BLOCK; + +extern ESL_ORF_BLOCK *esl_gencode_OrfBlockCreate (int count); +extern void esl_gencode_OrfBlockReuse (ESL_ORF_BLOCK *block); +extern void esl_gencode_OrfBlockDestroy(ESL_ORF_BLOCK *block); + /* struct esl_gencode_workstate_s * keeps state in DNA sequence , allowing us to process a sequence * either in a single gulp (using ReadSeq) or in overlapping windows @@ -42,7 +71,7 @@ typedef struct esl_gencode_workstate_s { int is_revcomp; // TRUE|FALSE: TRUE if we're doing reverse complement strand int orfcount; // >=0: How many ORFs we've processed so far - ESL_SQ_BLOCK *orf_block; // sequences block to hold ORFs + ESL_ORF_BLOCK *orf_block; // block to hold ORFs /* one-time configuration information (from options) */ int do_watson; // TRUE|FALSE: TRUE if we translate the top strand