« back to all changes in this revision
Viewing changes to algo/blast/core/blast_hits.c
- browse files at revision 8
- compare with another revision
- download diff
- download tarball
- view history from revision 8
- 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
- files added:
- algo/blast/api/blast_message_api.c
- files removed:
- debian/libvibrant6.install
- files modified:
- VERSION
added
removed
algo/blast/core/blast_hits.c
1
/* $Id: blast_hits.c,v 1.173 2005/11/16 14:27:03 madden Exp $
1
/* $Id: blast_hits.c,v 1.188 2006/04/05 18:47:42 camacho Exp $
2
2
* ===========================================================================
3
3
*
4
4
* PUBLIC DOMAIN NOTICE
32
32
33
33
#ifndef SKIP_DOXYGEN_PROCESSING
34
34
static char const rcsid[] =
35
"$Id: blast_hits.c,v 1.173 2005/11/16 14:27:03 madden Exp $";
35
"$Id: blast_hits.c,v 1.188 2006/04/05 18:47:42 camacho Exp $";
36
36
#endif /* SKIP_DOXYGEN_PROCESSING */
37
37
38
38
#include <algo/blast/core/blast_options.h>
39
39
#include <algo/blast/core/blast_extend.h>
40
40
#include <algo/blast/core/blast_hits.h>
41
41
#include <algo/blast/core/blast_util.h>
42
#include <algo/blast/core/blast_hspstream.h>
42
43
#include "blast_hits_priv.h"
43
44
#include "blast_itree.h"
44
45
68
69
else
69
70
(*retval)->prelim_hitlist_size = hit_options->hitlist_size;
70
71
71
(*retval)->options = hit_options;
72
(*retval)->hsp_num_max = BlastHspNumMax(scoring_options->gapped_calculation, hit_options);
72
73
73
74
return 0;
74
75
}
78
79
{
79
80
if (param)
80
81
{
81
param->options = NULL;
82
82
sfree(param);
83
83
}
84
84
return NULL;
85
85
}
86
86
87
88
87
89
/********************************************************************************
88
90
Functions manipulating BlastHSP's
89
91
********************************************************************************/
148
150
return 0;
149
151
}
150
152
153
Int4 BlastHspNumMax(Boolean gapped_calculation, const BlastHitSavingOptions* options)
154
{
155
Int4 retval=0;
156
157
if (options->hsp_num_max <= 0)
158
{
159
if (!gapped_calculation || options->program_number == eBlastTypeTblastx)
160
retval = kUngappedHSPNumMax;
161
else
162
retval = INT4_MAX;
163
}
164
else
165
retval = options->hsp_num_max;
166
167
return retval;
168
}
169
151
170
/** Copies all contents of a BlastHSP structure. Used in PHI BLAST for splitting
152
171
* results corresponding to different pattern occurrences in query.
153
172
* @param hsp Original HSP [in]
168
187
new_hsp->num = hsp->num;
169
188
new_hsp->num_ident = hsp->num_ident;
170
189
new_hsp->bit_score = hsp->bit_score;
190
new_hsp->comp_adjustment_method = hsp->comp_adjustment_method;
171
191
if (hsp->gap_info) {
172
GapEditScript* edit_script = hsp->gap_info;
173
GapEditScript* new_edit_script;
174
new_hsp->gap_info = new_edit_script =
175
(GapEditScript*) BlastMemDup(edit_script, sizeof(GapEditScript));
176
/* Copy the linked list of edit scripts. */
177
for (edit_script = edit_script->next; edit_script;
178
edit_script = edit_script->next) {
179
new_edit_script->next =
180
(GapEditScript*) BlastMemDup(edit_script, sizeof(GapEditScript));
181
new_edit_script = new_edit_script->next;
182
}
192
new_hsp->gap_info = GapEditScriptDup(hsp->gap_info);
183
193
}
184
194
185
195
if (hsp->pat_info) {
228
238
* @param best_q_end Pointer to end of the new alignment in query [in]
229
239
* @param best_s_start Pointer to start of the new alignment in subject [in]
230
240
* @param best_s_end Pointer to end of the new alignment in subject [in]
231
* @param best_start_esp Link in the edit script chain where the new alignment
241
* @param best_start_esp_index index of the edit script array where the new alignment
232
242
* starts. [in]
233
* @param best_prev_esp Link in the edit script chain immediately preceding
234
* best_start_esp. [in]
235
* @param best_end_esp Link in the edit script chain where the new alignment
243
* @param best_end_esp_index index in the edit script array where the new alignment
236
244
* ends. [in]
237
245
* @param best_end_esp_num Number of edit operations in the last edit script,
238
246
* that are included in the alignment. [in]
244
252
Int4 score, Uint1* query_start, Uint1* subject_start,
245
253
Uint1* best_q_start, Uint1* best_q_end,
246
254
Uint1* best_s_start, Uint1* best_s_end,
247
GapEditScript* best_start_esp,
248
GapEditScript* best_prev_esp, GapEditScript* best_end_esp,
249
Int4 best_end_esp_num)
255
int best_start_esp_index,
256
int best_end_esp_index,
257
int best_end_esp_num)
250
258
{
251
259
Boolean delete_hsp = TRUE;
252
260
253
261
hsp->score = score;
254
262
255
/* Make corrections in edit block and free any parts that are no longer
256
needed */
257
if (gapped && best_start_esp && best_start_esp != hsp->gap_info) {
258
/* best_prev_esp is the link in the chain exactly preceding the starting
259
edit script of the best part of the alignment. If best alignment was
260
found, but it does not start from the original start, best_prev_esp
261
cannot be NULL. */
262
ASSERT (best_prev_esp);
263
/* Unlink the good part of the alignment from the previous
264
(negative-scoring) part that is being deleted. */
265
best_prev_esp->next = NULL;
266
GapEditScriptDelete(hsp->gap_info);
267
hsp->gap_info = best_start_esp;
268
}
269
270
263
if (hsp->score >= cutoff_score) {
271
264
/* Update all HSP offsets. */
272
265
hsp->query.offset = best_q_start - query_start;
275
268
hsp->subject.end = hsp->subject.offset + best_s_end - best_s_start;
276
269
277
270
if (gapped) {
278
if (best_end_esp->next != NULL) {
279
GapEditScriptDelete(best_end_esp->next);
280
best_end_esp->next = NULL;
271
int last_num=hsp->gap_info->size - 1;
272
if (best_end_esp_index != last_num|| best_start_esp_index > 0)
273
{
274
GapEditScript* esp_temp = GapEditScriptNew(best_end_esp_index-best_start_esp_index+1);
275
GapEditScriptPartialCopy(esp_temp, 0, hsp->gap_info, best_start_esp_index, best_end_esp_index);
276
hsp->gap_info = GapEditScriptDelete(hsp->gap_info);
277
hsp->gap_info = esp_temp;
281
278
}
282
best_end_esp->num = best_end_esp_num;
279
last_num = hsp->gap_info->size - 1;
280
hsp->gap_info->num[last_num] = best_end_esp_num;
281
ASSERT(best_end_esp_num >= 0);
283
282
}
284
283
delete_hsp = FALSE;
285
284
}
294
293
BlastScoreBlk* sbp)
295
294
{
296
295
Int4 sum, score, gap_open, gap_extend;
296
Int4 index; /* loop index */
297
297
298
GapEditScript* best_start_esp; /* Starting edit script for the best scoring
299
piece of the alignment. */
300
GapEditScript* best_end_esp; /* Ending edit script for the best scoring piece
301
of the alignment. */
302
GapEditScript* best_prev_esp; /* Previous link in the edit script chain
303
before best_end_esp. */
304
Int4 best_end_esp_num; /* Number of operations inside the ending edit script
305
for the best scoring piece. */
298
int best_start_esp_index = 0;
299
int best_end_esp_index = 0;
300
int current_start_esp_index = 0;
301
int best_end_esp_num = 0;
302
GapEditScript* esp; /* Used to hold GapEditScript of hsp->gap_info */
306
303
307
304
Uint1* best_q_start; /* Start of the best scoring part in query. */
308
305
Uint1* best_s_start; /* Start of the best scoring part in subject. */
309
306
Uint1* best_q_end; /* End of the best scoring part in query. */
310
307
Uint1* best_s_end; /* End of the best scoring part in subject. */
311
308
312
GapEditScript* current_start_esp; /* Starting edit script for the current
313
part of the alignment. */
314
GapEditScript* current_esp; /* Ending edit script for the current part
315
of the alignment. */
316
GapEditScript* current_prev_esp; /* Previous link in the edit script chain
317
before current_esp. */
318
GapEditScript* current_start_prev_esp; /* Previous link in the edit script
319
chain before current_start_esp. */
320
309
321
310
Uint1* current_q_start; /* Start of the current part of the alignment in
322
311
query. */
323
312
Uint1* current_s_start; /* Start of the current part of the alignment in
324
313
subject. */
325
Int4 op_index; /* Index of an operation within a single edit script. */
326
314
327
315
Uint1* query,* subject;
328
316
Int4** matrix;
352
340
/* Point all pointers to the beginning of the alignment. */
353
341
best_q_start = best_q_end = current_q_start = query;
354
342
best_s_start = best_s_end = current_s_start = subject;
355
best_start_esp = best_end_esp = current_start_esp = current_esp =
356
hsp->gap_info;
357
343
/* There are no previous edit scripts at the beginning. */
358
best_prev_esp = current_prev_esp = current_start_prev_esp = NULL;
359
best_end_esp_num = 0;
360
op_index = 0;
361
344
362
while (current_esp) {
363
/* Process substitutions one operation at a time, full gaps in one
364
step. */
365
if (current_esp->op_type == eGapAlignSub) {
366
sum += factor*matrix[*query & kResidueMask][*subject];
367
query++;
368
subject++;
369
op_index++;
370
} else if (current_esp->op_type == eGapAlignDel) {
371
sum -= gap_open + gap_extend * current_esp->num;
372
subject += current_esp->num;
373
op_index += current_esp->num;
374
} else if (current_esp->op_type == eGapAlignIns) {
375
sum -= gap_open + gap_extend * current_esp->num;
376
query += current_esp->num;
377
op_index += current_esp->num;
378
}
345
best_end_esp_num = -1;
346
esp = hsp->gap_info;
347
for (index=0; index<esp->size; index++)
348
{
349
int op_index = 0; /* Index of an operation within a single edit script. */
350
for (op_index=0; op_index<esp->num[index]; )
351
{
352
/* Process substitutions one operation at a time, full gaps in one step. */
353
if (esp->op_type[index] == eGapAlignSub) {
354
sum += factor*matrix[*query & kResidueMask][*subject];
355
query++;
356
subject++;
357
op_index++;
358
} else if (esp->op_type[index] == eGapAlignDel) {
359
sum -= gap_open + gap_extend * esp->num[index];
360
subject += esp->num[index];
361
op_index += esp->num[index];
362
} else if (esp->op_type[index] == eGapAlignIns) {
363
sum -= gap_open + gap_extend * esp->num[index];
364
query += esp->num[index];
365
op_index += esp->num[index];
366
}
379
367
380
if (sum < 0) {
368
if (sum < 0) {
381
369
/* Point current edit script chain start to the new place.
382
370
If we are in the middle of an edit script, reduce its length and
383
371
point operation index to the beginning of a modified edit script;
384
372
if we are at the end, move to the next edit script. */
385
if (op_index < current_esp->num) {
386
current_esp->num -= op_index;
387
current_start_esp = current_esp;
388
current_start_prev_esp = current_prev_esp;
389
op_index = 0;
390
} else {
391
current_start_esp = current_esp->next;
392
current_start_prev_esp = current_esp;
393
}
394
/* Set sum to 0, to start a fresh count. */
395
sum = 0;
396
/* Set current starting positions in sequences to the new start. */
397
current_q_start = query;
398
current_s_start = subject;
373
if (op_index < esp->num[index]) {
374
esp->num[index] -= op_index;
375
current_start_esp_index = index;
376
op_index = 0;
377
} else {
378
current_start_esp_index = index + 1;
379
}
380
/* Set sum to 0, to start a fresh count. */
381
sum = 0;
382
/* Set current starting positions in sequences to the new start. */
383
current_q_start = query;
384
current_s_start = subject;
399
385
400
/* If score has passed the cutoff at some point, leave the best score
401
and edit scripts positions information untouched, otherwise reset
402
the best score to 0 and point the best edit script positions to
403
the new start. */
404
if (score < hit_params->cutoff_score) {
405
/* Start from new offset; discard all previous information. */
406
best_q_start = query;
407
best_s_start = subject;
408
score = 0;
386
/* If score has passed the cutoff at some point, leave the best score
387
and edit scripts positions information untouched, otherwise reset
388
the best score to 0 and point the best edit script positions to
389
the new start. */
390
if (score < hit_params->cutoff_score) {
391
/* Start from new offset; discard all previous information. */
392
best_q_start = query;
393
best_s_start = subject;
394
score = 0;
409
395
410
/* Set best start and end edit script pointers to new start. */
411
best_start_esp = current_start_esp;
412
best_end_esp = current_start_esp;
413
best_prev_esp = current_prev_esp;
414
}
415
} else if (sum > score) {
416
/* Remember this point as the best scoring end point, and the current
396
/* Set best start and end edit script pointers to new start. */
397
best_start_esp_index = current_start_esp_index;
398
best_end_esp_index = current_start_esp_index;
399
}
400
/* break; */ /* start on next GapEditScript. */
401
} else if (sum > score) {
402
/* Remember this point as the best scoring end point, and the current
417
403
start of the alignment as the start of the best alignment. */
418
score = sum;
419
420
best_q_start = current_q_start;
421
best_s_start = current_s_start;
422
best_q_end = query;
423
best_s_end = subject;
424
425
best_start_esp = current_start_esp;
426
best_end_esp = current_esp;
427
best_prev_esp = current_start_prev_esp;
428
best_end_esp_num = op_index;
429
}
430
/* If operation index has reached the end of the current edit script,
431
move on to the next link in the edit script chain. */
432
if (op_index >= current_esp->num) {
433
op_index = 0;
434
current_prev_esp = current_esp;
435
current_esp = current_esp->next;
404
score = sum;
405
406
best_q_start = current_q_start;
407
best_s_start = current_s_start;
408
best_q_end = query;
409
best_s_end = subject;
410
411
best_start_esp_index = current_start_esp_index;
412
best_end_esp_index = index;
413
best_end_esp_num = op_index;
414
}
436
415
}
437
416
} /* loop on edit scripts */
438
417
443
422
s_UpdateReevaluatedHSP(hsp, TRUE, hit_params->cutoff_score, score,
444
423
query_start, subject_start, best_q_start,
445
424
best_q_end, best_s_start, best_s_end,
446
best_start_esp, best_prev_esp, best_end_esp,
425
best_start_esp_index, best_end_esp_index,
447
426
best_end_esp_num);
448
427
}
449
428
470
449
return
471
450
s_UpdateReevaluatedHSP(hsp, FALSE, cutoff_score, score, query_start,
472
451
subject_start, best_q_start, best_q_end,
473
best_s_start, best_s_end, NULL, NULL, NULL, 0);
452
best_s_start, best_s_end, 0, 0, 0);
474
453
}
475
454
476
455
Boolean
540
519
{
541
520
Int4 i, num_ident, align_length, q_off, s_off;
542
521
Uint1* q,* s;
543
GapEditScript* esp;
544
522
Int4 q_length = hsp->query.end - hsp->query.offset;
545
523
Int4 s_length = hsp->subject.end - hsp->subject.offset;
546
524
568
546
num_ident++;
569
547
}
570
548
} else {
571
for (esp = hsp->gap_info; esp; esp = esp->next) {
572
align_length += esp->num;
573
switch (esp->op_type) {
549
Int4 index;
550
GapEditScript* esp = hsp->gap_info;
551
for (index=0; index<esp->size; index++)
552
{
553
align_length += esp->num[index];
554
switch (esp->op_type[index]) {
574
555
case eGapAlignSub:
575
for (i=0; i<esp->num; i++) {
556
for (i=0; i<esp->num[index]; i++) {
576
557
if (*q++ == *s++)
577
558
num_ident++;
578
559
}
579
560
break;
580
561
case eGapAlignDel:
581
s += esp->num;
562
s += esp->num[index];
582
563
break;
583
564
case eGapAlignIns:
584
q += esp->num;
565
q += esp->num[index];
585
566
break;
586
567
default:
587
s += esp->num;
588
q += esp->num;
568
s += esp->num[index];
569
q += esp->num[index];
589
570
break;
590
571
}
591
572
}
601
582
BlastHSP* hsp, EBlastProgramType program,
602
583
Int4* num_ident_ptr, Int4* align_length_ptr)
603
584
{
604
Int4 i, num_ident, align_length;
585
Int4 num_ident, align_length;
586
Int4 index;
605
587
Uint1* q,* s;
606
588
GapEditScript* esp;
607
589
620
602
num_ident = 0;
621
603
align_length = 0;
622
604
623
for (esp = hsp->gap_info; esp; esp = esp->next) {
624
switch (esp->op_type) {
605
esp = hsp->gap_info;
606
for (index=0; index<esp->size; index++)
607
{
608
int i;
609
switch (esp->op_type[index]) {
625
610
case eGapAlignSub: /* Substitution */
626
align_length += esp->num;
627
for (i=0; i<esp->num; i++) {
611
align_length += esp->num[index];
612
for (i=0; i<esp->num[index]; i++) {
628
613
if (*q == *s)
629
614
num_ident++;
630
615
++q;
632
617
}
633
618
break;
634
619
case eGapAlignIns: /* Insertion */
635
align_length += esp->num;
636
s += esp->num * CODON_LENGTH;
620
align_length += esp->num[index];
621
s += esp->num[index] * CODON_LENGTH;
637
622
break;
638
623
case eGapAlignDel: /* Deletion */
639
align_length += esp->num;
640
q += esp->num;
624
align_length += esp->num[index];
625
q += esp->num[index];
641
626
break;
642
627
case eGapAlignDel2: /* Gap of two nucleotides. */
643
628
s -= 2;
652
637
s += 2;
653
638
break;
654
639
default:
655
s += esp->num * CODON_LENGTH;
656
q += esp->num;
640
s += esp->num[index] * CODON_LENGTH;
641
q += esp->num[index];
657
642
break;
658
643
}
659
644
}
700
685
701
686
/* Check whether this HSP passes the percent identity and minimal hit
702
687
length criteria, and delete it if it does not. */
703
if ((hsp->num_ident * 100 <
688
if ((hsp->num_ident * 100.0 <
704
689
align_length * hit_options->percent_identity) ||
705
690
align_length < hit_options->min_hit_length) {
706
691
delete_hsp = TRUE;
719
704
720
705
if (hsp->gap_info) {
721
706
GapEditScript* esp = hsp->gap_info;
722
for ( ; esp; esp = esp->next) {
723
if (esp->op_type == eGapAlignDel) {
724
length += esp->num;
725
gaps += esp->num;
726
++gap_opens;
727
} else if (esp->op_type == eGapAlignIns) {
728
++gap_opens;
729
gaps += esp->num;
707
Int4 index;
708
for (index=0; index<esp->size; index++) {
709
if (esp->op_type[index] == eGapAlignDel) {
710
length += esp->num[index];
711
gaps += esp->num[index];
712
++gap_opens;
713
} else if (esp->op_type[index] == eGapAlignIns) {
714
++gap_opens;
715
gaps += esp->num[index];
730
716
}
731
717
}
732
718
} else if (s_length > length) {
1059
1045
(h2->query.offset - h2->subject.offset);
1060
1046
}
1061
1047
1048
1062
1049
/** An auxiliary structure used for merging HSPs */
1063
1050
typedef struct BlastHSPSegment {
1064
Int4 q_start, /**< Start of segment in query. */
1065
q_end; /**< End of segment in query. */
1066
Int4 s_start, /**< Start of segment in subject. */
1067
s_end; /**< End of segment in subject. */
1068
1051
struct BlastHSPSegment* next; /**< Next link in the chain. */
1052
Int4 q_start, /**< Query start of segment. */
1053
s_start; /**< Subject start of segment. */
1054
Int4 length; /**< length of segments. */
1055
Int4 diagonal; /**< diagonal (q_start - s_start). */
1056
GapEditScript* esp; /**< pointer to edit script, not owned! */
1057
Int4 esp_offset; /**< element of above esp this node refers to */
1069
1058
} BlastHSPSegment;
1070
1059
1071
/** Maximal diagonal distance between HSP starting offsets, within which HSPs
1072
* from search of different chunks of subject sequence are considered for
1073
* merging.
1074
*/
1075
#define OVERLAP_DIAG_CLOSE 10
1076
1077
/** Merge the two HSPs if they intersect.
1078
* @param hsp1 The first HSP; also contains the result of merge. [in] [out]
1079
* @param hsp2 The second HSP [in]
1080
* @param start The starting offset of the subject coordinates where the
1081
* intersection is possible [in]
1082
*/
1083
static Boolean
1084
s_BlastHSPsMerge(BlastHSP* hsp1, BlastHSP* hsp2, Int4 start)
1085
{
1086
BlastHSPSegment* segments1,* segments2,* new_segment1,* new_segment2;
1087
GapEditScript* esp1,* esp2,* esp;
1088
Int4 end = start + DBSEQ_CHUNK_OVERLAP - 1;
1089
Int4 min_diag, max_diag, num1, num2, dist = 0, next_dist = 0;
1090
Int4 diag1_start, diag1_end, diag2_start, diag2_end;
1091
Int4 index;
1092
Uint1 intersection_found;
1093
EGapAlignOpType op_type = eGapAlignSub;
1094
1095
if (!hsp1->gap_info || !hsp2->gap_info) {
1060
1061
/** function to create and populate a BlastHSPSegement.
1062
* The new segemnt is added to an existing chain unless NULL
1063
*
1064
* @param segment object to be allocated and populated [in|out]
1065
* @param q_start start of match on query [in]
1066
* @param s_start start of match on subject [in]
1067
* @param esp GapEditScript for this alignment [in]
1068
* @param esp_offset relevant element of GapEditScript [in]
1069
*/
1070
static void
1071
s_AddHSPSegment(BlastHSPSegment** segment, int q_start, int s_start, GapEditScript* esp, int esp_offset)
1072
{
1073
BlastHSPSegment* new = malloc(sizeof(BlastHSPSegment));
1074
new->next = NULL;
1075
new->q_start = q_start;
1076
new->s_start = s_start;
1077
new->length = esp->num[esp_offset];
1078
new->diagonal = q_start - s_start;
1079
new->esp = esp;
1080
new->esp_offset = esp_offset;
1081
if (*segment)
1082
{
1083
BlastHSPSegment* var = *segment;
1084
while (var->next)
1085
var = var->next;
1086
var->next = new;
1087
}
1088
else
1089
*segment = new;
1090
1091
return;
1092
}
1093
1094
/** Deallocates all segments in chain.
1095
* Does not free underlying data
1096
* @param segments linked list of segments
1097
*/
1098
static BlastHSPSegment*
1099
s_DeleteAllHSPSegments(BlastHSPSegment* segments)
1100
{
1101
BlastHSPSegment* var = segments;
1102
BlastHSPSegment* next;
1103
while (var)
1104
{
1105
next = var->next;
1106
sfree(var);
1107
var = next;
1108
}
1109
return NULL;
1110
}
1111
1112
1113
1114
Boolean
1115
BlastMergeTwoHSPs(BlastHSP* hsp1, BlastHSP* hsp2, Int4 start)
1116
{
1117
BlastHSPSegment *segment1=NULL, *segment2=NULL;
1118
Boolean found = FALSE;
1119
GapEditScript* esp;
1120
Int4 q_start, s_start; /* start on query and subject sequences. */
1121
Int4 max_s_end = -1, max_q_end = -1; /* furthest extent of first HSP. */
1122
Int4 index; /* loop index. */
1123
1124
if (!hsp1->gap_info || !hsp2->gap_info)
1125
{
1096
1126
/* Assume that this is an ungapped alignment, hence simply compare
1097
1127
diagonals. Do not merge if they are on different diagonals */
1098
1128
if (s_DiagCompareHSPs(&hsp1, &hsp2) == 0 &&
1103
1133
} else
1104
1134
return FALSE;
1105
1135
}
1106
/* Find whether these HSPs have an intersection point */
1107
segments1 = (BlastHSPSegment*) calloc(1, sizeof(BlastHSPSegment));
1108
1109
esp1 = hsp1->gap_info;
1110
esp2 = hsp2->gap_info;
1111
1112
segments1->q_start = hsp1->query.offset;
1113
segments1->s_start = hsp1->subject.offset;
1114
while (segments1->s_start < start) {
1115
if (esp1->op_type == eGapAlignIns)
1116
segments1->q_start += esp1->num;
1117
else if (segments1->s_start + esp1->num < start) {
1118
if (esp1->op_type == eGapAlignSub) {
1119
segments1->s_start += esp1->num;
1120
segments1->q_start += esp1->num;
1121
} else if (esp1->op_type == eGapAlignDel)
1122
segments1->s_start += esp1->num;
1123
} else
1124
break;
1125
esp1 = esp1->next;
1126
}
1127
/* Current esp is the first segment within the overlap region */
1128
segments1->s_end = segments1->s_start + esp1->num - 1;
1129
if (esp1->op_type == eGapAlignSub)
1130
segments1->q_end = segments1->q_start + esp1->num - 1;
1131
else
1132
segments1->q_end = segments1->q_start;
1133
1134
new_segment1 = segments1;
1135
1136
for (esp = esp1->next; esp; esp = esp->next) {
1137
new_segment1->next = (BlastHSPSegment*)
1138
calloc(1, sizeof(BlastHSPSegment));
1139
new_segment1->next->q_start = new_segment1->q_end + 1;
1140
new_segment1->next->s_start = new_segment1->s_end + 1;
1141
new_segment1 = new_segment1->next;
1142
if (esp->op_type == eGapAlignSub) {
1143
new_segment1->q_end += esp->num - 1;
1144
new_segment1->s_end += esp->num - 1;
1145
} else if (esp->op_type == eGapAlignIns) {
1146
new_segment1->q_end += esp->num - 1;
1147
new_segment1->s_end = new_segment1->s_start;
1148
} else {
1149
new_segment1->s_end += esp->num - 1;
1150
new_segment1->q_end = new_segment1->q_start;
1151
}
1152
}
1153
1154
/* Now create the second segments list */
1155
1156
segments2 = (BlastHSPSegment*) calloc(1, sizeof(BlastHSPSegment));
1157
segments2->q_start = hsp2->query.offset;
1158
segments2->s_start = hsp2->subject.offset;
1159
segments2->q_end = segments2->q_start + esp2->num - 1;
1160
segments2->s_end = segments2->s_start + esp2->num - 1;
1161
1162
new_segment2 = segments2;
1163
1164
for (esp = esp2->next; esp && new_segment2->s_end < end;
1165
esp = esp->next) {
1166
new_segment2->next = (BlastHSPSegment*)
1167
calloc(1, sizeof(BlastHSPSegment));
1168
new_segment2->next->q_start = new_segment2->q_end + 1;
1169
new_segment2->next->s_start = new_segment2->s_end + 1;
1170
new_segment2 = new_segment2->next;
1171
if (esp->op_type == eGapAlignIns) {
1172
new_segment2->s_end = new_segment2->s_start;
1173
new_segment2->q_end = new_segment2->q_start + esp->num - 1;
1174
} else if (esp->op_type == eGapAlignDel) {
1175
new_segment2->s_end = new_segment2->s_start + esp->num - 1;
1176
new_segment2->q_end = new_segment2->q_start;
1177
} else if (esp->op_type == eGapAlignSub) {
1178
new_segment2->s_end = new_segment2->s_start + esp->num - 1;
1179
new_segment2->q_end = new_segment2->q_start + esp->num - 1;
1180
}
1181
}
1182
1183
new_segment1 = segments1;
1184
new_segment2 = segments2;
1185
intersection_found = 0;
1186
num1 = num2 = 0;
1187
while (new_segment1 && new_segment2 && !intersection_found) {
1188
if (new_segment1->s_end < new_segment2->s_start ||
1189
new_segment1->q_end < new_segment2->q_start) {
1190
new_segment1 = new_segment1->next;
1191
num1++;
1192
continue;
1193
}
1194
if (new_segment2->s_end < new_segment1->s_start ||
1195
new_segment2->q_end < new_segment1->q_start) {
1196
new_segment2 = new_segment2->next;
1197
num2++;
1198
continue;
1199
}
1200
diag1_start = new_segment1->s_start - new_segment1->q_start;
1201
diag2_start = new_segment2->s_start - new_segment2->q_start;
1202
diag1_end = new_segment1->s_end - new_segment1->q_end;
1203
diag2_end = new_segment2->s_end - new_segment2->q_end;
1204
1205
if (diag1_start == diag1_end && diag2_start == diag2_end &&
1206
diag1_start == diag2_start) {
1207
/* Both segments substitutions, on same diagonal */
1208
intersection_found = 1;
1209
dist = new_segment2->s_end - new_segment1->s_start + 1;
1210
break;
1211
} else if (diag1_start != diag1_end && diag2_start != diag2_end) {
1212
/* Both segments gaps - must intersect */
1213
intersection_found = 3;
1214
1215
op_type = eGapAlignIns;
1216
dist = new_segment2->s_end - new_segment1->s_start + 1;
1217
next_dist = new_segment2->q_end - new_segment1->q_start - dist + 1;
1218
if (new_segment2->q_end - new_segment1->q_start < dist) {
1219
dist = new_segment2->q_end - new_segment1->q_start + 1;
1220
op_type = eGapAlignDel;
1221
next_dist = new_segment2->s_end - new_segment1->s_start - dist + 1;
1222
}
1223
break;
1224
} else if (diag1_start != diag1_end) {
1225
max_diag = MAX(diag1_start, diag1_end);
1226
min_diag = MIN(diag1_start, diag1_end);
1227
if (diag2_start >= min_diag && diag2_start <= max_diag) {
1228
intersection_found = 2;
1229
dist = diag2_start - min_diag + 1;
1230
if (new_segment1->s_end == new_segment1->s_start)
1231
next_dist = new_segment2->s_end - new_segment1->s_end + 1;
1232
else
1233
next_dist = new_segment2->q_end - new_segment1->q_end + 1;
1234
break;
1235
}
1236
} else if (diag2_start != diag2_end) {
1237
max_diag = MAX(diag2_start, diag2_end);
1238
min_diag = MIN(diag2_start, diag2_end);
1239
if (diag1_start >= min_diag && diag1_start <= max_diag) {
1240
intersection_found = 2;
1241
next_dist = max_diag - diag1_start + 1;
1242
if (new_segment2->s_end == new_segment2->s_start)
1243
dist = new_segment2->s_start - new_segment1->s_start + 1;
1244
else
1245
dist = new_segment2->q_start - new_segment1->q_start + 1;
1246
break;
1247
}
1248
}
1249
if (new_segment1->s_end <= new_segment2->s_end) {
1250
new_segment1 = new_segment1->next;
1251
num1++;
1252
} else {
1253
new_segment2 = new_segment2->next;
1254
num2++;
1255
}
1256
}
1257
1258
/* Free the segments linked lists that are no longer needed. */
1259
for ( ; segments1; segments1 = new_segment1) {
1260
new_segment1 = segments1->next;
1261
sfree(segments1);
1262
}
1263
1264
for ( ; segments2; segments2 = new_segment2) {
1265
new_segment2 = segments2->next;
1266
sfree(segments2);
1267
}
1268
1269
if (intersection_found) {
1270
esp = NULL;
1271
for (index = 0; index < num1-1; index++)
1272
esp1 = esp1->next;
1273
for (index = 0; index < num2-1; index++) {
1274
esp = esp2;
1275
esp2 = esp2->next;
1276
}
1277
if (intersection_found < 3) {
1278
if (num1 > 0)
1279
esp1 = esp1->next;
1280
if (num2 > 0) {
1281
esp = esp2;
1282
esp2 = esp2->next;
1283
}
1284
}
1285
switch (intersection_found) {
1286
case 1:
1287
esp1->num = dist;
1288
esp1->next = esp2->next;
1289
esp2->next = NULL;
1290
break;
1291
case 2:
1292
esp1->num = dist;
1293
esp2->num = next_dist;
1294
esp1->next = esp2;
1295
if (esp)
1296
esp->next = NULL;
1297
break;
1298
case 3:
1299
esp1->num += dist;
1300
esp2->op_type = op_type;
1301
esp2->num = next_dist;
1302
esp1->next = esp2;
1303
if (esp)
1304
esp->next = NULL;
1305
break;
1306
default: break;
1307
}
1308
hsp1->query.end = hsp2->query.end;
1309
hsp1->subject.end = hsp2->subject.end;
1310
}
1311
1312
return (Boolean) intersection_found;
1136
1137
esp = hsp1->gap_info;
1138
q_start = hsp1->query.offset;
1139
s_start = hsp1->subject.offset;
1140
for (index=0; index<esp->size; index++)
1141
{
1142
if (esp->op_type[index] == eGapAlignSub)
1143
{
1144
if (s_start+(esp->num[index]) > start) /* End of segment within overlap region. */
1145
{
1146
s_AddHSPSegment(&segment1, q_start, s_start, esp, index);
1147
max_q_end = MAX(max_q_end, q_start+esp->num[index]);
1148
max_s_end = MAX(max_s_end, s_start+esp->num[index]);
1149
}
1150
q_start += esp->num[index];
1151
s_start += esp->num[index];
1152
}
1153
else if (esp->op_type[index] == eGapAlignIns)
1154
q_start += esp->num[index];
1155
else if (esp->op_type[index] == eGapAlignDel)
1156
s_start += esp->num[index];
1157
}
1158
1159
esp = hsp2->gap_info;
1160
q_start = hsp2->query.offset;
1161
s_start = hsp2->subject.offset;
1162
for (index=0; index<esp->size; index++)
1163
{
1164
if (esp->op_type[index] == eGapAlignSub)
1165
{
1166
if (max_s_end > s_start)
1167
s_AddHSPSegment(&segment2, q_start, s_start, esp, index);
1168
else
1169
break;
1170
1171
q_start += esp->num[index];
1172
s_start += esp->num[index];
1173
}
1174
else if (esp->op_type[index] == eGapAlignIns)
1175
q_start += esp->num[index];
1176
else if (esp->op_type[index] == eGapAlignDel)
1177
s_start += esp->num[index];
1178
}
1179
1180
if (segment1 && segment2)
1181
{
1182
BlastHSPSegment *segment1_var=segment1, *segment2_var=segment2;
1183
while (segment1_var)
1184
{
1185
BlastHSPSegment* segment2_var = segment2;
1186
while (segment2_var)
1187
{
1188
if (segment2_var->diagonal == segment1_var->diagonal)
1189
{
1190
if (segment1_var->q_start <= segment2_var->q_start
1191
&& (segment1_var->q_start+segment1_var->length) >= segment2_var->q_start)
1192
found = TRUE; /* start of 2 in extent of 1. */
1193
else if (segment1_var->q_start >= segment2_var->q_start
1194
&& (segment1_var->q_start) <= segment2_var->q_start+segment2_var->length)
1195
found = TRUE; /* start of 1 in extent of 2. */
1196
}
1197
if (found)
1198
break;
1199
segment2_var = segment2_var->next;
1200
}
1201
if (found)
1202
break;
1203
segment1_var = segment1_var->next;
1204
}
1205
1206
if (found == TRUE)
1207
{
1208
int end_1 = segment1_var->q_start + segment1_var->length;
1209
int end_2 = segment2_var->q_start + segment2_var->length;
1210
GapEditScript* esp_temp_1 = segment1_var->esp;
1211
GapEditScript* esp_temp_2 = segment2_var->esp;
1212
GapEditScript* new_esp = NULL;
1213
int new_size = 0;
1214
1215
/* we use the GapEditScript element pointed to be segment1_var, ignoring the one
1216
pointed to by segment2_var as it overlaps with the first. */
1217
esp_temp_1->num[segment1_var->esp_offset] += (end_2 - end_1);
1218
new_size = segment1_var->esp_offset + esp_temp_2->size - segment2_var->esp_offset;
1219
new_esp = GapEditScriptNew(new_size);
1220
GapEditScriptPartialCopy(new_esp, 0, esp_temp_1, 0, segment1_var->esp_offset);
1221
GapEditScriptPartialCopy(new_esp, 1+segment1_var->esp_offset, esp_temp_2,
1222
segment2_var->esp_offset+1, esp_temp_2->size-1);
1223
hsp1->gap_info = GapEditScriptDelete(hsp1->gap_info);
1224
hsp1->gap_info = new_esp;
1225
hsp1->query.end = hsp2->query.end;
1226
hsp1->subject.end = hsp2->subject.end;
1227
hsp1->score = MAX(hsp1->score, hsp2->score); /* Neither score may be correct, so we use the highest one. */
1228
}
1229
}
1230
1231
segment1 = s_DeleteAllHSPSegments(segment1);
1232
segment2 = s_DeleteAllHSPSegments(segment2);
1233
1234
return found;
1313
1235
}
1314
1236
1237
/** Maximal diagonal distance between HSP starting offsets, within which HSPs
1238
* from search of different chunks of subject sequence are considered for
1239
* merging.
1240
*/
1241
#define OVERLAP_DIAG_CLOSE 10
1315
1242
/********************************************************************************
1316
1243
Functions manipulating BlastHSPList's
1317
1244
********************************************************************************/
1351
1278
return hsp_list;
1352
1279
}
1353
1280
1354
/** Duplicate HSPList's contents, copying only pointers to individual HSPs,
1355
* effectively transferring the ownership of the pointers to the return value.
1356
* The caller is still responsible for calling Blast_HSPListFree on the first
1357
* argument to this function.
1358
* @param hsp_list source Blast_HSPList to copy HSP pointers from [in|out]
1359
* @return duplicate of hsp_list or NULL if out of memory
1360
*/
1361
static BlastHSPList*
1362
s_BlastHSPListReleaseHSPs(BlastHSPList* hsp_list)
1363
{
1364
BlastHSPList* new_hsp_list = (BlastHSPList*)
1365
BlastMemDup(hsp_list, sizeof(BlastHSPList));
1366
if ( !new_hsp_list ) {
1367
return NULL;
1368
}
1369
1370
new_hsp_list->hsp_array = (BlastHSP**)
1371
BlastMemDup(hsp_list->hsp_array, hsp_list->allocated*sizeof(BlastHSP*));
1372
if (!new_hsp_list->hsp_array) {
1373
sfree(new_hsp_list);
1374
return NULL;
1375
}
1376
new_hsp_list->allocated = hsp_list->allocated;
1377
ASSERT(hsp_list->hspcnt <= hsp_list->allocated);
1378
1379
/* hsp_list->hsp_array elements past hsp_list->hspcnt should be NULL, so it
1380
* doesn't hurt to zero hsp_list->allocated */
1381
memset(hsp_list->hsp_array, 0, hsp_list->allocated*sizeof(BlastHSP*));
1382
hsp_list->hspcnt = 0;
1383
1384
return new_hsp_list;
1385
}
1386
1387
1281
/** This is a copy of a static function from ncbimisc.c.
1388
1282
* Turns array into a heap with respect to a given comparison function.
1389
1283
*/
1758
1652
return 0;
1759
1653
}
1760
1654
1761
/** Callback for sorting HSPs by starting offset in query. The sorting criteria
1762
* in order of priority: context, starting offset in query, starting offset in
1763
* subject. Null HSPs are moved to the end of the array.
1655
/** Callback for sorting HSPs by starting offset in query. Sorting is by
1656
* increasing context, then increasing query start offset, then increasing
1657
* subject start offset, then decreasing score, then increasing query end
1658
* offset, then increasing subject end offset. Null HSPs are moved to the
1659
* end of the array.
1764
1660
* @param v1 pointer to first HSP [in]
1765
1661
* @param v2 pointer to second HSP [in]
1766
1662
* @return Result of comparison.
1768
1664
static int
1769
1665
s_QueryOffsetCompareHSPs(const void* v1, const void* v2)
1770
1666
{
1771
BlastHSP* h1,* h2;
1772
BlastHSP** hp1,** hp2;
1667
BlastHSP* h1,* h2;
1668
BlastHSP** hp1,** hp2;
1773
1669
1774
hp1 = (BlastHSP**) v1;
1775
hp2 = (BlastHSP**) v2;
1776
h1 = *hp1;
1777
h2 = *hp2;
1670
hp1 = (BlastHSP**) v1;
1671
hp2 = (BlastHSP**) v2;
1672
h1 = *hp1;
1673
h2 = *hp2;
1778
1674
1779
1675
if (!h1 && !h2)
1780
1676
return 0;
1789
1685
if (h1->context > h2->context)
1790
1686
return 1;
1791
1687
1792
if (h1->query.offset < h2->query.offset)
1793
return -1;
1794
if (h1->query.offset > h2->query.offset)
1795
return 1;
1796
1797
if (h1->subject.offset < h2->subject.offset)
1798
return -1;
1799
if (h1->subject.offset > h2->subject.offset)
1800
return 1;
1801
1802
return 0;
1688
if (h1->query.offset < h2->query.offset)
1689
return -1;
1690
if (h1->query.offset > h2->query.offset)
1691
return 1;
1692
1693
if (h1->subject.offset < h2->subject.offset)
1694
return -1;
1695
if (h1->subject.offset > h2->subject.offset)
1696
return 1;
1697
1698
/* tie breakers: sort by decreasing score, then
1699
by increasing size of query range, then by
1700
increasing subject range. */
1701
1702
if (h1->score < h2->score)
1703
return 1;
1704
if (h1->score > h2->score)
1705
return -1;
1706
1707
if (h1->query.end < h2->query.end)
1708
return -1;
1709
if (h1->query.end > h2->query.end)
1710
return 1;
1711
1712
if (h1->subject.end < h2->subject.end)
1713
return -1;
1714
if (h1->subject.end > h2->subject.end)
1715
return 1;
1716
1717
return 0;
1803
1718
}
1804
1719
1805
/** Callback for sorting HSPs by ending offset in query. The sorting criteria
1806
* in order of priority: context, ending offset in query, ending offset in
1807
* subject. Null HSPs are moved to the end of the array.
1720
/** Callback for sorting HSPs by ending offset in query. Sorting is by
1721
* increasing context, then increasing query end offset, then increasing
1722
* subject end offset, then decreasing score, then decreasing query start
1723
* offset, then decreasing subject start offset. Null HSPs are moved to the
1724
* end of the array.
1808
1725
* @param v1 pointer to first HSP [in]
1809
1726
* @param v2 pointer to second HSP [in]
1810
1727
* @return Result of comparison.
1812
1729
static int
1813
1730
s_QueryEndCompareHSPs(const void* v1, const void* v2)
1814
1731
{
1815
BlastHSP* h1,* h2;
1816
BlastHSP** hp1,** hp2;
1732
BlastHSP* h1,* h2;
1733
BlastHSP** hp1,** hp2;
1817
1734
1818
hp1 = (BlastHSP**) v1;
1819
hp2 = (BlastHSP**) v2;
1820
h1 = *hp1;
1821
h2 = *hp2;
1735
hp1 = (BlastHSP**) v1;
1736
hp2 = (BlastHSP**) v2;
1737
h1 = *hp1;
1738
h2 = *hp2;
1822
1739
1823
1740
if (!h1 && !h2)
1824
1741
return 0;
1833
1750
if (h1->context > h2->context)
1834
1751
return 1;
1835
1752
1836
if (h1->query.end < h2->query.end)
1837
return -1;
1838
if (h1->query.end > h2->query.end)
1839
return 1;
1840
1841
if (h1->subject.end < h2->subject.end)
1842
return -1;
1843
if (h1->subject.end > h2->subject.end)
1844
return 1;
1845
1846
return 0;
1753
if (h1->query.end < h2->query.end)
1754
return -1;
1755
if (h1->query.end > h2->query.end)
1756
return 1;
1757
1758
if (h1->subject.end < h2->subject.end)
1759
return -1;
1760
if (h1->subject.end > h2->subject.end)
1761
return 1;
1762
1763
/* tie breakers: sort by decreasing score, then
1764
by increasing size of query range, then by
1765
increasing size of subject range. The shortest range
1766
means the *largest* sequence offset must come
1767
first */
1768
if (h1->score < h2->score)
1769
return 1;
1770
if (h1->score > h2->score)
1771
return -1;
1772
1773
if (h1->query.offset < h2->query.offset)
1774
return 1;
1775
if (h1->query.offset > h2->query.offset)
1776
return -1;
1777
1778
if (h1->subject.offset < h2->subject.offset)
1779
return 1;
1780
if (h1->subject.offset > h2->subject.offset)
1781
return -1;
1782
1783
return 0;
1847
1784
}
1848
1785
1849
1786
Int4
1852
1789
1853
1790
{
1854
1791
BlastHSP** hsp_array; /* hsp_array to purge. */
1855
Int4 index = 0; /* loop index. */
1856
Int4 increment = 1;
1857
Int2 retval = 0;
1792
Int4 i, j;
1793
Int2 retval;
1858
1794
Int4 hsp_count;
1859
1795
1860
1796
/* If HSP list is empty, return immediately. */
1870
1806
hsp_count = hsp_list->hspcnt;
1871
1807
1872
1808
qsort(hsp_array, hsp_count, sizeof(BlastHSP*), s_QueryOffsetCompareHSPs);
1873
while (index < hsp_count-increment)
1874
{
1875
if (hsp_array[index+increment] == NULL)
1876
{
1877
increment++;
1878
continue;
1879
}
1880
1881
if (hsp_array[index] && hsp_array[index]->query.offset == hsp_array[index+increment]->query.offset &&
1882
hsp_array[index]->subject.offset == hsp_array[index+increment]->subject.offset &&
1883
hsp_array[index]->context == hsp_array[index+increment]->context)
1884
{
1885
if (hsp_array[index]->score > hsp_array[index+increment]->score) {
1886
hsp_array[index+increment] =
1887
Blast_HSPFree(hsp_array[index+increment]);
1888
increment++;
1889
} else {
1890
hsp_array[index] = Blast_HSPFree(hsp_array[index]);
1891
index++;
1892
increment = 1;
1893
}
1894
} else {
1895
index++;
1896
increment = 1;
1897
}
1809
i = 0;
1810
while (i < hsp_count) {
1811
j = 1;
1812
while (i+j < hsp_count &&
1813
hsp_array[i] && hsp_array[i+j] &&
1814
hsp_array[i]->context == hsp_array[i+j]->context &&
1815
hsp_array[i]->query.offset == hsp_array[i+j]->query.offset &&
1816
hsp_array[i]->subject.offset == hsp_array[i+j]->subject.offset) {
1817
hsp_array[i+j] = Blast_HSPFree(hsp_array[i+j]);
1818
j++;
1819
}
1820
i += j;
1898
1821
}
1899
1822
1900
1823
qsort(hsp_array, hsp_count, sizeof(BlastHSP*), s_QueryEndCompareHSPs);
1901
index=0;
1902
increment=1;
1903
while (index < hsp_count-increment)
1904
{
1905
if (hsp_array[index+increment] == NULL)
1906
{
1907
increment++;
1908
continue;
1909
}
1910
1911
if (hsp_array[index] &&
1912
hsp_array[index]->query.end == hsp_array[index+increment]->query.end &&
1913
hsp_array[index]->subject.end == hsp_array[index+increment]->subject.end &&
1914
hsp_array[index]->context == hsp_array[index+increment]->context)
1915
{
1916
if (hsp_array[index]->score > hsp_array[index+increment]->score) {
1917
hsp_array[index+increment] =
1918
Blast_HSPFree(hsp_array[index+increment]);
1919
increment++;
1920
} else {
1921
hsp_array[index] = Blast_HSPFree(hsp_array[index]);
1922
index++;
1923
increment = 1;
1924
}
1925
} else {
1926
index++;
1927
increment = 1;
1928
}
1824
i = 0;
1825
while (i < hsp_count) {
1826
j = 1;
1827
while (i+j < hsp_count &&
1828
hsp_array[i] && hsp_array[i+j] &&
1829
hsp_array[i]->context == hsp_array[i+j]->context &&
1830
hsp_array[i]->query.end == hsp_array[i+j]->query.end &&
1831
hsp_array[i]->subject.end == hsp_array[i+j]->subject.end) {
1832
hsp_array[i+j] = Blast_HSPFree(hsp_array[i+j]);
1833
j++;
1834
}
1835
i += j;
1929
1836
}
1930
1837
1931
1838
retval = Blast_HSPListPurgeNullHSPs(hsp_list);
1935
1842
return hsp_list->hspcnt;
1936
1843
}
1937
1844
1938
/** Status values returned by an HSP inclusion test function. */
1939
typedef enum EHSPInclusionStatus {
1940
eEqual = 0, /**< Identical */
1941
eFirstInSecond, /**< First included in rectangle formed by second */
1942
eSecondInFirst, /**< Second included in rectangle formed by first */
1943
eDiagNear, /**< Diagonals are near, but neither HSP is included in
1944
the other. */
1945
eDiagDistant /**< Diagonals are far apart, or different contexts */
1946
} EHSPInclusionStatus;
1947
1948
1845
/** Diagonal distance between HSPs, outside of which one HSP cannot be
1949
1846
* considered included in the other.
1950
1847
*/
1951
1848
#define MIN_DIAG_DIST 60
1952
1849
1953
/** HSP inclusion criterion for megablast: one HSP must be included in a
1954
* diagonal strip of a certain width around the other, and also in a rectangle
1955
* formed by the other HSP's endpoints.
1956
*/
1957
static EHSPInclusionStatus
1958
s_BlastHSPInclusionTest(BlastHSP* hsp1, BlastHSP* hsp2)
1959
{
1960
if (hsp1->context != hsp2->context ||
1961
!MB_HSP_CLOSE(hsp1->query.offset, hsp2->query.offset,
1962
hsp1->subject.offset, hsp2->subject.offset,
1963
MIN_DIAG_DIST))
1964
return eDiagDistant;
1965
1966
if (hsp1->query.offset == hsp2->query.offset &&
1967
hsp1->query.end == hsp2->query.end &&
1968
hsp1->subject.offset == hsp2->subject.offset &&
1969
hsp1->subject.end == hsp2->subject.end &&
1970
hsp1->score == hsp2->score) {
1971
return eEqual;
1972
} else if (hsp1->query.offset >= hsp2->query.offset &&
1973
hsp1->query.end <= hsp2->query.end &&
1974
hsp1->subject.offset >= hsp2->subject.offset &&
1975
hsp1->subject.end <= hsp2->subject.end &&
1976
hsp1->score <= hsp2->score) {
1977
return eFirstInSecond;
1978
} else if (hsp1->query.offset <= hsp2->query.offset &&
1979
hsp1->query.end >= hsp2->query.end &&
1980
hsp1->subject.offset <= hsp2->subject.offset &&
1981
hsp1->subject.end >= hsp2->subject.end &&
1982
hsp1->score >= hsp2->score) {
1983
return eSecondInFirst;
1984
}
1985
return eDiagNear;
1986
}
1987
1988
/** How many HSPs to check for inclusion for each new HSP? */
1989
#define MAX_NUM_CHECK_INCLUSION 20
1990
1991
1850
/** Remove redundant HSPs in an HSP list based on a diagonal inclusion test: if
1992
1851
* an HSP is within a certain diagonal distance of another HSP, and its endpoints
1993
1852
* are contained in a rectangle formed by another HSP, then it is removed.
1994
1853
* Performed only after a single-phase greedy gapped extension, when there is no
1995
1854
* extra traceback stage that could fix the inclusions.
1996
* @param hsp_list HSP list to check [in] [out]
1855
* @param query_blk Query sequence for the HSPs
1856
* @param subject_blk Subject sequence for the HSPs
1857
* @param query_info Used to map HSPs uniquely onto the complete
1858
* set of query sequences
1859
* @param hsp_list HSP list to check (must be in standard sorted order) [in/out]
1997
1860
*/
1998
static Int2
1999
s_BlastHSPListCheckDiagonalInclusion(BlastHSPList* hsp_list)
1861
static void
1862
s_BlastHSPListCheckDiagonalInclusion(BLAST_SequenceBlk* query_blk,
1863
BLAST_SequenceBlk* subject_blk,
1864
const BlastQueryInfo* query_info,
1865
BlastHSPList* hsp_list)
2000
1866
{
2001
Int4 index, new_hspcnt, index1, index2;
1867
Int4 index;
2002
1868
BlastHSP** hsp_array = hsp_list->hsp_array;
2003
Boolean shift_needed = FALSE;
2004
EHSPInclusionStatus inclusion_status = eDiagNear;
1869
BlastIntervalTree *tree;
2005
1870
2006
1871
if (hsp_list->hspcnt <= 1)
2007
return 0;
2008
2009
qsort(hsp_array, hsp_list->hspcnt, sizeof(BlastHSP*),
2010
s_DiagCompareHSPs);
2011
2012
for (index=1, new_hspcnt=0; index<hsp_list->hspcnt; index++) {
2013
if (!hsp_array[index])
2014
break;
2015
inclusion_status = eDiagNear;
2016
for (index1 = new_hspcnt; inclusion_status != eDiagDistant &&
2017
index1 >= 0 && new_hspcnt-index1 < MAX_NUM_CHECK_INCLUSION;
2018
index1--) {
2019
inclusion_status =
2020
s_BlastHSPInclusionTest(hsp_array[index], hsp_array[index1]);
2021
if (inclusion_status == eFirstInSecond ||
2022
inclusion_status == eEqual) {
2023
/* Free the new HSP and break out of the inclusion test loop */
2024
hsp_array[index] = Blast_HSPFree(hsp_array[index]);
2025
break;
2026
} else if (inclusion_status == eSecondInFirst) {
2027
hsp_array[index1] = Blast_HSPFree(hsp_array[index1]);
2028
shift_needed = TRUE;
2029
}
2030
}
2031
2032
/* If some lower indexed HSPs have been removed, shift the subsequent
2033
HSPs */
2034
if (shift_needed) {
2035
/* Find the first non-NULL HSP, going backwards */
2036
while (index1 >= 0 && !hsp_array[index1])
2037
index1--;
2038
/* Go forward, and shift any non-NULL HSPs */
2039
for (index2 = ++index1; index1 <= new_hspcnt; index1++) {
2040
if (hsp_array[index1])
2041
hsp_array[index2++] = hsp_array[index1];
2042
}
2043
new_hspcnt = index2 - 1;
2044
shift_needed = FALSE;
2045
}
2046
if (hsp_array[index] != NULL)
2047
hsp_array[++new_hspcnt] = hsp_array[index];
1872
return;
1873
1874
/* Remove any HSPs that are contained within other HSPs.
1875
Since the list is sorted by score already, any HSP
1876
contained by a previous HSP is guaranteed to have a
1877
lower score, and may be purged. */
1878
1879
tree = Blast_IntervalTreeInit(0, query_blk->length + 1,
1880
0, subject_blk->length + 1);
1881
1882
for (index = 0; index < hsp_list->hspcnt; index++) {
1883
BlastHSP *hsp = hsp_array[index];
1884
1885
if (BlastIntervalTreeContainsHSP(tree, hsp, query_info,
1886
MIN_DIAG_DIST)) {
1887
hsp_array[index] = Blast_HSPFree(hsp);
1888
}
1889
else {
1890
BlastIntervalTreeAddHSP(hsp, tree, query_info,
1891
eQueryAndSubject);
1892
}
2048
1893
}
2049
/* Set all HSP pointers that will not be used any more to NULL */
2050
memset(&hsp_array[new_hspcnt+1], 0,
2051
(hsp_list->hspcnt - new_hspcnt - 1)*sizeof(BlastHSP*));
2052
hsp_list->hspcnt = new_hspcnt + 1;
2053
2054
/* Sort the HSP array by score again. */
2055
Blast_HSPListSortByScore(hsp_list);
2056
2057
return 0;
1894
tree = Blast_IntervalTreeFree(tree);
1895
Blast_HSPListPurgeNullHSPs(hsp_list);
2058
1896
}
2059
1897
2060
1898
Int2
2103
1941
/* Return the packed sequence to the database */
2104
1942
BlastSeqSrcReleaseSequence(seq_src, (void*) &seq_arg);
2105
1943
/* Get the unpacked sequence */
2106
BlastSeqSrcGetSequence(seq_src, (void*) &seq_arg);
1944
if ((status=BlastSeqSrcGetSequence(seq_src, (void*) &seq_arg)))
1945
return status;
2107
1946
}
2108
1947
2109
1948
if (kTranslateSubject) {
2182
2021
Blast_HSPListPurgeNullHSPs(hsp_list);
2183
2022
}
2184
2023
2024
/* If greedy extension with traceback was used, remove
2025
any HSPs that share a starting or ending diagonal
2026
with a higher-scoring HSP. */
2027
if (gapped)
2028
Blast_HSPListPurgeHSPsWithCommonEndpoints(program, hsp_list);
2029
2030
/* Sort the HSP array by score (scores may have changed!) */
2031
Blast_HSPListSortByScore(hsp_list);
2032
2185
2033
/* If greedy extension with traceback was used, check for
2186
2034
HSP inclusion in a diagonal strip around another HSP. */
2187
2035
if (gapped)
2188
status = s_BlastHSPListCheckDiagonalInclusion(hsp_list);
2189
2190
/* Sort the HSP array by score (scores may have changed!) */
2191
Blast_HSPListSortByScore(hsp_list);
2036
s_BlastHSPListCheckDiagonalInclusion(query_blk, subject_blk,
2037
query_info, hsp_list);
2192
2038
2193
2039
Blast_HSPListAdjustOddBlastnScores(hsp_list, gapped, sbp);
2194
2040
2195
return status;
2041
return 0;
2196
2042
}
2197
2043
2198
2044
/** Combine two HSP lists, without altering the individual HSPs, and without
2265
2111
hsp_list->hspcnt = 0;
2266
2112
}
2267
2113
2268
Int2 Blast_HSPListAppend(BlastHSPList* hsp_list,
2114
Int2 Blast_HSPListAppend(BlastHSPList** old_hsp_list_ptr,
2269
2115
BlastHSPList** combined_hsp_list_ptr, Int4 hsp_num_max)
2270
2116
{
2271
2117
BlastHSPList* combined_hsp_list = *combined_hsp_list_ptr;
2272
BlastHSP** new_hsp_array;
2118
BlastHSPList* hsp_list = *old_hsp_list_ptr;
2273
2119
Int4 new_hspcnt;
2274
2120
2275
2121
if (!hsp_list || hsp_list->hspcnt == 0)
2276
2122
return 0;
2277
2123
2278
/* If no previous HSP list, just return a copy of the new one. */
2124
/* If no previous HSP list, return a pointer to the old one */
2279
2125
if (!combined_hsp_list) {
2280
if (!(combined_hsp_list = s_BlastHSPListReleaseHSPs(hsp_list)))
2281
return 1;
2282
*combined_hsp_list_ptr = combined_hsp_list;
2283
2126
*combined_hsp_list_ptr = hsp_list;
2127
*old_hsp_list_ptr = NULL;
2284
2128
return 0;
2285
2129
}
2130
2286
2131
/* Just append new list to the end of the old list, in case of
2287
2132
multiple frames of the subject sequence */
2288
2133
new_hspcnt = MIN(combined_hsp_list->hspcnt + hsp_list->hspcnt,
2290
2135
if (new_hspcnt > combined_hsp_list->allocated &&
2291
2136
!combined_hsp_list->do_not_reallocate) {
2292
2137
Int4 new_allocated = MIN(2*new_hspcnt, hsp_num_max);
2138
BlastHSP** new_hsp_array;
2293
2139
new_hsp_array = (BlastHSP**)
2294
2140
realloc(combined_hsp_list->hsp_array,
2295
2141
new_allocated*sizeof(BlastHSP*));
2307
2153
2308
2154
s_BlastHSPListsCombineByScore(hsp_list, combined_hsp_list, new_hspcnt);
2309
2155
2156
hsp_list = Blast_HSPListFree(hsp_list);
2157
*old_hsp_list_ptr = NULL;
2158
2310
2159
return 0;
2311
2160
}
2312
2161
2313
Int2 Blast_HSPListsMerge(BlastHSPList* hsp_list,
2162
Int2 Blast_HSPListsMerge(BlastHSPList** hsp_list_ptr,
2314
2163
BlastHSPList** combined_hsp_list_ptr,
2315
2164
Int4 hsp_num_max, Int4 start, Boolean merge_hsps)
2316
2165
{
2317
2166
BlastHSPList* combined_hsp_list = *combined_hsp_list_ptr;
2167
BlastHSPList* hsp_list = *hsp_list_ptr;
2318
2168
BlastHSP* hsp, *hsp_var;
2319
2169
BlastHSP** hspp1,** hspp2;
2320
Int4 index, index1;
2170
Int4 index, index1, last_index1;
2321
2171
Int4 hspcnt1, hspcnt2, new_hspcnt = 0;
2322
2172
BlastHSP** new_hsp_array;
2323
2173
2326
2176
2327
2177
/* If no previous HSP list, just return a copy of the new one. */
2328
2178
if (!combined_hsp_list) {
2329
*combined_hsp_list_ptr = s_BlastHSPListReleaseHSPs(hsp_list);
2179
*combined_hsp_list_ptr = hsp_list;
2180
*hsp_list_ptr = NULL;
2330
2181
return 0;
2331
2182
}
2332
2183
2362
2213
qsort(hspp1, hspcnt1, sizeof(BlastHSP*), s_DiagCompareHSPs);
2363
2214
qsort(hspp2, hspcnt2, sizeof(BlastHSP*), s_DiagCompareHSPs);
2364
2215
2365
for (index=0; index<hspcnt1; index++) {
2366
for (index1 = 0; index1<hspcnt2; index1++) {
2216
for (index=last_index1=0; index<hspcnt1; index++) {
2217
2218
/* scan through hspp2 until an HSP is found whose
2219
starting diagonal is less than OVERLAP_DIAG_CLOSE
2220
diagonals ahead of the current HSP from hspp1 */
2221
2222
for (index1 = last_index1; index1<hspcnt2; index1++,last_index1++) {
2223
/* Skip already deleted HSPs */
2224
if (!hspp2[index1])
2225
continue;
2226
if (s_DiagCompareHSPs(&hspp1[index], &hspp2[index1]) <
2227
OVERLAP_DIAG_CLOSE) {
2228
break;
2229
}
2230
}
2231
2232
/* attempt to merge the HSPs in hspp2 until their
2233
diagonals occur too far away */
2234
2235
for (; index1<hspcnt2; index1++) {
2367
2236
/* Skip already deleted HSPs */
2368
2237
if (!hspp2[index1])
2369
2238
continue;
2371
2240
ABS(s_DiagCompareHSPs(&hspp1[index], &hspp2[index1])) <
2372
2241
OVERLAP_DIAG_CLOSE) {
2373
2242
if (merge_hsps) {
2374
if (s_BlastHSPsMerge(hspp1[index], hspp2[index1], start)) {
2243
if (BlastMergeTwoHSPs(hspp1[index], hspp2[index1], start)) {
2375
2244
/* Free the second HSP. */
2376
2245
hspp2[index1] = Blast_HSPFree(hspp2[index1]);
2377
2246
}
2427
2296
2428
2297
s_BlastHSPListsCombineByScore(hsp_list, combined_hsp_list, new_hspcnt);
2429
2298
2299
hsp_list = Blast_HSPListFree(hsp_list);
2300
*hsp_list_ptr = NULL;
2301
2430
2302
return 0;
2431
2303
}
2432
2304
2433
2305
void Blast_HSPListAdjustOffsets(BlastHSPList* hsp_list, Int4 offset)
2434
2306
{
2435
2307
Int4 index;
2436
BlastHSP* hsp;
2308
2309
if (offset == 0) {
2310
return;
2311
}
2437
2312
2438
2313
for (index=0; index<hsp_list->hspcnt; index++) {
2439
hsp = hsp_list->hsp_array[index];
2314
BlastHSP* hsp = hsp_list->hsp_array[index];
2440
2315
hsp->subject.offset += offset;
2441
2316
hsp->subject.end += offset;
2442
2317
hsp->subject.gapped_start += offset;
3025
2900
3026
2901
if (!(tmp_hsp_list = hsp_list_array[query_index])) {
3027
2902
hsp_list_array[query_index] = tmp_hsp_list =
3028
Blast_HSPListNew(blasthit_params->options->hsp_num_max);
2903
Blast_HSPListNew(blasthit_params->hsp_num_max);
3029
2904
if (tmp_hsp_list == NULL)
3030
2905
{
3031
2906
sfree(hsp_list_array);
3175
3050
3176
3051
return phi_results;
3177
3052
}
3053
3054
BlastHSPResults*
3055
Blast_HSPResultsFromHSPStream(BlastHSPStream* hsp_stream,
3056
size_t num_queries,
3057
const BlastHitSavingOptions* hit_options,
3058
const BlastExtensionOptions* ext_options,
3059
const BlastScoringOptions* scoring_options)
3060
{
3061
BlastHSPResults* retval = NULL;
3062
SBlastHitsParameters* bhp = NULL;
3063
BlastHSPList* hsp_list = NULL;
3064
3065
SBlastHitsParametersNew(hit_options, ext_options, scoring_options, &bhp);
3066
3067
retval = Blast_HSPResultsNew((Int4) num_queries);
3068
3069
while (BlastHSPStreamRead(hsp_stream, &hsp_list) != kBlastHSPStream_Eof) {
3070
Blast_HSPResultsInsertHSPList(retval, hsp_list,
3071
bhp->prelim_hitlist_size);
3072
}
3073
bhp = SBlastHitsParametersFree(bhp);
3074
return retval;
3075
}
3076
3077
/** Comparison function for sorting HSP lists in increasing order of the
3078
* number of HSPs in a hit. Needed for s_TrimResultsByTotalHSPLimit below.
3079
* @param v1 Pointer to the first HSP list [in]
3080
* @param v2 Pointer to the second HSP list [in]
3081
*/
3082
static int
3083
s_CompareHsplistHspcnt(const void* v1, const void* v2)
3084
{
3085
BlastHSPList* r1 = *((BlastHSPList**) v1);
3086
BlastHSPList* r2 = *((BlastHSPList**) v2);
3087
3088
if (r1->hspcnt < r2->hspcnt)
3089
return -1;
3090
else if (r1->hspcnt > r2->hspcnt)
3091
return 1;
3092
else
3093
return 0;
3094
}
3095
3096
/** Removes extra results if a limit is imposed on the total number of HSPs
3097
* returned. If the search involves multiple query sequences, the total HSP
3098
* limit is applied separately to each query.
3099
* The trimming algorithm makes sure that at least 1 HSP is returned for each
3100
* database sequence hit. Suppose results for a given query consist of HSP
3101
* lists for N database sequences, and the limit is T. HSP lists are sorted in
3102
* order of increasing number of HSPs in each list. Then the algorithm proceeds
3103
* by leaving at most i*T/N HSPs for the first i HSP lists, for every
3104
* i = 1, 2, ..., N.
3105
* @param results Results after preliminary stage of a BLAST search [in|out]
3106
* @param total_hsp_limit Limit on total number of HSPs [in]
3107
* @return TRUE if HSP limit has been exceeded, FALSE otherwise.
3108
*/
3109
static Boolean
3110
s_TrimResultsByTotalHSPLimit(BlastHSPResults* results, Uint4 total_hsp_limit)
3111
{
3112
int query_index;
3113
Boolean hsp_limit_exceeded = FALSE;
3114
3115
if (total_hsp_limit == 0) {
3116
return hsp_limit_exceeded;
3117
}
3118
3119
for (query_index = 0; query_index < results->num_queries; ++query_index) {
3120
BlastHitList* hit_list = NULL;
3121
BlastHSPList** hsplist_array = NULL;
3122
Int4 hsplist_count = 0;
3123
int subj_index;
3124
3125
if ( !(hit_list = results->hitlist_array[query_index]) )
3126
continue;
3127
/* The count of HSPs is separate for each query. */
3128
hsplist_count = hit_list->hsplist_count;
3129
hsplist_array = (BlastHSPList**)
3130
malloc(hsplist_count*sizeof(BlastHSPList*));
3131
3132
for (subj_index = 0; subj_index < hsplist_count; ++subj_index) {
3133
hsplist_array[subj_index] = hit_list->hsplist_array[subj_index];
3134
}
3135
3136
qsort((void*)hsplist_array, hsplist_count,
3137
sizeof(BlastHSPList*), s_CompareHsplistHspcnt);
3138
3139
{
3140
Int4 tot_hsps = 0; /* total number of HSPs */
3141
Uint4 hsp_per_seq = MAX(1, total_hsp_limit/hsplist_count);
3142
for (subj_index = 0; subj_index < hsplist_count; ++subj_index) {
3143
Int4 allowed_hsp_num = ((subj_index+1)*hsp_per_seq) - tot_hsps;
3144
BlastHSPList* hsp_list = hsplist_array[subj_index];
3145
if (hsp_list->hspcnt > allowed_hsp_num) {
3146
/* Free the extra HSPs */
3147
int hsp_index;
3148
for (hsp_index = allowed_hsp_num;
3149
hsp_index < hsp_list->hspcnt; ++hsp_index) {
3150
Blast_HSPFree(hsp_list->hsp_array[hsp_index]);
3151
}
3152
hsp_list->hspcnt = allowed_hsp_num;
3153
hsp_limit_exceeded = TRUE;
3154
}
3155
tot_hsps += hsp_list->hspcnt;
3156
}
3157
}
3158
sfree(hsplist_array);
3159
}
3160
3161
return hsp_limit_exceeded;
3162
}
3163
3164
BlastHSPResults*
3165
Blast_HSPResultsFromHSPStreamWithLimit(BlastHSPStream* hsp_stream,
3166
Uint4 num_queries,
3167
const BlastHitSavingOptions* hit_options,
3168
const BlastExtensionOptions* ext_options,
3169
const BlastScoringOptions* scoring_options,
3170
Uint4 max_num_hsps,
3171
Boolean* removed_hsps)
3172
{
3173
Boolean rm_hsps = FALSE;
3174
BlastHSPResults* retval = Blast_HSPResultsFromHSPStream(hsp_stream,
3175
num_queries,
3176
hit_options,
3177
ext_options,
3178
scoring_options);
3179
3180
rm_hsps = s_TrimResultsByTotalHSPLimit(retval, max_num_hsps);
3181
if (removed_hsps) {
3182
*removed_hsps = rm_hsps;
3183
}
3184
return retval;
3185
}
3186
Loggerhead is a web-based interface for Breezy
Version: 2.0.1