← Back to branch summary

~ubuntu-branches/ubuntu/edgy/ncbi-tools6/edgy

« back to all changes in this revision

Viewing changes to algo/blast/core/blast_parameters.c

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-19 23:28:07 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20060719232807-et3cdmcjgmnyleyx
Tags: 6.1.20060507-3ubuntu1
Re-merge with Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/core/blast_parameters.c
1
 
/* $Id: blast_parameters.c,v 1.12 2005/11/16 14:27:03 madden Exp $
 
1
/* $Id: blast_parameters.c,v 1.16 2006/01/12 20:34:32 camacho Exp $
2
2
 * ===========================================================================
3
3
 *
4
4
 *                            PUBLIC DOMAIN NOTICE
30
30
 
31
31
#ifndef SKIP_DOXYGEN_PROCESSING
32
32
static char const rcsid[] = 
33
 
    "$Id: blast_parameters.c,v 1.12 2005/11/16 14:27:03 madden Exp $";
 
33
    "$Id: blast_parameters.c,v 1.16 2006/01/12 20:34:32 camacho Exp $";
34
34
#endif /* SKIP_DOXYGEN_PROCESSING */
35
35
 
36
36
#include <algo/blast/core/blast_parameters.h>
48
48
 */
49
49
static Boolean s_BlastKarlinBlkIsValid(const Blast_KarlinBlk* kbp)
50
50
{
51
 
    return (kbp->Lambda > 0 && kbp->K > 0 && kbp->H > 0);
 
51
    if ( !kbp ) {
 
52
        return FALSE;
 
53
    } else {
 
54
        return (kbp->Lambda > 0 && kbp->K > 0 && kbp->H > 0);
 
55
    }
52
56
}
53
57
 
54
58
/** Returns the first valid Karlin-Altchul block from the list of blocks.
69
73
    ASSERT(kbp_in && query_info && kbp_ret);
70
74
 
71
75
    for (i=query_info->first_context; i<=query_info->last_context; i++) {
 
76
         ASSERT(s_BlastKarlinBlkIsValid(kbp_in[i]) ==
 
77
                query_info->contexts[i].is_valid);
72
78
         if (s_BlastKarlinBlkIsValid(kbp_in[i])) {
73
79
              *kbp_ret = kbp_in[i];
74
80
              status = 0;
97
103
    ASSERT(kbp_in && query_info);
98
104
 
99
105
    for (i=query_info->first_context; i<=query_info->last_context; i++) {
 
106
        ASSERT(s_BlastKarlinBlkIsValid(kbp_in[i]) ==
 
107
               query_info->contexts[i].is_valid);
100
108
        if (s_BlastKarlinBlkIsValid(kbp_in[i])) {
101
109
            if (min_lambda > kbp_in[i]->Lambda)
102
110
            {
141
149
               retval = eRight;
142
150
         break;
143
151
     case MB_LOOKUP_TABLE:
144
 
         if (((BlastMBLookupTable*)lookup_wrap->lut)->ag_scanning_mode == TRUE)
145
 
               retval = eRightAndLeft;
146
 
         else if (((BlastMBLookupTable*)lookup_wrap->lut)->template_length > 0)
 
152
         if (((BlastMBLookupTable*)lookup_wrap->lut)->template_length > 0)
147
153
               retval = eUpdateDiag;   /* Used for discontiguous megablast. */
148
154
         else
149
 
               retval = eRight;
 
155
               retval = eRightAndLeft;
150
156
         break;
151
157
   }
152
158
   ASSERT(retval != eMaxSeedExtensionMethod);
202
208
{
203
209
   Blast_KarlinBlk* kbp_std;
204
210
   Int2 status = 0;
205
 
   const int kQueryLenForStacks = 50;  /* Use stacks rather than diag for 
206
 
                                     any query longer than this.  Only for blastn. */
 
211
   const int kQueryLenForStacks = 8000;  /* For blastn, use stacks rather 
 
212
                                            than diags for any query longer 
 
213
                                            than this */
207
214
 
208
215
   /* If parameters pointer is NULL, there is nothing to fill, 
209
216
      so don't do anything */
238
245
               hit_params, sbp, query_info,
239
246
               subject_length, *parameters);
240
247
 
 
248
   if (program_number == eBlastTypeBlastn) {
 
249
      Int4 i;
 
250
      Int4 reward = sbp->reward;
 
251
      Int4 penalty = sbp->penalty;
 
252
      Int4 *table = (*parameters)->nucl_score_table;
 
253
 
 
254
      /* nucleotide ungapped extensions are first computed
 
255
         approximately, and then recomputed exactly if the
 
256
         approximate score is high enough. The approximate
 
257
         computation considers nucleotides in batches of 4,
 
258
         so a table is needed that contains the combined scores
 
259
         of all combinations of 4 matches/mismatches */
 
260
 
 
261
      for (i = 0; i < 256; i++) {
 
262
         /* break the bit pattern of i into four 2-bit groups.
 
263
            If bits in a group are set, that group represents 
 
264
            a mismatch */
 
265
 
 
266
         Int4 score = 0;
 
267
         if (i & 3) score += penalty; else score += reward;
 
268
         if ((i >> 2) & 3) score += penalty; else score += reward;
 
269
         if ((i >> 4) & 3) score += penalty; else score += reward;
 
270
         if (i >> 6) score += penalty; else score += reward;
 
271
         table[i] = score;
 
272
      }
 
273
   }
 
