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

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

生物技术

开发平台：

C/C++

blast_setup.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_setup.c,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:07:47  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.87
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: blast_setup.c,v 1000.3 2004/06/01 18:07:47 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 official 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.  *
  35.  * Author: Tom Madden
  36.  *
  37.  */
  39. /** @file blast_setup.c
  40.  * Utilities initialize/setup BLAST.
  41.  */
  44. static char const rcsid[] =
  45.     "$Id: blast_setup.c,v 1000.3 2004/06/01 18:07:47 gouriano Exp $";
  47. #include <algo/blast/core/blast_setup.h>
  48. #include <algo/blast/core/blast_util.h>
  49. #include <algo/blast/core/blast_filter.h>
  50. #include <algo/blast/core/blast_encoding.h>
  52. /* See description in blast_setup.h */
  53. Int2
  54. BlastScoreBlkGappedFill(BlastScoreBlk * sbp,
  55.                         const BlastScoringOptions * scoring_options,
  56.                         Uint1 program, BlastQueryInfo * query_info)
  57. {
  58.     Int2 tmp_index;
  60.     if (sbp == NULL || scoring_options == NULL)
  61.         return 1;
  63.     /* At this stage query sequences are nucleotide only for blastn */
  65.     if (program == blast_type_blastn) {
  67.         /* for blastn, duplicate the ungapped Karlin 
  68.            structures for use in gapped alignments */
  70.         for (tmp_index = query_info->first_context;
  71.              tmp_index <= query_info->last_context; tmp_index++) {
  72.             if (sbp->kbp_std[tmp_index] != NULL) {
  73.                 Blast_KarlinBlk *kbp_gap;
  74.                 Blast_KarlinBlk *kbp;
  75.  
  76.                 sbp->kbp_gap_std[tmp_index] = Blast_KarlinBlkCreate();
  77.                 kbp_gap = sbp->kbp_gap_std[tmp_index];
  78.                 kbp     = sbp->kbp_std[tmp_index];
  80.                 kbp_gap->Lambda = kbp->Lambda;
  81.                 kbp_gap->K      = kbp->K;
  82.                 kbp_gap->logK   = kbp->logK;
  83.                 kbp_gap->H      = kbp->H;
  84.                 kbp_gap->paramC = kbp->paramC;
  85.             }
  86.         }
  88.     } else {
  89.         for (tmp_index = query_info->first_context;
  90.              tmp_index <= query_info->last_context; tmp_index++) {
  91.             sbp->kbp_gap_std[tmp_index] = Blast_KarlinBlkCreate();
  92.             Blast_KarlinBlkGappedCalc(sbp->kbp_gap_std[tmp_index],
  93.                                       scoring_options->gap_open,
  94.                                       scoring_options->gap_extend,
  95.                                       scoring_options->decline_align,
  96.                                       sbp->name, NULL);
  97.         }
  98.     }
  100.     sbp->kbp_gap = sbp->kbp_gap_std;
  101.     return 0;
  102. }
  104. static Int2
  105. PHIScoreBlkFill(BlastScoreBlk* sbp, const BlastScoringOptions* options,
  106.    Blast_Message** blast_message)
  107. {
  108.    Blast_KarlinBlk* kbp;
  109.    char buffer[1024];
  110.    Int2 status = 0;
  112.    sbp->read_in_matrix = TRUE;
  113.    sbp->name = strdup(options->matrix);
  114.    if ((status = BLAST_ScoreBlkMatFill(sbp, options->matrix_path)) != 0)
  115.       return status;
  116.    kbp = sbp->kbp_gap_std[0] = Blast_KarlinBlkCreate();
  117.    sbp->kbp_std = sbp->kbp_gap_std;
  119.    if (0 == strcmp("BLOSUM62", options->matrix)) {
  120.       kbp->paramC = 0.50;
  121.       if ((11 == options->gap_open) && (1 == options->gap_extend)) {
  122.          kbp->Lambda = 0.270;
  123.          kbp->K = 0.047;
  124.          return status;
  125.       }
  126.       if ((9 == options->gap_open) && (2 == options->gap_extend)) {
  127.          kbp->Lambda = 0.285;
  128.          kbp->K = 0.075;
  129.          return status;
  130.       }
  131.       if ((8 == options->gap_open) && (2 == options->gap_extend)) {
  132.          kbp->Lambda = 0.265;
  133.          kbp->K = 0.046;
  134.          return status;
  135.       }
  136.       if ((7 == options->gap_open) && (2 == options->gap_extend)) {
  137.          kbp->Lambda = 0.243;
  138.          kbp->K = 0.032;
  139.          return status;
  140.       }
  141.       if ((12 == options->gap_open) && (1 == options->gap_extend)) {
  142.          kbp->Lambda = 0.281;
  143.          kbp->K = 0.057;
  144.          return status;
  145.       }
  146.       if ((10 == options->gap_open) && (1 == options->gap_extend)) {
  147.          kbp->Lambda = 0.250;
  148.          kbp->K = 0.033;
  149.          return status;
  150.       }
  151.       sprintf(buffer, "The combination %d for gap opening cost and %d for gap extension is not supported in PHI-BLAST with matrix %sn", options->gap_open, options->gap_extend, options->matrix);
  152.       Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  153.    }
  154.    else {
  155.       if (0 == strcmp("PAM30", options->matrix)) { 
  156.          kbp->paramC = 0.30;
  157.          if ((9 == options->gap_open) && (1 == options->gap_extend)) {
  158.             kbp->Lambda = 0.295;
  159.             kbp->K = 0.13;
  160.             return status;
  161.          }
  162.          if ((7 == options->gap_open) && (2 == options->gap_extend)) {
  163.             kbp->Lambda = 0.306;
  164.             kbp->K = 0.15;
  165.             return status;
  166.          }
  167.          if ((6 == options->gap_open) && (2 == options->gap_extend)) {
  168.             kbp->Lambda = 0.292;
  169.             kbp->K = 0.13;
  170.             return status;
  171.          }
  172.          if ((5 == options->gap_open) && (2 == options->gap_extend)) {
  173.             kbp->Lambda = 0.263;
  174.             kbp->K = 0.077;
  175.             return status;
  176.          }
  177.          if ((10 == options->gap_open) && (1 == options->gap_extend)) {
  178.             kbp->Lambda = 0.309;
  179.             kbp->K = 0.15;
  180.             return status;
  181.          }
  182.          if ((8 == options->gap_open) && (1 == options->gap_extend)) {
  183.             kbp->Lambda = 0.270;
  184.             kbp->K = 0.070;
  185.             return status;
  186.          }
  187.          sprintf(buffer, "The combination %d for gap opening cost and %d for gap extension is not supported in PHI-BLAST with matrix %sn", options->gap_open, options->gap_extend, options->matrix);
  188.          Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  189.       }
  190.       else {
  191.          if (0 == strcmp("PAM70", options->matrix)) { 
  192.             kbp->paramC = 0.35;
  193.             if ((10 == options->gap_open) && (1 == options->gap_extend)) {
  194.                kbp->Lambda = 0.291;
  195.                kbp->K = 0.089;
  196.                return status;
  197.             }
  198.             if ((8 == options->gap_open) && (2 == options->gap_extend)) {
  199.                kbp->Lambda = 0.303;
  200.                kbp->K = 0.13;
  201.                return status;
  202.             }
  203.             if ((7 == options->gap_open) && (2 == options->gap_extend)) {
  204.                kbp->Lambda = 0.287;
  205.                kbp->K = 0.095;
  206.                return status;
  207.             }
  208.             if ((6 == options->gap_open) && (2 == options->gap_extend)) {
  209.                kbp->Lambda = 0.269;
  210.                kbp->K = 0.079;
  211.                return status;
  212.             }
  213.             if ((11 == options->gap_open) && (1 == options->gap_extend)) {
  214.                kbp->Lambda = 0.307;
  215.                kbp->K = 0.13;
  216.                return status;
  217.             }
  218.             if ((9 == options->gap_open) && (1 == options->gap_extend)) {
  219.                kbp->Lambda = 0.269;
  220.                kbp->K = 0.058;
  221.                return status;
  222.             }
  223.             sprintf(buffer, "The combination %d for gap opening cost and %d for gap extension is not supported in PHI-BLAST with matrix %sn", options->gap_open, options->gap_extend, options->matrix);
  224.             Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  225.          }
  226.          else {
  227.             if (0 == strcmp("BLOSUM80", options->matrix)) { 
  228.                kbp->paramC = 0.40;
  229.                if ((10 == options->gap_open) && (1 == options->gap_extend)) {
  230.                   kbp->Lambda = 0.300;
  231.                   kbp->K = 0.072;
  232.                   return status;
  233.                }
  234.                if ((8 == options->gap_open) && (2 == options->gap_extend)) {
  235.                   kbp->Lambda = 0.308;
  236.                   kbp->K = 0.089;
  237.                   return status;
  238.                }
  239.                if ((7 == options->gap_open) && (2 == options->gap_extend)) {
  240.                   kbp->Lambda = 0.295;
  241.                   kbp->K = 0.077;
  242.                   return status;
  243.                }
  244.                if ((6 == options->gap_open) && (2 == options->gap_extend)) {
  245.                   kbp->Lambda = 0.271;
  246.                   kbp->K = 0.051;
  247.                   return status;
  248.                }
  249.                if ((11 == options->gap_open) && (1 == options->gap_extend)) {
  250.                   kbp->Lambda = 0.314;
  251.                   kbp->K = 0.096;
  252.                   return status;
  253.                }
  254.                if ((9 == options->gap_open) && (1 == options->gap_extend)) {
  255.                   kbp->Lambda = 0.277;
  256.                   kbp->K = 0.046;
  257.                   return status;
  258.                }
  259.                sprintf(buffer, "The combination %d for gap opening cost and %d for gap extension is not supported in PHI-BLAST with matrix %sn", options->gap_open, options->gap_extend, options->matrix);
  260.                Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  261.             }
  262.             else {
  263.                if (0 == strcmp("BLOSUM45", options->matrix)) { 
  264.                   kbp->paramC = 0.60;
  265.                   if ((14 == options->gap_open) && (2 == options->gap_extend)) {
  266.                      kbp->Lambda = 0.199;
  267.                      kbp->K = 0.040;
  268.                      return status;
  269.                   }
  270.                   if ((13 == options->gap_open) && (3 == options->gap_extend)) {
  271.                      kbp->Lambda = 0.209;
  272.                      kbp->K = 0.057;
  273.                      return status;
  274.                   }
  275.                   if ((12 == options->gap_open) && (3 == options->gap_extend)) {
  276.                      kbp->Lambda = 0.203;
  277.                      kbp->K = 0.049;
  278.                      return status;
  279.                   }
  280.                   if ((11 == options->gap_open) && (3 == options->gap_extend)) {
  281.                      kbp->Lambda = 0.193;
  282.                      kbp->K = 0.037;
  283.                      return status;
  284.                   }
  285.                   if ((10 == options->gap_open) && (3 == options->gap_extend)) {
  286.                      kbp->Lambda = 0.182;
  287.                      kbp->K = 0.029;
  288.                      return status;
  289.                   }
  290.                   if ((15 == options->gap_open) && (2 == options->gap_extend)) {
  291.                      kbp->Lambda = 0.206;
  292.                      kbp->K = 0.049;
  293.                      return status;
  294.                   }
  295.                   if ((13 == options->gap_open) && (2 == options->gap_extend)) {
  296.                      kbp->Lambda = 0.190;
  297.                      kbp->K = 0.032;
  298.                      return status;
  299.                   }
  300.                   if ((12 == options->gap_open) && (2 == options->gap_extend)) {
  301.                      kbp->Lambda = 0.177;
  302.                      kbp->K = 0.023;
  303.                      return status;
  304.                   }
  305.                   if ((19 == options->gap_open) && (1 == options->gap_extend)) {
  306.                      kbp->Lambda = 0.209;
  307.                      kbp->K = 0.049;
  308.                      return status;
  309.                   }
  310.                   if ((18 == options->gap_open) && (1 == options->gap_extend)) {
  311.                      kbp->Lambda = 0.202;
  312.                      kbp->K = 0.041;
  313.                      return status;
  314.                   }
  315.                   if ((17 == options->gap_open) && (1 == options->gap_extend)) {
  316.                      kbp->Lambda = 0.195;
  317.                      kbp->K = 0.034;
  318.                      return status;
  319.                   }
  320.                   if ((16 == options->gap_open) && (1 == options->gap_extend)) {
  321.                      kbp->Lambda = 0.183;
  322.                      kbp->K = 0.024;
  323.                      return status;
  324.                   }
  325.                   sprintf(buffer, "The combination %d for gap opening cost and %d for gap extension is not supported in PHI-BLAST with matrix %sn", options->gap_open, options->gap_extend, options->matrix);
  326.                   Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  327.                }
  328.                else {
  329.                   sprintf(buffer, "Matrix %s not allowed in PHI-BLASTn", options->matrix);
  330.                   Blast_MessageWrite(blast_message, BLAST_SEV_WARNING, 2, 1, buffer);
  331.                }
  332.             }
  333.          }
  334.       }
  335.    }
  336. return status;
  337. }
  339. Int2
  340. BlastScoreBlkMatrixInit(Uint1 program_number, 
  341.                   const BlastScoringOptions* scoring_options,
  342.                   BlastScoreBlk* sbp)
  343. {
  344.    if (!sbp || !scoring_options)
  345.       return 1;
  347.    if (program_number == blast_type_blastn) {
  349.       BLAST_ScoreSetAmbigRes(sbp, 'N');
  350.       sbp->penalty = scoring_options->penalty;
  351.       sbp->reward = scoring_options->reward;
  352.       if (scoring_options->matrix_path &&
  353.           *scoring_options->matrix_path != NULLB &&
  354.           scoring_options->matrix && *scoring_options->matrix != NULLB) {
  356.           sbp->read_in_matrix = TRUE;
  357.           sbp->name = strdup(scoring_options->matrix);
  359.        } else {
  360.           char buffer[50];
  361.           sbp->read_in_matrix = FALSE;
  362.           sprintf(buffer, "blastn matrix:%ld %ld",
  363.                   (long) sbp->reward, (long) sbp->penalty);
  364.           sbp->name = strdup(buffer);
  365.        }
  367.     } else {
  368.        char* p = NULL;
  370.        sbp->read_in_matrix = TRUE;
  371.        BLAST_ScoreSetAmbigRes(sbp, 'X');
  372.        sbp->name = strdup(scoring_options->matrix);
  373.        /* protein matrices are in all caps by convention */
  374.        for (p = sbp->name; *p != NULLB; p++)
  375.           *p = toupper(*p);
  376.    }
  378.    return BLAST_ScoreBlkMatFill(sbp, scoring_options->matrix_path);
  379. }
  382. Int2 
  383. BlastSetup_GetScoreBlock(BLAST_SequenceBlk* query_blk, 
  384.                          BlastQueryInfo* query_info, 
  385.                          const BlastScoringOptions* scoring_options, 
  386.                          Uint1 program_number, 
  387.                          Boolean phi_align, 
  388.                          BlastScoreBlk* *sbpp, 
  389.                          double scale_factor, 
  390.                          Blast_Message* *blast_message)
  391. {
  392.     BlastScoreBlk* sbp;
  393.     Int2 status=0;      /* return value. */
  394.     Int4 context; /* loop variable. */
  396.     if (sbpp == NULL)
  397.        return 1;
  399.     if (program_number == blast_type_blastn)
  400.        sbp = BlastScoreBlkNew(BLASTNA_SEQ_CODE, query_info->last_context + 1);
  401.     else
  402.        sbp = BlastScoreBlkNew(BLASTAA_SEQ_CODE, query_info->last_context + 1);
  404.     if (!sbp)
  405.        return 1;
  407.     *sbpp = sbp;
  408.     sbp->scale_factor = scale_factor;
  410.     status = BlastScoreBlkMatrixInit(program_number, scoring_options, sbp);
  411.     if (status != 0)
  412.        return status;
  414.     for (context = query_info->first_context;
  415.          context <= query_info->last_context; ++context) {
  416.       
  417.         Int4 context_offset;
  418.         Int4 query_length;
  419.         Uint1 *buffer;              /* holds sequence */
  421.         /* For each query, check if forward strand is present */
  422.         if ((query_length = BLAST_GetQueryLength(query_info, context)) < 0)
  423.             continue;
  425.         context_offset = query_info->context_offsets[context];
  426.         buffer = &query_blk->sequence[context_offset];
  428.         if (!phi_align &&
  429.            (status = BLAST_ScoreBlkFill(sbp, (char *) buffer,
  430.                                         query_length, context))) {
  431.             Blast_MessageWrite(blast_message, BLAST_SEV_ERROR, 2, 1, 
  432.                    "Query completely filtered; nothing left to search");
  433.             return status;
  434.         }
  435.     }
  438.     /* Fills in block for gapped blast. */
  439.     if (phi_align) {
  440.        PHIScoreBlkFill(sbp, scoring_options, blast_message);
  441.     } else if (scoring_options->gapped_calculation) {
  442.        status = BlastScoreBlkGappedFill(sbp, scoring_options, 
  443.                                         program_number, query_info);
  444.        if (status) {
  445.           Blast_MessageWrite(blast_message, BLAST_SEV_ERROR, 2, 1, 
  446.                              "Unable to initialize scoring block");
  447.           return status;
  448.        }
  449.     }
  450.    
  451.     /* Get "ideal" values if the calculated Karlin-Altschul params bad. */
  452.     if (program_number == blast_type_blastx ||
  453.         program_number == blast_type_tblastx ||
  454.         program_number == blast_type_rpstblastn) {
  455.         /* Adjust the ungapped Karlin parameters */
  456.         sbp->kbp = sbp->kbp_std;
  457.         Blast_KarlinBlkStandardCalc(sbp, query_info->first_context,
  458.                                     query_info->last_context);
  459.     }
  461.     if (program_number == blast_type_blastx ||
  462.         program_number == blast_type_tblastx) {
  463.         /* Adjust the gapped Karlin parameters, if it is a gapped search */
  464.         if (scoring_options->gapped_calculation) {
  465.            sbp->kbp = sbp->kbp_gap_std;
  466.           Blast_KarlinBlkStandardCalc(sbp, query_info->first_context,
  467.                                        query_info->last_context);
  468.         }
  469.     }
  471.     /* Why are there so many Karlin block names?? */
  472.     sbp->kbp = sbp->kbp_std;
  473.     if (scoring_options->gapped_calculation)
  474.        sbp->kbp_gap = sbp->kbp_gap_std;
  476.     return 0;
  477. }
  479. Int2 BLAST_MainSetUp(Uint1 program_number,
  480.     const QuerySetUpOptions *qsup_options,
  481.     const BlastScoringOptions *scoring_options,
  482.     const BlastHitSavingOptions *hit_options,
  483.     BLAST_SequenceBlk *query_blk,
  484.     BlastQueryInfo *query_info,
  485.     double scale_factor,
  486.     BlastSeqLoc **lookup_segments, 
  487.     BlastMaskLoc **filter_out,
  488.     BlastScoreBlk **sbpp, 
  489.     Blast_Message **blast_message)
  490. {
  491.     Boolean mask_at_hash = FALSE; /* mask only for making lookup table? */
  492.     Int2 status = 0;            /* return value */
  493.     BlastMaskLoc *filter_maskloc = NULL;   /* Local variable for mask locs. */
  496.     status = BlastSetUp_GetFilteringLocations(query_blk, 
  497.                                               query_info, 
  498.                                               program_number, 
  499.                                               qsup_options->filter_string, 
  500.                                               &filter_maskloc, 
  501.                                               &mask_at_hash, 
  502.                                               blast_message);
  503.     if (status) {
  504.         return status;
  505.     } 
  507.     if (!mask_at_hash)
  508.     {
  509.         status = BlastSetUp_MaskQuery(query_blk, query_info, filter_maskloc, 
  510.                                       program_number);
  511.         if (status != 0) {
  512.             return status;
  513.         }
  514.     }
  516.     /* If there was a lower case mask, its contents have now been moved to 
  517.      * filter_maskloc and are no longer needed in the query block.
  518.     */
  519.     query_blk->lcase_mask = NULL;
  521.     if (filter_out)
  522.        *filter_out = filter_maskloc;
  523.     else 
  524.         filter_maskloc = BlastMaskLocFree(filter_maskloc);
  526.     if (program_number == blast_type_blastx && scoring_options->is_ooframe) {
  527.         BLAST_InitDNAPSequence(query_blk, query_info);
  528.     }
  530.     BLAST_ComplementMaskLocations(program_number, query_info, filter_maskloc, 
  531.                                   lookup_segments);
  534.     status = BlastSetup_GetScoreBlock(query_blk, query_info, scoring_options, 
  535.                                       program_number, hit_options->phi_align, 
  536.                                       sbpp, scale_factor, blast_message);
  537.     if (status > 0) {
  538.         return status;
  539.     }
  541.     return 0;
  542. }
  545. Int2 BLAST_CalcEffLengths (Uint1 program_number, 
  546.    const BlastScoringOptions* scoring_options,
  547.    const BlastEffectiveLengthsParameters* eff_len_params, 
  548.    const BlastScoreBlk* sbp, BlastQueryInfo* query_info)
  549. {
  550.    double alpha=0, beta=0; /*alpha and beta for new scoring system */
  551.    Int4 length_adjustment = 0;  /* length adjustment for current iteration. */
  552.    Int4 index; /* loop index. */
  553.    Int4 db_num_seqs; /* number of sequences in database. */
  554.    Int8 db_length; /* total length of database. */
  555.    Blast_KarlinBlk* *kbp_ptr; /* Array of Karlin block pointers */
  556.    Int4 query_length;   /* length of an individual query sequence */
  557.    Int8 effective_search_space = 0; /* Effective search space for a given 
  558.                                    sequence/strand/frame */
  559.    const BlastEffectiveLengthsOptions* eff_len_options = eff_len_params->options;
  561.    if (sbp == NULL)
  562.       return 1;
  564.    /* use overriding value from effective lengths options or the real value
  565.       from effective lengths parameters. */
  566.    if (eff_len_options->db_length > 0)
  567.       db_length = eff_len_options->db_length;
  568.    else
  569.       db_length = eff_len_params->real_db_length;
  571.    /* If database (subject) length is not available at this stage, 
  572.       do nothing */
  573.    if (db_length == 0)
  574.       return 0;
  576.    if (program_number == blast_type_tblastn || 
  577.        program_number == blast_type_tblastx)
  578.       db_length = db_length/3;
  580.    if (eff_len_options->dbseq_num > 0)
  581.       db_num_seqs = eff_len_options->dbseq_num;
  582.    else
  583.       db_num_seqs = eff_len_params->real_num_seqs;
  584.    
  585.    if (program_number != blast_type_blastn) {
  586.       if (scoring_options->gapped_calculation) {
  587.          BLAST_GetAlphaBeta(sbp->name,&alpha,&beta,TRUE, 
  588.             scoring_options->gap_open, scoring_options->gap_extend);
  589.       }
  590.    }
  591.    
  592.    if (scoring_options->gapped_calculation) 
  593.       kbp_ptr = sbp->kbp_gap_std;
  594.    else
  595.       kbp_ptr = sbp->kbp_std; 
  596.    
  597.    for (index = query_info->first_context;
  598.         index <= query_info->last_context;
  599.         index++) {
  600.        Blast_KarlinBlk * kbp;   /* statistical parameters for the
  601.                                    current context */
  602.        kbp = kbp_ptr[index];
  604.        if ( (query_length = BLAST_GetQueryLength(query_info, index)) > 0) {
  605.           /* Use the correct Karlin block. For blastn, two identical Karlin
  606.              blocks are allocated for each sequence (one per strand), but we
  607.              only need one of them.
  608.           */
  609.           if (program_number != blast_type_blastn &&
  610.               scoring_options->gapped_calculation) {
  611.              BLAST_ComputeLengthAdjustment(kbp->K, kbp->logK,
  612.                                            alpha/kbp->Lambda, beta,
  613.                                            query_length, db_length,
  614.                                            db_num_seqs, &length_adjustment);
  615.           } else {
  616.              BLAST_ComputeLengthAdjustment(kbp->K, kbp->logK, 1/kbp->H, 0,
  617.                                            query_length, db_length,
  618.                                            db_num_seqs, &length_adjustment);
  619.           }        
  620.           /* If overriding search space value is provided in options, 
  621.              do not calculate it. */
  622.           if (eff_len_options->searchsp_eff) {
  623.              effective_search_space = eff_len_options->searchsp_eff;
  624.           } else {
  625.              effective_search_space =
  626.                 (query_length - length_adjustment) * 
  627.                 (db_length - db_num_seqs*length_adjustment);
  629.              /* For translated RPS blast, the DB size is left unchanged
  630.                 and the query size is divided by 3 (for conversion to 
  631.                 a protein sequence) and multiplied by 6 (for 6 frames) */
  632.              if (program_number == blast_type_rpstblastn)
  633.                 effective_search_space *= (Int8)(NUM_FRAMES / CODON_LENGTH);
  634.           }
  635.        }
  636.        query_info->eff_searchsp_array[index] = effective_search_space;
  637.        query_info->length_adjustments[index] = length_adjustment;
  638.    }
  640.    return 0;
  641. }
  643. Int2 
  644. BLAST_GapAlignSetUp(Uint1 program_number,
  645.     const BlastSeqSrc* seq_src,
  646.     const BlastScoringOptions* scoring_options,
  647.     const BlastEffectiveLengthsOptions* eff_len_options,
  648.     const BlastExtensionOptions* ext_options,
  649.     const BlastHitSavingOptions* hit_options,
  650.     BlastQueryInfo* query_info, 
  651.     BlastScoreBlk* sbp, 
  652.     BlastScoringParameters** score_params,
  653.     BlastExtensionParameters** ext_params,
  654.     BlastHitSavingParameters** hit_params,
  655.     BlastEffectiveLengthsParameters** eff_len_params,
  656.     BlastGapAlignStruct** gap_align)
  657. {
  658.    Int2 status = 0;
  659.    Uint4 max_subject_length;
  660.    Int8 total_length;
  661.    Int4 num_seqs;
  663.    total_length = BLASTSeqSrcGetTotLen(seq_src);
  664.    
  665.    if (total_length > 0) {
  666.       num_seqs = BLASTSeqSrcGetNumSeqs(seq_src);
  667.    } else {
  668.       /* Not a database search; each subject sequence is considered
  669.          individually */
  670.       num_seqs = 1;
  671.    }
  673.    /* Initialize the effective length parameters with real values of
  674.       database length and number of sequences */
  675.    BlastEffectiveLengthsParametersNew(eff_len_options, total_length, num_seqs, 
  676.                                       eff_len_params);
  677.    if ((status = BLAST_CalcEffLengths(program_number, scoring_options, 
  678.                     *eff_len_params, sbp, query_info)) != 0)
  679.       return status;
  681.    BlastScoringParametersNew(scoring_options, sbp, score_params);
  683.    BlastExtensionParametersNew(program_number, ext_options, sbp, 
  684.                                query_info, ext_params);
  686.    BlastHitSavingParametersNew(program_number, hit_options, *ext_params, 
  687.                                sbp, query_info, hit_params);
  689.    /* To initialize the gapped alignment structure, we need to know the 
  690.       maximal subject sequence length */
  691.    max_subject_length = BLASTSeqSrcGetMaxSeqLen(seq_src);
  693.    if ((status = BLAST_GapAlignStructNew(*score_params, *ext_params, 
  694.                     max_subject_length, query_info->max_length, sbp, 
  695.                     gap_align)) != 0) {
  696.       return status;
  697.    }
  699.    return status;
  700. }
  702. Int2 BLAST_OneSubjectUpdateParameters(Uint1 program_number,
  703.                     Uint4 subject_length,
  704.                     const BlastScoringOptions* scoring_options,
  705.                     BlastQueryInfo* query_info, 
  706.                     BlastScoreBlk* sbp, 
  707.                     const BlastExtensionParameters* ext_params,
  708.                     BlastHitSavingParameters* hit_params,
  709.                     BlastInitialWordParameters* word_params,
  710.                     BlastEffectiveLengthsParameters* eff_len_params)
  711. {
  712.    Int2 status = 0;
  713.    eff_len_params->real_db_length = subject_length;
  714.    if ((status = BLAST_CalcEffLengths(program_number, scoring_options, 
  715.                                       eff_len_params, sbp, query_info)) != 0)
  716.       return status;
  717.    /* Update cutoff scores in hit saving parameters */
  718.    BlastHitSavingParametersUpdate(program_number, ext_params,
  719.                                   sbp, query_info, 
  720.                                   hit_params);
  721.    
  722.    if (word_params) {
  723.       /* Update cutoff scores in initial word parameters */
  724.       BlastInitialWordParametersUpdate(program_number, hit_params, ext_params,
  725.          sbp, query_info, subject_length, word_params);
  726.    }
  727.    return status;
  728. }