← 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_gapalign.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_gapalign.c
1
 
/* $Id: blast_gapalign.c,v 1.163 2005/11/30 18:29:14 papadopo Exp $
 
1
/* $Id: blast_gapalign.c,v 1.171 2006/05/01 17:28:08 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_gapalign.c,v 1.163 2005/11/30 18:29:14 papadopo Exp $";
 
35
    "$Id: blast_gapalign.c,v 1.171 2006/05/01 17:28:08 camacho Exp $";
36
36
#endif /* SKIP_DOXYGEN_PROCESSING */
37
37
 
38
38
#include <algo/blast/core/blast_options.h>
190
190
   if (score_params->reward % 2 == 1) {
191
191
      reward = 2*score_params->reward;
192
192
      penalty = -2*score_params->penalty;
193
 
      Xdrop = 2*ext_params->gap_x_dropoff;
 
193
      Xdrop = 2*MAX(ext_params->gap_x_dropoff,
 
194
                    ext_params->gap_x_dropoff_final);
194
195
      gap_open = 2*score_params->gap_open;
195
196
      gap_extend = 2*score_params->gap_extend;
196
197
   } else {
197
198
      reward = score_params->reward;
198
199
      penalty = -score_params->penalty;
199
 
      Xdrop = ext_params->gap_x_dropoff;
 
200
      Xdrop = MAX(ext_params->gap_x_dropoff,
 
201
                  ext_params->gap_x_dropoff_final);
200
202
      gap_open = score_params->gap_open;
201
203
      gap_extend = score_params->gap_extend;
202
204
   }
2021
2023
    return best_score;
2022
2024
}
2023
2025
 
 
2026
/** Simple wrapper around binary search function for BlastQueryInfo structure
 
2027
 * to obtain the context in which this initial ungapped HSP is found
 
2028
 * @param query_info Query information structure [in]
 
2029
 * @param init_hsp Initial HSP [in]
 
2030
 */
 
2031
static NCBI_INLINE Int4
 
2032
s_GetUngappedHSPContext(const BlastQueryInfo* query_info, 
 
2033
                        const BlastInitHSP* init_hsp)
 
2034
{
 
2035
    return BSearchContextInfo(init_hsp->offsets.qs_offsets.q_off, query_info);
 
2036
}
 
2037
 
 
2038
/** Adjust the HSP offsets in the initial ungapped HSP structure given the
 
2039
 * query start 
 
2040
 * @param init_hsp Initial HSP [in|out]
 
2041
 * @param query_start Starting offset of the query sequence in the query info
 
2042
 * structure [in]
 
2043
 */
 
2044
static NCBI_INLINE void
 
2045
s_AdjustInitialHSPOffsets(BlastInitHSP* init_hsp, Int4 query_start)
 
2046
{
 
2047
    init_hsp->offsets.qs_offsets.q_off -= query_start;
 
2048
    if (init_hsp->ungapped_data) {
 
2049
        init_hsp->ungapped_data->q_start -= query_start;
 
2050
    }
 
2051
}
 
2052
 
 
2053
/** Set up a BLAST_SequenceBlk structure for a single query sequence
 
2054
 * @param concatenated_query Query sequence block for concatenated query [in]
 
2055
 * @param query_info Query information structure for concatenated query [in]
 
2056
 * @param context Context in which the HSP ocurred [in]
 
2057
 * @param single_query query Output sequence block which will contain adjusted
 
2058
 * sequence pointer and sequence length [out]
 
2059
 * @param query_start Starting offset of the single query sequence in the
 
2060
 * concatenated query [out]
 
2061
  */
 
2062
static void
 
2063
s_SetUpLocalBlastSequenceBlk(const BLAST_SequenceBlk* concatenated_query,
 
2064
                             const BlastQueryInfo* query_info, Int4 context,
 
2065
                             BLAST_SequenceBlk* single_query,
 
2066
                             Int4* query_start)
 
2067
{
 
2068
    Int4 query_length = 0;
 
2069
    if (concatenated_query->oof_sequence) {
 
2070
        /* Out-of-frame blastx case: all frames of the same parity are mixed
 
2071
           together in a special sequence. */
 
2072
        Int4 query_end_pt = 0;
 
2073
        Int4 mixed_frame_context = context - context % CODON_LENGTH;
 
2074
        
 
2075
        *query_start = query_info->contexts[mixed_frame_context].query_offset;
 
2076
        query_end_pt =
 
2077
            query_info->contexts[mixed_frame_context+CODON_LENGTH-1].query_offset +
 
2078
            query_info->contexts[mixed_frame_context+CODON_LENGTH-1].query_length;
 
2079
        query_length = query_end_pt - *query_start;
 
2080
        single_query->sequence = NULL;
 
2081
        single_query->oof_sequence = concatenated_query->oof_sequence + *query_start;
 
2082
    } else {
 
2083
        *query_start  = query_info->contexts[context].query_offset;
 
2084
        query_length = query_info->contexts[context].query_length;
 
2085
        single_query->sequence = concatenated_query->sequence + *query_start;
 
2086
        single_query->oof_sequence = NULL;
 
2087
    }
 
2088
    single_query->length = query_length;
 
2089
}
 
2090
 
 
2091
 
