← Back to branch summary

~ubuntu-branches/ubuntu/feisty/ncbi-tools6/feisty

« back to all changes in this revision

Viewing changes to algo/blast/core/blast_engine.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_engine.c
1
 
/* $Id: blast_engine.c,v 1.203 2005/11/22 13:44:13 coulouri Exp $
 
1
/* $Id: blast_engine.c,v 1.218 2006/05/05 13:44:35 coulouri Exp $
2
2
 * ===========================================================================
3
3
 *
4
4
 *                            PUBLIC DOMAIN NOTICE
55
55
 
56
56
#ifndef SKIP_DOXYGEN_PROCESSING
57
57
static char const rcsid[] = 
58
 
    "$Id: blast_engine.c,v 1.203 2005/11/22 13:44:13 coulouri Exp $";
 
58
    "$Id: blast_engine.c,v 1.218 2006/05/05 13:44:35 coulouri Exp $";
59
59
#endif /* SKIP_DOXYGEN_PROCESSING */
60
60
 
61
61
#include <algo/blast/core/blast_engine.h>
73
73
 
74
74
NCBI_XBLAST_EXPORT const int   kBlastMajorVersion = 2;
75
75
NCBI_XBLAST_EXPORT const int   kBlastMinorVersion = 2;
76
 
NCBI_XBLAST_EXPORT const int   kBlastPatchVersion = 13;
77
 
NCBI_XBLAST_EXPORT const char* kBlastReleaseDate = "Nov-27-2005";
 
76
NCBI_XBLAST_EXPORT const int   kBlastPatchVersion = 14;
 
77
NCBI_XBLAST_EXPORT const char* kBlastReleaseDate = "May-07-2006";
78
78
 
79
79
/** Structure to be passed to s_BlastSearchEngineCore, containing pointers 
80
80
    to various preallocated structures and arrays. */
87
87
                                              pointer */
88
88
   BlastInitHitList* init_hitlist; /**< Placeholder for HSPs after 
89
89
                                        ungapped extension */
90
 
   BlastHSPList* hsp_list; /**< Placeholder for HSPs after gapped 
91
 
                                extension */
92
90
   BlastOffsetPair* offset_pairs; /**< Array of offset pairs for initial seeds. */
93
91
   Uint1* translation_buffer; /**< Placeholder for translated subject
94
92
                                   sequences */
