blast_gapalign.c

ncbi源码

        *a /= g;
        *b /= g;
        *c /= g;
    }
    return g;
}

/** Deallocate the memory for greedy gapped alignment */
static SGreedyAlignMem* BLAST_GreedyAlignsfree(SGreedyAlignMem* gamp)
{
   if (gamp->flast_d) {
      sfree(gamp->flast_d[0]);
      sfree(gamp->flast_d);
   } else {
      if (gamp->flast_d_affine) {
         sfree(gamp->flast_d_affine[0]);
         sfree(gamp->flast_d_affine);
      }
      sfree(gamp->uplow_free);
   }
   sfree(gamp->max_row_free);
   if (gamp->space)
      MBSpaceFree(gamp->space);
   sfree(gamp);
   return NULL;
}

/** Allocate memory for the greedy gapped alignment algorithm
 * @param score_params Parameters related to scoring [in]
 * @param ext_params Options and parameters related to the extension [in]
 * @param max_dbseq_length The length of the longest sequence in the 
 *        database [in]
 * @return The allocated SGreedyAlignMem structure
 */
static SGreedyAlignMem* 
BLAST_GreedyAlignMemAlloc(const BlastScoringParameters* score_params,
		       const BlastExtensionParameters* ext_params,
		       Int4 max_dbseq_length)
{
#define ERROR_FRACTION 2    /* N.B.: This value should match the value of
                               ERROR_FRACTION in the anonymous enum in
                               greedy_align.c */
#define ICEIL(x,y) ((((x)-1)/(y))+1)    /* FIXME: duplicated from greedy_align.c */
   SGreedyAlignMem* gamp;
   Int4 max_d, max_d_1, Xdrop, d_diff, max_cost, gd, i;
   Int4 reward, penalty, gap_open, gap_extend;
   Int4 Mis_cost, GE_cost;
   Boolean do_traceback;
   
   if (score_params == NULL || ext_params == NULL) 
      return NULL;
   
   do_traceback = 
      (ext_params->options->ePrelimGapExt != eGreedyExt);

   if (score_params->reward % 2 == 1) {
      reward = 2*score_params->reward;
      penalty = -2*score_params->penalty;
      Xdrop = 2*ext_params->gap_x_dropoff;
      gap_open = 2*score_params->gap_open;
      gap_extend = 2*score_params->gap_extend;
   } else {
      reward = score_params->reward;
      penalty = -score_params->penalty;
      Xdrop = ext_params->gap_x_dropoff;
      gap_open = score_params->gap_open;
      gap_extend = score_params->gap_extend;
   }

   if (gap_open == 0 && gap_extend == 0)
      gap_extend = reward / 2 + penalty;

   max_d = (Int4) (max_dbseq_length / ERROR_FRACTION + 1);

   gamp = (SGreedyAlignMem*) calloc(1, sizeof(SGreedyAlignMem));

   if (score_params->gap_open==0 && score_params->gap_extend==0) {
      d_diff = ICEIL(Xdrop+reward/2, penalty+reward);
   
      gamp->flast_d = (Int4**) malloc((max_d + 2) * sizeof(Int4*));
      if (gamp->flast_d == NULL) {
         sfree(gamp);
         return NULL;
      }
      gamp->flast_d[0] = 
         (Int4*) malloc((max_d + max_d + 6) * sizeof(Int4) * 2);
      if (gamp->flast_d[0] == NULL) {
#ifdef ERR_POST_EX_DEFINED
	 ErrPostEx(SEV_WARNING, 0, 0, 
              "Failed to allocate %ld bytes for greedy alignment", 
              (max_d + max_d + 6) * sizeof(Int4) * 2);
#endif
         BLAST_GreedyAlignsfree(gamp);
         return NULL;
      }

      gamp->flast_d[1] = gamp->flast_d[0] + max_d + max_d + 6;
      gamp->flast_d_affine = NULL;
      gamp->uplow_free = NULL;
   } else {
      gamp->flast_d = NULL;
      Mis_cost = reward + penalty;
      GE_cost = gap_extend+reward/2;
      max_d_1 = max_d;
      max_d *= GE_cost;
      max_cost = MAX(Mis_cost, gap_open+GE_cost);
      gd = gdb3(&Mis_cost, &gap_open, &GE_cost);
      d_diff = ICEIL(Xdrop+reward/2, gd);
      gamp->uplow_free = (Int4*) calloc(2*(max_d+1+max_cost), sizeof(Int4));
      gamp->flast_d_affine = (SThreeVal**) 
	 malloc((MAX(max_d, max_cost) + 2) * sizeof(SThreeVal*));
      if (!gamp->uplow_free || !gamp->flast_d_affine) {
         BLAST_GreedyAlignsfree(gamp);
         return NULL;
      }
      gamp->flast_d_affine[0] = (SThreeVal*)
	 calloc((2*max_d_1 + 6) , sizeof(SThreeVal) * (max_cost+1));
      for (i = 1; i <= max_cost; i++)
	 gamp->flast_d_affine[i] = 
	    gamp->flast_d_affine[i-1] + 2*max_d_1 + 6;
      if (!gamp->flast_d_affine || !gamp->flast_d_affine[0])
         BLAST_GreedyAlignsfree(gamp);
   }
   gamp->max_row_free = (Int4*) malloc(sizeof(Int4) * (max_d + 1 + d_diff));

   if (do_traceback)
      gamp->space = MBSpaceNew();
   if (!gamp->max_row_free || (do_traceback && !gamp->space))
      /* Failure in one of the memory allocations */
      BLAST_GreedyAlignsfree(gamp);

   return gamp;
}