2024
2092
/** Find the HSP offsets relative to the individual query sequence instead of
2025
 
 * the concatenated sequence.
 
2093
 * the concatenated sequence and retrieve relevant portions of the query
 
2094
 * sequence data.
2026
2095
 * @param query Query sequence block [in]
2027
2096
 * @param query_info Query information structure, including context offsets 
2028
2097
 *                   array [in]
2029
 
 * @param init_hsp Initial HSP [in] [out]
2030
 
 * @param query_out Optional: query sequence block with modified sequence 
 
2098
 * @param init_hsp Initial HSP, this will be modified to hold the adjusted
 
2099
 *                 offsets [in] [out]
 
2100
 * @param query_out query sequence block with modified sequence 
2031
2101
 *                  pointer and sequence length [out]
2032
 
 * @param init_hsp_out Optional: Pointer to initial HSP structure to hold
2033
 
 *                     adjusted offsets; the input init_hsp will be modified if
2034
 
 *                     this parameter is NULL [out]
2035
 
 * @param context_out Which context this HSP belongs to? [out]
 
2102
 * @param context Which context this HSP belongs to? [out]
2036
2103
 */
2037
2104
static void 
2038
 
s_GetRelativeCoordinates(const BLAST_SequenceBlk* query, 
2039
 
                       BlastQueryInfo* query_info, 
2040
 
                       BlastInitHSP* init_hsp, BLAST_SequenceBlk* query_out,
2041
 
                       BlastInitHSP* init_hsp_out, Int4* context_out)
2042
 
{
2043
 
   Int4 context = 0;
2044
 
   Int4 query_start = 0, query_length = 0;
2045
 
 
2046
 
   context = BSearchContextInfo(init_hsp->offsets.qs_offsets.q_off, query_info);
2047
 
 
2048
 
   if (query && query->oof_sequence) {
2049
 
       /* Out-of-frame blastx case: all frames of the same parity are mixed
2050
 
          together in a special sequence. */
2051
 
       Int4 query_end_pt = 0;
2052
 
       Int4 mixed_frame_context = context - context % CODON_LENGTH;
2053
 
       
2054
 
       query_start = query_info->contexts[mixed_frame_context].query_offset;
2055
 
       query_end_pt =
2056
 
           query_info->contexts[mixed_frame_context+CODON_LENGTH-1].query_offset +
2057
 
           query_info->contexts[mixed_frame_context+CODON_LENGTH-1].query_length;
2058
 
       query_length =
2059
 
           query_end_pt - query_start;
2060
 
   } else {
2061
 
       query_start  = query_info->contexts[context].query_offset;
2062
 
       query_length = query_info->contexts[context].query_length;
2063
 
   }
2064
 
   
2065
 
   if (query && query_out) {
2066
 
      if (query->oof_sequence) {
2067
 
         query_out->sequence = NULL;
2068
 
         query_out->oof_sequence = query->oof_sequence + query_start;
2069
 
      } else {
2070
 
         query_out->sequence = query->sequence + query_start;
2071
 
         query_out->oof_sequence = NULL;
2072
 
      }
2073
 
      query_out->length = query_length;
2074
 
   }
2075
 
 
2076
 
   if (init_hsp_out) {
2077
 
      init_hsp_out->offsets.qs_offsets.q_off = 
2078
 
          init_hsp->offsets.qs_offsets.q_off - query_start;
2079
 
      init_hsp_out->offsets.qs_offsets.s_off = 
2080
 
          init_hsp->offsets.qs_offsets.s_off;
2081
 
      if (init_hsp->ungapped_data) {
2082
 
         if (!init_hsp_out->ungapped_data) {
2083
 
            init_hsp_out->ungapped_data = (BlastUngappedData*) 
2084
 
               BlastMemDup(init_hsp->ungapped_data, sizeof(BlastUngappedData));
2085
 
         } else {
2086
 
            memcpy(init_hsp_out->ungapped_data, init_hsp->ungapped_data,
2087
 
                   sizeof(BlastUngappedData));
2088
 
         }
2089
 
         init_hsp_out->ungapped_data->q_start = 
2090
 
            init_hsp->ungapped_data->q_start - query_start;
2091
 
      }
2092
 
   } else {
2093
 
      init_hsp->offsets.qs_offsets.q_off -= query_start;
2094
 
      if (init_hsp->ungapped_data)
2095
 
         init_hsp->ungapped_data->q_start -= query_start;
2096
 
   }
2097
 
 
2098
 
   *context_out = context;
2099
 
}
2100
 
 
2101
 
Int2 BLAST_MbGetGappedScore(EBlastProgramType program_number, 
2102
 
             BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
2103
 
                            BLAST_SequenceBlk* subject,
2104
 
                            BlastGapAlignStruct* gap_align,
2105
 
                            const BlastScoringParameters* score_params, 
2106
 
                            const BlastExtensionParameters* ext_params,
2107
 
                            const BlastHitSavingParameters* hit_params,
2108
 
                            BlastInitHitList* init_hitlist,
2109
 
                            BlastHSPList** hsp_list_ptr, BlastGappedStats* gapped_stats)
2110
 
