开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
blast_options.c
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：53k
源码类别：

生物技术

开发平台：

C/C++

blast_options.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_options.c,v $
  4.  * PRODUCTION Revision 1000.6  2004/06/01 18:07:31  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.111
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: blast_options.c,v 1000.6 2004/06/01 18:07:31 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's offical duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  */
  36. /** @file blast_options.c
  37.  * @todo FIXME needs file description & doxygen comments
  38.  */
  40. static char const rcsid[] = 
  41.     "$Id: blast_options.c,v 1000.6 2004/06/01 18:07:31 gouriano Exp $";
  43. #include <algo/blast/core/blast_options.h>
  44. #include <algo/blast/core/blast_filter.h>
  45. #include <algo/blast/core/blast_encoding.h>
  47. QuerySetUpOptions*
  48. BlastQuerySetUpOptionsFree(QuerySetUpOptions* options)
  50. {
  51.    sfree(options->filter_string);
  53.    sfree(options);
  54.    return NULL;
  55. }
  57. Int2
  58. BlastQuerySetUpOptionsNew(QuerySetUpOptions* *options)
  59. {
  60.    *options = (QuerySetUpOptions*) calloc(1, sizeof(QuerySetUpOptions));
  61.    
  62.    if (*options == NULL)
  63.       return 1;
  64.    
  65.    (*options)->genetic_code = BLAST_GENETIC_CODE;
  67.    return 0;
  68. }
  70. Int2 BLAST_FillQuerySetUpOptions(QuerySetUpOptions* options,
  71.         Uint1 program, const char *filter_string, Uint1 strand_option)
  72. {
  73.    if (options == NULL)
  74.       return 1;
  75.    
  76.    if (strand_option && 
  77.        (program == blast_type_blastn || program == blast_type_blastx ||
  78.        program == blast_type_tblastx)) {
  79.       options->strand_option = strand_option;
  80.    }
  82.    /* "L" indicates low-complexity (seg for proteins, 
  83.       dust for nucleotides). */
  84.    if (!filter_string || !strcasecmp(filter_string, "T")) {
  85.       if (program == blast_type_blastn)
  86.          options->filter_string = strdup("D");
  87.       else
  88.          options->filter_string = strdup("S");
  89.    } else {
  90.       options->filter_string = strdup(filter_string); 
  91.    }
  93.    return 0;
  94. }
  96. BlastInitialWordOptions*
  97. BlastInitialWordOptionsFree(BlastInitialWordOptions* options)
  99. {
  101. sfree(options);
  103. return NULL;
  104. }
  107. Int2
  108. BlastInitialWordOptionsNew(Uint1 program, 
  109.    BlastInitialWordOptions* *options)
  110. {
  111.    *options = 
  112.       (BlastInitialWordOptions*) calloc(1, sizeof(BlastInitialWordOptions));
  113.    if (*options == NULL)
  114.       return 1;
  116.    if (program != blast_type_blastn) { /* protein-protein options. */
  117.       (*options)->window_size = BLAST_WINDOW_SIZE_PROT;
  118.       (*options)->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_PROT;
  119.    } else {
  120.       (*options)->window_size = BLAST_WINDOW_SIZE_NUCL;
  121.    }
  123.    return 0;
  124. }
  126. Int2
  127. BLAST_FillInitialWordOptions(BlastInitialWordOptions* options, 
  128.    Uint1 program, Boolean greedy, Int4 window_size, 
  129.    Boolean variable_wordsize, Boolean ag_blast, Boolean mb_lookup,
  130.    double xdrop_ungapped)
  131. {
  132.    if (!options)
  133.       return 1;
  135.    /* Ungapped extension is performed in all cases except when greedy
  136.       gapped extension is used */
  137.    if (program != blast_type_blastn) {
  138.       options->ungapped_extension = TRUE;
  139.       options->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_PROT;
  140.    } else if (!greedy) {
  141.       options->ungapped_extension = TRUE;
  142.       options->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_NUCL;
  143.    } else {
  144.       options->ungapped_extension = FALSE;
  145.    }
  147.    if (window_size != 0)
  148.       options->window_size = window_size;
  149.    if (xdrop_ungapped != 0)
  150.       options->x_dropoff = xdrop_ungapped;
  152.    options->variable_wordsize = variable_wordsize;
  154.    if (ag_blast) {
  155.       options->extension_method = eRightAndLeft;
  156.    } else {
  157.       options->extension_method = eRight;
  158.       if (mb_lookup)
  159.          options->container_type = eMbStacks;
  160.       else
  161.          options->container_type = eDiagArray;
  162.    }
  164.    return 0;
  165. }
  169. BlastInitialWordParameters*
  170. BlastInitialWordParametersFree(BlastInitialWordParameters* parameters)
  172. {
  173. sfree(parameters);
  174. return NULL;
  176. }
  178. /** Compute the default cutoff expect value for ungapped extensions
  179.  * @param program The blast program type
  180.  * @return The default per-program expect value
  181.  */
  182. static double GetUngappedCutoff(Uint1 program)
  183. {
  184.    switch(program) {
  185.    case blast_type_blastn:
  186.       return UNGAPPED_CUTOFF_E_BLASTN;
  187.    case blast_type_blastp: 
  188.    case blast_type_rpsblast: 
  189.       return UNGAPPED_CUTOFF_E_BLASTP;
  190.    case blast_type_blastx: 
  191.       return UNGAPPED_CUTOFF_E_BLASTX;
  192.    case blast_type_tblastn:
  193.    case blast_type_rpstblastn:
  194.       return UNGAPPED_CUTOFF_E_TBLASTN;
  195.    case blast_type_tblastx:
  196.       return UNGAPPED_CUTOFF_E_TBLASTX;
  197.    }
  198.    return 0;
  199. }
  201. Int2
  202. BlastInitialWordParametersNew(Uint1 program_number, 
  203.    const BlastInitialWordOptions* word_options, 
  204.    const BlastHitSavingParameters* hit_params, 
  205.    const BlastExtensionParameters* ext_params, BlastScoreBlk* sbp, 
  206.    BlastQueryInfo* query_info, 
  207.    Uint4 subject_length,
  208.    BlastInitialWordParameters* *parameters)
  209. {
  210.    Int2 status = 0;
  212.    /* If parameters pointer is NULL, there is nothing to fill, 
  213.       so don't do anything */
  214.    if (!parameters)
  215.       return 0;
  217.    ASSERT(sbp);
  218.    ASSERT(word_options);
  220.    *parameters = (BlastInitialWordParameters*) 
  221.       calloc(1, sizeof(BlastInitialWordParameters));
  223.    (*parameters)->options = (BlastInitialWordOptions *) word_options;
  225.    ASSERT(sbp->kbp_std[query_info->first_context]);
  227.    (*parameters)->x_dropoff = (Int4)
  228.       ceil(sbp->scale_factor * word_options->x_dropoff * NCBIMATH_LN2/
  229.            sbp->kbp_std[query_info->first_context]->Lambda);
  231.    status = BlastInitialWordParametersUpdate(program_number,
  232.                hit_params, ext_params, sbp, query_info,
  233.                subject_length, *parameters);
  235.    return status;
  236. }
  238. Int2
  239. BlastInitialWordParametersUpdate(Uint1 program_number, 
  240.    const BlastHitSavingParameters* hit_params, 
  241.    const BlastExtensionParameters* ext_params, BlastScoreBlk* sbp, 
  242.    BlastQueryInfo* query_info, Uint4 subj_length,
  243.    BlastInitialWordParameters* parameters)
  244. {
  245.    Blast_KarlinBlk* kbp;
  246.    Boolean gapped_calculation = TRUE;
  248.    ASSERT(sbp);
  249.    ASSERT(hit_params);
  250.    ASSERT(ext_params);
  251.    ASSERT(query_info);
  253.    /* kbp_gap is only non-NULL for gapped searches! */
  254.    if (sbp->kbp_gap) {
  255.       kbp = sbp->kbp_gap[query_info->first_context];
  256.    } else {
  257.       kbp = sbp->kbp_std[query_info->first_context];
  258.       gapped_calculation = FALSE;
  259.    }
  261.    ASSERT(kbp);
  262.    /* For non-blastn programs cutoff score should not be larger than 
  263.       gap trigger. */
  264.    if (gapped_calculation && program_number != blast_type_blastn) {
  265.       parameters->cutoff_score = 
  266.          MIN(ext_params->gap_trigger, hit_params->cutoff_score);
  267.    } else {
  268.       Int4 s2 = 0;
  269.       double e2;
  270.       Int4 qlen;
  271.       /* Calculate score cutoff corresponding to a fixed e-value (0.05);
  272.          If it is smaller, then use this one */
  273.       qlen = query_info->context_offsets[query_info->last_context+1] - 1;
  274.       
  275.       e2 = GetUngappedCutoff(program_number);
  277.       BLAST_Cutoffs(&s2, &e2, kbp, MIN(subj_length, (Uint4) qlen)*subj_length, TRUE, 
  278.                     hit_params->gap_decay_rate);
  279.       s2 *= (Int4)sbp->scale_factor;
  280.       parameters->cutoff_score = MIN(hit_params->cutoff_score, s2);
  282.       if (parameters->x_dropoff != 0 && parameters->cutoff_score != 0) {
  283.          parameters->x_dropoff = 
  284.             MIN(parameters->x_dropoff, parameters->cutoff_score);
  285.       } else if (parameters->cutoff_score != 0) {
  286.          parameters->x_dropoff = parameters->cutoff_score;
  287.       }
  288.    }
  290.    return 0;
  291. }
  293. BlastExtensionOptions*
  294. BlastExtensionOptionsFree(BlastExtensionOptions* options)
  296. {
  298. sfree(options);
  300. return NULL;
  301. }
  305. Int2
  306. BlastExtensionOptionsNew(Uint1 program, BlastExtensionOptions* *options)
  308. {
  309. *options = (BlastExtensionOptions*) 
  310.            calloc(1, sizeof(BlastExtensionOptions));
  312. if (*options == NULL)
  313. return 1;
  315. if (program != blast_type_blastn) /* protein-protein options. */
  316. {
  317. (*options)->gap_x_dropoff = BLAST_GAP_X_DROPOFF_PROT;
  318. (*options)->gap_x_dropoff_final = 
  319.                    BLAST_GAP_X_DROPOFF_FINAL_PROT;
  320. (*options)->gap_trigger = BLAST_GAP_TRIGGER_PROT;
  321. (*options)->ePrelimGapExt = eDynProgExt;
  322. }
  323. else
  324. {
  325. (*options)->gap_trigger = BLAST_GAP_TRIGGER_NUCL;
  326. }
  328. return 0;
  329. }
  331. Int2
  332. BLAST_FillExtensionOptions(BlastExtensionOptions* options, 
  333.    Uint1 program, Boolean greedy, double x_dropoff, double x_dropoff_final)
  334. {
  335.    if (!options)
  336.       return 1;
  338.    if (program == blast_type_blastn) {
  339.       if (greedy) {
  340.          options->gap_x_dropoff = BLAST_GAP_X_DROPOFF_GREEDY;
  341.          options->ePrelimGapExt = eGreedyWithTracebackExt;
  342.          options->eTbackExt = eGreedyTbck;
  343.       } else {
  344.          options->gap_x_dropoff = BLAST_GAP_X_DROPOFF_NUCL;
  345.          options->gap_x_dropoff_final = BLAST_GAP_X_DROPOFF_FINAL_NUCL;
  346.          options->ePrelimGapExt = eDynProgExt;
  347.          options->eTbackExt = eDynProgTbck;
  348.       }
  349.    }
  351.    if (x_dropoff)
  352.       options->gap_x_dropoff = x_dropoff;
  353.    if (x_dropoff_final) {
  354.       options->gap_x_dropoff_final = x_dropoff_final;
  355.    } else {
  356.       /* Final X-dropoff can't be smaller than preliminary X-dropoff */
  357.       options->gap_x_dropoff_final = 
  358.          MAX(options->gap_x_dropoff_final, x_dropoff);
  359.    }
  361.    return 0;
  363. }
  365. Int2 
  366. BlastExtensionOptionsValidate(Uint1 program_number, 
  367.    const BlastExtensionOptions* options, Blast_Message* *blast_msg)
  369. {
  370. if (options == NULL)
  371. return 1;
  373. if (program_number != blast_type_blastn)
  374. {
  375. if (options->ePrelimGapExt == eGreedyWithTracebackExt || 
  376.              options->ePrelimGapExt == eGreedyExt)
  377. {
  378. Int4 code=2;
  379. Int4 subcode=1;
  380. Blast_MessageWrite(blast_msg, BLAST_SEV_WARNING, code, subcode, 
  381.                             "Greedy extension only supported for BLASTN");
  382. return (Int2) code;
  383. }
  384. }
  386. return 0;
  387. }
  389. Int2 BlastExtensionParametersNew(Uint1 program_number, 
  390.         const BlastExtensionOptions* options, BlastScoreBlk* sbp,
  391.         BlastQueryInfo* query_info, BlastExtensionParameters* *parameters)
  392. {
  393.    Blast_KarlinBlk* kbp,* kbp_gap;
  394.    BlastExtensionParameters* params;
  396.    /* If parameters pointer is NULL, there is nothing to fill, 
  397.       so don't do anything */
  398.    if (!parameters)
  399.       return 0;
  401.    if (sbp->kbp) {
  402.       kbp = sbp->kbp[query_info->first_context];
  403.    } else {
  404.       /* The Karlin block is not found, can't do any calculations */
  405.       *parameters = NULL;
  406.       return -1;
  407.    }
  409.    if (sbp->kbp_gap) {
  410.       kbp_gap = sbp->kbp_gap[query_info->first_context];
  411.    } else {
  412.       kbp_gap = kbp;
  413.    }
  414.    
  415.    *parameters = params = (BlastExtensionParameters*) 
  416.       malloc(sizeof(BlastExtensionParameters));
  418.    params->options = (BlastExtensionOptions *) options;
  419.    params->gap_x_dropoff = 
  420.       (Int4) (options->gap_x_dropoff*NCBIMATH_LN2 / kbp_gap->Lambda);
  421.    params->gap_x_dropoff_final = (Int4) 
  422.       (options->gap_x_dropoff_final*NCBIMATH_LN2 / kbp_gap->Lambda);
  424.    params->gap_trigger = (Int4)
  425.       ((options->gap_trigger*NCBIMATH_LN2 + kbp->logK) / kbp->Lambda);
  427.    if (sbp->scale_factor > 1.0) {
  428.       params->gap_trigger *= (Int4)sbp->scale_factor;
  429.       params->gap_x_dropoff *= (Int4)sbp->scale_factor;
  430.       params->gap_x_dropoff_final *= (Int4)sbp->scale_factor;
  431.    }
  432.    return 0;
  433. }
  435. BlastExtensionParameters*
  436. BlastExtensionParametersFree(BlastExtensionParameters* parameters)
  437. {
  438.   sfree(parameters);
  439.   return NULL;
  440. }
  443. BlastScoringOptions*
  444. BlastScoringOptionsFree(BlastScoringOptions* options)
  446. {
  447. if (options == NULL)
  448. return NULL;
  450. sfree(options->matrix);
  451.    sfree(options->matrix_path);
  452. sfree(options);
  454. return NULL;
  455. }
  457. Int2 
  458. BlastScoringOptionsNew(Uint1 program_number, BlastScoringOptions* *options)
  459. {
  460.    *options = (BlastScoringOptions*) calloc(1, sizeof(BlastScoringOptions));
  462.    if (*options == NULL)
  463.       return 1;
  464.    
  465.    if (program_number != blast_type_blastn) { /* protein-protein options. */
  466.       (*options)->shift_pen = INT2_MAX;
  467.       (*options)->is_ooframe = FALSE;
  468.       (*options)->gap_open = BLAST_GAP_OPEN_PROT;
  469.       (*options)->gap_extend = BLAST_GAP_EXTN_PROT;
  470.    } else { /* nucleotide-nucleotide options. */
  471.       (*options)->penalty = BLAST_PENALTY;
  472.       (*options)->reward = BLAST_REWARD;
  473.    }
  474.    (*options)->decline_align = INT2_MAX;
  475.    (*options)->gapped_calculation = TRUE;
  476.    
  477.    return 0;
  478. }
  480. Int2 
  481. BLAST_FillScoringOptions(BlastScoringOptions* options, 
  482.    Uint1 program_number, Boolean greedy_extension, Int4 penalty, Int4 reward, 
  483.    const char *matrix, Int4 gap_open, Int4 gap_extend)
  484. {
  485.    if (!options)
  486.       return 1;
  488.    if (program_number != blast_type_blastn) { /* protein-protein options. */
  489.       if (matrix) {
  490.          unsigned int i;
  491.          options->matrix = strdup(matrix);
  492.          /* Make it all upper case */
  493.          for (i=0; i<strlen(options->matrix); ++i)
  494.             options->matrix[i] = toupper(options->matrix[i]);
  495.       } else {
  496.          options->matrix = strdup("BLOSUM62");
  497.       }
  498.    } else { /* nucleotide-nucleotide options. */
  499.       if (penalty)
  500.          options->penalty = penalty;
  501.       if (reward)
  502.          options->reward = reward;
  504.       if (greedy_extension) {
  505.          options->gap_open = BLAST_GAP_OPEN_MEGABLAST;
  506.          options->gap_extend = BLAST_GAP_EXTN_MEGABLAST;
  507.       } else {
  508.          options->gap_open = BLAST_GAP_OPEN_NUCL;
  509.          options->gap_extend = BLAST_GAP_EXTN_NUCL;
  510.       }
  511.    }
  512.    if (gap_open)
  513.       options->gap_open = gap_open;
  514.    if (gap_extend)
  515.       options->gap_extend = gap_extend;
  517.    return 0;
  518. }
  520. Int2 
  521. BlastScoringOptionsValidate(Uint1 program_number, 
  522.    const BlastScoringOptions* options, Blast_Message* *blast_msg)
  524. {
  525. if (options == NULL)
  526. return 1;
  528.    if (program_number == blast_type_tblastx && 
  529.               options->gapped_calculation)
  530.    {
  531. Int4 code=2;
  532. Int4 subcode=1;
  533.       Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  534.          "Gapped search is not allowed for tblastx");
  535. return (Int2) code;
  536.    }
  538. if (program_number == blast_type_blastn)
  539. {
  540. if (options->penalty >= 0)
  541. {
  542. Int4 code=2;
  543. Int4 subcode=1;
  544. Blast_MessageWrite(blast_msg, BLAST_SEV_WARNING, code, subcode, 
  545.                             "BLASTN penalty must be negative");
  546. return (Int2) code;
  547. }
  548.                 if (options->gap_open > 0 && options->gap_extend == 0) 
  549.                 {
  550.                         Int4 code=2;
  551.                         Int4 subcode=1;
  552.                         Blast_MessageWrite(blast_msg, BLAST_SEV_WARNING, 
  553.                            code, subcode, 
  554.                            "BLASTN gap extension penalty cannot be 0");
  555.                         return (Int2) code;
  556.                 }
  557. }
  558. else
  559. {
  560. Int2 status=0;
  562. if ((status=Blast_KarlinkGapBlkFill(NULL, options->gap_open, 
  563.                      options->gap_extend, options->decline_align, 
  564.                      options->matrix)) != 0)
  565. {
  566. if (status == 1)
  567. {
  568. char* buffer;
  569. Int4 code=2;
  570. Int4 subcode=1;
  572. buffer = BLAST_PrintMatrixMessage(options->matrix); 
  573.             Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR,
  574.                                code, subcode, buffer);
  575. sfree(buffer);
  576. return (Int2) code;
  578. }
  579. else if (status == 2)
  580. {
  581. char* buffer;
  582. Int4 code=2;
  583. Int4 subcode=1;
  585. buffer = BLAST_PrintAllowedValues(options->matrix, 
  586.                         options->gap_open, options->gap_extend, 
  587.                         options->decline_align); 
  588.             Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  589.                                buffer);
  590. sfree(buffer);
  591. return (Int2) code;
  592. }
  593. }
  595. }
  597. if (program_number != blast_type_blastx && 
  598.        program_number != blast_type_tblastn && options->is_ooframe)
  599. {
  600.       Int4 code=2;
  601.       Int4 subcode=1;
  602.       Blast_MessageWrite(blast_msg, BLAST_SEV_WARNING, code, subcode, 
  603.          "Out-of-frame only permitted for blastx and tblastn");
  604.       return (Int2) code;
  605. }
  607. return 0;
  609. }
  611. Int2 
  612. BlastScoringOptionsDup(BlastScoringOptions* *new_opt, const BlastScoringOptions* old_opt)
  613. {
  614.     if (old_opt == NULL || new_opt == NULL)
  615.        return -1;
  617.     *new_opt = (BlastScoringOptions*) BlastMemDup(old_opt, sizeof(BlastScoringOptions));
  618.     if (*new_opt == NULL)
  619.        return -1;
  621.     if (old_opt->matrix)
  622.        (*new_opt)->matrix = strdup(old_opt->matrix);
  624.     if (old_opt->matrix_path)
  625.        (*new_opt)->matrix_path = strdup(old_opt->matrix_path);
  627.     return 0;
  628. }
  630. BlastScoringParameters*
  631. BlastScoringParametersFree(BlastScoringParameters* parameters)
  632. {
  633. sfree(parameters);
  634. return NULL;
  635. }
  638. Int2
  639. BlastScoringParametersNew(const BlastScoringOptions* score_options, 
  640.                           BlastScoreBlk* sbp, 
  641.                           BlastScoringParameters* *parameters)
  642. {
  643.    BlastScoringParameters *params;
  644.    double scale_factor;
  646.    if (score_options == NULL)
  647.       return 1;
  649.    *parameters = params = (BlastScoringParameters*) 
  650.                         calloc(1, sizeof(BlastScoringParameters));
  651.    if (params == NULL)
  652.       return 2;
  654.    params->options = (BlastScoringOptions *)score_options;
  655.    scale_factor = sbp->scale_factor;
  656.    params->scale_factor = scale_factor;
  657.    params->reward = score_options->reward;
  658.    params->penalty = score_options->penalty;
  659.    params->gap_open = score_options->gap_open * (Int4)scale_factor;
  660.    params->gap_extend = score_options->gap_extend * (Int4)scale_factor;
  661.    params->decline_align = score_options->decline_align * (Int4)scale_factor;
  662.    params->shift_pen = score_options->shift_pen * (Int4)scale_factor;
  663.    return 0;
  664. }
  667. BlastEffectiveLengthsOptions*
  668. BlastEffectiveLengthsOptionsFree(BlastEffectiveLengthsOptions* options)
  670. {
  671. sfree(options);
  673. return NULL;
  674. }
  677. Int2 
  678. BlastEffectiveLengthsOptionsNew(BlastEffectiveLengthsOptions* *options)
  680. {
  681.    *options = (BlastEffectiveLengthsOptions*)
  682.       calloc(1, sizeof(BlastEffectiveLengthsOptions));
  684.    if (*options == NULL)
  685.       return 1;
  686.    
  687.    return 0;
  688. }
  690. BlastEffectiveLengthsParameters*
  691. BlastEffectiveLengthsParametersFree(BlastEffectiveLengthsParameters* parameters)
  693. {
  694. sfree(parameters);
  696. return NULL;
  697. }
  699. Int2 
  700. BlastEffectiveLengthsParametersNew(const BlastEffectiveLengthsOptions* options, 
  701.                                Int8 db_length, Int4 num_seqs,
  702.                                BlastEffectiveLengthsParameters* *parameters)
  703. {
  704.    *parameters = (BlastEffectiveLengthsParameters*) 
  705.       calloc(1, sizeof(BlastEffectiveLengthsParameters));
  706.    (*parameters)->options = (BlastEffectiveLengthsOptions*) options;
  707.    (*parameters)->real_db_length = db_length;
  708.    (*parameters)->real_num_seqs = num_seqs;
  709.    return 0;
  710. }
  712. Int2 
  713. BLAST_FillEffectiveLengthsOptions(BlastEffectiveLengthsOptions* options, 
  714.    Int4 dbseq_num, Int8 db_length, Int8 searchsp_eff)
  715. {
  716.    if (!options)
  717.       return 1;
  719.    if (searchsp_eff) {
  720.       /* dbnum_seq and dblen are used to calculate effective search space, so 
  721.          if it is already set don't bother with those. */
  722.       options->searchsp_eff = searchsp_eff;
  723.       return 0;
  724.    }
  726.    options->dbseq_num = dbseq_num;
  727.    options->db_length = db_length;
  729.    return 0;
  730. }
  732. LookupTableOptions*
  733. LookupTableOptionsFree(LookupTableOptions* options)
  735. {
  736.       sfree(options->phi_pattern);
  737.    
  738. sfree(options);
  739. return NULL;
  740. }
  742. Int2 
  743. LookupTableOptionsNew(Uint1 program_number, LookupTableOptions* *options)
  744. {
  745.    *options = (LookupTableOptions*) calloc(1, sizeof(LookupTableOptions));
  746.    
  747.    if (*options == NULL)
  748.       return 1;
  749.    
  750.    if (program_number != blast_type_blastn) {
  751.       (*options)->word_size = BLAST_WORDSIZE_PROT;
  752.       (*options)->alphabet_size = BLASTAA_SIZE;
  753.       (*options)->lut_type = AA_LOOKUP_TABLE;
  754.       
  755.       if (program_number == blast_type_blastp ||
  756.           program_number == blast_type_rpsblast)
  757.          (*options)->threshold = BLAST_WORD_THRESHOLD_BLASTP;
  758.       else if (program_number == blast_type_blastx)
  759.          (*options)->threshold = BLAST_WORD_THRESHOLD_BLASTX;
  760.       else if (program_number == blast_type_tblastn ||
  761.                program_number == blast_type_rpstblastn)
  762.          (*options)->threshold = BLAST_WORD_THRESHOLD_TBLASTN;
  763.       else if (program_number == blast_type_tblastx)
  764.          (*options)->threshold = BLAST_WORD_THRESHOLD_TBLASTX;
  765.       
  766.    } else {
  767.       (*options)->alphabet_size = BLASTNA_SIZE;
  768.    }
  770.    return 0;
  771. }
  773. Int4 CalculateBestStride(Int4 word_size, Boolean var_words, Int4 lut_type)
  774. {
  775.    Int4 lut_width;
  776.    Int4 extra = 1;
  777.    Uint1 remainder;
  778.    Int4 stride;
  780.    if (lut_type == MB_LOOKUP_TABLE)
  781.       lut_width = 12;
  782.    else if (word_size >= 8)
  783.       lut_width = 8;
  784.    else
  785.       lut_width = 4;
  787.    remainder = word_size % COMPRESSION_RATIO;
  789.    if (var_words && (remainder == 0) )
  790.       extra = COMPRESSION_RATIO;
  792.    stride = word_size - lut_width + extra;
  794.    remainder = stride % 4;
  796.    if (stride > 8 || (stride > 4 && remainder == 1) ) 
  797.       stride -= remainder;
  798.    return stride;
  799. }
  801. Int2 
  802. BLAST_FillLookupTableOptions(LookupTableOptions* options, 
  803.    Uint1 program_number, Boolean is_megablast, Int4 threshold,
  804.    Int2 word_size, Boolean ag_blast, Boolean variable_wordsize,
  805.    Boolean use_pssm)
  806. {
  807.    if (!options)
  808.       return 1;
  810.    if (program_number == blast_type_blastn) {
  811.       if (is_megablast) {
  812.          options->word_size = BLAST_WORDSIZE_MEGABLAST;
  813.          options->lut_type = MB_LOOKUP_TABLE;
  814.          options->max_positions = INT4_MAX;
  815.       } else {
  816.          options->lut_type = NA_LOOKUP_TABLE;
  817.          options->word_size = BLAST_WORDSIZE_NUCL;
  818.       }
  819.    } else {
  820.       options->lut_type = AA_LOOKUP_TABLE;
  821.    }
  823.    /* if the supplied threshold is -1, disable neighboring words */
  824.    if (threshold == -1)
  825.       options->threshold = 0;
  827.    /* if the supplied threshold is > 0, use it */
  828.    if (threshold > 0)
  829.       options->threshold = threshold;
  831.    /* otherwise, use the default */
  833.    if (use_pssm)
  834.       options->use_pssm = use_pssm;
  835.    if (program_number == blast_type_rpsblast ||
  836.        program_number == blast_type_rpstblastn)
  837.       options->lut_type = RPS_LOOKUP_TABLE;
  838.    if (word_size)
  839.       options->word_size = word_size;
  840.    if (program_number == blast_type_blastn) {
  841.       if (!ag_blast) {
  842.          options->scan_step = COMPRESSION_RATIO;
  843.       } else {
  844.          options->scan_step = CalculateBestStride(options->word_size, variable_wordsize,
  845.                                                   options->lut_type);
  846.       }
  847.    }
  848.    return 0;
  849. }
  851. /** Validate options for the discontiguous word megablast
  852.  * Word size must be 11 or 12; template length 16, 18 or 21; 
  853.  * template type 0, 1 or 2.
  854.  * @param word_size Word size option [in]
  855.  * @param template_length Discontiguous template length [in]
  856.  * @param template_type Discontiguous template type [in]
  857.  * @return TRUE if options combination valid.
  858.  */
  859. static Boolean 
  860. DiscWordOptionsValidate(Int2 word_size, Uint1 template_length,
  861.                         Uint1 template_type)
  862. {
  863.    if (template_length == 0)
  864.       return TRUE;
  866.    if (word_size != 11 && word_size != 12)
  867.       return FALSE;
  868.    if (template_length != 16 && template_length != 18 && 
  869.        template_length != 21)
  870.       return FALSE;
  871.    if (template_type > 2)
  872.       return FALSE;
  874.    return TRUE;
  875. }
  877. Int2 
  878. LookupTableOptionsValidate(Uint1 program_number, 
  879.    const LookupTableOptions* options, Blast_Message* *blast_msg)
  881. {
  882.    Int4 code=2;
  883.    Int4 subcode=1;
  885. if (options == NULL)
  886. return 1;
  888.         if (program_number == blast_type_rpsblast ||
  889.             program_number == blast_type_rpstblastn)
  890.                 return 0;
  892. if (program_number != blast_type_blastn && options->threshold <= 0)
  893. {
  894. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  895.                          "Non-zero threshold required");
  896. return (Int2) code;
  897. }
  899. if (options->word_size <= 0)
  900. {
  901. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  902.                          "Word-size must be greater than zero");
  903. return (Int2) code;
  904. } else if (program_number == blast_type_blastn && options->word_size < 7)
  905. {
  906. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  907.                          "Word-size must be 7" 
  908.                          "or greater for nucleotide comparison");
  909. return (Int2) code;
  910. } else if (program_number != blast_type_blastn && options->word_size > 5)
  911. {
  912. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  913.                          "Word-size must be less"
  914.                          "than 6 for protein comparison");
  915. return (Int2) code;
  916. }
  919. if (program_number != blast_type_blastn && 
  920.        options->lut_type == MB_LOOKUP_TABLE)
  921. {
  922. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  923.                          "Megablast lookup table only supported with blastn");
  924. return (Int2) code;
  925. }
  927.    if (options->lut_type == MB_LOOKUP_TABLE && options->word_size < 12 && 
  928.        options->mb_template_length == 0) {
  929.       Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  930.                          "Word size must be 12 or greater with megablast"
  931.                          " lookup table");
  932.       return (Int2) code;
  933.    }
  935.    if (program_number == blast_type_blastn && 
  936.        options->mb_template_length > 0) {
  937.       if (!DiscWordOptionsValidate(options->word_size,
  938.             options->mb_template_length, options->mb_template_type)) {
  939.          Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  940.             "Invalid discontiguous template parameters");
  941.          return (Int2) code;
  942.       } else if (options->lut_type != MB_LOOKUP_TABLE) {
  943.          Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  944.             "Invalid lookup table type for discontiguous Mega BLAST");
  945.          return (Int2) code;
  946.       } else if (options->scan_step != 1 && 
  947.                  options->scan_step != COMPRESSION_RATIO) {
  948.          Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  949.             "Invalid scanning stride for discontiguous Mega BLAST");
  950.          return (Int2) code;
  951.       }
  952.    }
  954. return 0;
  955. }
  957. BlastHitSavingOptions*
  958. BlastHitSavingOptionsFree(BlastHitSavingOptions* options)
  960. {
  961.   sfree(options);
  962.   return NULL;
  963. }
  966. Int2 BlastHitSavingOptionsNew(Uint1 program_number, 
  967.         BlastHitSavingOptions* *options)
  968. {
  969.    *options = (BlastHitSavingOptions*) calloc(1, sizeof(BlastHitSavingOptions));
  970.    
  971.    if (*options == NULL)
  972.       return 1;
  974.    (*options)->hitlist_size = 500;
  975.    (*options)->prelim_hitlist_size = 500;
  976.    (*options)->expect_value = BLAST_EXPECT_VALUE;
  978.    /* other stuff?? */
  979.    
  980.    return 0;
  982. }
  984. Int2
  985. BLAST_FillHitSavingOptions(BlastHitSavingOptions* options, 
  986.                            double evalue, Int4 hitlist_size)
  987. {
  988.    if (!options)
  989.       return 1;
  991.    if (hitlist_size)
  992.       options->hitlist_size = options->prelim_hitlist_size = hitlist_size;
  993.    if (evalue)
  994.       options->expect_value = evalue;
  996.    return 0;
  998. }
  1000. Int2
  1001. BlastHitSavingOptionsValidate(Uint1 program_number,
  1002.    const BlastHitSavingOptions* options, Blast_Message* *blast_msg)
  1003. {
  1004. if (options == NULL)
  1005. return 1;
  1007. if (options->hitlist_size < 1 || options->prelim_hitlist_size < 1)
  1008. {
  1009. Int4 code=1;
  1010. Int4 subcode=1;
  1011. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  1012.                          "No hits are being saved");
  1013. return (Int2) code;
  1014. }
  1016. if (options->expect_value <= 0.0 && options->cutoff_score <= 0)
  1017. {
  1018. Int4 code=2;
  1019. Int4 subcode=1;
  1020. Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
  1021.          "expect value or cutoff score must be greater than zero");
  1022. return (Int2) code;
  1023. }
  1025. return 0;
  1026. }
  1029. BlastHitSavingParameters*
  1030. BlastHitSavingParametersFree(BlastHitSavingParameters* parmameters)
  1032. {
  1033.   sfree(parmameters);
  1034.   return NULL;
  1035. }
  1038. Int2
  1039. BlastHitSavingParametersNew(Uint1 program_number, 
  1040.    const BlastHitSavingOptions* options, 
  1041.    const BlastExtensionParameters* ext_params, 
  1042.    BlastScoreBlk* sbp, BlastQueryInfo* query_info, 
  1043.    BlastHitSavingParameters* *parameters)
  1044. {
  1045.    Boolean gapped_calculation = TRUE;
  1046.    Int2 status = 0;
  1047.    BlastHitSavingParameters* params;
  1049.    /* If parameters pointer is NULL, there is nothing to fill, 
  1050.       so don't do anything */
  1051.    if (!parameters)
  1052.       return 0;
  1054.    ASSERT(options);
  1055.    ASSERT(sbp);
  1056.    
  1057.    if (!sbp->kbp_gap)
  1058.       gapped_calculation = FALSE;
  1060.    /* If parameters have not yet been created, allocate and fill all
  1061.       parameters that are constant throughout the search */
  1062.    *parameters = params = (BlastHitSavingParameters*) 
  1063.       calloc(1, sizeof(BlastHitSavingParameters));
  1065.    if (params == NULL)
  1066.       return 1;
  1068.    params->options = (BlastHitSavingOptions *) options;
  1070.    /* If sum statistics use is forced by the options, 
  1071.       set it in the paramters */
  1072.    params->do_sum_stats = options->do_sum_stats;
  1073.    /* Sum statistics is used anyway for all ungapped searches and all 
  1074.       translated gapped searches (except RPS translated searches) */
  1075.    if (!gapped_calculation || 
  1076.        (program_number != blast_type_blastn && 
  1077.         program_number != blast_type_blastp &&
  1078.         program_number != blast_type_rpsblast &&
  1079.         program_number != blast_type_rpstblastn))
  1080.       params->do_sum_stats = TRUE;
  1081.    if (program_number == blast_type_blastn || !gapped_calculation) {
  1082.       params->gap_prob = BLAST_GAP_PROB;
  1083.       params->gap_decay_rate = BLAST_GAP_DECAY_RATE;
  1084.    } else {
  1085.       params->gap_prob = BLAST_GAP_PROB_GAPPED;
  1086.       params->gap_decay_rate = BLAST_GAP_DECAY_RATE_GAPPED;
  1087.    }
  1088.    params->gap_size = BLAST_GAP_SIZE;
  1089.    params->cutoff_big_gap = 0;
  1091.    status = BlastHitSavingParametersUpdate(program_number, 
  1092.                ext_params, sbp, query_info, params);
  1093.    return status;
  1094. }
  1096. Int2
  1097. BlastHitSavingParametersUpdate(Uint1 program_number, 
  1098.    const BlastExtensionParameters* ext_params, 
  1099.    BlastScoreBlk* sbp, BlastQueryInfo* query_info, 
  1100.    BlastHitSavingParameters* params)
  1101. {
  1102.    BlastHitSavingOptions* options;
  1103.    Blast_KarlinBlk* kbp;
  1104.    double evalue;
  1105.    double scale_factor = sbp->scale_factor;
  1107.    ASSERT(params);
  1108.    ASSERT(query_info);
  1109.    ASSERT(ext_params);
  1111.    options = params->options;
  1112.    evalue = options->expect_value;
  1114.    /* Scoring options are not available here, but we can determine whether
  1115.       this is a gapped or ungapped search by checking whether gapped
  1116.       Karlin blocks have been set. */
  1117.    if (sbp->kbp_gap) {
  1118.       kbp = sbp->kbp_gap[query_info->first_context];
  1119.    } else {
  1120.       kbp = sbp->kbp[query_info->first_context];
  1121.    }
  1123.    /* Calculate cutoffs based on the current effective lengths information */
  1124.    if (options->cutoff_score > 0) {
  1125.       params->cutoff_score = options->cutoff_score * (Int4) sbp->scale_factor;
  1126.    } else if (!options->phi_align) {
  1127.       Int4 context = query_info->first_context;
  1128.       Int8 searchsp = query_info->eff_searchsp_array[context];
  1130.       /* translated RPS searches must scale the search space down */
  1131.       if (program_number == blast_type_rpstblastn)
  1132.          searchsp = searchsp / NUM_FRAMES;
  1134.       params->cutoff_score = 0;
  1135.       BLAST_Cutoffs(&params->cutoff_score, &evalue, kbp, searchsp, FALSE, 0);
  1136.       params->cutoff_score *= (Int4) scale_factor;
  1138.       /* When sum statistics is used, all HSPs above the gap trigger 
  1139.          cutoff are saved until the sum statistics is applied to potentially
  1140.          link them with other HSPs and improve their e-values. 
  1141.          However this does not apply to the ungapped search! */
  1142.       if (params->do_sum_stats) {
  1143.          params->cutoff_score = 
  1144.             MIN(params->cutoff_score, ext_params->gap_trigger);
  1145.       }
  1146.    } else {
  1147.       params->cutoff_score = 0;
  1148.    }
  1149.    
  1150.    params->cutoff_small_gap = 
  1151.       MIN(params->cutoff_score, ext_params->gap_trigger);
  1152.       
  1153.    return 0;
  1154. }
  1156. Int2 PSIBlastOptionsNew(PSIBlastOptions** psi_options)
  1157. {
  1158.    PSIBlastOptions* options;
  1159.    if (!psi_options)
  1160.       return 0;
  1161.    options = (PSIBlastOptions*)calloc(1, sizeof(PSIBlastOptions));
  1162.    *psi_options = options;
  1163.    options->inclusion_ethresh = PSI_INCLUSION_ETHRESH;
  1164.    options->pseudo_count = PSI_PSEUDO_COUNT_CONST;
  1165.    
  1166.    return 0;
  1167. }
  1169. PSIBlastOptions* PSIBlastOptionsFree(PSIBlastOptions* psi_options)
  1170. {
  1171.    sfree(psi_options);
  1172.    return NULL;
  1173. }
  1175. Int2 BlastDatabaseOptionsNew(BlastDatabaseOptions** db_options)
  1176. {
  1177.    BlastDatabaseOptions* options = (BlastDatabaseOptions*)
  1178.       calloc(1, sizeof(BlastDatabaseOptions));
  1180.    options->genetic_code = BLAST_GENETIC_CODE;
  1181.    *db_options = options;
  1183.    return 0;
  1184. }
  1186. BlastDatabaseOptions* 
  1187. BlastDatabaseOptionsFree(BlastDatabaseOptions* db_options)
  1188. {
  1189.    sfree(db_options->gen_code_string);
  1190.    sfree(db_options);
  1191.    return NULL;
  1192. }
  1194. Int2 BLAST_InitDefaultOptions(Uint1 program_number,
  1195.    LookupTableOptions** lookup_options,
  1196.    QuerySetUpOptions** query_setup_options, 
  1197.    BlastInitialWordOptions** word_options,
  1198.    BlastExtensionOptions** ext_options,
  1199.    BlastHitSavingOptions** hit_options,
  1200.    BlastScoringOptions** score_options,
  1201.    BlastEffectiveLengthsOptions** eff_len_options,
  1202.    PSIBlastOptions** psi_options,
  1203.    BlastDatabaseOptions** db_options)
  1204. {
  1205.    Int2 status;
  1207.    if ((status = LookupTableOptionsNew(program_number, lookup_options)))
  1208.       return status;
  1210.    if ((status=BlastQuerySetUpOptionsNew(query_setup_options)))
  1211.       return status;
  1213.    if ((status=BlastInitialWordOptionsNew(program_number, word_options)))
  1214.       return status;
  1216.    if ((status = BlastExtensionOptionsNew(program_number, ext_options)))
  1217.       return status;
  1219.    if ((status=BlastHitSavingOptionsNew(program_number, hit_options)))
  1220.       return status;
  1222.    if ((status=BlastScoringOptionsNew(program_number, score_options)))
  1223.       return status;
  1225.    if ((status=BlastEffectiveLengthsOptionsNew(eff_len_options)))
  1226.       return status;
  1227.    
  1228.    if ((status=PSIBlastOptionsNew(psi_options)))
  1229.       return status;
  1231.    if ((status=BlastDatabaseOptionsNew(db_options)))
  1232.       return status;
  1234.    return 0;
  1236. }
  1238. Int2 BLAST_ValidateOptions(Uint1 program_number,
  1239.                            const BlastExtensionOptions* ext_options,
  1240.                            const BlastScoringOptions* score_options, 
  1241.                            const LookupTableOptions* lookup_options, 
  1242.                            const BlastHitSavingOptions* hit_options,
  1243.                            Blast_Message* *blast_msg)
  1244. {
  1245.    Int2 status = 0;
  1247.    if ((status = BlastExtensionOptionsValidate(program_number, ext_options,
  1248.                                                blast_msg)) != 0)
  1249.        return status;
  1250.    if ((status = BlastScoringOptionsValidate(program_number, score_options,
  1251.                                                blast_msg)) != 0)
  1252.        return status;
  1253.    if ((status = LookupTableOptionsValidate(program_number, 
  1254.                     lookup_options, blast_msg)) != 0)   
  1255.        return status;
  1256.    if ((status = BlastHitSavingOptionsValidate(program_number, hit_options,
  1257.                                                blast_msg)) != 0)
  1258.        return status;
  1259.    return status;
  1260. }
  1262. /** machine epsilon assumed by CalculateLinkHSPCutoffs */
  1263. #define MY_EPS 1.0e-9
  1265. void
  1266. CalculateLinkHSPCutoffs(Uint1 program, BlastQueryInfo* query_info, 
  1267.    BlastScoreBlk* sbp, BlastHitSavingParameters* hit_params, 
  1268.    BlastExtensionParameters* ext_params,
  1269.    Int8 db_length, Int4 subject_length)
  1270. {
  1271. double gap_prob, gap_decay_rate, x_variable, y_variable;
  1272. Blast_KarlinBlk* kbp;
  1273. Int4 expected_length, gap_size, query_length;
  1274. Int8 search_sp;
  1275.         Boolean translated_subject = (program == blast_type_tblastn || 
  1276.                                  program == blast_type_rpstblastn || 
  1277.                                  program == blast_type_tblastx);
  1279. /* Do this for the first context, should this be changed?? */
  1280. kbp = sbp->kbp[query_info->first_context];
  1281. gap_size = hit_params->gap_size;
  1282. gap_prob = hit_params->gap_prob;
  1283. gap_decay_rate = hit_params->gap_decay_rate;
  1284.    /* Use average query length */
  1285. query_length = query_info->context_offsets[query_info->last_context+1] /
  1286.       (query_info->last_context + 1);
  1288.    if (translated_subject) {
  1289.       /* Lengths in subsequent calculations should be on the protein scale */
  1290.       subject_length /= CODON_LENGTH;
  1291.       db_length /= CODON_LENGTH;
  1292.    }
  1294. /* Subtract off the expected score. */
  1295.    expected_length = BLAST_Nint(log(kbp->K*((double) query_length)*
  1296.                                     ((double) subject_length))/(kbp->H));
  1297.    query_length = query_length - expected_length;
  1299.    subject_length = subject_length - expected_length;
  1300.    query_length = MAX(query_length, 1);
  1301.    subject_length = MAX(subject_length, 1);
  1303.    /* If this is a database search, use database length, else the single 
  1304.       subject sequence length */
  1305.    if (db_length > subject_length) {
  1306.       y_variable = log((double) (db_length)/(double) subject_length)*(kbp->K)/
  1307.          (gap_decay_rate);
  1308.    } else {
  1309.       y_variable = log((double) (subject_length + expected_length)/
  1310.                        (double) subject_length)*(kbp->K)/(gap_decay_rate);
  1311.    }
  1312.    search_sp = ((Int8) query_length)* ((Int8) subject_length);
  1313.    x_variable = 0.25*y_variable*((double) search_sp);
  1315.    /* To use "small" gaps the query and subject must be "large" compared to
  1316.       the gap size. If small gaps may be used, then the cutoff values must be
  1317.       adjusted for the "bayesian" possibility that both large and small gaps 
  1318.       are being checked for. */
  1320.    if (search_sp > 8*gap_size*gap_size) {
  1321.       x_variable /= (1.0 - gap_prob + MY_EPS);
  1322.       hit_params->cutoff_big_gap = 
  1323.          (Int4) floor((log(x_variable)/kbp->Lambda)) + 1;
  1324.       x_variable = y_variable*(gap_size*gap_size);
  1325.       x_variable /= (gap_prob + MY_EPS);
  1326.       hit_params->cutoff_small_gap = 
  1327.          MAX(ext_params->gap_trigger, 
  1328.              (Int4) floor((log(x_variable)/kbp->Lambda)) + 1);
  1330.       hit_params->ignore_small_gaps = FALSE;
  1331.    } else {
  1332.       hit_params->cutoff_big_gap = 
  1333.          (Int4) floor((log(x_variable)/kbp->Lambda)) + 1;
  1334.       hit_params->cutoff_small_gap = hit_params->cutoff_big_gap;
  1335.       hit_params->ignore_small_gaps = TRUE;
  1336.    }
  1338.    hit_params->cutoff_big_gap *= (Int4)sbp->scale_factor;
  1339.    hit_params->cutoff_small_gap *= (Int4)sbp->scale_factor;
  1340. }
  1343. /*
  1344.  * ===========================================================================
  1345.  *
  1346.  * $Log: blast_options.c,v $
  1347.  * Revision 1000.6  2004/06/01 18:07:31  gouriano
  1348.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.111
  1349.  *
  1350.  * Revision 1.111  2004/05/26 16:04:54  papadopo
  1351.  * fix doxygen errors
  1352.  *
  1353.  * Revision 1.110  2004/05/24 17:26:21  camacho
  1354.  * Fix PC warning
  1355.  *
  1356.  * Revision 1.109  2004/05/20 16:29:30  madden
  1357.  * Make searchsp an Int8 consistent with rest of blast
  1358.  *
  1359.  * Revision 1.108  2004/05/19 14:52:02  camacho
  1360.  * 1. Added doxygen tags to enable doxygen processing of algo/blast/core
  1361.  * 2. Standardized copyright, CVS $Id string, $Log and rcsid formatting and i
  1362.  *    location
  1363.  * 3. Added use of @todo doxygen keyword
  1364.  *
  1365.  * Revision 1.107  2004/05/17 15:30:20  madden
  1366.  * Int algorithm_type replaced with enum EBlastPrelimGapExt, removed include for blast_gapalign.h
  1367.  *
  1368.  * Revision 1.106  2004/05/14 17:11:03  dondosha
  1369.  * Minor correction in setting X-dropoffs
  1370.  *
  1371.  * Revision 1.105  2004/05/14 13:14:15  camacho
  1372.  * Use correct definition for inclusion threshold
  1373.  *
  1374.  * Revision 1.104  2004/05/12 12:18:06  madden
  1375.  * Clean out PSIBlast options, add fields to ExtensionOptions to support smith-waterman and composition-based stats
  1376.  *
  1377.  * Revision 1.103  2004/05/10 14:27:23  madden
  1378.  * Correction to CalculateLinkHSPCutoffs to use gap_trigger in calculation of small cutoff
  1379.  *
  1380.  * Revision 1.102  2004/05/07 15:22:15  papadopo
  1381.  * 1. add functions to allocate and free BlastScoringParameters structures
  1382.  * 2. apply a scaling factor to all cutoffs generated in HitSavingParameters
  1383.  *    or ExtentionParameters structures
  1384.  *
  1385.  * Revision 1.101  2004/04/29 17:41:05  papadopo
  1386.  * Scale down the search space when calculating the S2 cutoff score for a translated RPS search
  1387.  *
  1388.  * Revision 1.100  2004/04/29 15:08:43  madden
  1389.  * Add BlastScoringOptionsDup
  1390.  *
  1391.  * Revision 1.99  2004/04/23 14:02:25  papadopo
  1392.  * ignore validation of LookupTableOptions if performing an RPS search
  1393.  *
  1394.  * Revision 1.98  2004/04/22 22:18:03  dondosha
  1395.  * Set lookup table type correctly in BLAST_FillLookupTableOptions - needed for C driver only
  1396.  *
  1397.  * Revision 1.97  2004/04/21 17:00:59  madden
  1398.  * Removed set but not read variable
  1399.  *
  1400.  * Revision 1.96  2004/04/19 12:58:44  madden
  1401.  * Changed BLAST_KarlinBlk to Blast_KarlinBlk to avoid conflict with blastkar.h structure, renamed some functions to start with Blast_Karlin, made Blast_KarlinBlkDestruct public
  1402.  *
  1403.  * Revision 1.95  2004/04/16 14:17:06  papadopo
  1404.  * add use of RPS-specific defines, remove RPS argument to FillLookupTableOptions
  1405.  *
  1406.  * Revision 1.94  2004/04/07 03:06:16  camacho
  1407.  * Added blast_encoding.[hc], refactoring blast_stat.[hc]
  1408.  *
  1409.  * Revision 1.93  2004/03/26 20:46:00  dondosha
  1410.  * Made gap_trigger parameter an integer, as in the old code
  1411.  *
  1412.  * Revision 1.92  2004/03/22 20:11:37  dondosha
  1413.  * Do not allow small gaps cutoff to be less than gap trigger
  1414.  *
  1415.  * Revision 1.91  2004/03/17 15:19:10  camacho
  1416.  * Add missing casts
  1417.  *
  1418.  * Revision 1.90  2004/03/11 23:58:10  dondosha
  1419.  * Set cutoff_score to 0 before calling BLAST_Cutoffs, so it knows what to calculate
  1420.  *
  1421.  * Revision 1.89  2004/03/11 20:41:49  camacho
  1422.  * Remove dead code
  1423.  *
  1424.  * Revision 1.88  2004/03/10 17:33:10  papadopo
  1425.  * Make a separate lookup table type for RPS blast
  1426.  *
  1427.  * Revision 1.87  2004/03/09 22:37:26  dondosha
  1428.  * Added const qualifiers to parameter arguments wherever relevant
  1429.  *
  1430.  * Revision 1.86  2004/03/09 18:46:24  dondosha
  1431.  * Corrected how cutoffs are calculated
  1432.  *
  1433.  * Revision 1.85  2004/03/04 21:07:48  papadopo
  1434.  * add RPS BLAST functionality
  1435.  *
  1436.  * Revision 1.84  2004/02/27 15:56:33  papadopo
  1437.  * Mike Gertz' modifications to unify handling of gapped Karlin blocks for protein and nucleotide searches. Also modified BLAST_MainSetUp to allocate gapped Karlin blocks last
  1438.  *
  1439.  * Revision 1.83  2004/02/24 17:57:14  dondosha
  1440.  * Added function to combine all options validation functions for the C engine
  1441.  *
  1442.  * Revision 1.82  2004/02/19 21:16:48  dondosha
  1443.  * Use enum type for severity argument in Blast_MessageWrite
  1444.  *
  1445.  * Revision 1.81  2004/02/17 22:10:30  dondosha
  1446.  * Set preliminary hitlist size in options initialization
  1447.  *
  1448.  * Revision 1.80  2004/02/07 15:48:30  ucko
  1449.  * PSIBlastOptionsNew: rearrange slightly so that declarations come first.
  1450.  *
  1451.  * Revision 1.79  2004/02/06 22:49:30  dondosha
  1452.  * Check for NULL pointer in PSIBlastOptionsNew
  1453.  *
  1454.  * Revision 1.78  2004/02/03 18:33:39  dondosha
  1455.  * Correction to previous change: word size can be 11 if discontiguous words
  1456.  *
  1457.  * Revision 1.77  2004/02/03 16:17:33  dondosha
  1458.  * Require word size to be >= 12 with megablast lookup table
  1459.  *
  1460.  * Revision 1.76  2004/02/02 18:49:32  dondosha
  1461.  * Fixes for minor compiler warnings
  1462.  *
  1463.  * Revision 1.75  2003/12/31 20:04:47  dondosha
  1464.  * Round best stride to a number divisible by 4 for all values except 6 and 7
  1465.  *
  1466.  * Revision 1.74  2003/12/31 16:04:37  coulouri
  1467.  * use -1 to disable protein neighboring words
  1468.  *
  1469.  * Revision 1.73  2003/12/08 16:03:05  coulouri
  1470.  * Propagate protein neighboring threshold even if it is zero
  1471.  *
  1472.  * Revision 1.72  2003/11/24 23:18:32  dondosha
  1473.  * Added gap_decay_rate argument to BLAST_Cutoffs; removed BLAST_Cutoffs_simple
  1474.  *
  1475.  * Revision 1.71  2003/11/12 18:17:46  dondosha
  1476.  * Correction in calculating scanning stride
  1477.  *
  1478.  * Revision 1.70  2003/11/04 23:22:47  dondosha
  1479.  * Do not calculate hit saving cutoff score for PHI BLAST
  1480.  *
  1481.  * Revision 1.69  2003/10/30 19:34:01  dondosha
  1482.  * Removed gapped_calculation from BlastHitSavingOptions structure
  1483.  *
  1484.  * Revision 1.68  2003/10/24 20:55:10  camacho
  1485.  * Rename GetDefaultStride
  1486.  *
  1487.  * Revision 1.67  2003/10/22 16:44:33  dondosha
  1488.  * Added function to calculate default stride value for AG method
  1489.  *
  1490.  * Revision 1.66  2003/10/21 22:15:34  camacho
  1491.  * Rearranging of C options structures, fix seed extension method
  1492.  *
  1493.  * Revision 1.65  2003/10/17 18:20:20  dondosha
  1494.  * Use separate variables for different initial word extension options
  1495.  *
  1496.  * Revision 1.64  2003/10/15 16:59:43  coulouri
  1497.  * type correctness fixes
  1498.  *
  1499.  * Revision 1.63  2003/10/07 17:26:11  dondosha
  1500.  * Lower case mask moved from options to the sequence block
  1501.  *
  1502.  * Revision 1.62  2003/10/02 22:08:34  dondosha
  1503.  * Corrections for one-strand translated searches
  1504.  *
  1505.  * Revision 1.61  2003/10/01 22:36:52  dondosha
  1506.  * Correction of setting of e2 in revision 1.57 was wrong
  1507.  *
  1508.  * Revision 1.60  2003/09/24 19:28:20  dondosha
  1509.  * Correction in setting extend word method: unset options that are set by default but overridden
  1510.  *
  1511.  * Revision 1.59  2003/09/12 17:26:01  dondosha
  1512.  * Added check that gap extension option cannot be 0 when gap open is not 0
  1513.  *
  1514.  * Revision 1.58  2003/09/10 19:48:08  dondosha
  1515.  * Removed dependency on mb_lookup.h
  1516.  *
  1517.  * Revision 1.57  2003/09/09 22:12:02  dondosha
  1518.  * Minor correction for ungapped cutoff calculation; added freeing of PHI pattern
  1519.  *
  1520.  * Revision 1.56  2003/09/08 12:55:57  madden
  1521.  * Allow use of PSSM to construct lookup table
  1522.  *
  1523.  * Revision 1.55  2003/08/27 15:05:37  camacho
  1524.  * Use symbolic name for alphabet sizes
  1525.  *
  1526.  * Revision 1.54  2003/08/26 21:53:33  madden
  1527.  * Protein alphabet is 26 chars, not 25
  1528.  *
  1529.  * Revision 1.53  2003/08/11 15:01:59  dondosha
  1530.  * Added algo/blast/core to all #included headers
  1531.  *
  1532.  * Revision 1.52  2003/08/01 17:26:19  dondosha
  1533.  * Use renamed versions of functions from local blastkar.h
  1534.  *
  1535.  * Revision 1.51  2003/07/31 17:45:17  dondosha
  1536.  * Made use of const qualifier consistent throughout the library
  1537.  *
  1538.  * Revision 1.50  2003/07/31 14:31:41  camacho
  1539.  * Replaced Char for char
  1540.  *
  1541.  * Revision 1.49  2003/07/31 14:19:28  camacho
  1542.  * Replaced FloatHi for double
  1543.  *
  1544.  * Revision 1.48  2003/07/31 00:32:37  camacho
  1545.  * Eliminated Ptr notation
  1546.  *
  1547.  * Revision 1.47  2003/07/30 22:06:25  dondosha
  1548.  * Convert matrix name to upper case when filling scoring options
  1549.  *
  1550.  * Revision 1.46  2003/07/30 19:39:14  camacho
  1551.  * Remove PNTRs
  1552.  *
  1553.  * Revision 1.45  2003/07/30 18:58:10  dondosha
  1554.  * Removed unused member matrixname from lookup table options
  1555.  *
  1556.  * Revision 1.44  2003/07/30 17:15:00  dondosha
  1557.  * Minor fixes for very strict compiler warnings
  1558.  *
  1559.  * Revision 1.43  2003/07/30 16:32:02  madden
  1560.  * Use ansi functions when possible
  1561.  *
  1562.  * Revision 1.42  2003/07/29 14:42:31  coulouri
  1563.  * use strdup() instead of StringSave()
  1564.  *
  1565.  * Revision 1.41  2003/07/28 19:04:15  camacho
  1566.  * Replaced all MemNews for calloc
  1567.  *
  1568.  * Revision 1.40  2003/07/25 21:12:28  coulouri
  1569.  * remove constructions of the form "return sfree();" and "a=sfree(a);"
  1570.  *
  1571.  * Revision 1.39  2003/07/25 17:25:43  coulouri
  1572.  * in progres:
  1573.  *  * use malloc/calloc/realloc instead of Malloc/Calloc/Realloc
  1574.  *  * add sfree() macro and __sfree() helper function to util.[ch]
  1575.  *  * use sfree() instead of MemFree()
  1576.  *
  1577.  * Revision 1.38  2003/07/23 17:31:10  camacho
  1578.  * BlastDatabaseParameters struct is deprecated
  1579.  *
  1580.  * Revision 1.37  2003/07/23 16:42:01  dondosha
  1581.  * Formatting options moved from blast_options.c to blast_format.c
  1582.  *
  1583.  * Revision 1.36  2003/07/22 20:26:16  dondosha
  1584.  * Initialize BlastDatabaseParameters structure outside engine
  1585.  *
  1586.  * Revision 1.35  2003/07/22 15:32:55  dondosha
  1587.  * Removed dependence on readdb API
  1588.  *
  1589.  * Revision 1.34  2003/07/21 20:31:47  dondosha
  1590.  * Added BlastDatabaseParameters structure with genetic code string
  1591.  *
  1592.  * Revision 1.33  2003/06/26 21:38:05  dondosha
  1593.  * Program number is removed from options structures, and passed explicitly as a parameter to functions that need it
  1594.  *
  1595.  * Revision 1.32  2003/06/26 20:24:06  camacho
  1596.  * Do not free options structure in BlastExtensionParametersFree
  1597.  *
  1598.  * Revision 1.31  2003/06/23 21:49:11  dondosha
  1599.  * Possibility of linking HSPs for tblastn activated
  1600.  *
  1601.  * Revision 1.30  2003/06/20 21:40:21  dondosha
  1602.  * Added parameters for linking HSPs
  1603.  *
  1604.  * Revision 1.29  2003/06/20 15:20:21  dondosha
  1605.  * Memory leak fixes
  1606.  *
  1607.  * Revision 1.28  2003/06/18 12:21:01  camacho
  1608.  * Added proper return value
  1609.  *
  1610.  * Revision 1.27  2003/06/17 20:42:43  camacho
  1611.  * Moved comments to header file, fixed includes
  1612.  *
  1613.  * Revision 1.26  2003/06/11 16:14:53  dondosha
  1614.  * Added initialization of PSI-BLAST and database options
  1615.  *
  1616.  * Revision 1.25  2003/06/09 20:13:17  dondosha
  1617.  * Minor type casting compiler warnings fixes
  1618.  *
  1619.  * Revision 1.24  2003/06/06 17:02:30  dondosha
  1620.  * Typo fix
  1621.  *
  1622.  * Revision 1.23  2003/06/04 20:16:51  coulouri
  1623.  * make prototypes and definitions agree
  1624.  *
  1625.  * Revision 1.22  2003/06/03 15:50:39  coulouri
  1626.  * correct function pointer argument
  1627.  *
  1628.  * Revision 1.21  2003/05/30 15:52:11  coulouri
  1629.  * various lint-induced cleanups
  1630.  *
  1631.  * Revision 1.20  2003/05/21 22:31:53  dondosha
  1632.  * Added forcing of ungapped search for tblastx to option validation
  1633.  *
  1634.  * Revision 1.19  2003/05/18 21:57:37  camacho
  1635.  * Use Uint1 for program name whenever possible
  1636.  *
  1637.  * Revision 1.18  2003/05/15 22:01:22  coulouri
  1638.  * add rcsid string to sources
  1639.  *
  1640.  * Revision 1.17  2003/05/13 20:41:48  dondosha
  1641.  * Correction in assigning of number of db sequences for 2 sequence case
  1642.  *
  1643.  * Revision 1.16  2003/05/13 15:11:34  dondosha
  1644.  * Changed some char * arguments to const char *
  1645.  *
  1646.  * Revision 1.15  2003/05/07 17:44:31  dondosha
  1647.  * Assign ungapped xdropoff default correctly for protein programs
  1648.  *
  1649.  * Revision 1.14  2003/05/06 20:29:57  dondosha
  1650.  * Fix in filling effective length options
  1651.  *
  1652.  * Revision 1.13  2003/05/06 14:34:51  dondosha
  1653.  * Fix in comment
  1654.  *
  1655.  * Revision 1.12  2003/05/01 16:56:30  dondosha
  1656.  * Fixed strict compiler warnings
  1657.  *
  1658.  * Revision 1.11  2003/05/01 15:33:39  dondosha
  1659.  * Reorganized the setup of BLAST search
  1660.  *
  1661.  * Revision 1.10  2003/04/24 14:27:35  dondosha
  1662.  * Correction for latest changes
  1663.  *
  1664.  * Revision 1.9  2003/04/23 20:04:49  dondosha
  1665.  * Added a function BLAST_InitAllDefaultOptions to initialize all various options structures with only default values
  1666.  *
  1667.  * Revision 1.8  2003/04/17 21:14:41  dondosha
  1668.  * Added cutoff score hit parameters that is calculated from e-value
  1669.  *
  1670.  * Revision 1.7  2003/04/16 22:25:37  dondosha
  1671.  * Correction to previous change
  1672.  *
  1673.  * Revision 1.6  2003/04/16 22:20:24  dondosha
  1674.  * Correction in calculation of cutoff score for ungapped extensions
  1675.  *
  1676.  * Revision 1.5  2003/04/11 22:35:48  dondosha
  1677.  * Minor corrections for blastn
  1678.  *
  1679.  * Revision 1.4  2003/04/03 22:57:50  dondosha
  1680.  * Uninitialized variable fix
  1681.  *
  1682.  * Revision 1.3  2003/04/02 17:20:41  dondosha
  1683.  * Added calculation of ungapped cutoff score in correct place
  1684.  *
  1685.  * Revision 1.2  2003/04/01 17:42:33  dondosha
  1686.  * Added arguments to BlastExtensionParametersNew
  1687.  *
  1688.  * Revision 1.1  2003/03/31 18:22:30  camacho
  1689.  * Moved from parent directory
  1690.  *
  1691.  * Revision 1.30  2003/03/28 23:12:34  dondosha
  1692.  * Added program argument to BlastFormattingOptionsNew
  1693.  *
  1694.  * Revision 1.29  2003/03/27 20:54:19  dondosha
  1695.  * Moved ungapped cutoff from hit options to word options
  1696.  *
  1697.  * Revision 1.28  2003/03/25 16:30:25  dondosha
  1698.  * Strict compiler warning fixes
  1699.  *
  1700.  * Revision 1.27  2003/03/24 20:39:17  dondosha
  1701.  * Added BlastExtensionParameters structure to hold raw gapped X-dropoff values
  1702.  *
  1703.  * Revision 1.26  2003/03/19 19:52:42  dondosha
  1704.  * 1. Added strand option argument to BlastQuerySetUpOptionsNew
  1705.  * 2. Added check of discontiguous template parameters in LookupTableOptionsValidate
  1706.  *
  1707.  * Revision 1.25  2003/03/14 19:08:53  dondosha
  1708.  * Added arguments to various OptionsNew functions, so all initialization can be done inside
  1709.  *
  1710.  * Revision 1.24  2003/03/12 17:03:41  dondosha
  1711.  * Set believe_query in formatting options to FALSE by default
  1712.  *
  1713.  * Revision 1.23  2003/03/11 20:40:32  dondosha
  1714.  * Correction in assigning gap_x_dropoff_final
  1715.  *
  1716.  * Revision 1.22  2003/03/10 16:44:42  dondosha
  1717.  * Added functions for initialization and freeing of formatting options structure
  1718.  *
  1719.  * Revision 1.21  2003/03/07 20:41:08  dondosha
  1720.  * Small corrections in option initialization functions
  1721.  *
  1722.  * Revision 1.20  2003/03/06 19:25:52  madden
  1723.  * Include blast_util.h
  1724.  *
  1725.  * Revision 1.19  2003/03/05 21:19:09  coulouri
  1726.  * set NA_LOOKUP_TABLE flag
  1727.  *
  1728.  * Revision 1.18  2003/03/05 20:58:50  dondosha
  1729.  * Corrections for handling effective search space for multiple queries
  1730.  *
  1731.  * Revision 1.17  2003/03/05 15:36:34  madden
  1732.  * Moved BlastNumber2Program and BlastProgram2Number from blast_options to blast_util
  1733.  *
  1734.  * Revision 1.16  2003/03/03 14:43:21  madden
  1735.  * Use BlastKarlinkGapBlkFill, PrintMatrixMessage, and PrintAllowedValuesMessage
  1736.  *
  1737.  * Revision 1.15  2003/02/26 15:42:50  madden
  1738.  * const charPtr becomes const char *, add BlastExtensionOptionsValidate
  1739.  *
  1740.  * Revision 1.14  2003/02/14 16:30:19  dondosha
  1741.  * Get rid of a compiler warning for type mismatch
  1742.  *
  1743.  * Revision 1.13  2003/02/13 21:42:25  madden
  1744.  * Added validation functions
  1745.  *
  1746.  * Revision 1.12  2003/02/04 13:14:36  dondosha
  1747.  * Changed the macro definitions for 
  1748.  *
  1749.  * Revision 1.11  2003/01/31 17:00:32  dondosha
  1750.  * Do not set the scan step in LookupTableOptionsNew
  1751.  *
  1752.  * Revision 1.10  2003/01/28 15:13:25  madden
  1753.  * Added functions and structures for parameters
  1754.  *
  1755.  * Revision 1.9  2003/01/22 20:49:31  dondosha
  1756.  * Set decline_align for blastn too
  1757.  *
  1758.  * Revision 1.8  2003/01/22 15:09:55  dondosha
  1759.  * Correction for default penalty assignment
  1760.  *
  1761.  * Revision 1.7  2003/01/17 22:10:45  madden
  1762.  * Added functions for BlastExtensionOptions, BlastInitialWordOptions as well as defines for default values
  1763.  *
  1764.  * Revision 1.6  2003/01/10 18:36:40  madden
  1765.  * Change call to BlastEffectiveLengthsOptionsNew
  1766.  *
  1767.  * Revision 1.5  2003/01/02 17:09:35  dondosha
  1768.  * Fill alphabet size when creating lookup table options structure
  1769.  *
  1770.  * Revision 1.4  2002/12/24 14:49:00  madden
  1771.  * Set defaults for LookupTableOptions for protein-protein searches
  1772.  *
  1773.  * Revision 1.3  2002/12/04 13:38:21  madden
  1774.  * Add function LookupTableOptionsNew
  1775.  *
  1776.  * Revision 1.2  2002/10/17 15:45:17  madden
  1777.  * Make BLOSUM62 default
  1778.  *
  1779.  * Revision 1.1  2002/10/07 21:05:12  madden
  1780.  * Sets default option values
  1781.  *
  1782.  * ===========================================================================
  1783.  */