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..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*/ @@ -1393,7 +1394,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 +1460,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 +1577,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/impl_avx/impl_avx.h b/src/impl_avx/impl_avx.h index 6dc3835e..ff0130c9 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 */ @@ -784,11 +785,18 @@ 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 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..07156a64 100644 --- a/src/impl_avx/ssvfilter.c +++ b/src/impl_avx/ssvfilter.c @@ -42,3 +42,82 @@ 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 there is no SIMD + * 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 +#ifdef eslENABLE_SSE + p7_SSVFilter_OrfBlock_sse(om, orf, n, xE); + return eslOK; +#else + return eslENORESULT; +#endif +} 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 */ 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; +} + + 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 80bb2ab4..97b41b99 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; @@ -105,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 @@ -282,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); } @@ -292,137 +295,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) { @@ -459,7 +331,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; @@ -475,7 +347,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++) { @@ -632,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); @@ -780,7 +654,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; @@ -1274,7 +1148,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; @@ -1443,7 +1317,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; @@ -1464,7 +1338,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; @@ -1473,9 +1347,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 +1385,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 +1468,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 @@ -1638,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. @@ -1695,7 +1601,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; @@ -1713,9 +1619,12 @@ 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; + 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 @@ -1727,16 +1636,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); @@ -1746,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]); @@ -1757,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); @@ -1766,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; @@ -1907,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); } @@ -1918,6 +1867,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); } @@ -1925,11 +1875,13 @@ 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); 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); } @@ -2032,6 +1984,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; @@ -2070,7 +2023,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); 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;