103
101
{
104
102
   BlastExtendWordFree(aux_struct->ewp);
105
103
   BLAST_InitHitListFree(aux_struct->init_hitlist);
106
 
   Blast_HSPListFree(aux_struct->hsp_list);
107
104
   sfree(aux_struct->offset_pairs);
108
105
   
109
106
   sfree(aux_struct);
123
120
static void 
124
121
s_TranslateHSPsToDNAPCoord(EBlastProgramType program, 
125
122
                           BlastInitHitList* init_hitlist, 
126
 
                           BlastQueryInfo* query_info,
 
123
                           const BlastQueryInfo* query_info,
127
124
                           Int2 subject_frame, Int4 subject_length, 
128
125
                           Int4 offset)
129
126
{
196
193
   OffsetArrayToContextOffsets(info, new_offsets, kProgram);
197
194
}
198
195
 
199
 
/** The core of the BLAST search: comparison between the (concatenated)
200
 
 * query against one subject sequence. Translation of the subject sequence
201
 
 * into 6 frames is done inside, if necessary. If subject sequence is 
202
 
 * too long, it can be split into several chunks. 
 
196
/** Searches only one context of a database sequence, but does all chunks if it is split.
203
197
 * @param program_number BLAST program type [in]
204
198
 * @param query Query sequence structure [in]
205
 
 * @param query_info_in Query information [in]
 
199
 * @param query_info Query information [in]
206
200
 * @param subject Subject sequence structure [in]
 
201
 * @param orig_length original length of query before translation [in]
207
202
 * @param lookup Lookup table [in]
208
203
 * @param gap_align Structure for gapped alignment information [in]
209
204
 * @param score_params Scoring parameters [in]
211
206
 *                    parameters [in]
212
207
 * @param ext_params Gapped extension parameters [in]
213
208
 * @param hit_params Hit saving parameters [in]
214
 
 * @param db_options Database options [in]
215
209
 * @param diagnostics Hit counts and other diagnostics [in] [out]
216
210
 * @param aux_struct Structure containing different auxiliary data and memory
217
211
 *                   for the preliminary stage of the BLAST search [in]
218
 
 * @param hsp_list_out List of HSPs found for a given subject sequence [in]
 
212
 * @param hsp_list_out_ptr List of HSPs found for a given subject sequence [out]
 
213
 * @param interrupt_search function callback to allow interruption of BLAST
 
214
 *                   search [in, optional]
 
215
 * @param progress_info contains information about the progress of the current
 
216
 *                   BLAST search [in|out]
219
217
 */
 
218
 
220
219
static Int2
221
 
s_BlastSearchEngineCore(EBlastProgramType program_number, BLAST_SequenceBlk* query, 
222
 
   BlastQueryInfo* query_info_in, BLAST_SequenceBlk* subject, 
223
 
   LookupTableWrap* lookup, BlastGapAlignStruct* gap_align, 
 
220
s_BlastSearchEngineOneContext(EBlastProgramType program_number, 
 
221
   BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
 
222
   BLAST_SequenceBlk* subject, Int4 orig_length, LookupTableWrap* lookup, 
 
223
   BlastGapAlignStruct* gap_align, 
224
224
   const BlastScoringParameters* score_params, 
225
225
   const BlastInitialWordParameters* word_params, 
226
226
   const BlastExtensionParameters* ext_params, 
227
227
   const BlastHitSavingParameters* hit_params, 
228
 
   const BlastDatabaseOptions* db_options,
229
228
   BlastDiagnostics* diagnostics,
230
229
   BlastCoreAuxStruct* aux_struct,
231
 
   BlastHSPList** hsp_list_out)
 
230
   BlastHSPList** hsp_list_out_ptr,
 
231
   TInterruptFnPtr interrupt_search, 
 
232
   SBlastProgress* progress_info)
232
233
{
233
 
   BlastInitHitList* init_hitlist = aux_struct->init_hitlist;
234
 
   BlastHSPList* hsp_list = aux_struct->hsp_list;
235
 
   Uint1* translation_buffer = NULL;
236
 
   Int4* frame_offsets   = NULL;
237
 
   Int4* frame_offsets_a = NULL; /* Will be freed if non-null */
238
 
   BlastHitSavingOptions* hit_options = hit_params->options;
239
 
   BlastScoringOptions* score_options = score_params->options;
240
 
   BlastHSPList* combined_hsp_list = NULL;
241
 
   Int2 status = 0;
242
 
   Uint4 context, first_context, last_context;
243
 
   Int4 orig_length = subject->length;
244
 
   Uint1* orig_sequence = subject->sequence;
245
 
   Int4 **matrix;
246
 
   Int4 hsp_num_max;
247
 
   BlastUngappedStats* ungapped_stats = NULL;
248
 
   BlastGappedStats* gapped_stats = NULL;
249
 
   BlastQueryInfo* query_info = query_info_in;
250
 
   Boolean prelim_traceback = 
251
 
      (ext_params->options->ePrelimGapExt == eGreedyWithTracebackExt);
252
 
 
253
 
   const Boolean kTranslatedSubject = 
254
 
        (Blast_SubjectIsTranslated(program_number) || program_number == eBlastTypeRpsTblastn);
255
 
   const Boolean kNucleotide = (program_number == eBlastTypeBlastn ||
256
 
                                program_number == eBlastTypePhiBlastn);
257
 
 
258
 
   if (kTranslatedSubject) {
259
 
      first_context = 0;
260
 
      last_context = 5;
261
 
      if (score_options->is_ooframe) {
262
 
         BLAST_GetAllTranslations(orig_sequence, eBlastEncodingNcbi2na,
263
 
            orig_length, db_options->gen_code_string, &translation_buffer,
264
 
            &frame_offsets, &subject->oof_sequence);
265
 
         subject->oof_sequence_allocated = TRUE;
266
 
         frame_offsets_a = frame_offsets;
267
 
      } else if (program_number == eBlastTypeRpsTblastn ) {
268
 
          /* For RPS tblastn, subject is actually query, which has already 
269
 
             been translated during the setup stage. */
270
 
          translation_buffer = orig_sequence - 1;
271
 
          frame_offsets_a = frame_offsets =
272
 
              ContextOffsetsToOffsetArray(query_info_in);
273
 
      } else {
274
 
         BLAST_GetAllTranslations(orig_sequence, eBlastEncodingNcbi2na,
275
 
            orig_length, db_options->gen_code_string, &translation_buffer,
276
 
            &frame_offsets, NULL);
277
 
         frame_offsets_a = frame_offsets;
278
 
      }
279
 
   } else if (kNucleotide) {
280
 
      first_context = 1;
281
 
      last_context = 1;
282
 
   } else {
283
 
      first_context = 0;
284
 
      last_context = 0;
285
 
   }
286
 
 
287
 
   *hsp_list_out = NULL;
288
 
 
289
 
   if (gap_align->positionBased)
290
 
      matrix = gap_align->sbp->psi_matrix->pssm->data;
291
 
   else
292
 
      matrix = gap_align->sbp->matrix->data;
293
 
 
294
 
   hsp_num_max = (hit_options->hsp_num_max ? hit_options->hsp_num_max : INT4_MAX);
295
 
 
296
 
   if (diagnostics) {
297
 
      ungapped_stats = diagnostics->ungapped_stat;
298
 
      gapped_stats = diagnostics->gapped_stat;
299
 
   }
300
 
 
301
 
   /* Substitute query info by concatenated database info for RPS BLAST search */
302
 
   if (Blast_ProgramIsRpsBlast(program_number)) {
303
 
      BlastRPSLookupTable* lut = (BlastRPSLookupTable*) lookup->lut;
304
 
      query_info = (BlastQueryInfo*) calloc(1, sizeof(BlastQueryInfo));
305
 
      query_info->num_queries = lut->num_profiles;
306
 
      query_info->last_context = lut->num_profiles - 1;
307
 
      /* Since this will really be "subject info", not "query info",
308
 
         pass program argument such that all frames will be set to 0. */
309
 
      s_RPSOffsetArrayToContextOffsets(query_info,
310
 
                                       lut->rps_seq_offsets);
311
 
   }
312
 
 
313
 
   /* Loop over frames of the subject sequence */
314
 
   for (context=first_context; context<=last_context; context++) {
 
234
      Int2 status = 0; /* return value */
315
235
      Int4 chunk; /* loop variable below. */
316
236
      Int4 num_chunks; /* loop variable below. */
317
237
      Int4 offset = 0; /* Used as offset into subject sequence (if chunked) */
318
238
      Int4 total_subject_length; /* Length of subject sequence used when split. */
319
 
 
320
 
      if (kTranslatedSubject) {
321
 
         subject->frame =
322
 
             BLAST_ContextToFrame(eBlastTypeBlastx, context);
323
 
         subject->sequence = 
324
 
            translation_buffer + frame_offsets[context] + 1;
325
 
         subject->length = 
326
 
           frame_offsets[context+1] - frame_offsets[context] - 1;
327
 
      } else {
328
 
         subject->frame = context;
 
239
      BlastHSPList* combined_hsp_list = NULL;
 
240
      BlastHSPList* hsp_list = NULL;
 
241
      BlastInitHitList* init_hitlist = aux_struct->init_hitlist;
 
242
      BlastScoringOptions* score_options = score_params->options;
 
243
      BlastUngappedStats* ungapped_stats = NULL;
 
244
      BlastGappedStats* gapped_stats = NULL;
 
245
      Int4 **matrix;
 
246
      const Boolean kPrelimTraceback = 
 
247
         (ext_params->options->ePrelimGapExt == eGreedyWithTracebackExt);
 
248
      const Boolean kTranslatedSubject = 
 
249
        (Blast_SubjectIsTranslated(program_number) || program_number == eBlastTypeRpsTblastn);
 
250
      const Boolean kNucleotide = (program_number == eBlastTypeBlastn ||
 
251
                                program_number == eBlastTypePhiBlastn);
 
252
      const int kHspNumMax = BlastHspNumMax(score_options->gapped_calculation, hit_params->options);
 
253
     
 
254
      if (diagnostics) {
 
255
         ungapped_stats = diagnostics->ungapped_stat;
 
256
         gapped_stats = diagnostics->gapped_stat;
329
257
      }
330
 
     
 
258
 
 
259
      if (gap_align->positionBased)
 
260
         matrix = gap_align->sbp->psi_matrix->pssm->data;
 
261
      else
 
262
         matrix = gap_align->sbp->matrix->data;
 
263
 
331
264
      /* Split subject sequence into chunks if it is too long */
332
265
      num_chunks = (subject->length - DBSEQ_CHUNK_OVERLAP) / 
333
266
         (MAX_DBSEQ_LEN - DBSEQ_CHUNK_OVERLAP) + 1;
334
267
      total_subject_length = subject->length;
335
268
      
336
269
      for (chunk = 0; chunk < num_chunks; ++chunk) {
 
270
         /* Delete if not done in last loop iteration to prevent memory leak. */
 
271
         hsp_list = Blast_HSPListFree(hsp_list);  /* In case this was not freed in above loop. */
337
272
         if (chunk > 0) {
338
273
            offset += subject->length - DBSEQ_CHUNK_OVERLAP;
339
274
            if (kNucleotide) {
391
326
               subject->length = prot_length;
392
327
         } else {
393
328
            BLAST_GetUngappedHSPList(init_hitlist, query_info, subject, 
394
 
                                     hit_options, &hsp_list);
 
329
                                     hit_params->options, &hsp_list);
395
330
         }
396
331
 
397
332
         if (hsp_list->hspcnt == 0)
400
335
         /* The subject ordinal id is not yet filled in this HSP list */
401
336
         hsp_list->oid = subject->oid;
402
337
 
 
338
         /* check for interrupt */
 
339
         if (interrupt_search && (*interrupt_search)(progress_info) == TRUE) {
 
340
            combined_hsp_list = Blast_HSPListFree(combined_hsp_list);
 
341
            BlastInitHitListReset(init_hitlist);
 
342
            status = BLASTERR_INTERRUPTED;
 
343
            break;
 
344
         }
 
345
 
403
346
         Blast_HSPListAdjustOffsets(hsp_list, offset);
404
347
         /* Allow merging of HSPs either if traceback is already 
405
348
            available, or if it is an ungapped search */
406
 
         Blast_HSPListsMerge(hsp_list, &combined_hsp_list, hsp_num_max, offset,
407
 
            (Boolean)(prelim_traceback || !score_options->gapped_calculation));
 
349
         status = Blast_HSPListsMerge(&hsp_list, &combined_hsp_list,  kHspNumMax, offset,
 
350
            (Boolean)(kPrelimTraceback || !score_options->gapped_calculation));
408
351
      } /* End loop on chunks of subject sequence */
409
352
 
410
 
      if (Blast_HSPListAppend(combined_hsp_list, hsp_list_out, hsp_num_max)) {
 
353
      hsp_list = Blast_HSPListFree(hsp_list);  /* In case this was not freed in above loop. */
 
354
 
 
355
      *hsp_list_out_ptr = combined_hsp_list;
 
356
 
 
357
      return status;
 
358
}
 
359
 
 
360
/** Clean up function for s_BlastSearchEngineCore
 
361
 * @param program_number BLAST program type [in]
 
362
 * @param query_info Query information structure local to
 
363
 * s_BlastSearchEngineCore, which may or may not be deallocated [in]
 
364
 * @param query_info_in Query information [in]
 
365
 * @param translation_buffer buffer containing translated sequence data [in]
 
366
 * @param frame_offsets_a FIXME
 
367
 */
 
368
static void
 
369
s_BlastSearchEngineCoreCleanUp(EBlastProgramType program_number,
 
370
             BlastQueryInfo* query_info,
 
371
             const BlastQueryInfo* query_info_in,
 
372
             Uint1* translation_buffer,
 
373
             Int4* frame_offsets_a)
 
374
{
 
375
   /* Free the local query info structure when needed (in RPS BLAST). */
 
376
   if (query_info != query_info_in)
 
377
      BlastQueryInfoFree(query_info);
 
378
 
 
379
   /* Free translation buffer and frame offsets, except for RPS tblastn,
 
380
    * where they are taken from different structures, and hence shouldn't 
 
381
    * be freed here. 
 
382
    */
 
383
   if (program_number != eBlastTypeRpsTblastn) {
 
384
      if (translation_buffer) {
 
385
         sfree(translation_buffer);
 
386
      }
 
387
   }
 
388
   
 
389
   if (frame_offsets_a) {
 
390
       sfree(frame_offsets_a);
 
391
   }
 
392
}
 
393
 
 
394
/** The core of the BLAST search: comparison between the (concatenated)
 
395
 * query against one subject sequence. Translation of the subject sequence
 
396
 * into 6 frames is done inside, if necessary. If subject sequence is 
 
397
 * too long, it can be split into several chunks. 
 
398
 * @param program_number BLAST program type [in]
 
399
 * @param query Query sequence structure [in]
 
400
 * @param query_info_in Query information [in]
 
401
 * @param subject Subject sequence structure [in]
 
402
 * @param lookup Lookup table [in]
 
403
 * @param gap_align Structure for gapped alignment information [in]
 
404
 * @param score_params Scoring parameters [in]
 
405
 * @param word_params Initial word finding and ungapped extension 
 
406
 *                    parameters [in]
 
407
 * @param ext_params Gapped extension parameters [in]
 
408
 * @param hit_params Hit saving parameters [in]
 
409
 * @param db_options Database options [in]
 
410
 * @param diagnostics Hit counts and other diagnostics [in] [out]
 
411
 * @param aux_struct Structure containing different auxiliary data and memory
 
412
 *                   for the preliminary stage of the BLAST search [in]
 
413
 * @param hsp_list_out_ptr List of HSPs found for a given subject sequence [in]
 
414
 * @param interrupt_search function callback to allow interruption of BLAST
 
415
 *                   search [in, optional]
 
416
 * @param progress_info contains information about the progress of the current
 
417
 *                   BLAST search [in|out]
 
418
 */
 
419
static Int2
 
420
s_BlastSearchEngineCore(EBlastProgramType program_number, BLAST_SequenceBlk* query, 
 
421
   BlastQueryInfo* query_info_in, BLAST_SequenceBlk* subject, 
 
422
   LookupTableWrap* lookup, BlastGapAlignStruct* gap_align, 
 
423
   const BlastScoringParameters* score_params, 
 
424
   const BlastInitialWordParameters* word_params, 
 
425
   const BlastExtensionParameters* ext_params, 
 
426
   const BlastHitSavingParameters* hit_params, 
 
427
   const BlastDatabaseOptions* db_options,
 
428
   BlastDiagnostics* diagnostics,
 
429
   BlastCoreAuxStruct* aux_struct,
 
430
   BlastHSPList** hsp_list_out_ptr,
 
431
   TInterruptFnPtr interrupt_search, 
 
432
   SBlastProgress* progress_info)
 
433
{
 
434
   BlastHSPList* hsp_list_out=NULL;
 
435
   Uint1* translation_buffer = NULL;
 
436
   Int4* frame_offsets   = NULL;
 
437
   Int4* frame_offsets_a = NULL; /* Will be freed if non-null */
 
438
   BlastHitSavingOptions* hit_options = hit_params->options;
 
439
   BlastScoringOptions* score_options = score_params->options;
 
440
   Int2 status = 0;
 
441
   Uint4 context, first_context, last_context;
 
442
   Int4 orig_length = subject->length;
 
443
   Uint1* orig_sequence = subject->sequence;
 
444
   BlastQueryInfo* query_info = query_info_in;
 
445
 
 
446
   const Boolean kTranslatedSubject = 
 
447
        (Blast_SubjectIsTranslated(program_number) || program_number == eBlastTypeRpsTblastn);
 
448
   const Boolean kNucleotide = (program_number == eBlastTypeBlastn ||
 
449
                                program_number == eBlastTypePhiBlastn);
 
450
   const int kHspNumMax = BlastHspNumMax(score_options->gapped_calculation, hit_options);
 
451
 
 
452
   *hsp_list_out_ptr = NULL;
 
453
 
 
454
   if (kTranslatedSubject) {
 
455
      first_context = 0;
 
456
      last_context = 5;
 
457
      if (score_options->is_ooframe) {
 
458
         BLAST_GetAllTranslations(orig_sequence, eBlastEncodingNcbi2na,
 
459
            orig_length, db_options->gen_code_string, &translation_buffer,
 
460
            &frame_offsets, &subject->oof_sequence);
 
461
         subject->oof_sequence_allocated = TRUE;
 
462
         frame_offsets_a = frame_offsets;
 
463
      } else if (program_number == eBlastTypeRpsTblastn ) {
 
464
          /* For RPS tblastn, subject is actually query, which has already 
 
465
             been translated during the setup stage. */
 
466
          translation_buffer = orig_sequence - 1;
 
467
          frame_offsets_a = frame_offsets =
 
468
              ContextOffsetsToOffsetArray(query_info_in);
 
469
      } else {
 
470
         BLAST_GetAllTranslations(orig_sequence, eBlastEncodingNcbi2na,
 
471
            orig_length, db_options->gen_code_string, &translation_buffer,
 
472
            &frame_offsets, NULL);
 
473
         frame_offsets_a = frame_offsets;
 
474
      }
 
475
   } else if (kNucleotide) {
 
476
      first_context = 1;
 
477
      last_context = 1;
 
478
   } else {
 
479
      first_context = 0;
 
480
      last_context = 0;
 
481
   }
 
482
 
 
483
 
 
484
   /* Substitute query info by concatenated database info for RPS BLAST search */
 
485
   if (Blast_ProgramIsRpsBlast(program_number)) {
 
486
      BlastRPSLookupTable* lut = (BlastRPSLookupTable*) lookup->lut;
 
487
      query_info = BlastQueryInfoNew(eBlastTypeRpsBlast, lut->num_profiles);
 
488
      /* Since this will really be "subject info", not "query info",
 
489
         pass program argument such that all frames will be set to 0. */
 
490
      s_RPSOffsetArrayToContextOffsets(query_info, lut->rps_seq_offsets);
 
491
   }
 
492
 
 
493
   /* Loop over frames of the subject sequence */
 
494
   for (context=first_context; context<=last_context; context++) {
 
495
      BlastHSPList* hsp_list_for_chunks = NULL;
 
496
      if (kTranslatedSubject) {
 
497
         subject->frame =
 
498
             BLAST_ContextToFrame(eBlastTypeBlastx, context);
 
499
         subject->sequence = 
 
500
            translation_buffer + frame_offsets[context] + 1;
 
501
         subject->length = 
 
502
           frame_offsets[context+1] - frame_offsets[context] - 1;
 
503
      } else {
 
504
         subject->frame = context;
 
505
      }
 
506
      status = s_BlastSearchEngineOneContext(program_number, query, query_info, 
 
507
                                             subject, orig_length, lookup, 
 
508
                                             gap_align, score_params, 
 
509
                                             word_params, ext_params, 
 
510
                                             hit_params, diagnostics, 
 
511
                                             aux_struct, &hsp_list_for_chunks,
 
512
                                             interrupt_search, progress_info);
 
513
      if (status != 0) {
 
514
          break;
 
515
      }
 
516
     
 
517
      if (Blast_HSPListAppend(&hsp_list_for_chunks, &hsp_list_out, kHspNumMax)) {
411
518
         status = 1;
412
519
         break;
413
520
      }
414
 
      combined_hsp_list = Blast_HSPListFree(combined_hsp_list);
 
521
      
 
522
      /* if searching was interrupted, delete accumulated results
 
523
         but continue execution so temporary structures get freed */
 
524
      if (interrupt_search && (*interrupt_search)(progress_info) == TRUE) {
 
525
         status = BLASTERR_INTERRUPTED;
 
526
         break;
 
527
      }
415
528
   } /* End loop on frames */
416
529
 
417
530
   /* Restore the original contents of the subject block */
418
531
   subject->length = orig_length;
419
532
   subject->sequence = orig_sequence;
420
533
 
421
 
   if (status) 
422
 
      return status;
423
 
 
424
 
   hsp_list = *hsp_list_out;
 
534
   if (status) {
 
535
       hsp_list_out = Blast_HSPListFree(hsp_list_out);
 
536
       s_BlastSearchEngineCoreCleanUp(program_number, query_info, 
 
537
                                      query_info_in, translation_buffer, 
 
538
                                      frame_offsets_a);
 
539
       return status;
 
540
   }
425
541
 
426
542
   if (hit_params->link_hsp_params) {
427
 
      status = BLAST_LinkHsps(program_number, hsp_list, query_info,
 
543
      status = BLAST_LinkHsps(program_number, hsp_list_out, query_info,
428
544
                  subject->length, gap_align->sbp, hit_params->link_hsp_params, 
429
545
                  score_options->gapped_calculation);
430
546
   } else if (!Blast_ProgramIsPhiBlast(program_number) &&
432
548
       /* Calculate e-values for all HSPs. Skip this step
433
549
          for PHI and RPS BLAST, since calculating the E values 
434
550
          requires precomputation that has not been done yet */
435
 
       status = Blast_HSPListGetEvalues(query_info, hsp_list, 
 
551
       status = Blast_HSPListGetEvalues(query_info, hsp_list_out, 
436
552
                                        score_options->gapped_calculation, 
437
553
                                        gap_align->sbp, 0, 1.0);
438
554
   }
439
555
   
440
 
   /* Free the local query info structure when needed (in RPS BLAST). */
441
 
   if (query_info != query_info_in)
442
 
      BlastQueryInfoFree(query_info);
443
 
 
444
556
   /* Discard HSPs that don't pass the e-value test. */
445
 
   status = Blast_HSPListReapByEvalue(hsp_list, hit_options);
 
557
   status = Blast_HSPListReapByEvalue(hsp_list_out, hit_options);
446
558
 
447
559
   /* If there are no HSPs left, destroy the HSP list too. */
448
 
   if (hsp_list && hsp_list->hspcnt == 0)
449
 
      *hsp_list_out = hsp_list = Blast_HSPListFree(hsp_list);
 
560
   if (hsp_list_out && hsp_list_out->hspcnt == 0)
 
561
      *hsp_list_out_ptr = hsp_list_out = Blast_HSPListFree(hsp_list_out);
450
562
 
451
 
   if (gapped_stats && hsp_list && hsp_list->hspcnt > 0) {
 
563
   if (diagnostics && diagnostics->gapped_stat && hsp_list_out && hsp_list_out->hspcnt > 0) {
 
564
      BlastGappedStats* gapped_stats = diagnostics->gapped_stat;
452
565
      ++gapped_stats->num_seqs_passed;
453
 
      gapped_stats->good_extensions += hsp_list->hspcnt;
454
 
   }
455
 
 
456
 
   /* Free translation buffer and frame offsets, except for RPS tblastn,
457
 
    * where they are taken from different structures, and hence shouldn't 
458
 
    * be freed here. 
459
 
    */
460
 
   if (program_number != eBlastTypeRpsTblastn) {
461
 
      if (translation_buffer) {
462
 
         sfree(translation_buffer);
463
 
      }
464
 
   }
465
 
   
466
 
   if (frame_offsets_a) {
467
 
       sfree(frame_offsets_a);
468
 
   }
469
 
   
 
566
      gapped_stats->good_extensions += hsp_list_out->hspcnt;
 
567
   }
 
568
 
 
569
   s_BlastSearchEngineCoreCleanUp(program_number, query_info, query_info_in,
 
570
                                  translation_buffer, frame_offsets_a);
 
571
   
 
572
   *hsp_list_out_ptr = hsp_list_out;
 
573
 
470
574
   return status;
471
575
}
472
576
 
564
668
   /* Pick which gapped alignment algorithm to use. */
565
669
   if (phi_lookup)
566
670
      aux_struct->GetGappedScore = PHIGetGappedScore;
567
 
   else if (ext_options->ePrelimGapExt == eDynProgExt)
 
671
   else 
568
672
      aux_struct->GetGappedScore = BLAST_GetGappedScore;
569
 
   else
570
 
      aux_struct->GetGappedScore = BLAST_MbGetGappedScore;
571
673
 
572
 
   aux_struct->hsp_list = Blast_HSPListNew(hit_options->hsp_num_max);
573
674
   return status;
574
675
}
575
676
 
592
693
 * @param diagnostics Return statistics containing numbers of hits on 
593
694
 *                    different stages of the search. Statistics saved only 
594
695
 *                    for the allocated parts of the structure. [in] [out]
 
696
 * @param interrupt_search function callback to allow interruption of BLAST
 
697
 *                   search [in, optional]
 
698
 * @param progress_info contains information about the progress of the current
 
699
 *                   BLAST search [in|out]
595
700
 */
596
701
static Int2 
597
702
s_RPSPreliminarySearchEngine(EBlastProgramType program_number, 
603
708
   const BlastExtensionParameters* ext_params, 
604
709
   BlastGapAlignStruct* gap_align,
605
710
   const BlastHitSavingParameters* hit_params,
606
 
   BlastHSPStream* hsp_stream, BlastDiagnostics* diagnostics)
 
711
   BlastHSPStream* hsp_stream, BlastDiagnostics* diagnostics,
 
712
   TInterruptFnPtr interrupt_search, SBlastProgress* progress_info)
607
713
{
608
714
   BlastHSPList* hsp_list = NULL;
609
715
   Int2 status = 0;
673
779
          s_BlastSearchEngineCore(program_number, &concat_db, one_query_info, 
674
780
             one_query, lookup_wrap, gap_align, score_params, 
675
781
             word_params, ext_params, hit_params, NULL, 
676
 
             diagnostics, aux_struct, &hsp_list);
 
782
             diagnostics, aux_struct, &hsp_list, interrupt_search, 
 
783
             progress_info);
 
784
 
 
785
       if (interrupt_search && (*interrupt_search)(progress_info) == TRUE) {
 
786
           hsp_list = Blast_HSPListFree(hsp_list);
 
787
           status = BLASTERR_INTERRUPTED;
 
788
           break;
 
789
       }
677
790
 
678
791
       /* Save the resulting list of HSPs. 'query' and 'subject' are
679
792
          still reversed */
712
825
   BlastEffectiveLengthsParameters* eff_len_params,
713
826
   const PSIBlastOptions* psi_options, 
714
827
   const BlastDatabaseOptions* db_options,
715
 
   BlastHSPStream* hsp_stream, BlastDiagnostics* diagnostics)
 
828
   BlastHSPStream* hsp_stream, BlastDiagnostics* diagnostics,
 
829
   TInterruptFnPtr interrupt_search, SBlastProgress* progress_info)
716
830
{
717
831
   BlastCoreAuxStruct* aux_struct = NULL;
718
832
   BlastHSPList* hsp_list = NULL; 
739
853
          ext_options, hit_options, query, &aux_struct)) != 0)
