1
/* $Id: blast_traceback.c,v 1.209 2009/01/12 14:19:33 kazimird Exp $
1
/* $Id: blast_traceback.c,v 1.212 2009/07/09 15:49:31 kazimird Exp $
2
2
* ===========================================================================
4
4
* PUBLIC DOMAIN NOTICE
51
51
#ifndef SKIP_DOXYGEN_PROCESSING
52
52
static char const rcsid[] =
53
"$Id: blast_traceback.c,v 1.209 2009/01/12 14:19:33 kazimird Exp $";
53
"$Id: blast_traceback.c,v 1.212 2009/07/09 15:49:31 kazimird Exp $";
54
54
#endif /* SKIP_DOXYGEN_PROCESSING */
56
56
#include <algo/blast/core/blast_traceback.h>
359
359
const Boolean kGreedyTraceback = (ext_options->eTbackExt == eGreedyTbck);
360
360
const Boolean kTranslateSubject =
361
361
(Blast_SubjectIsTranslated(program_number) || program_number == eBlastTypeRpsTblastn);
362
const Boolean kFullTranslation = (fence_hit && *fence_hit);
362
363
const Boolean kSmithWaterman = (ext_options->eTbackExt ==
363
364
eSmithWatermanTbckFull);
364
365
BlastQueryInfo* query_info = query_info_in;
479
481
subject = subject_blk->oof_sequence + CODON_LENGTH;
480
482
if (hsp->subject.frame < 0)
481
483
subject += subject_length + 1;
484
} else if (kFullTranslation) {
485
/* previous traceback has hit the fence, will do full translation */
486
BlastHSP* hsp_new = Blast_HSPNew();
487
hsp_new->subject.frame = hsp->subject.frame;
488
hsp_new->subject.offset = -1;
489
subject = Blast_HSPGetTargetTranslation(target_t, hsp_new, &subject_length);
490
hsp_new = Blast_HSPFree(hsp_new);
483
subject = Blast_HSPGetTargetTranslation(target_t, hsp, &subject_length);
492
subject = Blast_HSPGetTargetTranslation(target_t, hsp, &subject_length);
486
495
if (!kIsOutOfFrame && (((hsp->query.gapped_start == 0 &&
560
569
adjusted_subject, gap_align, score_params, q_start, s_start,
561
570
query_length, adjusted_s_length,
572
ASSERT(!(kFullTranslation && *fence_hit));
565
575
fence_error = (fence_hit && *fence_hit);
1150
1163
interrupt_search, progress_info);
1151
1164
} else if (ext_params->options->compositionBasedStats > 0 ||
1152
1165
ext_params->options->eTbackExt == eSmithWatermanTbck) {
1166
/* FIXME partial sequence fetching/translation could lead to fence hit
1154
1169
Blast_RedoAlignmentCore(program_number, query, query_info, sbp,
1155
1170
hsp_stream, seq_src, default_db_genetic_code,
1163
1178
const Boolean kPhiBlast = Blast_ProgramIsPhiBlast(program_number);
1164
1179
BlastHSPStreamResultBatch *batch =
1165
1180
Blast_HSPStreamResultBatchInit(query_info->num_queries);
1166
const Boolean kSubjectRanges = !(hit_params->restricted_align ||
1167
score_params->options->is_ooframe ||
1168
ext_params->options->eTbackExt ==
1169
eSmithWatermanTbck);
1171
1182
memset((void*) &seq_arg, 0, sizeof(seq_arg));
1187
1198
if (perform_traceback) {
1188
1199
seq_arg.oid = batch->hsplist_array[0]->oid;
1189
1200
seq_arg.encoding = encoding;
1190
seq_arg.enable_ranges = kSubjectRanges;
1191
1201
seq_arg.check_oid_exclusion = TRUE;
1202
seq_arg.reset_ranges = FALSE;
1193
1204
BlastSequenceBlkClean(seq_arg.seq);
1194
1205
if (BlastSeqSrcGetSequence(seq_src, &seq_arg) < 0) {
1242
1253
query_info, pattern_blk);
1244
1255
Boolean fence_hit = FALSE;
1245
Boolean * fence_hit_ptr = NULL;
1247
/* prepare to deal with partial subject sequence
1248
ranges, if these are configured and a previous
1249
HSP list has not turned them off */
1250
if (seq_arg.enable_ranges == TRUE)
1251
fence_hit_ptr = &fence_hit;
1253
1256
Blast_TracebackFromHSPList(program_number, hsp_list, query,
1254
1257
seq_arg.seq, query_info,
1255
1258
gap_align, sbp, score_params,
1256
1259
ext_params->options, hit_params,
1257
1260
seq_arg.seq->gen_code_string,
1260
1263
if (fence_hit) {
1261
1264
/* Disable range support and refetch the
1262
1265
(whole) subject sequence */
1264
seq_arg.enable_ranges = FALSE;
1267
seq_arg.reset_ranges = TRUE;
1265
1268
BlastSeqSrcReleaseSequence(seq_src, &seq_arg);
1266
1269
BlastSeqSrcGetSequence(seq_src, &seq_arg);
1268
/* Retry the alignment */
1271
/* The C toolkit will erase genetic_code, so do it again */
1272
if (Blast_SubjectIsTranslated(program_number) &&
1273
seq_arg.seq->gen_code_string == NULL) {
1274
seq_arg.seq->gen_code_string =
1275
GenCodeSingletonFind(default_db_genetic_code);
1276
ASSERT(seq_arg.seq->gen_code_string);
1279
/* Retry the alignment with fence_hit set*/
1270
1280
Blast_TracebackFromHSPList(program_number, hsp_list,
1271
1281
query, seq_arg.seq,
1272
1282
query_info, gap_align,
1305
1316
} /* loop over all batches */
1307
1318
batch = Blast_HSPStreamResultBatchFree(batch);
1319
/* post-traceback pipes */
1320
BlastHSPStreamTBackClose(hsp_stream, results);
1308
1322
BlastSequenceBlkFree(seq_arg.seq);
1311
if (hit_params->options->culling_limit > 0)
1312
Blast_HSPResultsPerformCulling(results, query_info,
1313
hit_params->options->culling_limit,
1316
1325
/* Re-sort the hit lists according to their best e-values, because
1317
1326
they could have changed. Only do this for a database search. */
1318
1327
if (BlastSeqSrcGetTotLen(seq_src) > 0)