{
2111
 
   const BlastExtensionOptions* ext_options = ext_params->options;
2112
 
   Int4 index;
2113
 
   BlastInitHSP* init_hsp;
2114
 
   BlastHSPList* hsp_list;
2115
 
   const BlastHitSavingOptions* hit_options = hit_params->options;
2116
 
   BLAST_SequenceBlk query_tmp;
2117
 
   Int4 context;
2118
 
   BlastIntervalTree *tree;
2119
 
   
2120
 
   if (*hsp_list_ptr == NULL)
2121
 
      *hsp_list_ptr = hsp_list = Blast_HSPListNew(hit_options->hsp_num_max);
2122
 
   else 
2123
 
      hsp_list = *hsp_list_ptr;
2124
 
 
2125
 
   tree = Blast_IntervalTreeInit(0, query->length + 1,
2126
 
                                 0, subject->length + 1);
2127
 
 
2128
 
   for (index=0; index<init_hitlist->total; index++) {
2129
 
      BlastHSP tmp_hsp;
2130
 
      Int4 q_start, q_end, s_start, s_end, score;
2131
 
 
2132
 
      init_hsp = &init_hitlist->init_hsp_array[index];
2133
 
 
2134
 
      /* Change query coordinates to relative in the initial HSP right here */
2135
 
      s_GetRelativeCoordinates(query, query_info, init_hsp, &query_tmp, NULL, 
2136
 
                             &context);
2137
 
 
2138
 
      if (!init_hsp->ungapped_data) {
2139
 
         q_start = q_end = init_hsp->offsets.qs_offsets.q_off;
2140
 
         s_start = s_end = init_hsp->offsets.qs_offsets.s_off;
2141
 
         score = INT4_MIN;
2142
 
      } else {
2143
 
         q_start = init_hsp->ungapped_data->q_start;
2144
 
         q_end = q_start + init_hsp->ungapped_data->length;
2145
 
         s_start = init_hsp->ungapped_data->s_start;
2146
 
         s_end = s_start + init_hsp->ungapped_data->length;
2147
 
         score = init_hsp->ungapped_data->score;
2148
 
      }
2149
 
 
2150
 
      tmp_hsp.score = score;
2151
 
      tmp_hsp.context = context;
2152
 
      tmp_hsp.query.offset = q_start;
2153
 
      tmp_hsp.query.end = q_end;
2154
 
      tmp_hsp.query.frame = query_info->contexts[context].frame;
2155
 
      tmp_hsp.subject.offset = s_start;
2156
 
      tmp_hsp.subject.end = s_end;
2157
 
      tmp_hsp.subject.frame = 1;
2158
 
 
2159
 
      if (!BlastIntervalTreeContainsHSP(tree, &tmp_hsp, query_info,
2160
 
                                        MB_DIAG_CLOSE))
2161
 
      {
2162
 
         if (gapped_stats)
2163
 
            ++gapped_stats->extensions;
2164
 
 
2165
 
         BLAST_GreedyGappedAlignment(query_tmp.sequence, 
2166
 
            subject->sequence, query_tmp.length, subject->length, gap_align, 
2167
 
            score_params, init_hsp->offsets.qs_offsets.q_off, 
2168
 
            init_hsp->offsets.qs_offsets.s_off, (Boolean) TRUE, 
2169
 
            (Boolean) (ext_options->ePrelimGapExt == eGreedyWithTracebackExt));
2170
 
 
2171
 
         if (gap_align->score >= hit_options->cutoff_score) {
2172
 
            /* gap_align contains alignment endpoints; init_hsp contains 
2173
 
               the offsets to start the alignment from, if traceback is to 
2174
 
               be performed later */
2175
 
            BlastHSP* new_hsp;
2176
 
            Blast_HSPInit(gap_align->query_start, gap_align->query_stop, 
2177
 
                          gap_align->subject_start, gap_align->subject_stop,
2178
 
                          init_hsp->offsets.qs_offsets.q_off, 
2179
 
                          init_hsp->offsets.qs_offsets.s_off, context, 
2180
 
                          query_info->contexts[context].frame, 1, 
2181
 
                          gap_align->score, &(gap_align->edit_script), 
2182
 
                          &new_hsp);
2183
 
            Blast_HSPListSaveHSP(hsp_list, new_hsp);
2184
 
            BlastIntervalTreeAddHSP(new_hsp, tree, query_info, 
2185
 
                                    eQueryAndSubject);
2186
 
         }
2187
 
         else
2188
 
         {
2189
 
            gap_align->edit_script = GapEditScriptDelete(gap_align->edit_script);
2190
 
         }
2191
 
      }
2192
 
   }
2193
 
 
2194
 
   tree = Blast_IntervalTreeFree(tree);
2195
 
 
2196
 
   return 0;
 
2105
s_AdjustHspOffsetsAndGetQueryData(const BLAST_SequenceBlk* query, 
 
2106
                                  const BlastQueryInfo* query_info, 
 
2107
                                  BlastInitHSP* init_hsp, 
 
2108
                                  BLAST_SequenceBlk* query_out, 
 
2109
                                  Int4* context)
 
2110
{
 
2111
    Int4 query_start = 0;
 
2112
 
 
2113
    ASSERT(query);
 
2114
    ASSERT(query_info);
 
2115
    ASSERT(query_out);
 
2116
    ASSERT(init_hsp);
 
2117
    ASSERT(context);
 
2118
 
 
2119
    *context = s_GetUngappedHSPContext(query_info, init_hsp);
 
2120
    s_SetUpLocalBlastSequenceBlk(query, query_info, *context,
 
2121
                                 query_out, &query_start);
 
2122
    s_AdjustInitialHSPOffsets(init_hsp, query_start);
2197
2123
}
2198
2124
 
2199
2125
 
2210
2136
Blast_PrelimEditBlockToGapEditScript (GapPrelimEditBlock* rev_prelim_tback,
2211
2137
                                      GapPrelimEditBlock* fwd_prelim_tback)
2212
2138
{
2213
 
   GapEditScript* esp_start = NULL,* esp;
 
2139
   Boolean merge_ops = FALSE;
 
2140
   GapEditScript* esp;
2214
2141
   GapPrelimEditScript *op;
2215
2142
   Int4 i;
 
2143
   Int4 index=0;
 
2144
   Int4 size = 0;
2216
2145
 
2217
2146
   if (rev_prelim_tback == NULL || fwd_prelim_tback == NULL)
2218
2147
      return NULL;
2219
2148
 
2220
 
   /* turn the right extension into a linked list (built
2221
 
      last to first). Reverse the edit script in the process */
2222
 
 
2223
 
   for (i=0; i < fwd_prelim_tback->num_ops; i++) {
2224
 
      esp = (GapEditScript*) malloc(sizeof(GapEditScript));
 
2149
   /* The fwd_prelim_tback script will get reversed here as the traceback started from the highest scoring point
 
2150
     and worked backwards. The rev_prelim_tback script does NOT get reversed.  Since it was reversed when the 
 
2151
     traceback was produced it's already "forward" */
 
2152
 
 
2153
   if (fwd_prelim_tback->num_ops > 0 && rev_prelim_tback->num_ops > 0 &&
 
2154
       fwd_prelim_tback->edit_ops[(fwd_prelim_tback->num_ops)-1].op_type == 
 
2155
         rev_prelim_tback->edit_ops[(rev_prelim_tback->num_ops)-1].op_type)
 
2156
     merge_ops = TRUE;
 
2157
 
 
2158
   size = fwd_prelim_tback->num_ops+rev_prelim_tback->num_ops;
 
2159
   if (merge_ops)
 
2160
     size--;
 
2161
 
 
2162
   esp = GapEditScriptNew(size);
 
2163
 
 
2164
   index = 0;
 
2165
   for (i=0; i < rev_prelim_tback->num_ops; i++) {
 
2166
      op = rev_prelim_tback->edit_ops + i;
 
2167
      esp->op_type[index] = op->op_type;
 
2168
      esp->num[index] = op->num;
 
2169
      index++;
 
2170
   }
 
2171
 
 
2172
   if (fwd_prelim_tback->num_ops == 0)
 
2173
      return esp;
 
2174
 
 
2175
   if (merge_ops)
 
2176
       esp->num[index-1] += fwd_prelim_tback->edit_ops[(fwd_prelim_tback->num_ops)-1].num;
 
2177
 
 
2178
   /* If we merge, then we skip the first one. */
 
2179
   if (merge_ops)
 
2180
      i = fwd_prelim_tback->num_ops - 2;
 
2181
   else
 
2182
      i = fwd_prelim_tback->num_ops - 1;
 
2183
 
 
2184
   for (; i >= 0; i--) {
2225
2185
      op = fwd_prelim_tback->edit_ops + i;
2226
 
      esp->op_type = op->op_type;
2227
 
      esp->num = op->num;
2228
 
      esp->next = esp_start;
2229
 
      esp_start = esp;
2230
 
   }
2231
 
 
2232
 
   /* if the first operation in the forward script matches the
2233
 
      last operation in the reverse script, merge the two 
2234
 
      operations */
2235
 
 
2236
 
   if (rev_prelim_tback->num_ops == 0)
2237
 
       return esp_start;
2238
 
 
2239
 
   i = rev_prelim_tback->num_ops - 1;
2240
 
   op = rev_prelim_tback->edit_ops + i;
2241
 
   if (esp_start && esp_start->op_type == op->op_type) {
2242
 
       esp_start->num += op->num;
2243
 
       i--;
2244
 
   }
2245
 
 
2246
 
   /* prepend the reverse script to the list of actions */
2247
 
 
2248
 
   for (; i >= 0; i--) {
2249
 
      esp = (GapEditScript*) malloc(sizeof(GapEditScript));
2250
 
      op = rev_prelim_tback->edit_ops + i;
2251
 
      esp->op_type = op->op_type;
2252
 
      esp->num = op->num;
2253
 
      esp->next = esp_start;
2254
 
      esp_start = esp;
2255
 
   }
2256
 
 
2257
 
   return esp_start;
 
2186
      esp->op_type[index] = op->op_type;
 
2187
      esp->num[index] = op->num;
 
2188
      index++;
 
2189
   }
 
2190
 
 
2191
   return esp;
2258
2192
}
2259
2193
 
2260
2194
/** Fills the BlastGapAlignStruct structure with the results of a gapped 
2723
2657
    return max_offset;
2724
2658
}
2725
2659
 
2726
 
/** Test for an array of ungapped HSPs to decide whether gapped alignment 
2727
 
 * should be performed on all of them. If the best HSP already has score above
2728
 
 * the cutoff for gapped alignments, then just return. Otherwise attempt 
2729
 
 * gapped alignment for all HSPs with scores above gap trigger, and see any of
2730
 
 * them receives score above the cutoff.
2731
 
 * @param program_number Type of BLAST program [in]
2732
 
 * @param query Query sequence [in]
2733
 
 * @param query_info Additional query information [in]
2734
 
 * @param subject Subject sequence [in]
2735
 
 * @param gap_align Gapped alignment structure [in]
2736
 
 * @param score_params Scoring parameters [in]
2737
 
 * @param hit_params Hit saving parameters [in]
2738
 
 * @param init_hitlist Initial hit list obtained by ungapped alignment [in]
2739
 
 * @param gapped_stats Gapped extension stats [in]
2740
 
 * @return TRUE if the set of initial HSPs has passed the test, so gapped 
2741
 
 *         alignment needs to be performed.
2742
 
 */
2743
 
static Boolean 
2744
 
s_BlastGappedScorePrelimTest(EBlastProgramType program_number, 
2745
 
        BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
2746
 
        BLAST_SequenceBlk* subject, 
2747
 
        BlastGapAlignStruct* gap_align,
2748
 
        const BlastScoringParameters* score_params,
2749
 
        const BlastHitSavingParameters* hit_params,
2750
 
        BlastInitHitList* init_hitlist,
2751
 
        BlastGappedStats* gapped_stats)
2752
 
{
2753
 
    BlastInitHSP* init_hsp = NULL;
2754
 
    BlastInitHSP* init_hsp_array;
2755
 
    Int4 init_hsp_count;
2756
 
    Int4 index;
2757
 
    BLAST_SequenceBlk query_tmp;
2758
 
    Int4 context;
2759
 
    Int4 **rpsblast_pssms = NULL;   /* Pointer to concatenated PSSMs in
2760
 
                                       RPS-BLAST database */
2761
 
    Boolean further_process = FALSE;
2762
 
    Int4 cutoff_score;
2763
 
    Boolean is_prot;
2764
 
    Int4 max_offset;
2765
 
    Int2 status = 0;
2766
 
 
2767
 
    cutoff_score = hit_params->cutoff_score;
2768
 
    is_prot = (program_number != eBlastTypeBlastn &&
2769
 
               program_number != eBlastTypePhiBlastn);
2770
 
 
2771
 
    if (program_number == eBlastTypeRpsTblastn || 
2772
 
        program_number == eBlastTypeRpsBlast) {
2773
 
        rpsblast_pssms = gap_align->sbp->psi_matrix->pssm->data;
2774
 
    }
2775
 
 
2776
 
    ASSERT(Blast_InitHitListIsSortedByScore(init_hitlist));
2777
 
 
2778
 
    init_hsp_array = init_hitlist->init_hsp_array;
2779
 
    init_hsp_count = init_hitlist->total;
2780
 
 
2781
 
    /* If no initial HSP passes the score threshold corresponding to the e-value
2782
 
       cutoff so far, check if any would do after gapped alignment, and exit if
2783
 
       none are found. 
2784
 
       Only attempt to extend initial HSPs whose scores are already above 
2785
 
       gap trigger */
2786
 
    
2787
 
    if (init_hsp_array[0].ungapped_data && 
2788
 
        init_hsp_array[0].ungapped_data->score < cutoff_score) {
2789
 
        BlastInitHSP init_hsp_tmp;
2790
 
        init_hsp_tmp.ungapped_data = NULL;
2791
 
        for (index=0; index<init_hsp_count; index++) {
2792
 
            init_hsp = &init_hsp_array[index];
2793
 
 
2794
 
            if (gapped_stats) {
2795
 
                ++gapped_stats->extra_extensions;
2796
 
                ++gapped_stats->extensions;
2797
 
            }
2798
 
 
2799
 
            /* Don't modify initial HSP's coordinates here, because it will be 
2800
 
               done again if further processing is required */
2801
 
            s_GetRelativeCoordinates(query, query_info, init_hsp, &query_tmp, 
2802
 
                                   &init_hsp_tmp, &context);
2803
 
            if (rpsblast_pssms)
2804
 
                gap_align->sbp->psi_matrix->pssm->data = rpsblast_pssms + 
2805
 
                    query_info->contexts[context].query_offset;
2806
 
 
2807
 
            if(is_prot && !score_params->options->is_ooframe) {
2808
 
                max_offset = 
2809
 
                    BlastGetStartForGappedAlignment(query_tmp.sequence, 
2810
 
                                                    subject->sequence, gap_align->sbp,
2811
 
                                                    init_hsp_tmp.ungapped_data->q_start,
2812
 
                                                    init_hsp_tmp.ungapped_data->length,
2813
 
                                                    init_hsp_tmp.ungapped_data->s_start,
2814
 
                                                    init_hsp_tmp.ungapped_data->length);
2815
 
                init_hsp_tmp.offsets.qs_offsets.s_off += 
2816
 
                    max_offset - init_hsp_tmp.offsets.qs_offsets.q_off;
2817
 
                init_hsp_tmp.offsets.qs_offsets.q_off = max_offset;
2818
 
            }
2819
 
 
2820
 
            if (is_prot) {
2821
 
                status =  
2822
 
                    s_BlastProtGappedAlignment(program_number, &query_tmp, 
2823
 
                                              subject, gap_align, score_params, &init_hsp_tmp);
2824
 
            } else {
2825
 
                status = 
2826
 
                    s_BlastDynProgNtGappedAlignment(&query_tmp, subject, 
2827
 
                                                   gap_align, score_params, &init_hsp_tmp);
2828
 
            }
2829
 
            if (status) {
2830
 
                further_process = FALSE;
2831
 
                break;
2832
 
            }
2833
 
            if (gap_align->score >= cutoff_score) {
2834
 
                further_process = TRUE;
2835
 
                break;
2836
 
            }
2837
 
        }
2838
 
        sfree(init_hsp_tmp.ungapped_data);
2839
 
    } else {
2840
 
        index = 0;
2841
 
        further_process = TRUE;
2842
 
        if (gapped_stats)
2843
 
            ++gapped_stats->seqs_ungapped_passed;
2844
 
    }
2845
 
 
2846
 
    if (rpsblast_pssms) {
2847
 
        gap_align->sbp->psi_matrix->pssm->data = rpsblast_pssms;
2848
 
    }
2849
 
 
2850
 
    return further_process;
2851
 
}
2852
 
 
2853
2660
Int2 BLAST_GetGappedScore (EBlastProgramType program_number, 
2854
2661
        BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
2855
2662
        BLAST_SequenceBlk* subject, 
2866
2673
   BlastInitHSP* init_hsp_array;
2867
2674
   Int4 q_start, s_start, q_end, s_end;
2868
2675
   Boolean is_prot;
2869
 
   Int4 max_offset;
 
2676
   Boolean is_greedy;
2870
2677
   Int2 status = 0;
2871
2678
   BlastHSPList* hsp_list = NULL;
2872
2679
   const BlastHitSavingOptions* hit_options = hit_params->options;
2876
2683
   Int4 score;
2877
2684
   Int4 **rpsblast_pssms = NULL;   /* Pointer to concatenated PSSMs in
2878
2685
                                       RPS-BLAST database */
 
2686
   const int kHspNumMax = BlastHspNumMax(TRUE, hit_options);
2879
2687
 
2880
2688
   if (!query || !subject || !gap_align || !score_params || !ext_params ||
2881
2689
       !hit_params || !init_hitlist || !hsp_list_ptr)
2886
2694
 
2887
2695
   is_prot = (program_number != eBlastTypeBlastn &&
2888
2696
              program_number != eBlastTypePhiBlastn);
 
2697
   is_greedy = (ext_params->options->ePrelimGapExt != eDynProgExt);
2889
2698
 
2890
2699
   if (program_number == eBlastTypeRpsTblastn || 
2891
2700
       program_number == eBlastTypeRpsBlast) {
2892
2701
       rpsblast_pssms = gap_align->sbp->psi_matrix->pssm->data;
2893
2702
   }
2894
2703
 
 
2704
   ASSERT(Blast_InitHitListIsSortedByScore(init_hitlist));
 
2705
 
2895
2706
   if (*hsp_list_ptr == NULL)
2896
 
      *hsp_list_ptr = hsp_list = Blast_HSPListNew(hit_options->hsp_num_max);
 
2707
      *hsp_list_ptr = hsp_list = Blast_HSPListNew(kHspNumMax);
2897
2708
   else 
2898
2709
      hsp_list = *hsp_list_ptr;
2899
2710
 
2900
 
   if (!s_BlastGappedScorePrelimTest(program_number, query, query_info, 
2901
 
           subject, gap_align, score_params, hit_params, 
2902
 
           init_hitlist, gapped_stats))
2903
 
      return 0;
2904
 
 
2905
2711
   init_hsp_array = init_hitlist->init_hsp_array;
2906
2712
 
2907
2713
   /* Initialize the interval tree with the maximum possible
2928
2734
      BlastHSP tmp_hsp;
2929
2735
      init_hsp = &init_hsp_array[index];
2930
2736
 
2931
 
      /* Now adjust the initial HSP's coordinates. */
2932
 
      s_GetRelativeCoordinates(query, query_info, init_hsp, &query_tmp, 
2933
 
                             NULL, &context);
 
2737
      s_AdjustHspOffsetsAndGetQueryData(query, query_info, init_hsp, 
 
2738
                                        &query_tmp, &context);
2934
2739
 
2935
2740
      if (rpsblast_pssms)
2936
2741
         gap_align->sbp->psi_matrix->pssm->data = rpsblast_pssms + 
2967
2772
         }
2968
2773
 
2969
2774
         if(is_prot && !score_params->options->is_ooframe) {
2970
 
            max_offset = 
 
2775
            Int4 offset =
2971
2776
               BlastGetStartForGappedAlignment(query_tmp.sequence, 
2972
2777
                  subject->sequence, gap_align->sbp,
2973
2778
                  init_hsp->ungapped_data->q_start,
2975
2780
                  init_hsp->ungapped_data->s_start,
2976
2781
                  init_hsp->ungapped_data->length);
2977
2782
            init_hsp->offsets.qs_offsets.s_off += 
2978
 
                max_offset - init_hsp->offsets.qs_offsets.q_off;
2979
 
            init_hsp->offsets.qs_offsets.q_off = max_offset;
 
2783
                offset - init_hsp->offsets.qs_offsets.q_off;
 
2784
            init_hsp->offsets.qs_offsets.q_off = offset;
2980
2785
         }
2981
2786
 
2982
2787
         if (is_prot) {
2983
2788
            status =  s_BlastProtGappedAlignment(program_number, &query_tmp, 
2984
2789
                         subject, gap_align, score_params, init_hsp);
 
2790
         } else if (is_greedy) {
 
2791
            status = BLAST_GreedyGappedAlignment(
 
2792
                         query_tmp.sequence, subject->sequence, 
 
2793
                         query_tmp.length, subject->length, 
 
2794
                         gap_align, score_params, 
 
2795
                         init_hsp->offsets.qs_offsets.q_off, 
 
2796
                         init_hsp->offsets.qs_offsets.s_off, 
 
2797
                         (Boolean) TRUE, 
 
2798
                         (Boolean) (ext_params->options->ePrelimGapExt == 
 
2799
                                    eGreedyWithTracebackExt));
2985
2800
         } else {
2986
2801
            /*  Start the gapped alignment on the fourth character of the
2987
2802
             *  eight character word that seeded the alignment; the start
3012
2827
                 query_frame = query_info->contexts[context].frame;
3013
2828
             }
3014
2829
 
3015
 
             Blast_HSPInit(gap_align->query_start, 
3016
 
                           gap_align->query_stop, gap_align->subject_start, 
3017
 
                           gap_align->subject_stop, 
 
2830
             Blast_HSPInit(gap_align->query_start, gap_align->query_stop, 
 
2831
                           gap_align->subject_start, gap_align->subject_stop,
3018
2832
                           init_hsp->offsets.qs_offsets.q_off, 
3019
2833
                           init_hsp->offsets.qs_offsets.s_off, context, 
3020
2834
                           query_frame, subject->frame, gap_align->score, 
3023
2837
             BlastIntervalTreeAddHSP(new_hsp, tree, query_info, 
3024
2838
                                     eQueryAndSubject);
3025
2839
         }
 
2840
         else {
 
2841
            /* a greedy alignment may have traceback associated with it;
 
2842
               free that traceback if the alignment will not be used */
 
2843
            if (is_greedy) {
 
2844
               gap_align->edit_script = GapEditScriptDelete(
 
2845
                                             gap_align->edit_script);
 
2846
            }
 
2847
         }
3026
2848
      }