740
854
      return status;
741
855
 
 
856
   /* remember the current search state */
 
857
   if (progress_info)
 
858
       progress_info->stage = ePrelimSearch;
 
859
 
742
860
   /* For RPS BLAST, there is no loop over subject sequences, so the preliminary
743
861
      search engine is done in a separate function. */
744
862
   if (Blast_ProgramIsRpsBlast(program_number)) {
745
863
      status =         
746
864
         s_RPSPreliminarySearchEngine(program_number, query, query_info, 
747
865
            seq_src, score_params, lookup_wrap, aux_struct, word_params, 
748
 
            ext_params, gap_align, hit_params, hsp_stream, diagnostics);
 
866
            ext_params, gap_align, hit_params, hsp_stream, diagnostics,
 
867
            interrupt_search, progress_info);
749
868
      word_params = BlastInitialWordParametersFree(word_params);
750
869
      s_BlastCoreAuxStructFree(aux_struct);
751
870
      return status;
764
883
 
765
884
   db_length = BlastSeqSrcGetTotLen(seq_src);
766
885
 
767
 
   itr = BlastSeqSrcIteratorNew();
 
886
   itr = BlastSeqSrcIteratorNewEx(MAX(BlastSeqSrcGetNumSeqs(seq_src)/100,1));
768
887
 