274
 
241
275
   return status;
242
276
}
243
277
 
292
326
         Blast_KarlinBlk* kbp_ungap = sbp->kbp_std[index];
293
327
         const BlastInitialWordOptions* kOptions = parameters->options;
294
328
 
 
329
         ASSERT(s_BlastKarlinBlkIsValid(kbp_ungap) ==
 
330
                query_info->contexts[index].is_valid);
295
331
         if (s_BlastKarlinBlkIsValid(kbp_ungap)) {
296
332
            gap_trigger = (Int4) ((kOptions->gap_trigger * NCBIMATH_LN2 + 
297
333
                                   kbp_ungap->logK) / kbp_ungap->Lambda);
302
338
          continue;
303
339
 
304
340
      kbp = kbp_array[index];
305
 
      if (!s_BlastKarlinBlkIsValid(kbp))  /* skip invalid Karlin blocks */
 
341
      if (!s_BlastKarlinBlkIsValid(kbp)) { /* skip invalid Karlin blocks */
 
342
          ASSERT(query_info->contexts[index].is_valid == FALSE);
306
343
          continue;
 
344
      }
307
345
 
308
346
      if (!gapped_calculation || program_number == eBlastTypeBlastn) {
309
347
         double cutoff_e = s_GetCutoffEvalue(program_number);
337
375
 
338
376
   parameters->cutoff_score = MIN(hit_params->cutoff_score_max, cutoff_s);
339
377
   
 
378
   /* Nucleotide searches first compute an approximate ungapped
 
379
      alignment and compare it to a reduced ungapped cutoff score */
 
380
   if (program_number == eBlastTypeBlastn) {
 
381
       parameters->reduced_nucl_cutoff_score = 
 
382
                          (Int4)(0.6 * parameters->cutoff_score);
 
383
   }
 
384
 
340
385
   /* Note that x_dropoff_init stays constant throughout the search.
341
386
      The cutoff_score and x_dropoff parameters may be updated multiple times, 
342
387
      if every subject sequence is treated individually. Hence we need to know 
387
432
         when rescaling and composition based statistics is applied, as we
388
433
         lose precision. Therefore this is redone in Kappa_RedoAlignmentCore */
389
434
      params->gap_x_dropoff_final = (Int4) 
390
 
          (options->gap_x_dropoff_final*NCBIMATH_LN2 / min_lambda);
 
435
          MAX(options->gap_x_dropoff_final*NCBIMATH_LN2 / min_lambda, params->gap_x_dropoff);
391
436
   }
392
437
   
393
438
   if (sbp->scale_factor > 1.0) {
678
723
         double evalue = options->expect_value;
679
724
 
680
725
         kbp = kbp_array[context];
681
 
         if (!s_BlastKarlinBlkIsValid(kbp))  /* skip invalid Karlin blocks */
 
726
         if (!s_BlastKarlinBlkIsValid(kbp)) { /* skip invalid Karlin blocks */
 
727
             ASSERT(query_info->contexts[context].is_valid == FALSE);
682
728
             continue;
 
729
         }
683
730
         searchsp = query_info->contexts[context].eff_searchsp;
684
731
         if (searchsp == 0)         /* skip invalid contexts */
685
732
            continue;
715
762
            Int4 new_cutoff = 0;
716
763
 
717
764
            kbp = kbp_array[context];
718
 
            if (!s_BlastKarlinBlkIsValid(kbp))  /* skip invalid Karlin blocks */
 
765
            if (!s_BlastKarlinBlkIsValid(kbp)) {/* skip invalid Karlin blocks */
 
766
                ASSERT(query_info->contexts[context].is_valid == FALSE);
719
767
                continue;
 
768
            }
720
769
            BLAST_Cutoffs(&new_cutoff, &evalue_hsp, kbp, searchsp,
721
770
                       TRUE, params->link_hsp_params->gap_decay_rate);
722
771
            /* Update the computed cutoff if new_cutoff is smaller */
826
875
 * ===========================================================================
827
876
 *
828
877
 * $Log: blast_parameters.c,v $
 
878
 * Revision 1.16  2006/01/12 20:34:32  camacho
 
879
 * + assertions for validity of context
 
880
 *
 
881
 * Revision 1.15  2006/01/03 17:53:20  papadopo
 
882
 * 1. increase the cutoff query size for using stacks
 
883
 * 2. initialize fields for approximate nucleotide ungapped alignment
 
884
 *
 
885
 * Revision 1.14  2006/01/03 14:18:42  madden
 
886
 * In BlastExtensionParametersNew raise gap_x_dropoff_final to gap_x_dropoff if it is lower
 
887
 *
 
888
 * Revision 1.13  2005/12/19 16:12:30  papadopo
 
889
 * remove the possibility of specifying eRight for megablast extension method
 
890
 *
829
891
 * Revision 1.12  2005/11/16 14:27:03  madden
830
892
 * Fix spelling in CRN
831
893
 *