From ee7a0ebfc53a4a6d4dda86fa568c133e7195f8fc Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Thu, 1 Oct 2026 14:47:51 -0500 Subject: [PATCH 1/7] pipeline: keep each ORF's window index in the pipeline's workspace p7_pli_Frameshift() recorded which DNA window an ORF belongs to in orfsq->idx, the only thing the pipeline wrote into the ORF block it is given. Keep it in the per-call workspace instead (pli_tmp->orf_win[], starting at -1). The ORF block is then read-only to the pipeline, and an ORF needs no idx field. Co-Authored-By: Claude Opus 5.5 --- src/p7_pipeline.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/p7_pipeline.c b/src/p7_pipeline.c index 80bb2ab4..b4a86562 100644 --- a/src/p7_pipeline.c +++ b/src/p7_pipeline.c @@ -29,6 +29,7 @@ typedef struct { P7_OMX **oxf_holder; // - a temporary list of forward parser matricies for ORFs float *fwdsc; double *P_orf; // - a temporary list or forwrad P values for ORFs + int *orf_win; // - window index of each ORF; kept here so the ORF block stays read-only } P7_PIPELINE_OBJS; @@ -1473,9 +1474,6 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF /* Build windows from ORF's that pass F4 */ p7_pli_BuildDNAWindows(pli, orf_block, dnasq, om, bg, data, &fwd_windowlist, 0., pli_tmp, hit_windows, hit_windows_start, complementarity); - /* An ORF belongs to window w only once the loop below finds it inside w */ - for(i = 0; i < orf_block->count; i++) orf_block->list[i].idx = -1; - for(w = 0; w < fwd_windowlist.count; w++) { dna_window = &(fwd_windowlist.windows[w]); @@ -1514,7 +1512,7 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF /* Only process ORF if it in inside the current window */ if(orf_start >= window_start && orf_end <= window_end) { - orfsq->idx = w; + pli_tmp->orf_win[i] = w; P_min = ESL_MIN(P_min, pli_tmp->P_orf[i]); tot_orfsc = p7_FLogsum(tot_orfsc, pli_tmp->fwdsc[i]); orf_cnt++; @@ -1597,7 +1595,7 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF orfsq = &(orf_block->list[i]); - if(orfsq->idx != w) continue; // This ORF does not overlap with this window + if(pli_tmp->orf_win[i] != w) continue; // This ORF does not overlap with this window if(pli_tmp->P_orf[i] > pli->F3) continue; // This ORF did not pass Forward if(pli_tmp->oxf_holder[i] == NULL) continue; // This ORF has already been aligned @@ -1727,16 +1725,19 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI pli_tmp->tmpseq = NULL; pli_tmp->oxf_holder = NULL; pli_tmp->P_orf = NULL; + pli_tmp->orf_win = NULL; pli_tmp->fwdsc = NULL; ESL_ALLOC(pli_tmp->fwdsc, sizeof(float) * orf_block->count); ESL_ALLOC(pli_tmp->P_orf, sizeof(double) * orf_block->count); + ESL_ALLOC(pli_tmp->orf_win, sizeof(int) * orf_block->count); ESL_ALLOC(pli_tmp->oxf_holder, sizeof(P7_OMX *) * orf_block->count); for(i = 0; i < orf_block->count; i++) { pli_tmp->oxf_holder[i] = NULL; pli_tmp->fwdsc[i] = -eslINFINITY; pli_tmp->P_orf[i] = 1.0; + pli_tmp->orf_win[i] = -1; } pli_tmp->tmpseq = esl_sq_CreateDigital(dnasq->abc); @@ -1918,6 +1919,7 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI if (pli_tmp->tmpseq != NULL) esl_sq_Destroy(pli_tmp->tmpseq); if (pli_tmp->oxf_holder != NULL) free(pli_tmp->oxf_holder); if (pli_tmp->P_orf != NULL) free(pli_tmp->P_orf); + if (pli_tmp->orf_win != NULL) free(pli_tmp->orf_win); if (pli_tmp->fwdsc != NULL) free(pli_tmp->fwdsc); free(pli_tmp); } @@ -1930,6 +1932,7 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI if (pli_tmp->tmpseq != NULL) esl_sq_Destroy(pli_tmp->tmpseq); if (pli_tmp->oxf_holder != NULL) free(pli_tmp->oxf_holder); if (pli_tmp->P_orf != NULL) free(pli_tmp->P_orf); + if (pli_tmp->orf_win != NULL) free(pli_tmp->orf_win); if (pli_tmp->fwdsc != NULL) free(pli_tmp->fwdsc); free(pli_tmp); } From 5d80965e278139e741049d1fe18784f93d8697ef Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sat, 3 Oct 2026 21:28:32 -0500 Subject: [PATCH 2/7] pipeline: remove p7_pli_ExtendAndMergeWindows_BATH() Nothing calls it, and with the window index kept in the pipeline's workspace it was the last code that wrote an ORF's idx. Co-Authored-By: Claude Opus 5.5 --- src/p7_pipeline.c | 131 ---------------------------------------------- 1 file changed, 131 deletions(-) diff --git a/src/p7_pipeline.c b/src/p7_pipeline.c index b4a86562..6e77864b 100644 --- a/src/p7_pipeline.c +++ b/src/p7_pipeline.c @@ -293,137 +293,6 @@ p7_pipeline_Destroy_BATH(P7_PIPELINE *pli) * 2. The pipeline API. *****************************************************************/ -/* Function: p7_pli_ExtendAndMergeWindows_BATH - * Synopsis: Turns a list of ORF ssv diagonals into DNA windows, and - * merges overlapping windows. - * - * Purpose: Accepts a of and creates a set of SSV - * diagonals. For each ORFs best (higest scoreing) - * diagonal, extends those windows based on a combination - * of the max_length value from and the prefix and - * suffix lengths stored in , and converts them to - * DNA coordinates. Then merges (in place) windows that - * overlap by more than percent, ensuring - * that windows stay within the bounds of 1... - * - * Returns: - */ -int -p7_pli_ExtendAndMergeWindows_BATH(P7_PIPELINE *pli, ESL_SQ_BLOCK *orf_block, ESL_SQ *dnasq, P7_OPROFILE *om, P7_BG *bg, const P7_SCOREDATA *data, P7_HMM_WINDOWLIST *windowlist, float pct_overlap, int complementarity) -{ - - int i, f; - int best_window; - float best_score; - P7_HMM_WINDOWLIST tmp_windowlist; - P7_HMM_WINDOW *prev_window = NULL; - P7_HMM_WINDOW *curr_window = NULL; - int64_t window_start; - int64_t window_end; - int32_t window_len; - int64_t overlap_start; - int64_t overlap_end; - int32_t overlap_len; - int new_hit_cnt = 0; - ESL_SQ *curr_orf; - - p7_hmmwindow_init(&tmp_windowlist); - - for(f = 0; f < orf_block->count; f++) - { - curr_orf = &(orf_block->list[f]); - - p7_oprofile_ReconfigLength(om, curr_orf->n); - p7_omx_GrowTo(pli->oxf, om->M, 0, curr_orf->n); /* expand the one-row omx if needed */ - - p7_SSVFilter_BATH(curr_orf->dsq, curr_orf->n, om, pli->oxf, data, bg, pli->F1, &tmp_windowlist); - - /* If the orf fails to produce a window use the full ORF coords aligned to the center of the model. */ - if(tmp_windowlist.count == 0) { - if(curr_orf->n >= om->M) { - window_start = (curr_orf->n - om->M) / 2 + 1; - window_end = om->M; - window_len = om->M; - } - else { - window_start = 1; - window_end = om->M - ((om->M - curr_orf->n) / 2); - window_len = curr_orf->n; - } - p7_hmmwindow_new(&tmp_windowlist, 0, window_start, window_end, window_len, 0.0, 0, curr_orf->n); - } - - best_window = 0; - best_score = tmp_windowlist.windows[0].score; - for(i = 1; i < tmp_windowlist.count; i++) { - if(tmp_windowlist.windows[i].score > best_score) { - best_window = i; - best_score = tmp_windowlist.windows[i].score; - } - } - - curr_window = tmp_windowlist.windows+best_window; - - /* Extend ORF coords */ - window_start = curr_window->n - (om->max_length * (0.1 + data->prefix_lengths[curr_window->k - curr_window->length + 1])) + 1; - window_end = curr_window->n + curr_window->length + (om->max_length * (0.1 + data->suffix_lengths[curr_window->k])) - 2; - - window_start = ESL_MIN(0, window_start); //move start to at least the begining of the ORF - window_end = ESL_MAX(curr_orf->n, window_end); // move end to at least the end of the ORF - - /* Convert to DNA coords */ - if(complementarity) { - window_start = ESL_MAX(1, (dnasq->n - curr_orf->start + 1) + (window_start * 3)); - window_end = ESL_MIN(dnasq->n, (dnasq->n - curr_orf->start + 1) + (window_end * 3)); - } - else { - window_start = ESL_MAX(1, curr_orf->start + (window_start * 3)); - window_end = ESL_MIN(dnasq->n, curr_orf->start + (window_end * 3)); - } - - p7_hmmwindow_new(windowlist, 0, window_start, curr_window->k, window_end-window_start+1, 0.0, complementarity, dnasq->n); - curr_orf->idx = windowlist->count - 1; // keep track of which window ORFs belong to - tmp_windowlist.count = 0; - } - - if (tmp_windowlist.windows != NULL) free (tmp_windowlist.windows); - - if( windowlist->count == 0) return eslOK; - - p7_hmmwindow_SortByStart(windowlist); - new_hit_cnt = 0; - - /* merge overlapping windows, compressing list in place. */ - for (i=1; icount; i++) { - prev_window = windowlist->windows+new_hit_cnt; - curr_window = windowlist->windows+i; - - overlap_start = ESL_MAX(prev_window->n, curr_window->n); - overlap_end = ESL_MIN(prev_window->n+prev_window->length-1, curr_window->n+curr_window->length-1); - overlap_len = overlap_end - overlap_start + 1; - - window_start = ESL_MIN(prev_window->n, curr_window->n); - window_end = ESL_MAX(prev_window->n+prev_window->length-1, curr_window->n+curr_window->length-1); - window_len = window_end - window_start + 1; - - - if(((float)(overlap_len)/ESL_MIN(prev_window->length, curr_window->length) > pct_overlap) && - window_len < ( 2 * (om->max_length * 3))) - { - prev_window->n = window_start; - prev_window->length = window_len; - } - else { - new_hit_cnt++; - windowlist->windows[new_hit_cnt] = windowlist->windows[i]; - } - orf_block->list[i].idx = new_hit_cnt; - } - windowlist->count = new_hit_cnt+1; - - return eslOK; -} - static int p7_pli_ComputeLocalCompo(const P7_SCOREDATA *data, const P7_OPROFILE *om, const P7_BG *bg, int k_start, int k_end, float *compo) { From a9c7c9604573c967c8586df4a58c22e3f85930f9 Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sat, 3 Oct 2026 21:28:56 -0500 Subject: [PATCH 3/7] pipeline: take translated ORFs as a compact ESL_ORF_BLOCK The pipeline reads only the residues, length and coordinates of an ORF, so it now takes easel's ESL_ORF_BLOCK in place of a block of ESL_SQ objects, about 75 bytes per ORF where an ESL_SQ took about 1 kB. bathsearch creates and empties its ORF blocks with the matching easel calls, and the seed-window scoring in p7_hmmwindow.c and the pipeline unit test pass their translated sequence as an ESL_ORF. This needs the easel commit that adds ESL_ORF_BLOCK. Output is unchanged. For one profile against a 352 Mb genome with 16 threads, peak memory falls from 133 MB to 66 MB and CPU time by 17%. Co-Authored-By: Claude Opus 5.5 --- src/bathsearch.c | 12 ++++++------ src/hmmer.h | 8 ++++---- src/p7_alidisplay.c | 2 +- src/p7_domaindef.c | 6 +++--- src/p7_hmmwindow.c | 9 ++++++++- src/p7_pipeline.c | 23 ++++++++++++++--------- 6 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/bathsearch.c b/src/bathsearch.c index 21e92485..f905a185 100644 --- a/src/bathsearch.c +++ b/src/bathsearch.c @@ -830,7 +830,7 @@ serial_master(ESL_GETOPTS *go, struct cfg_s *cfg) info[i].hw = p7_hmmwindow_CreateList(); info[i].gcode = gcode; info[i].wrk = esl_gencode_WorkstateCreate(go, gcode); - info[i].wrk->orf_block = esl_sq_CreateDigitalBlock(BLOCK_SIZE, abcAA); + info[i].wrk->orf_block = esl_gencode_OrfBlockCreate(BLOCK_SIZE); info[i].th = p7_tophits_Create(); info[i].om = p7_oprofile_Clone(om); info[i].gm = p7_profile_Clone(gm); @@ -923,7 +923,7 @@ serial_master(ESL_GETOPTS *go, struct cfg_s *cfg) p7_hmmwindow_DestroyList(info[i].hw); if(info[i].wrk->orf_block != NULL) { - esl_sq_DestroyBlock(info[i].wrk->orf_block); + esl_gencode_OrfBlockDestroy(info[i].wrk->orf_block); info[i].wrk->orf_block = NULL; esl_gencode_WorkstateDestroy(info[i].wrk); } @@ -1098,7 +1098,7 @@ serial_loop(WORKER_INFO *info, ID_LENGTH_LIST *id_length_list, ESL_SQFILE *dbfp, p7_Pipeline_BATH(info->pli, info->om, info->gm, info->om_fs3, info->om_fs5, info->gm_fs5, info->scoredata, info->bg, info->th, info->pli->nseqs, dbsq_dna, info->wrk->orf_block, info->gcode, info->hw, p7_NOCOMPLEMENT); p7_pipeline_Reuse_BATH(info->pli); // prepare for next search - esl_sq_ReuseBlock(info->wrk->orf_block); + esl_gencode_OrfBlockReuse(info->wrk->orf_block); } if (info->pli->strands != p7_STRAND_TOPONLY) @@ -1112,7 +1112,7 @@ serial_loop(WORKER_INFO *info, ID_LENGTH_LIST *id_length_list, ESL_SQFILE *dbfp, p7_Pipeline_BATH(info->pli, info->om, info->gm, info->om_fs3, info->om_fs5, info->gm_fs5, info->scoredata, info->bg, info->th, info->pli->nseqs, dbsq_dna, info->wrk->orf_block, info->gcode, info->hw, p7_COMPLEMENT); p7_pipeline_Reuse_BATH(info->pli); // prepare for next search - esl_sq_ReuseBlock(info->wrk->orf_block); + esl_gencode_OrfBlockReuse(info->wrk->orf_block); /* Reverse sequence back to original */ esl_sq_ReverseComplement(dbsq_dna); @@ -1283,7 +1283,7 @@ pipeline_thread(void *arg) p7_Pipeline_BATH(info->pli, info->om, info->gm, info->om_fs3, info->om_fs5, info->gm_fs5, info->scoredata, info->bg, info->th, block->first_seqidx + i, dnaSeq, info->wrk->orf_block, info->gcode, info->hw, p7_NOCOMPLEMENT); p7_pipeline_Reuse_BATH(info->pli); // prepare for next search - esl_sq_ReuseBlock(info->wrk->orf_block); + esl_gencode_OrfBlockReuse(info->wrk->orf_block); } if (info->pli->strands != p7_STRAND_TOPONLY) { @@ -1294,7 +1294,7 @@ pipeline_thread(void *arg) p7_Pipeline_BATH(info->pli, info->om, info->gm, info->om_fs3, info->om_fs5, info->gm_fs5, info->scoredata, info->bg, info->th, block->first_seqidx + i, dnaSeq, info->wrk->orf_block, info->gcode, info->hw, p7_COMPLEMENT); p7_pipeline_Reuse_BATH(info->pli); // prepare for next search - esl_sq_ReuseBlock(info->wrk->orf_block); + esl_gencode_OrfBlockReuse(info->wrk->orf_block); esl_sq_ReverseComplement(dnaSeq); } } diff --git a/src/hmmer.h b/src/hmmer.h index 09c2e24e..792d6479 100644 --- a/src/hmmer.h +++ b/src/hmmer.h @@ -1393,7 +1393,7 @@ extern int p7_tracealign_getMSAandStats(P7_HMM *hmm, ESL_SQ **sq, int N, ESL_MS /* p7_alidisplay.c */ extern P7_ALIDISPLAY *p7_alidisplay_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *sq, const ESL_SQ *ntsq); extern P7_ALIDISPLAY *p7_alidisplay_fs_Create(const P7_TRACE *tr, int which, const P7_FS_PROFILE *gm_fs5, const ESL_SQ *sq, int show_cigar); -extern P7_ALIDISPLAY *p7_alidisplay_nonfs_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *sq, const ESL_SQ *orfsq, int orf_pos, int show_cigar); +extern P7_ALIDISPLAY *p7_alidisplay_nonfs_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *sq, const ESL_ORF *orfsq, int orf_pos, int show_cigar); extern P7_ALIDISPLAY *p7_alidisplay_splice_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *target_seq, const ESL_SQ *amino_sq, int amino_pos, int splice_cnt, int show_cigar); extern P7_ALIDISPLAY *p7_alidisplay_splice_fs_Create(const P7_TRACE *tr, int which, const P7_FS_PROFILE *gm_fs, const ESL_SQ *sq, ESL_DSQ *nuc_dsq, int64_t *nuc_index, int nuc_pos, int splice_cnt, int show_cigar); extern P7_ALIDISPLAY *p7_alidisplay_Create_empty(); @@ -1459,7 +1459,7 @@ extern int p7_domaindef_DumpPosteriors(FILE *ofp, P7_DOMAINDEF *ddef); extern void p7_domaindef_Destroy_BATH(P7_DOMAINDEF *ddef); extern int p7_domaindef_ByPosteriorHeuristics_Frameshift_BATH(P7_PIPELINE *pli, ESL_SQ *windowsq, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_BG *bg, const ESL_GENCODE *gcode); -extern int p7_domaindef_ByPosteriorHeuristics_BATH(const ESL_SQ *orfsq, const ESL_SQ *sq, const int64_t ntsqlen, P7_OPROFILE *om, P7_PROFILE *gm, P7_OMX *oxf, P7_OMX *oxb, P7_OMX *fwd, P7_OMX *bck, P7_DOMAINDEF *ddef); +extern int p7_domaindef_ByPosteriorHeuristics_BATH(const ESL_ORF *orfsq, const ESL_SQ *sq, const int64_t ntsqlen, P7_OPROFILE *om, P7_PROFILE *gm, P7_OMX *oxf, P7_OMX *oxb, P7_OMX *fwd, P7_OMX *bck, P7_DOMAINDEF *ddef); /* p7_gmx.c */ extern P7_GMX *p7_gmx_Create (int allocM, int allocL, int allocLx, int nscells); @@ -1576,9 +1576,9 @@ extern int p7_pli_DomainIncludable (P7_PIPELINE *pli, float dom_score, double l extern int p7_pli_NewModel (P7_PIPELINE *pli, const P7_OPROFILE *om, P7_BG *bg); extern int p7_pli_NewModelThresholds(P7_PIPELINE *pli, const P7_OPROFILE *om); extern int p7_pli_NewSeq (P7_PIPELINE *pli, const ESL_SQ *sq); -extern int p7_pli_computeAliScores_BATH(P7_DOMAIN *dom, P7_TRACE *tr, const ESL_SQ *orfsq, const P7_PROFILE *gm); +extern int p7_pli_computeAliScores_BATH(P7_DOMAIN *dom, P7_TRACE *tr, const ESL_ORF *orfsq, const P7_PROFILE *gm); extern int p7_pli_computeAliScores_Frameshift_BATH(P7_DOMAIN *dom, P7_TRACE *tr, const ESL_SQ *seq, const P7_FS_PROFILE *gm_fs); -extern int p7_Pipeline_BATH (P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_SQ *dnasq, ESL_SQ_BLOCK *orf_block, ESL_GENCODE *gcode, P7_HMM_WINDOWLIST *splcing_windows, int complementarity); +extern int p7_Pipeline_BATH (P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_SQ *dnasq, ESL_ORF_BLOCK *orf_block, ESL_GENCODE *gcode, P7_HMM_WINDOWLIST *splcing_windows, int complementarity); extern int p7_pli_Statistics(FILE *ofp, P7_PIPELINE *pli, ESL_STOPWATCH *w); diff --git a/src/p7_alidisplay.c b/src/p7_alidisplay.c index daebde8f..d843d471 100644 --- a/src/p7_alidisplay.c +++ b/src/p7_alidisplay.c @@ -934,7 +934,7 @@ p7_alidisplay_fs_Create(const P7_TRACE *tr, int which, const P7_FS_PROFILE *gm_f * in the data. */ P7_ALIDISPLAY * -p7_alidisplay_nonfs_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *sq, const ESL_SQ *orfsq, int orf_pos, int show_cigar) +p7_alidisplay_nonfs_Create(const P7_TRACE *tr, int which, const P7_OPROFILE *om, const ESL_SQ *sq, const ESL_ORF *orfsq, int orf_pos, int show_cigar) { P7_ALIDISPLAY *ad = NULL; diff --git a/src/p7_domaindef.c b/src/p7_domaindef.c index 57927a6e..d880501d 100644 --- a/src/p7_domaindef.c +++ b/src/p7_domaindef.c @@ -23,7 +23,7 @@ static int is_multidomain_region_frameshift (P7_DOMAINDEF *ddef, int i, int j); static int region_trace_ensemble (P7_DOMAINDEF *ddef, const P7_OPROFILE *om, const ESL_DSQ *dsq, int ireg, int jreg, const P7_OMX *fwd, P7_OMX *wrk, int *ret_nc); static int region_trace_ensemble_frameshift (P7_DOMAINDEF *ddef, const P7_FS_OPROFILE *om_fs, const ESL_DSQ *dsq, const ESL_ALPHABET *abc, int ireg, int jreg, const P7_OMX *fwd, int *ret_nc); static int rescore_isolated_domain_frameshift(P7_DOMAINDEF *ddef, P7_PIPELINE *pli, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, ESL_SQ *windowsq, int i, int j, P7_BG *bg, const ESL_GENCODE *gcode); -static int rescore_isolated_domain_bath(P7_DOMAINDEF *ddef, P7_OPROFILE *om, P7_PROFILE *gm, const ESL_SQ *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OMX *ox1, P7_OMX *ox2, int i, int j, int null2_is_done); +static int rescore_isolated_domain_bath(P7_DOMAINDEF *ddef, P7_OPROFILE *om, P7_PROFILE *gm, const ESL_ORF *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OMX *ox1, P7_OMX *ox2, int i, int j, int null2_is_done); /***************************************************************** * 1. The P7_DOMAINDEF object: allocation, reuse, destruction @@ -496,7 +496,7 @@ p7_domaindef_ByPosteriorHeuristics_Frameshift_BATH(P7_PIPELINE *pli, ESL_SQ *win * models. */ int -p7_domaindef_ByPosteriorHeuristics_BATH(const ESL_SQ *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OPROFILE *om, P7_PROFILE *gm, P7_OMX *oxf, P7_OMX *oxb, P7_OMX *fwd, P7_OMX *bck, P7_DOMAINDEF *ddef) +p7_domaindef_ByPosteriorHeuristics_BATH(const ESL_ORF *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OPROFILE *om, P7_PROFILE *gm, P7_OMX *oxf, P7_OMX *oxb, P7_OMX *fwd, P7_OMX *bck, P7_DOMAINDEF *ddef) { int i, j; int triggered; @@ -1223,7 +1223,7 @@ rescore_isolated_domain_frameshift(P7_DOMAINDEF *ddef, P7_PIPELINE *pli, P7_FS_O * */ static int -rescore_isolated_domain_bath(P7_DOMAINDEF *ddef, P7_OPROFILE *om, P7_PROFILE *gm, const ESL_SQ *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OMX *ox1, P7_OMX *ox2, int i, int j, int null2_is_done) +rescore_isolated_domain_bath(P7_DOMAINDEF *ddef, P7_OPROFILE *om, P7_PROFILE *gm, const ESL_ORF *orfsq, const ESL_SQ *windowsq, const int64_t ntsqlen, P7_OMX *ox1, P7_OMX *ox2, int i, int j, int null2_is_done) { P7_DOMAIN *dom = NULL; diff --git a/src/p7_hmmwindow.c b/src/p7_hmmwindow.c index 1a09af74..a80e7c0d 100644 --- a/src/p7_hmmwindow.c +++ b/src/p7_hmmwindow.c @@ -587,7 +587,14 @@ p7_hmmwindow_GetSeedHits(P7_HMM_WINDOWLIST *hw, const P7_TOPHITS *th, P7_HMM *hm hit->dcl->scores_per_pos = NULL; hit->dcl->k_per_pos = NULL; - p7_pli_computeAliScores_BATH(hit->dcl, hit->dcl->tr, orfsq, gm); + { + ESL_ORF orf; /* the translated window, as the ORF the scoring takes */ + orf.dsq = orfsq->dsq; + orf.n = n; + orf.start = 0; + orf.end = 0; + p7_pli_computeAliScores_BATH(hit->dcl, hit->dcl->tr, &orf, gm); + } esl_sq_Destroy(orfsq); diff --git a/src/p7_pipeline.c b/src/p7_pipeline.c index 6e77864b..722c5cc1 100644 --- a/src/p7_pipeline.c +++ b/src/p7_pipeline.c @@ -329,7 +329,7 @@ p7_pli_ComputeLocalCompo(const P7_SCOREDATA *data, const P7_OPROFILE *om, const int -p7_pli_BuildDNAWindows(P7_PIPELINE *pli, ESL_SQ_BLOCK *orf_block, ESL_SQ *dnasq, P7_OPROFILE *om, P7_BG *bg, const P7_SCOREDATA *data, P7_HMM_WINDOWLIST *windowlist, float pct_overlap, P7_PIPELINE_OBJS *pli_tmp, P7_HMM_WINDOWLIST *hit_windows, int hit_windows_start, int complementarity) +p7_pli_BuildDNAWindows(P7_PIPELINE *pli, ESL_ORF_BLOCK *orf_block, ESL_SQ *dnasq, P7_OPROFILE *om, P7_BG *bg, const P7_SCOREDATA *data, P7_HMM_WINDOWLIST *windowlist, float pct_overlap, P7_PIPELINE_OBJS *pli_tmp, P7_HMM_WINDOWLIST *hit_windows, int hit_windows_start, int complementarity) { int i, f, w; @@ -345,7 +345,7 @@ p7_pli_BuildDNAWindows(P7_PIPELINE *pli, ESL_SQ_BLOCK *orf_block, ESL_SQ *dnasq, int64_t overlap_end; int32_t overlap_len; int new_hit_cnt = 0; - ESL_SQ *curr_orf; + ESL_ORF *curr_orf; for(f = 0; f < orf_block->count; f++) { @@ -650,7 +650,7 @@ p7_pipeline_Merge(P7_PIPELINE *p1, P7_PIPELINE *p2) * Throws: on allocation failure. */ int -p7_pli_computeAliScores_BATH(P7_DOMAIN *dom, P7_TRACE *tr, const ESL_SQ *orfsq, const P7_PROFILE *gm) +p7_pli_computeAliScores_BATH(P7_DOMAIN *dom, P7_TRACE *tr, const ESL_ORF *orfsq, const P7_PROFILE *gm) { int i, k, n; @@ -1144,7 +1144,7 @@ p7_pli_postDomainDef_Frameshift_BATH(P7_PIPELINE *pli, P7_FS_PROFILE *gm_fs5, P7 */ static int -p7_pli_postDomainDef_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, int window_start, ESL_SQ *orfsq, ESL_SQ *dnasq, ESL_SQ *windowsq, int complementarity) +p7_pli_postDomainDef_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, int window_start, ESL_ORF *orfsq, ESL_SQ *dnasq, ESL_SQ *windowsq, int complementarity) { int d; @@ -1313,7 +1313,7 @@ return eslOK; * */ static int -p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_SQ_BLOCK *orf_block, ESL_SQ *dnasq, ESL_GENCODE *gcode, P7_PIPELINE_OBJS *pli_tmp, P7_HMM_WINDOWLIST *hit_windows, int hit_windows_start, int complementarity) +p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_ORF_BLOCK *orf_block, ESL_SQ *dnasq, ESL_GENCODE *gcode, P7_PIPELINE_OBJS *pli_tmp, P7_HMM_WINDOWLIST *hit_windows, int hit_windows_start, int complementarity) { int i, w, h; @@ -1334,7 +1334,7 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF double P_null; /* P-value of frameshift forward for window w/o bias adjustments*/ double P_tot; /* P-value of summed forward score for all ORFs */ double P_min; /* lowest p-value produced by an ORF */ - ESL_SQ *orfsq; + ESL_ORF *orfsq; P7_HMM_WINDOWLIST fwd_windowlist; P7_HMM_WINDOW *dna_window; @@ -1562,7 +1562,7 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF * Xref: J4/25. */ int -p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_SQ *dnasq, ESL_SQ_BLOCK *orf_block, ESL_GENCODE *gcode,P7_HMM_WINDOWLIST *hit_windows, int complementarity) +p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFILE *om_fs3, P7_FS_OPROFILE *om_fs5, P7_FS_PROFILE *gm_fs5, P7_SCOREDATA *data, P7_BG *bg, P7_TOPHITS *hitlist, int64_t seqidx, ESL_SQ *dnasq, ESL_ORF_BLOCK *orf_block, ESL_GENCODE *gcode,P7_HMM_WINDOWLIST *hit_windows, int complementarity) { int i, w; @@ -1580,7 +1580,7 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI float filtersc; /* bias and null score */ float local_filtersc; /* bias and null score using local_compo */ double P; /* p-value holder */ - ESL_SQ *orfsq; /* ORF sequence */ + ESL_ORF *orfsq; /* ORF sequence */ P7_HMM_WINDOW *window; P7_PIPELINE_OBJS *pli_tmp; @@ -1904,6 +1904,7 @@ utest_computeAliScores_matches_frameshift(ESL_RANDOMNESS *r, P7_HMM *hmm, P7_PRO P7_TRACE *tr_amino = NULL; P7_TRACE *tr_dna = NULL; P7_DOMAIN dom_new, dom_old; + ESL_ORF orf; float vsc; int idx, z, i, n; double tol = 1e-5; @@ -1942,7 +1943,11 @@ utest_computeAliScores_matches_frameshift(ESL_RANDOMNESS *r, P7_HMM *hmm, P7_PRO memset(&dom_new, 0, sizeof(dom_new)); memset(&dom_old, 0, sizeof(dom_old)); - if (p7_pli_computeAliScores_BATH(&dom_new, tr_amino, sq, gm) != eslOK) esl_fatal(msg); + orf.dsq = sq->dsq; /* the emitted protein, as the ORF the scoring takes */ + orf.n = sq->n; + orf.start = 0; + orf.end = 0; + if (p7_pli_computeAliScores_BATH(&dom_new, tr_amino, &orf, gm) != eslOK) esl_fatal(msg); if (p7_pli_computeAliScores_Frameshift_BATH(&dom_old, tr_dna, windowsq, gm_fs5) != eslOK) esl_fatal(msg); if (dom_new.per_pos_len != dom_old.per_pos_len) esl_fatal(msg); From 8e4d09e25c20bf09e401496e339acdc6137e555e Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sat, 3 Oct 2026 22:43:12 -0500 Subject: [PATCH 4/7] p7_trace_fs_Clone(): allocate for the trace's length The clone took the source's allocation sizes, nalloc and ndomalloc, and copied that much. The domain definition's working trace keeps the size of the longest trace it has held, so every hit's copy was that large: about 2,000 states allocated for an average of 122 used, 34 kB a hit. The clone now allocates and copies what the trace holds. Hits are unchanged. This is a prerequisite for the one-pass bathsearch driver, which holds the hits of every query until the pass ends: with it, a search of 5,222 profiles against a 352 Mb genome peaks at 5.9 GB where it peaked at 28.0 GB. Co-Authored-By: Claude Opus 5.5 --- src/p7_trace.c | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/p7_trace.c b/src/p7_trace.c index d3a83e17..163b6640 100644 --- a/src/p7_trace.c +++ b/src/p7_trace.c @@ -299,8 +299,8 @@ p7_trace_fs_Clone(const P7_TRACE *tr) dup->pp = NULL; dup->N = tr->N; - dup->nalloc = tr->nalloc; - dup->ndomalloc = tr->ndomalloc; + dup->nalloc = ESL_MAX(tr->N, 1); /* room for what holds, not for what it has allocated: a */ + dup->ndomalloc = ESL_MAX(tr->ndom, 1); /* working trace keeps the size of the longest one it has held */ dup->ndom = tr->ndom; dup->M = tr->M; dup->L = tr->L; @@ -312,48 +312,48 @@ p7_trace_fs_Clone(const P7_TRACE *tr) /* The trace data itself */ ESL_ALLOC(dup->st, sizeof(char) * dup->nalloc); - memcpy(dup->st, tr->st, sizeof(char) * dup->nalloc); + memcpy(dup->st, tr->st, sizeof(char) * dup->N); ESL_ALLOC(dup->k, sizeof(int) * dup->nalloc); - memcpy(dup->k, tr->k, sizeof(int) * dup->nalloc); + memcpy(dup->k, tr->k, sizeof(int) * dup->N); ESL_ALLOC(dup->i, sizeof(int) * dup->nalloc); - memcpy(dup->i, tr->i, sizeof(int) * dup->nalloc); + memcpy(dup->i, tr->i, sizeof(int) * dup->N); if (tr->c != NULL) { ESL_ALLOC(dup->c, sizeof(int) * dup->nalloc); - memcpy(dup->c, tr->c, sizeof(int) * dup->nalloc); + memcpy(dup->c, tr->c, sizeof(int) * dup->N); } if (tr->pp != NULL) { ESL_ALLOC(dup->pp, sizeof(float) * dup->nalloc); - memcpy(dup->pp, tr->pp, sizeof(float) * dup->nalloc); + memcpy(dup->pp, tr->pp, sizeof(float) * dup->N); } if (tr->sp != NULL) { ESL_ALLOC(dup->sp, sizeof(float) * dup->nalloc); - memcpy(dup->sp, tr->sp, sizeof(float) * dup->nalloc); + memcpy(dup->sp, tr->sp, sizeof(float) * dup->N); } /* The trace's index: table of domain start/stop coords */ ESL_ALLOC(dup->tfrom, sizeof(int) * dup->ndomalloc); - memcpy(dup->tfrom, tr->tfrom, sizeof(int) * dup->ndomalloc); + memcpy(dup->tfrom, tr->tfrom, sizeof(int) * dup->ndom); ESL_ALLOC(dup->tto, sizeof(int) * dup->ndomalloc); - memcpy(dup->tto, tr->tto, sizeof(int) * dup->ndomalloc); + memcpy(dup->tto, tr->tto, sizeof(int) * dup->ndom); ESL_ALLOC(dup->sqfrom, sizeof(int) * dup->ndomalloc); - memcpy(dup->sqfrom, tr->sqfrom, sizeof(int) * dup->ndomalloc); + memcpy(dup->sqfrom, tr->sqfrom, sizeof(int) * dup->ndom); ESL_ALLOC(dup->sqto, sizeof(int) * dup->ndomalloc); - memcpy(dup->sqto, tr->sqto, sizeof(int) * dup->ndomalloc); + memcpy(dup->sqto, tr->sqto, sizeof(int) * dup->ndom); ESL_ALLOC(dup->hmmfrom, sizeof(int) * dup->ndomalloc); - memcpy(dup->hmmfrom, tr->hmmfrom, sizeof(int) * dup->ndomalloc); + memcpy(dup->hmmfrom, tr->hmmfrom, sizeof(int) * dup->ndom); ESL_ALLOC(dup->hmmto, sizeof(int) * dup->ndomalloc); - memcpy(dup->hmmto, tr->hmmto, sizeof(int) * dup->ndomalloc); + memcpy(dup->hmmto, tr->hmmto, sizeof(int) * dup->ndom); return dup; From 862561b6f9fd2f0397e89572906a94ae2dff54be Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sun, 4 Oct 2026 00:46:50 -0500 Subject: [PATCH 5/7] ssvfilter: SSV maximum of every ORF in a block p7_SSVFilter_OrfBlock() gives the highest cell of the SSV matrix for each ORF of a block in one call, laying out the bands once for the model. p7_SSVFilter_FromXE() turns such a maximum into the status and score that p7_SSVFilter() returns; p7_SSVFilter_avx() now ends by calling it. AVX2 only: on other processors the block call returns eslENORESULT. Co-Authored-By: Claude Opus 5.5 --- src/impl_avx/impl_avx.h | 7 +++ src/impl_avx/ssvfilter.c | 74 ++++++++++++++++++++++++++++++ src/impl_avx/ssvfilter_avx.c | 89 +++++++++++++++++++----------------- 3 files changed, 129 insertions(+), 41 deletions(-) diff --git a/src/impl_avx/impl_avx.h b/src/impl_avx/impl_avx.h index 6dc3835e..7b505fce 100644 --- a/src/impl_avx/impl_avx.h +++ b/src/impl_avx/impl_avx.h @@ -15,6 +15,7 @@ #include "esl_alphabet.h" #include "esl_random.h" #include "esl_cpu.h" +#include "esl_gencode.h" /* ESL_ORF */ #ifdef eslENABLE_SSE #include /* SSE */ @@ -788,7 +789,13 @@ extern int p7_SSVFilter_sse(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, fl /* ssvfilter_avx.c — AVX2 implementation */ #ifdef eslENABLE_AVX extern int p7_SSVFilter_avx(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc); +extern void p7_SSVFilter_OrfBlock_avx(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE); #endif +/* SSV over a block of ORFs: the pipeline can take each ORF's SSV maximum from + * one call, and turn it into the SSV filter's status and score itself */ +#define p7_SSV_ORFBLOCK 1 +extern int p7_SSVFilter_OrfBlock(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE); +extern int p7_SSVFilter_FromXE(int xE, const P7_OPROFILE *om, float *ret_sc); /* stotrace.c — self-patching dispatcher */ extern int (*p7_StochasticTrace)(ESL_RANDOMNESS *rng, const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, const P7_OMX *ox, P7_TRACE *tr); diff --git a/src/impl_avx/ssvfilter.c b/src/impl_avx/ssvfilter.c index 39974b3e..decb136f 100644 --- a/src/impl_avx/ssvfilter.c +++ b/src/impl_avx/ssvfilter.c @@ -42,3 +42,77 @@ p7_SSVFilter_Dispatcher(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float return eslENORESULT; #endif } + + +/* Function: p7_SSVFilter_FromXE() + * Synopsis: The SSV filter's result, from the SSV matrix's highest cell. + * + * Purpose: is the highest cell of the SSV matrix of a sequence against + * , as p7_SSVFilter_OrfBlock() gives it. Return what + * p7_SSVFilter() returns for that sequence, and its score in + * <*ret_sc>; must be configured for the sequence's length. + * + * Returns: , or , as p7_SSVFilter(). + */ +int +p7_SSVFilter_FromXE(int xE_in, const P7_OPROFILE *om, float *ret_sc) +{ + uint16_t xE = (uint16_t) xE_in; + uint16_t xJ; + + if (om->tjb_b + om->tbm_b + om->tec_b + om->bias_b >= 127) + return eslENORESULT; + + /* Saturation floors every diagonal at the begin score (128), so a + * max of 128 means no diagonal scored above it and the true best + * may be lower; let the full MSV filter compute it. */ + if (xE <= 128) return eslENORESULT; + + if (xE >= 255 - om->bias_b) { + *ret_sc = eslINFINITY; + + if (om->base_b - om->tjb_b - om->tbm_b < 128) + return eslENORESULT; + + return eslERANGE; + } + + xE += om->base_b - om->tjb_b - om->tbm_b; + xE -= 128; + + if (xE >= 255 - om->bias_b) { + *ret_sc = eslINFINITY; + return eslERANGE; + } + + xJ = xE - om->tec_b; + + if (xJ > om->base_b) return eslENORESULT; + + *ret_sc = ((float)(xJ - om->tjb_b) - (float) om->base_b); + *ret_sc /= om->scale_b; + *ret_sc -= 3.0f; + + return eslOK; +} + + +/* Function: p7_SSVFilter_OrfBlock() + * Synopsis: SSV maximum of every ORF in a block. + * + * Purpose: For each of the ORFs in , store in the highest + * cell of its SSV matrix against . It does not depend on the + * length is configured for; p7_SSVFilter_FromXE() makes the + * filter's status and score from it. + * + * Returns: on success; if this processor has no + * implementation, and is unset. + */ +int +p7_SSVFilter_OrfBlock(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE) +{ +#ifdef eslENABLE_AVX + if (esl_cpu_has_avx()) { p7_SSVFilter_OrfBlock_avx(om, orf, n, xE); return eslOK; } +#endif + return eslENORESULT; +} diff --git a/src/impl_avx/ssvfilter_avx.c b/src/impl_avx/ssvfilter_avx.c index e2254274..0bb9428e 100644 --- a/src/impl_avx/ssvfilter_avx.c +++ b/src/impl_avx/ssvfilter_avx.c @@ -436,6 +436,18 @@ calc_band_18(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, int q, __m256i be * 2. p7_SSVFilter_avx() *****************************************************************/ +/* the band function for each band width */ +static __m256i (* const fs[MAX_BANDS + 1])(const ESL_DSQ *, int, const P7_OPROFILE *, int, register __m256i, __m256i) + = { NULL + , calc_band_1, calc_band_2, calc_band_3, calc_band_4, calc_band_5, calc_band_6 +#if MAX_BANDS > 6 + , calc_band_7, calc_band_8, calc_band_9, calc_band_10, calc_band_11, calc_band_12, calc_band_13, calc_band_14 +#endif +#if MAX_BANDS > 14 + , calc_band_15, calc_band_16, calc_band_17, calc_band_18 +#endif + }; + static uint8_t get_xE_avx(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om) { @@ -448,17 +460,6 @@ get_xE_avx(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om) int last_q = 0; int i; - __m256i (*fs[MAX_BANDS + 1])(const ESL_DSQ *, int, const P7_OPROFILE *, int, register __m256i, __m256i) - = { NULL - , calc_band_1, calc_band_2, calc_band_3, calc_band_4, calc_band_5, calc_band_6 -#if MAX_BANDS > 6 - , calc_band_7, calc_band_8, calc_band_9, calc_band_10, calc_band_11, calc_band_12, calc_band_13, calc_band_14 -#endif -#if MAX_BANDS > 14 - , calc_band_15, calc_band_16, calc_band_17, calc_band_18 -#endif - }; - beginv = _mm256_set1_epi8(-128); xEv = beginv; @@ -497,45 +498,51 @@ get_xE_avx(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om) int p7_SSVFilter_avx(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) { - uint16_t xE; - uint16_t xJ; - if (om->tjb_b + om->tbm_b + om->tec_b + om->bias_b >= 127) return eslENORESULT; - xE = get_xE_avx(dsq, L, om); - - /* Saturation floors every diagonal at the begin score (128), so a - * max of 128 means no diagonal scored above it and the true best - * may be lower; let the full MSV filter compute it. */ - if (xE <= 128) return eslENORESULT; - - if (xE >= 255 - om->bias_b) { - *ret_sc = eslINFINITY; + return p7_SSVFilter_FromXE(get_xE_avx(dsq, L, om), om, ret_sc); +} - if (om->base_b - om->tjb_b - om->tbm_b < 128) - return eslENORESULT; - return eslERANGE; +/* Function: p7_SSVFilter_OrfBlock_avx() + * Synopsis: SSV maximum of every ORF in a block, AVX2 path. + * + * Purpose: See p7_SSVFilter_OrfBlock(). is what get_xE_avx() gives + * for ORF . The bands depend only on the model, so they are laid + * out once for the block. + */ +#define SSV_BLOCK_BANDS 64 /* bands laid out up front; enough for a model of 28,672 positions */ +void +p7_SSVFilter_OrfBlock_avx(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE) +{ + __m256i beginv = _mm256_set1_epi8(-128); + __m256i xEv; + int Q = p7O_NQB_AVX(om->M); + int bands = (Q + MAX_BANDS - 1) / MAX_BANDS; + int q0[SSV_BLOCK_BANDS]; /* first vector of each band */ + int w[SSV_BLOCK_BANDS]; /* its width */ + int last_q = 0; + int q, b, i; + + if (bands > SSV_BLOCK_BANDS) { + for (i = 0; i < n; i++) xE[i] = get_xE_avx(orf[i].dsq, (int) orf[i].n, om); + return; } - xE += om->base_b - om->tjb_b - om->tbm_b; - xE -= 128; - - if (xE >= 255 - om->bias_b) { - *ret_sc = eslINFINITY; - return eslERANGE; + for (b = 0; b < bands; b++) { + q = (Q * (b + 1)) / bands; + q0[b] = last_q; + w[b] = q - last_q; + last_q = q; } - xJ = xE - om->tec_b; - - if (xJ > om->base_b) return eslENORESULT; - - *ret_sc = ((float)(xJ - om->tjb_b) - (float) om->base_b); - *ret_sc /= om->scale_b; - *ret_sc -= 3.0f; - - return eslOK; + for (i = 0; i < n; i++) { + xEv = beginv; + for (b = 0; b < bands; b++) + xEv = fs[w[b]](orf[i].dsq, (int) orf[i].n, om, q0[b], beginv, xEv); + xE[i] = esl_avx_hmax_epu8(xEv); + } } #endif /* eslENABLE_AVX */ From de3d6cbc03f48560164ce9861d6f62cafe8a52d8 Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sun, 4 Oct 2026 00:46:50 -0500 Subject: [PATCH 6/7] pipeline: reject ORFs on their SSV maximum For every ORF, p7_Pipeline_BATH() configured the model and the null model for its length and took logs and exps to turn its MSV score into a P-value, although about 99% of ORFs are rejected there by SSV alone. It now runs SSV over the whole block first, and keeps for each ORF length the lowest SSV maximum that is not such a rejection, found by taking the same steps once per length. An ORF below it is skipped. The rest take the MSV stage as before, and so does an ORF whose maximum is at the floor, where the SSV filter gives no result. The output is unchanged. Building with -Dp7_SSV_ORFBLOCK_CHECK runs both and stops if they disagree. CPU cycles for bathsearch against a 352 Mb genome fall 5-17% for one profile and 20% for 50 profiles (15% with --fs). Co-Authored-By: Claude Opus 5.5 --- src/hmmer.h | 1 + src/p7_pipeline.c | 80 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) diff --git a/src/hmmer.h b/src/hmmer.h index 792d6479..0dda2845 100644 --- a/src/hmmer.h +++ b/src/hmmer.h @@ -1130,6 +1130,7 @@ typedef struct p7_pipeline_s { enum p7_pipemodes_e mode; /* p7_SCAN_MODELS | p7_SEARCH_SEQS */ int spliced; /* TRUE if user uses --splice slaf to enable spliced alignments */ int fs_pipe; /* TRUE if bathsearch is allowed to use the frameshift aware pipeline branch (use --fs flag) */ + uint16_t *ssv_cut; /* ssv_cut[L]: lowest SSV maximum at which an ORF of length L is not rejected by SSV alone; 0 if not yet known */ int std_pipe; /* TRUE if bathsearch is allowed to use the standard translation pipeline (do not use --fsonly flag) */ int strands; /* p7_STRAND_TOPONLY | p7_STRAND_BOTTOMONLY | p7_STRAND_BOTH */ int W; /* window length for nhmmer scan - essentially maximum length of model that we expect to find*/ diff --git a/src/p7_pipeline.c b/src/p7_pipeline.c index 722c5cc1..97b41b99 100644 --- a/src/p7_pipeline.c +++ b/src/p7_pipeline.c @@ -106,6 +106,7 @@ p7_pipeline_Create_BATH(ESL_GETOPTS *go, int M_hint, int L_hint, enum p7_pipemod pli->spliced = (go ? esl_opt_IsUsed(go, "--splice") : 0); pli->fs_pipe = (go ? (esl_opt_IsUsed(go, "--fs") || esl_opt_IsUsed(go, "--fsonly")) : 0); pli->std_pipe = (go ? !esl_opt_IsUsed(go, "--fsonly") : 1); + pli->ssv_cut = NULL; /* Create sparce memeory forward and backward optimized matricies for use in the * non-frameshift pipeline branch @@ -283,6 +284,7 @@ p7_pipeline_Destroy_BATH(P7_PIPELINE *pli) p7_omx_Destroy(pli->bck); esl_randomness_Destroy(pli->r); p7_domaindef_Destroy_BATH(pli->ddef); + if (pli->ssv_cut != NULL) free(pli->ssv_cut); free(pli); } @@ -502,6 +504,8 @@ p7_pli_NewModel(P7_PIPELINE *pli, const P7_OPROFILE *om, P7_BG *bg) if (pli->do_biasfilter) p7_bg_SetFilter(bg, om->M, om->compo); + if (pli->ssv_cut != NULL) { free(pli->ssv_cut); pli->ssv_cut = NULL; } /* the cutoffs belong to the model */ + if (pli->mode == p7_SEARCH_SEQS) status = p7_pli_NewModelThresholds(pli, om); @@ -1505,6 +1509,41 @@ p7_pli_Frameshift(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROF +#ifdef p7_SSV_ORFBLOCK +#define p7_SSV_CUT_MAXL 4096 /* ORFs this long or longer always take the full MSV stage */ + +/* ssv_cutoff() + * The lowest SSV maximum (see p7_SSVFilter_OrfBlock()) at which an ORF of + * length is not rejected by the SSV filter alone: below it, the filter + * returns with a score whose P-value is above F1. Found by taking the + * steps of the MSV stage of p7_Pipeline_BATH() for each in turn, so an + * ORF below it is one that stage rejects. 128 is the lowest there is, and + * at 128 the SSV filter gives no result, so the search starts above it and the + * caller sends such an ORF to the MSV stage; 256 means no gets through. + * Leaves and configured for , as that stage does. + */ +static int +ssv_cutoff(const P7_PIPELINE *pli, P7_OPROFILE *om, P7_BG *bg, int L) +{ + float nullsc, usc, seqsc; + double P; + int xE; + + p7_bg_SetLength(bg, L); + p7_oprofile_ReconfigLength(om, L); + p7_bg_NullOne (bg, NULL, L, &nullsc); + + for (xE = 129; xE < 256; xE++) + { + if (p7_SSVFilter_FromXE(xE, om, &usc) != eslOK) break; + seqsc = (usc - nullsc) / eslCONST_LOG2; + P = esl_gumbel_surv( seqsc, om->evparam[p7_MMU], om->evparam[p7_MLAMBDA]); + if (! (P > pli->F1)) break; + } + return xE; +} +#endif /*p7_SSV_ORFBLOCK*/ + /* Function: p7_Pipeline_BATH() * Synopsis: Sequence to profile comparison pipeline for * frameshift aware translated search - bathsearch. @@ -1583,6 +1622,9 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI ESL_ORF *orfsq; /* ORF sequence */ P7_HMM_WINDOW *window; P7_PIPELINE_OBJS *pli_tmp; + uint8_t *ssv_xE = NULL; /* SSV maximum of each ORF, if SSV was run over the block */ + int ssv_len = 0; /* length of the last ORF rejected on its SSV maximum since , were configured */ + int ssv_rej; /* TRUE if this ORF's SSV maximum rejects it */ if (dnasq->n < 15) return eslOK; //DNA to short if (orf_block->count == 0) return eslOK; //No ORFS translated @@ -1616,6 +1658,20 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI if ((status = esl_sq_SetAccession(pli_tmp->tmpseq, dnasq->acc)) != eslOK) goto ERROR; if ((status = esl_sq_SetDesc (pli_tmp->tmpseq, dnasq->desc)) != eslOK) goto ERROR; +#ifdef p7_SSV_ORFBLOCK + /* Run SSV over the whole block. Most ORFs are rejected by SSV alone, and for + * those its maximum and the ORF's length decide it; see ssv_cutoff(). */ + if (pli->F1 < 1.0) + { + ESL_ALLOC(ssv_xE, sizeof(uint8_t) * orf_block->count); + if (pli->ssv_cut == NULL) { + ESL_ALLOC(pli->ssv_cut, sizeof(uint16_t) * p7_SSV_CUT_MAXL); + for (i = 0; i < p7_SSV_CUT_MAXL; i++) pli->ssv_cut[i] = 0; + } + if (p7_SSVFilter_OrfBlock(om, orf_block->list, orf_block->count, ssv_xE) != eslOK) { free(ssv_xE); ssv_xE = NULL; } + } +#endif + for (i = 0; i < orf_block->count; ++i) { orfsq = &(orf_block->list[i]); @@ -1627,6 +1683,19 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI { vfsc = -eslINFINITY; + ssv_rej = FALSE; +#ifdef p7_SSV_ORFBLOCK + if (ssv_xE != NULL && orfsq->n < p7_SSV_CUT_MAXL) + { + if (pli->ssv_cut[orfsq->n] == 0) { pli->ssv_cut[orfsq->n] = ssv_cutoff(pli, om, bg, orfsq->n); ssv_len = 0; } + if (ssv_xE[i] > 128 && ssv_xE[i] < pli->ssv_cut[orfsq->n]) ssv_rej = TRUE; + } +#ifndef p7_SSV_ORFBLOCK_CHECK + if (ssv_rej) { ssv_len = orfsq->n; continue; } +#endif +#endif + ssv_len = 0; + p7_bg_SetLength(bg, orfsq->n); p7_oprofile_ReconfigLength(om, orfsq->n); p7_bg_NullOne (bg, orfsq->dsq, orfsq->n, &nullsc); @@ -1636,6 +1705,9 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI p7_MSVFilter(orfsq->dsq, orfsq->n, om, pli->oxf, &usc); seqsc = (usc - nullsc) / eslCONST_LOG2; P = esl_gumbel_surv( seqsc, om->evparam[p7_MMU], om->evparam[p7_MLAMBDA]); +#ifdef p7_SSV_ORFBLOCK_CHECK + if (ssv_rej && ! (P > pli->F1)) p7_Die("SSV block: rejected an ORF the MSV stage passes (L %d, xE %d, cutoff %d, P %g)", (int) orfsq->n, ssv_xE[i], pli->ssv_cut[orfsq->n], P); +#endif if (P > pli->F1 ) continue; pli->pos_past_msv += orfsq->n * 3; @@ -1777,6 +1849,13 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI } } + /* Leave and configured for the last ORF looked at, as they are when every ORF takes the MSV stage */ + if (ssv_len > 0) { + p7_bg_SetLength(bg, ssv_len); + p7_oprofile_ReconfigLength(om, ssv_len); + } + if (ssv_xE != NULL) { free(ssv_xE); ssv_xE = NULL; } + if(pli->fs_pipe) { p7_pli_Frameshift(pli, om, gm, om_fs3, om_fs5, gm_fs5, data, bg, hitlist, seqidx, orf_block, dnasq, gcode, pli_tmp, hit_windows, hit_windows_start, complementarity); } @@ -1796,6 +1875,7 @@ p7_Pipeline_BATH(P7_PIPELINE *pli, P7_OPROFILE *om, P7_PROFILE *gm, P7_FS_OPROFI return eslOK; ERROR: + if (ssv_xE != NULL) free(ssv_xE); if (pli_tmp != NULL) { if (pli_tmp->tmpseq != NULL) esl_sq_Destroy(pli_tmp->tmpseq); From 9b5e41ba5e6233f06e25c6e951f92ffa4746d17d Mon Sep 17 00:00:00 2001 From: Kevin Moran <39282665+Thernn88@users.noreply.github.com> Date: Sun, 4 Oct 2026 00:57:03 -0500 Subject: [PATCH 7/7] ssvfilter: SSV over an ORF block for the SSE builds p7_SSVFilter_OrfBlock() and p7_SSVFilter_FromXE() are added to the SSE code of the AVX build and to the SSE-only build, each as a loop over its get_xE(). The NEON build has neither and keeps the per-ORF MSV stage. The two SSE builds give the same output as before for 23 profiles against 14 Mb in default, --fs and --splice modes. Co-Authored-By: Claude Opus 5.5 --- src/impl_avx/impl_avx.h | 1 + src/impl_avx/ssvfilter.c | 7 +++++- src/impl_avx/ssvfilter_sse.c | 46 ++++++++++-------------------------- src/impl_sse/impl_sse.h | 5 ++++ src/impl_sse/ssvfilter.c | 43 +++++++++++++++++++++++++++++---- 5 files changed, 64 insertions(+), 38 deletions(-) diff --git a/src/impl_avx/impl_avx.h b/src/impl_avx/impl_avx.h index 7b505fce..ff0130c9 100644 --- a/src/impl_avx/impl_avx.h +++ b/src/impl_avx/impl_avx.h @@ -785,6 +785,7 @@ extern int (*p7_SSVFilter)(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, flo /* ssvfilter_sse.c — SSE implementation */ #ifdef eslENABLE_SSE extern int p7_SSVFilter_sse(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc); +extern void p7_SSVFilter_OrfBlock_sse(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE); #endif /* ssvfilter_avx.c — AVX2 implementation */ #ifdef eslENABLE_AVX diff --git a/src/impl_avx/ssvfilter.c b/src/impl_avx/ssvfilter.c index decb136f..07156a64 100644 --- a/src/impl_avx/ssvfilter.c +++ b/src/impl_avx/ssvfilter.c @@ -105,7 +105,7 @@ p7_SSVFilter_FromXE(int xE_in, const P7_OPROFILE *om, float *ret_sc) * length is configured for; p7_SSVFilter_FromXE() makes the * filter's status and score from it. * - * Returns: on success; if this processor has no + * Returns: on success; if there is no SIMD * implementation, and is unset. */ int @@ -114,5 +114,10 @@ p7_SSVFilter_OrfBlock(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t #ifdef eslENABLE_AVX if (esl_cpu_has_avx()) { p7_SSVFilter_OrfBlock_avx(om, orf, n, xE); return eslOK; } #endif +#ifdef eslENABLE_SSE + p7_SSVFilter_OrfBlock_sse(om, orf, n, xE); + return eslOK; +#else return eslENORESULT; +#endif } diff --git a/src/impl_avx/ssvfilter_sse.c b/src/impl_avx/ssvfilter_sse.c index eafc426e..11a25028 100644 --- a/src/impl_avx/ssvfilter_sse.c +++ b/src/impl_avx/ssvfilter_sse.c @@ -500,45 +500,25 @@ get_xE_sse(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om) int p7_SSVFilter_sse(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) { - uint16_t xE; - uint16_t xJ; - if (om->tjb_b + om->tbm_b + om->tec_b + om->bias_b >= 127) return eslENORESULT; - xE = get_xE_sse(dsq, L, om); - - /* Saturation floors every diagonal at the begin score (128), so a - * max of 128 means no diagonal scored above it and the true best - * may be lower; let the full MSV filter compute it. */ - if (xE <= 128) return eslENORESULT; - - if (xE >= 255 - om->bias_b) { - *ret_sc = eslINFINITY; - - if (om->base_b - om->tjb_b - om->tbm_b < 128) - return eslENORESULT; - - return eslERANGE; - } - - xE += om->base_b - om->tjb_b - om->tbm_b; - xE -= 128; - - if (xE >= 255 - om->bias_b) { - *ret_sc = eslINFINITY; - return eslERANGE; - } - - xJ = xE - om->tec_b; + return p7_SSVFilter_FromXE(get_xE_sse(dsq, L, om), om, ret_sc); +} - if (xJ > om->base_b) return eslENORESULT; - *ret_sc = ((float)(xJ - om->tjb_b) - (float) om->base_b); - *ret_sc /= om->scale_b; - *ret_sc -= 3.0f; +/* Function: p7_SSVFilter_OrfBlock_sse() + * Synopsis: SSV maximum of every ORF in a block, SSE path. + * + * Purpose: See p7_SSVFilter_OrfBlock(). is what get_xE_sse() gives + * for ORF . + */ +void +p7_SSVFilter_OrfBlock_sse(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE) +{ + int i; - return eslOK; + for (i = 0; i < n; i++) xE[i] = get_xE_sse(orf[i].dsq, (int) orf[i].n, om); } #endif /* eslENABLE_SSE */ diff --git a/src/impl_sse/impl_sse.h b/src/impl_sse/impl_sse.h index abbbe9cc..133b58ae 100644 --- a/src/impl_sse/impl_sse.h +++ b/src/impl_sse/impl_sse.h @@ -502,6 +502,11 @@ extern void p7_oprofile_DestroyBlock(P7_OM_BLOCK *block); /* ssvfilter.c */ extern int p7_SSVFilter (const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc); +/* SSV over a block of ORFs: the pipeline can take each ORF's SSV maximum from + * one call, and turn it into the SSV filter's status and score itself */ +#define p7_SSV_ORFBLOCK 1 +extern int p7_SSVFilter_OrfBlock(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE); +extern int p7_SSVFilter_FromXE(int xE, const P7_OPROFILE *om, float *ret_sc); /* msvfilter.c */ extern int p7_MSVFilter (const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, P7_OMX *ox, float *ret_sc); diff --git a/src/impl_sse/ssvfilter.c b/src/impl_sse/ssvfilter.c index 4fb7f24b..2a59a6b8 100644 --- a/src/impl_sse/ssvfilter.c +++ b/src/impl_sse/ssvfilter.c @@ -872,11 +872,19 @@ get_xE(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om) } +/* Function: p7_SSVFilter_FromXE() + * Synopsis: The SSV filter's result, from the SSV matrix's highest cell. + * + * Purpose: is the highest cell of the SSV matrix of a sequence + * against , as get_xE() gives it. Return what p7_SSVFilter() + * returns for that sequence, and its score in <*ret_sc>; must + * be configured for the sequence's length. + */ int -p7_SSVFilter(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) +p7_SSVFilter_FromXE(int xE_in, const P7_OPROFILE *om, float *ret_sc) { /* Use 16 bit values to avoid overflow due to moved baseline */ - uint16_t xE; + uint16_t xE = (uint16_t) xE_in; uint16_t xJ; if (om->tjb_b + om->tbm_b + om->tec_b + om->bias_b >= 127) { @@ -885,8 +893,6 @@ p7_SSVFilter(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) return eslENORESULT; } - xE = get_xE(dsq, L, om); - /* Saturation floors every diagonal at the begin score (128), so a * max of 128 means no diagonal scored above it and the true best * may be lower; let the full MSV filter compute it. */ @@ -930,3 +936,32 @@ p7_SSVFilter(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) } +int +p7_SSVFilter(const ESL_DSQ *dsq, int L, const P7_OPROFILE *om, float *ret_sc) +{ + if (om->tjb_b + om->tbm_b + om->tec_b + om->bias_b >= 127) return eslENORESULT; /* see p7_SSVFilter_FromXE() */ + + return p7_SSVFilter_FromXE(get_xE(dsq, L, om), om, ret_sc); +} + + +/* Function: p7_SSVFilter_OrfBlock() + * Synopsis: SSV maximum of every ORF in a block. + * + * Purpose: For each of the ORFs in , store in the highest + * cell of its SSV matrix against . It does not depend on the + * length is configured for; p7_SSVFilter_FromXE() makes the + * filter's status and score from it. + * + * Returns: . + */ +int +p7_SSVFilter_OrfBlock(const P7_OPROFILE *om, const ESL_ORF *orf, int n, uint8_t *xE) +{ + int i; + + for (i = 0; i < n; i++) xE[i] = get_xE(orf[i].dsq, (int) orf[i].n, om); + return eslOK; +} + +