769
888
   /* iterate over all subject sequences */
770
889
   while ( (seq_arg.oid = BlastSeqSrcIteratorNext(seq_src, itr)) 
798
917
         s_BlastSearchEngineCore(program_number, query, query_info,
799
918
            seq_arg.seq, lookup_wrap, gap_align, score_params, word_params, 
800
919
            ext_params, hit_params, db_options, diagnostics, aux_struct, 
801
 
            &hsp_list);
802
 
      if (status)
803
 
         break;
 
920
            &hsp_list, interrupt_search, progress_info);
 
921
      if (status) {
 
922
          break;
 
923
      }
804
924
 
805
925
      if (hsp_list && hsp_list->hspcnt > 0) {
806
926
         if (!gapped_calculation || prelim_traceback) {
813
933
                     hsp_list, query, seq_arg.seq, word_params, hit_params, 
814
934
                     query_info, sbp, score_params, seq_src, 
815
935
                     (db_options ? db_options->gen_code_string : NULL));
 
936
               if (status) {
 
937
                  BlastSeqSrcReleaseSequence(seq_src, (void*) &seq_arg);
 
938
                  return status;
 
939
               }
816
940
               /* Relink HSPs if sum statistics is used, because scores might
817
941
                * have changed after reevaluation with ambiguities, and there
818
942
                * will be no traceback stage where relinking is done normally.
844
968
      }
845
969
      
846
970
      BlastSeqSrcReleaseSequence(seq_src, (void*) &seq_arg);
 
971
 
 
972
      /* check for interrupt */
 
973
      if (interrupt_search && (*interrupt_search)(progress_info) == TRUE) {
 
974
          status = BLASTERR_INTERRUPTED;
 
975
          break;
 
976
      }
847
977
   }