/* Documented in blast_gapalign.h */
BlastGapAlignStruct* 
BLAST_GapAlignStructFree(BlastGapAlignStruct* gap_align)
{
   GapEditBlockDelete(gap_align->edit_block);
   if (gap_align->greedy_align_mem)
      BLAST_GreedyAlignsfree(gap_align->greedy_align_mem);
   GapStateFree(gap_align->state_struct);

   sfree(gap_align);
   return NULL;
}

/* Documented in blast_gapalign.h */
Int2
BLAST_GapAlignStructNew(const BlastScoringParameters* score_params, 
   const BlastExtensionParameters* ext_params, 
   Uint4 max_subject_length, Int4 query_length, 
   BlastScoreBlk* sbp, BlastGapAlignStruct** gap_align_ptr)
{
   Int2 status = 0;
   BlastGapAlignStruct* gap_align;

   /* If pointer to output structure is NULL, just don't do anything */ 
   if (!gap_align_ptr)
      return 0;

   if (!gap_align_ptr || !sbp || !score_params || !ext_params)
      return -1;

   gap_align = (BlastGapAlignStruct*) calloc(1, sizeof(BlastGapAlignStruct));

   *gap_align_ptr = gap_align;

   gap_align->sbp = sbp;

   gap_align->gap_x_dropoff = ext_params->gap_x_dropoff;

   if (ext_params->options->ePrelimGapExt != eDynProgExt) {
      max_subject_length = MIN(max_subject_length, MAX_DBSEQ_LEN);
      gap_align->greedy_align_mem = 
         BLAST_GreedyAlignMemAlloc(score_params, ext_params, 
                                max_subject_length);
      if (!gap_align->greedy_align_mem)
         gap_align = BLAST_GapAlignStructFree(gap_align);
   }

   if (!gap_align)
      return -1;

   gap_align->positionBased = (sbp->posMatrix != NULL);

   return status;
}

/** Low level function to perform dynamic programming gapped extension 
 * with 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 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.
*/

#define SCRIPT_SUB      0x00
#define SCRIPT_INS      0x01
#define SCRIPT_DEL      0x02
#define SCRIPT_DECLINE  0x03

