blast_gapalign.c

ncbi源码

                }
            }

            if (score_gap_row < score_decline) {
                score_gap_row = score_decline;
                script_row = SCRIPT_DEL_ROW;
            }
            else {
                script_row = SCRIPT_INS_ROW;
                if (score < score_gap_row) {
                    script = SCRIPT_INS;
                    score = score_gap_row;
                }
            }

            if (best_score - score > x_dropoff) {

                if (first_b_index == b_index)
                    first_b_index++;
                else
                    score_array[b_index].best = MININT;
            }
            else {
                last_b_index = b_index;
                if (score > best_score) {
                    best_score = score;
                    *a_offset = a_index;
                    *b_offset = b_index;
                }

                score_gap_row -= gap_extend;
                score_gap_col -= gap_extend;
                if (score_gap_col < (score - gap_open_extend)) {
                    score_array[b_index].best_gap = score - gap_open_extend;
                }
                else {
                    score_array[b_index].best_gap = score_gap_col;
                    script += script_col;
                }

                if (score_gap_row < (score - gap_open_extend)) 
                    score_gap_row = score - gap_open_extend;
                else 
                    script += script_row;

                if (score_decline < (score - gap_open)) {
                    score_array[b_index].best_decline = score - gap_open - decline_penalty;
                }
                else {
                    score_array[b_index].best_decline = score_decline - decline_penalty;
                    script += SCRIPT_OPEN_GAP;
                }
                score_array[b_index].best = score;
            }

            score = next_score;
            score_decline = next_score_decline;
            edit_script_row[b_index] = script;
        }
  
        if (first_b_index == b_size) {
            a_index++;
            break;
        }

        if (last_b_index < b_size - 1) {
            b_size = last_b_index + 1;
        }
        else {
            while (score_gap_row >= (best_score - x_dropoff) && b_size <= N) {
                score_array[b_size].best = score_gap_row;
                score_array[b_size].best_gap = score_gap_row - gap_open_extend;
                score_array[b_size].best_decline = score_gap_row - gap_open -
                                                        decline_penalty;
                score_gap_row -= gap_extend;
                edit_script_row[b_size] = SCRIPT_INS;
                b_size++;
            }
        }

        if (b_size <= N) {
            score_array[b_size].best = MININT;
            score_array[b_size].best_gap = MININT;
            score_array[b_size].best_decline = MININT;
            b_size++;
        }
        state_struct->used += MAX(b_index, b_size) - orig_b_index + 1;
    }
    
    /* Pick the optimal path through the now complete
       edit_script[][]. This is equivalent to flattening 
       the 2-D array into a 1-D list of actions. */

    a_index = *a_offset;
    b_index = *b_offset;
    script = 0;
    edit_script_row = (Uint1 *)malloc(a_index + b_index);

    for (i = 0; a_index > 0 || b_index > 0; i++) {
        next_script = edit_script[a_index][b_index];

        switch(script) {
        case SCRIPT_INS:
            script = next_script & SCRIPT_OP_MASK;
            if (next_script & SCRIPT_INS_ROW)
                script = SCRIPT_INS;
            else if (next_script & SCRIPT_DEL_ROW)
                script = SCRIPT_DECLINE;
            break;

        case SCRIPT_DEL:
            script = next_script & SCRIPT_OP_MASK;
            if (next_script & SCRIPT_INS_COL)
                script = SCRIPT_DEL;
            else if (next_script & SCRIPT_DEL_COL)
                script = SCRIPT_DECLINE;
            break;

        case SCRIPT_DECLINE:
            script = next_script & SCRIPT_OP_MASK;
            if (next_script & SCRIPT_OPEN_GAP)
                script = SCRIPT_DECLINE;
            break;

        default:
            script = next_script & SCRIPT_OP_MASK;
            break;
        }

        if (script == SCRIPT_INS)
            b_index--;
        else if (script == SCRIPT_DEL)
            a_index--;
        else {
            a_index--;
            b_index--;
        }
        edit_script_row[i] = script;
    }

    /* Traceback proceeded backwards through edit_script, 
       so the output traceback information is written 
       in reverse order */

    align_len = i;
    for (i--; i >= 0; i--) {
        if (edit_script_row[i] == SCRIPT_SUB) 
            REP_
        else if (edit_script_row[i] == SCRIPT_INS) 
            INS_(1)
        else if (edit_script_row[i] == SCRIPT_DECLINE) 
            REPP_                     
        else 
            DEL_(1)      
    }
      
    sfree(edit_script_row);
    sfree(edit_script);
    sfree(score_array);
    align_len -= data.sapp - S;
    if (align_len > 0)
        memset(data.sapp, 0, align_len);
    *sapp = data.sapp;

    return best_score;
}