848
978
   
849
979
   BlastSequenceBlkFree(seq_arg.seq);
904
1034
                                      lookup_wrap, word_options, 
905
1035
                                      ext_params, hit_params, eff_len_params,
906
1036
                                      psi_options, db_options, hsp_stream, 
907
 
                                      local_diagnostics)) != 0)
 
1037
                                      local_diagnostics, 0, 0)) != 0)
908
1038
      return status;
909
1039
 
910
1040
   /* Do not destruct score block here */
923
1053
   return status;
924
1054
}
925
1055
 
 
1056
/** Function to deallocate data structures allocated in Blast_RunFullSearch */
 
1057
static void
 
1058
s_BlastRunFullSearchCleanUp(BlastGapAlignStruct* gap_align,
 
1059
                            BlastScoringParameters* score_params,
 
1060
                            BlastExtensionParameters* ext_params,
 
1061
                            BlastHitSavingParameters* hit_params,
 
1062
                            BlastEffectiveLengthsParameters* eff_len_params)
 
1063
{
 
1064
    /* Do not destruct score block here */
 
1065
    gap_align->sbp = NULL;
 
1066
    BLAST_GapAlignStructFree(gap_align);
 
1067
 
 
1068
    BlastScoringParametersFree(score_params);
 
1069
    BlastHitSavingParametersFree(hit_params);
 
1070
    BlastExtensionParametersFree(ext_params);
 
1071
    BlastEffectiveLengthsParametersFree(eff_len_params);
 
1072
}
 