#define SCRIPT_OP_MASK  0x03

#define SCRIPT_OPEN_GAP 0x04
#define SCRIPT_INS_ROW  0x10
#define SCRIPT_DEL_ROW  0x20
#define SCRIPT_INS_COL  0x40
#define SCRIPT_DEL_COL  0x80

Int4
ALIGN_EX(Uint1* A, Uint1* B, Int4 M, Int4 N, Int4* S, Int4* a_offset, 
	Int4* b_offset, Int4** sapp, BlastGapAlignStruct* gap_align, 
	const BlastScoringParameters* score_params, Int4 query_offset, 
        Boolean reversed, Boolean reverse_sequence)
	
{ 
    /* See SEMI_G_ALIGN_EX for more general comments on 
       what this code is doing; comments in this function
       only apply to the traceback computations */

    BlastGapDP* score_array;
    Int4 i; 
    Int4 a_index;
    Int4 b_index, b_size, first_b_index, last_b_index, b_increment;
    Uint1* b_ptr;
  
    Int4 gap_open;
    Int4 gap_extend;
    Int4 gap_open_extend;
    Int4 decline_penalty;
    Int4 x_dropoff;
    Int4 best_score;
  
    Int4* *matrix;
    Int4* matrix_row;
  
    Int4 score;
    Int4 score_gap_row;
    Int4 score_gap_col;
    Int4 score_decline;
    Int4 next_score;
    Int4 next_score_decline;
  
    GapData data;                       /* traceback variables */
    GapStateArrayStruct* state_struct;
    Uint1** edit_script;
    Uint1* edit_script_row;
    Int4 orig_b_index;
    Uint1 script, next_script, script_row, script_col;
    Int4 align_len;

    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;
  
    /* Initialize traceback information. edit_script[] is
       a 2-D array which is filled in row by row as the 
       dynamic programming takes place */

    data.sapp = S;
    data.last = 0;
    *sapp = S;
    GapPurgeState(gap_align->state_struct);
    edit_script = (Uint1**) malloc(sizeof(Uint1*)*(M+1));
    if (gap_extend > 0)
        state_struct = GapGetState(&gap_align->state_struct, 
                                   x_dropoff / gap_extend + 5);
    else
        state_struct = GapGetState(&gap_align->state_struct, N + 3);
    edit_script[0] = state_struct->state_array;
    edit_script_row = state_struct->state_array;

    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;
        edit_script_row[i] = SCRIPT_INS;
    }
    state_struct->used = i + 1;
  
    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++) {

        /* Set up the next row of the edit script; this involves
           allocating memory for the row, then pointing to it */

        if (gap_extend > 0)
            state_struct = GapGetState(&gap_align->state_struct, 
                          b_size - first_b_index + x_dropoff / gap_extend + 5);
        else
            state_struct = GapGetState(&gap_align->state_struct, 
                          N + 3 - first_b_index);
        edit_script[a_index] = state_struct->state_array + 
                                state_struct->used + 1;
        edit_script_row = edit_script[a_index];
        orig_b_index = first_b_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];

        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;

            /* script, script_row and script_col contain the
               actions specified by the dynamic programming.
               when the inner loop has finished, 'script' con-
               tains all of the actions to perform, and is
               written to edit_script[a_index][b_index]. Otherwise,
               this inner loop is exactly the same as the one
               in SEMI_G_ALIGN_EX() */

            if (score_decline > score) {
                script = SCRIPT_DECLINE;
                score = score_decline;
            }
            else {
                script = SCRIPT_SUB;
            }

            if (score_gap_col < score_decline) {
                score_gap_col = score_decline;
                script_col = SCRIPT_DEL_COL;
            }
            else {
                script_col = SCRIPT_INS_COL;
                if (score < score_gap_col) {
                    script = SCRIPT_DEL;
                    score = score_gap_col;