/** Low level function to perform gapped extension in one direction with 
 * or without traceback.
 * @param A The query sequence [in]
 * @param B The subject sequence [in]
 * @param M Maximal extension length in query [in]
 * @param N Maximal extension length in subject [in]
 * @param S The traceback information from previous extension [in]
 * @param a_offset Resulting starting offset in query [out]
 * @param b_offset Resulting starting offset in subject [out]
 * @param score_only Only find the score, without saving traceback [in]
 * @param sapp The traceback information [out]
 * @param gap_align Structure holding various information and allocated 
 *        memory for the gapped alignment [in]
 * @param score_params Parameters related to scoring [in]
 * @param query_offset The starting offset in query [in]
 * @param reversed Has the sequence been reversed? Used for psi-blast [in]
 * @param reverse_sequence Do reverse the sequence [in]
 * @return The best alignment score found.
 */
static Int4 SEMI_G_ALIGN_EX(Uint1* A, Uint1* B, Int4 M, Int4 N,
   Int4* S, Int4* a_offset, Int4* b_offset, Boolean score_only, Int4** sapp,
   BlastGapAlignStruct* gap_align, const BlastScoringParameters* score_params, 
   Int4 query_offset, Boolean reversed, Boolean reverse_sequence)
{
    BlastGapDP* score_array;            /* sequence pointers and indices */
    Int4 i; 
    Int4 a_index;
    Int4 b_index, b_size, first_b_index, last_b_index, b_increment;
    Uint1* b_ptr;
  
    Int4 gap_open;              /* alignment penalty variables */
    Int4 gap_extend;
    Int4 gap_open_extend;
    Int4 decline_penalty;
    Int4 x_dropoff;
  
    Int4* *matrix;              /* pointers to the score matrix */
    Int4* matrix_row;
  
    Int4 score;                 /* score tracking variables */
    Int4 score_gap_row;
    Int4 score_gap_col;
    Int4 score_decline;
    Int4 next_score;
    Int4 next_score_decline;
    Int4 best_score;
  
    if (!score_only) {
        return ALIGN_EX(A, B, M, N, S, a_offset, b_offset, sapp, gap_align, 
                      score_params, query_offset, reversed, reverse_sequence);
    }
    
    /* do initialization and sanity-checking */

    matrix = gap_align->sbp->matrix;
    *a_offset = 0;
    *b_offset = 0;
    gap_open = score_params->gap_open;
    gap_extend = score_params->gap_extend;
    gap_open_extend = gap_open + gap_extend;
    decline_penalty = score_params->decline_align;
    x_dropoff = gap_align->gap_x_dropoff;
  
    if (x_dropoff < gap_open_extend)
        x_dropoff = gap_open_extend;
  
    if(N <= 0 || M <= 0) 
        return 0;
  
    /* Allocate and fill in the auxiliary 
       bookeeping structures (one struct
       per letter of B) */

    score_array = (BlastGapDP*)malloc((N + 2) * sizeof(BlastGapDP));
    score = -gap_open_extend;
    score_array[0].best = 0;
    score_array[0].best_gap = -gap_open_extend;
    score_array[0].best_decline = -gap_open_extend - decline_penalty;
  
    for (i = 1; i <= N; i++) {
        if (score < -x_dropoff) 
            break;

        score_array[i].best = score;
        score_array[i].best_gap = score - gap_open_extend; 
        score_array[i].best_decline = score - gap_open_extend - decline_penalty;
        score -= gap_extend;
    }
  
    /* The inner loop below examines letters of B from 
       index 'first_b_index' to 'b_size' */

    b_size = i;
    best_score = 0;
    first_b_index = 0;
    if (reverse_sequence)
        b_increment = -1;
    else
        b_increment = 1;
  
    for (a_index = 1; a_index <= M; a_index++) {
        /* pick out the row of the score matrix 
           appropriate for A[a_index] */

        if (!(gap_align->positionBased)) {
            if(reverse_sequence)
                matrix_row = matrix[ A[ M - a_index ] ];
            else
                matrix_row = matrix[ A[ a_index ] ];
        }
        else {
            if(reversed || reverse_sequence)
                matrix_row = gap_align->sbp->posMatrix[M - a_index];
            else 
                matrix_row = gap_align->sbp->posMatrix[a_index + query_offset];
        }

        if(reverse_sequence)
            b_ptr = &B[N - first_b_index];
        else
            b_ptr = &B[first_b_index];

        /* initialize running-score variables */
        score = MININT;
        score_gap_row = MININT;
        score_decline = MININT;
        last_b_index = first_b_index;

        for (b_index = first_b_index; b_index < b_size; b_index++) {

            b_ptr += b_increment;
            score_gap_col = score_array[b_index].best_gap;
            next_score = score_array[b_index].best + matrix_row[ *b_ptr ];
            next_score_decline = score_array[b_index].best_decline;

            /* decline the alignment if that improves the score */

            score = MAX(score, score_decline);
            
            /* decline the best row score if that improves it;
               if not, make it the new high score if it's
               an improvement */

            if (score_gap_col < score_decline)
                score_gap_col = score_decline;
            else if (score < score_gap_col)
                score = score_gap_col;

            /* decline the best column score if that improves it;
               if not, make it the new high score if it's
               an improvement */

            if (score_gap_row < score_decline)
                score_gap_row = score_decline;
            else if (score < score_gap_row)
                score = score_gap_row;

            if (best_score - score > x_dropoff) {

                /* the current best score failed the X-dropoff
                   criterion. Note that this does not stop the
                   inner loop, only forces future iterations to
                   skip this column of B. 

                   Also, if the very first letter of B that was
                   tested failed the X dropoff criterion, make
                   sure future inner loops start one letter to 
                   the right */

                if (b_index == first_b_index)
                    first_b_index++;
                else
                    score_array[b_index].best = MININT;
            }
            else {
                last_b_index = b_index;
                if (score > best_score) {
                    best_score = score;
                    *a_offset = a_index;
                    *b_offset = b_index;
                }

                /* If starting a gap at this position will improve
                   the best row, column, or declined alignment score, 
                   update them to reflect that. */

                score_gap_row -= gap_extend;
                score_gap_col -= gap_extend;
                score_array[b_index].best_gap = MAX(score - gap_open_extend,
                                                    score_gap_col);
                score_gap_row = MAX(score - gap_open_extend, score_gap_row);

                score_array[b_index].best_decline = 
                        MAX(score_decline, score - gap_open) - decline_penalty;
                score_array[b_index].best = score;
            }

            score = next_score;
            score_decline = next_score_decline;
        }
  
        /* Finish aligning if the best scores for all positions
           of B will fail the X-dropoff test, i.e. the inner loop 
           bounds have converged to each other */

        if (first_b_index == b_size)
            break;

        if (last_b_index < b_size - 1) {
            /* This row failed the X-dropoff test earlier than
               the last row did; just shorten the loop bounds
               before doing the next row */

            b_size = last_b_index + 1;
        }
        else {
            /* The inner loop finished without failing the X-dropoff
               test; initialize extra bookkeeping structures until
               the X dropoff test fails or we run out of letters in B. 
               The next inner loop will have larger bounds */

            while (score_gap_row >= (best_score - x_dropoff) && b_size <= N) {
                score_array[b_size].best = score_gap_row;
                score_array[b_size].best_gap = score_gap_row - gap_open_extend;
                score_array[b_size].best_decline = score_gap_row - gap_open -