1073
 
926
1074
Int4 
927
1075
Blast_RunFullSearch(EBlastProgramType program_number, 
928
1076
   BLAST_SequenceBlk* query, BlastQueryInfo* query_info,
936
1084
   const PSIBlastOptions* psi_options, 
937
1085
   const BlastDatabaseOptions* db_options,
938
1086
   BlastHSPStream* hsp_stream, const BlastRPSInfo* rps_info, 
939
 
   BlastDiagnostics* diagnostics, BlastHSPResults** results)
 
1087
   BlastDiagnostics* diagnostics, BlastHSPResults** results,
 
1088
   TInterruptFnPtr interrupt_search,
 
1089
   SBlastProgress* progress_info)
940
1090
{
941
1091
   Int4 status = 0;
942
1092
   BlastScoringParameters* score_params = NULL;
943
1093
   BlastExtensionParameters* ext_params = NULL;
944
1094
   BlastHitSavingParameters* hit_params = NULL;
945
1095
   BlastEffectiveLengthsParameters* eff_len_params = NULL;
946
 
   BlastGapAlignStruct* gap_align;
 
1096
   BlastGapAlignStruct* gap_align = NULL;
947
1097
   SPHIPatternSearchBlk* pattern_blk = NULL;
948
1098
 
949
1099
   if ((status = 
950
1100
        BLAST_GapAlignSetUp(program_number, seq_src, score_options, 
951
1101
           eff_len_options, ext_options, hit_options, query_info, sbp, 
952
1102
           &score_params, &ext_params, &hit_params, &eff_len_params, 
953
 
           &gap_align)) != 0)
954
 
      return status;
 
1103
           &gap_align)) != 0) {
 
1104
       s_BlastRunFullSearchCleanUp(gap_align, score_params, ext_params, 
 
1105
                                   hit_params, eff_len_params);
 
1106
       return status;
 
1107
   }