3027
2849
   }   
3028
2850
 
3260
3082
    Int4 last_num;
3261
3083
    GapPrelimEditBlock *tmp_prelim_tback;
3262
3084
    Int4 i, num_nuc;
 
3085
    int extra_needed=0;
3263
3086
    
3264
3087
    /* prepend a substitution, since the input sequences were
3265
3088
       shifted prior to the OOF alignment */
3354
3177
    }
3355
3178
 
3356
3179
    /* form the final edit block */
3357
 
 
3358
 
    *edit_script_ptr = e_script =
3359
 
        Blast_PrelimEditBlockToGapEditScript(rev_prelim_tback,
3360
 
                                         fwd_prelim_tback);
 
3180
    e_script = Blast_PrelimEditBlockToGapEditScript(tmp_prelim_tback, fwd_prelim_tback);
 
3181
    GapPrelimEditBlockFree(tmp_prelim_tback);
3361
3182
 
3362
3183
    /* postprocess the edit script */
3363
 
 
3364
3184
    num_nuc = 0;
3365
 
    while (e_script != NULL) {
3366
 
 
 
3185
    for (i=0; i<e_script->size; i++)
 
3186
    {
 
3187
        int total_actions=0;
3367
3188
        /* Count the number of nucleotides in the next 
3368
3189
           traceback operation and delete any traceback operations
3369
3190
           that would make the alignment too long. This check is
3370
3191
           not needed for in-frame alignment because the
3371
3192
           traceback is never changed after it is computed */
3372
3193
 
3373
 
        last_op = e_script->op_type;
 
3194
        last_op = e_script->op_type[i];
3374
3195
 
3375
3196
        if (last_op == eGapAlignIns)
3376
3197
            last_op = eGapAlignSub;
3377
 
        i = last_op * e_script->num;
3378
 
 
3379
 
        if (num_nuc + i >= nucl_align_length) {
3380
 
            e_script->num = (nucl_align_length - num_nuc + 
 
3198
 
 
3199
        total_actions = last_op * e_script->num[i];
 
3200
 
 
3201
        if (num_nuc + total_actions >= nucl_align_length) {
 
3202
            e_script->num[i] = (nucl_align_length - num_nuc + 
3381
3203
                             last_op - 1) / last_op;
3382
 
            e_script->next = GapEditScriptDelete(e_script->next);
 
3204
            break; /* We delete the rest of the script. */
3383
3205
        }
3384
3206
        else {
3385
 
            num_nuc += i;
3386
 
        }
3387
 
 
3388
 
        /* Only one frame shift operation at a time is
 
3207
            num_nuc += total_actions;;
 
3208
        }
 
3209
    }
 
3210
    e_script->size = MIN(i+1, e_script->size);  /* If we broke out early then we truncate the edit script. */
 
3211
 
 
3212
    extra_needed = 0;
 
3213
    for (i=0; i<e_script->size; i++)
 
3214
    {
 
3215
        if (e_script->op_type[i] % 3 != 0 && e_script->num[i] > 1) {
 
3216
           extra_needed += e_script->num[i] - 1;
 
3217
        }
 
3218
    }
 
3219
 
 
3220
    if (extra_needed)
 
3221
    {
 
3222
        GapEditScript* new_esp = GapEditScriptNew(extra_needed+e_script->size);
 
3223
        int new_esp_i=0;
 
3224
        for (i=0; i<e_script->size; i++)
 
3225
        {
 
3226
           /* Only one frame shift operation at a time is
3389
3227
           allowed in a single edit script element. */
3390
 
        last_op = e_script->op_type;
3391
 
        if (last_op % 3 != 0 && e_script->num > 1) {
3392
 
            Int4 num_ops = e_script->num;
3393
 
            GapEditScript *last = e_script->next;
3394
 
            e_script->next = NULL;
3395
 
            e_script->num = 1;
3396
 
            for (i = 1; i < num_ops; i++) {
3397
 
                e_script = GapEditScriptNew(e_script);
3398
 
                e_script->op_type = last_op;
3399
 
                e_script->num = 1;
3400
 
            }
3401
 
            e_script->next = last;
3402
 
        }
3403
 
        e_script = e_script->next;
 
3228
           new_esp->num[new_esp_i] = e_script->num[i];
 
3229
           new_esp->op_type[new_esp_i] = e_script->op_type[i];
 
3230
           new_esp_i++;
 
3231
           last_op = e_script->op_type[i];
 
3232
           if (last_op % 3 != 0 && e_script->num[i] > 1) {
 
3233
               Int4 num_ops = e_script->num[i];
 
3234
               int esp_index=0;
 
3235
               new_esp->num[new_esp_i-1] = 1;
 
3236
               for (esp_index = 1; esp_index < num_ops; esp_index++) {
 
3237
                   new_esp->num[new_esp_i] = 1;
 
3238
                   new_esp->op_type[new_esp_i] = last_op;
 
3239
                   new_esp_i++;
 
3240
               }
 
3241
           }
 
3242
       }
 
3243
       e_script = GapEditScriptDelete(e_script);
 
3244
       e_script = new_esp;
3404
3245
    }
 
3246
    *edit_script_ptr = e_script;
3405
3247
 
3406
3248
    /* finally, add one to the size of any block of substitutions
3407
3249
       that follows a frame shift op */
3408
 
 
3409
 
    e_script = *edit_script_ptr;
3410
 
    last_op = e_script->op_type;
3411
 
    e_script = e_script->next;
3412
 
    while (e_script != NULL) {
3413
 
        if (e_script->op_type == eGapAlignSub && (last_op % 3) != 0)
3414
 
            e_script->num++;
3415
 
 
3416
 
        last_op = e_script->op_type;
3417
 
        e_script = e_script->next;
 
3250
    last_op = e_script->op_type[0];
 
3251
    for (i=1; i<e_script->size; i++)
 
3252
    {
 
3253
        if (e_script->op_type[i] == eGapAlignSub && (last_op % 3) != 0)
 
3254
            (e_script->num[i])++;
 
3255
 
 
3256
        last_op = e_script->op_type[i];
3418
3257
    }
3419
3258
 
3420
 
    GapPrelimEditBlockFree(tmp_prelim_tback);
3421
3259
    return 0;
3422
3260
}
3423
3261
 
3548
3386
   BlastHSPList* hsp_list = NULL;
3549
3387
   Int4 index;
3550
3388
   BlastInitHSP* init_hsp;
3551
 
   Int4 context;
 
3389
   const int kHspNumMax = BlastHspNumMax(FALSE, hit_options);
3552
3390
 
3553
3391
   /* The BlastHSPList structure can be allocated and passed from outside */
3554
3392
   if (*hsp_list_ptr != NULL)
3565
3403
   for (index = 0; index < init_hitlist->total; ++index) {
3566
3404
      BlastHSP* new_hsp;
3567
3405
      BlastUngappedData* ungapped_data=NULL;
 
3406
      Int4 context = 0;
3568
3407
      init_hsp = &init_hitlist->init_hsp_array[index];
3569
3408
      if (!init_hsp->ungapped_data) 
3570
3409
         continue;
3571
3410
 
 
3411
      if (!hsp_list) {
 
3412
         hsp_list = Blast_HSPListNew(kHspNumMax);
 
3413
         *hsp_list_ptr = hsp_list;
 
3414
      }
3572
3415
      /* Adjust the initial HSP's coordinates in case of concatenated 
3573
3416
         multiple queries/strands/frames */
3574
 
      s_GetRelativeCoordinates(NULL, query_info, init_hsp, NULL, 
3575
 
                             NULL, &context);
3576
 
      if (!hsp_list) {
3577
 
         hsp_list = Blast_HSPListNew(hit_options->hsp_num_max);
3578
 
         *hsp_list_ptr = hsp_list;
3579
 
      }
 
3417
      context = s_GetUngappedHSPContext(query_info, init_hsp);
 
3418
      s_AdjustInitialHSPOffsets(init_hsp,
 
3419
                                query_info->contexts[context].query_offset);
3580
3420
      ungapped_data = init_hsp->ungapped_data;
3581
3421
      Blast_HSPInit(ungapped_data->q_start, 
3582
3422
                    ungapped_data->length+ungapped_data->q_start,