955
1108
      
956
1109
   if ((status=
957
1110
        BLAST_PreliminarySearchEngine(program_number, query, query_info, 
958
1111
           seq_src, gap_align, score_params, lookup_wrap, word_options, 
959
1112
           ext_params, hit_params, eff_len_params, psi_options, 
960
 
           db_options, hsp_stream, diagnostics)) != 0)
961
 
      return status;
 
1113
           db_options, hsp_stream, diagnostics, interrupt_search, 
 
1114
           progress_info)) != 0) {
 
1115
       s_BlastRunFullSearchCleanUp(gap_align, score_params, ext_params, 
 
1116
                                   hit_params, eff_len_params);
 
1117
       return status;
 
1118
   }
962
1119
   
963
1120
   /* Prohibit any subsequent writing to the HSP stream. */
964
1121
   BlastHSPStreamClose(hsp_stream);
972
1129
        BLAST_ComputeTraceback(program_number, hsp_stream, query, query_info,
973
1130
                               seq_src, gap_align, score_params, ext_params, 
974
1131
                               hit_params, eff_len_params, db_options, 
975
 
                               psi_options, rps_info, pattern_blk, results))
976
 
       != 0)
 
1132
                               psi_options, rps_info, pattern_blk, results,
 
1133
                               interrupt_search, progress_info))
 
1134
       != 0) {
 
1135
       s_BlastRunFullSearchCleanUp(gap_align, score_params, ext_params, 
 
1136
                                   hit_params, eff_len_params);
977
1137
       return status;
978
 
 
979
 
   /* Do not destruct score block here */
980
 
   gap_align->sbp = NULL;
981
 
   BLAST_GapAlignStructFree(gap_align);
982
 
 
983
 
   score_params = BlastScoringParametersFree(score_params);
984
 
   hit_params = BlastHitSavingParametersFree(hit_params);
985
 
   ext_params = BlastExtensionParametersFree(ext_params);
986
 
   eff_len_params = BlastEffectiveLengthsParametersFree(eff_len_params);
 
1138
   }
 
1139
 
 
1140
   s_BlastRunFullSearchCleanUp(gap_align, score_params, ext_params, hit_params,
 
1141
                               eff_len_params);
987
1142
   return status;
988
1143
}