blast_gapalign.c

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: blast_gapalign.c,v $
 * PRODUCTION Revision 1000.5  2004/06/01 18:07:19  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.101
 * PRODUCTION
 * ===========================================================================
 */

/* $Id: blast_gapalign.c,v 1000.5 2004/06/01 18:07:19 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's offical duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Author: Ilya Dondoshansky
 *
 */

/** @file blast_gapalign.c
 * Functions to perform gapped alignment
 */

static char const rcsid[] = 
    "$Id: blast_gapalign.c,v 1000.5 2004/06/01 18:07:19 gouriano Exp $";

#include <algo/blast/core/blast_options.h>
#include <algo/blast/core/blast_def.h>
#include <algo/blast/core/blast_gapalign.h>
#include <algo/blast/core/blast_util.h> /* for NCBI2NA_UNPACK_BASE macros */
#include <algo/blast/core/blast_setup.h>
#include <algo/blast/core/greedy_align.h>

static Int2 BLAST_DynProgNtGappedAlignment(BLAST_SequenceBlk* query_blk, 
   BLAST_SequenceBlk* subject_blk, BlastGapAlignStruct* gap_align, 
   const BlastScoringParameters* score_params, BlastInitHSP* init_hsp);
static Int4 BLAST_AlignPackedNucl(Uint1* B, Uint1* A, Int4 N, Int4 M, 
   Int4* pej, Int4* pei, BlastGapAlignStruct* gap_align,
   const BlastScoringParameters* score_params, Boolean reverse_sequence);

static Int2 BLAST_ProtGappedAlignment(Uint1 program, 
   BLAST_SequenceBlk* query_in, BLAST_SequenceBlk* subject_in,
   BlastGapAlignStruct* gap_align,
   const BlastScoringParameters* score_params, BlastInitHSP* init_hsp);

/** Auxiliary structure for dynamic programming gapped extension */
typedef struct BlastGapDP {
  Int4 best;
  Int4 best_gap;
  Int4 best_decline;
} BlastGapDP;

typedef struct GapData {
  BlastGapDP* CD;
  Int4** v;
  Int4* sapp;                 /* Current script append ptr */
  Int4  last;
  Int4 h,  g;
} GapData;

/** Append "Delete k" op */
#define DEL_(k) \
data.last = (data.last < 0) ? (data.sapp[-1] -= (k)) : (*data.sapp++ = -(k));

/** Append "Insert k" op */
#define INS_(k) \
data.last = (data.last > 0) ? (data.sapp[-1] += (k)) : (*data.sapp++ = (k));

/** Append "Replace" op */
#define REP_ \
{data.last = *data.sapp++ = 0;}

/* Divide by two to prevent underflows. */
#define MININT INT4_MIN/2
#define REPP_ \
{*data.sapp++ = MININT; data.last = 0;}

/** Are the two HSPs within a given number of diagonals from each other? */
#define MB_HSP_CLOSE(q1, q2, s1, s2, c) \
(ABS((q1-s1) - (q2-s2)) < c)

/** Is one HSP contained in a diagonal strip around another? */
#define MB_HSP_CONTAINED(qo1,qo2,qe2,so1,so2,se2,c) \
(qo1>=qo2 && qo1<=qe2 && so1>=so2 && so1<=se2 && \
MB_HSP_CLOSE(qo1,qo2,so1,so2,c))

/** TRUE if c is between a and b; f between d and f. Determines if the
 * coordinates are already in an HSP that has been evaluated. 
 */
#define CONTAINED_IN_HSP(a,b,c,d,e,f) (((a <= c && b >= c) && (d <= f && e >= f)) ? TRUE : FALSE)

/** Callback for sorting HSPs by starting offset in query */ 
static int
query_offset_compare_hsps(const void* v1, const void* v2)
{
	BlastHSP* h1,* h2;
	BlastHSP** hp1,** hp2;

	hp1 = (BlastHSP**) v1;
	hp2 = (BlastHSP**) v2;
	h1 = *hp1;
	h2 = *hp2;

	if (h1 == NULL || h2 == NULL)
		return 0;

   /* If these are from different contexts, don't compare offsets */
   if (h1->context < h2->context) 
      return -1;
   if (h1->context > h2->context)
      return 1;

	if (h1->query.offset < h2->query.offset)
		return -1;
	if (h1->query.offset > h2->query.offset)
		return 1;

	return 0;
}

/** Callback for sorting HSPs by ending offset in query */
static int
query_end_compare_hsps(const void* v1, const void* v2)
{
	BlastHSP* h1,* h2;
	BlastHSP** hp1,** hp2;

	hp1 = (BlastHSP**) v1;
	hp2 = (BlastHSP**) v2;
	h1 = *hp1;
	h2 = *hp2;

	if (h1 == NULL || h2 == NULL)
		return 0;

   /* If these are from different contexts, don't compare offsets */
   if (h1->context < h2->context) 
      return -1;
   if (h1->context > h2->context)
      return 1;

	if (h1->query.end < h2->query.end)
		return -1;
	if (h1->query.end > h2->query.end)
		return 1;

	return 0;
}

/** Callback for sorting HSPs by score */
static int
score_compare_hsp(const void* v1, const void* v2)
{
	BlastHSP* h1,* h2;
	BlastHSP** hp1,** hp2;

	hp1 = (BlastHSP**) v1;
	hp2 = (BlastHSP**) v2;
	h1 = *hp1;
	h2 = *hp2;
        
        /* NULL HSPs are moved to the end */
	if (h1 == NULL)
           return 1;
        if (h2 == NULL)
           return -1;

	if (h1->score < h2->score) 
		return 1;
	if (h1->score > h2->score)
		return -1;

	return 0;
}

/** Check the gapped alignments for an overlap of two different alignments.
 * A sufficient overlap is when two alignments have the same start values
 * of have the same final values. 
 * @param hsp_array Pointer to an array of BlastHSP structures [in]
 * @param hsp_count The size of the hsp_array [in]
 * @param frame What frame these HSPs are from? [in]
 * @return The number of valid alignments remaining. 
*/
static Int4
CheckGappedAlignmentsForOverlap(BlastHSP* *hsp_array, Int4 hsp_count, Int2 frame)

{
   Int4 index1, index, increment;

   if (*hsp_array == NULL || hsp_count == 0)
      return 0;
   
   qsort(hsp_array, hsp_count, sizeof(BlastHSP*), query_offset_compare_hsps);
   index=0;
   increment=1;
   while (index < hsp_count-increment) {
      /* Check if both HSP's start on or end on the same digonal. */
      if (hsp_array[index+increment] == NULL) {
	 increment++;
	 continue;
      }
      
      if (frame != 0 && hsp_array[index+increment]->subject.frame != frame)
	 break;
      
      if (hsp_array[index] && hsp_array[index]->query.offset == hsp_array[index+increment]->query.offset &&
	  hsp_array[index]->subject.offset == hsp_array[index+increment]->subject.offset &&
	  SIGN(hsp_array[index]->query.frame) == SIGN(hsp_array[index+increment]->query.frame))
      {
	 if (hsp_array[index]->score > hsp_array[index+increment]->score) {
	   sfree(hsp_array[index+increment]);
	    increment++;
	 } else {
	    sfree(hsp_array[index]);
	    index++;
	    increment = 1;
	 }
      } else {
	 index++;
	 increment = 1;
      }
   }
   
   qsort(hsp_array, hsp_count, sizeof(BlastHSP*), query_end_compare_hsps);
   index=0;
   increment=1;
   while (index < hsp_count-increment)
   { /* Check if both HSP's start on or end on the same digonal. */
      if (hsp_array[index+increment] == NULL)
      {
	 increment++;
	 continue;
      }
      
      if (frame != 0 && hsp_array[index+increment]->subject.frame != frame)
	 break;
      
      if (hsp_array[index] &&
	  hsp_array[index]->query.end == hsp_array[index+increment]->query.end &&
	  hsp_array[index]->subject.end == hsp_array[index+increment]->subject.end &&
	  SIGN(hsp_array[index]->query.frame) == SIGN(hsp_array[index+increment]->query.frame))
      {
	 if (hsp_array[index]->score > hsp_array[index+increment]->score)
	 {
          sfree(hsp_array[index+increment]);
	    increment++;
	 } else	{
	    sfree(hsp_array[index]);
	    index++;
	    increment = 1;
	 }
      } else {
	 index++;
	 increment = 1;
      }
   }

   qsort(hsp_array,hsp_count,sizeof(BlastHSP*), score_compare_hsp);
   
   index1 = 0;
   for (index=0; index<hsp_count; index++)
   {
      if (hsp_array[index] != NULL)
	 index1++;
   }

   return index1;
}

/** Retrieve the state structure corresponding to a given length
 * @param head Pointer to the first element of the state structures 
 *        array [in]
 * @param length The length for which the state structure has to be 
 *        found [in]
 * @return The found or created state structure
 */
#define	CHUNKSIZE	2097152
static GapStateArrayStruct*
GapGetState(GapStateArrayStruct** head, Int4 length)

{
   GapStateArrayStruct*	retval,* var,* last;
   Int4	chunksize = MAX(CHUNKSIZE, length + length/3);

   length += length/3;	/* Add on about 30% so the end will get reused. */
   retval = NULL;
   if (*head == NULL) {
      retval = (GapStateArrayStruct*) 
         malloc(sizeof(GapStateArrayStruct));
      retval->state_array = (Uint1*) malloc(chunksize*sizeof(Uint1));
      retval->length = chunksize;
      retval->used = 0;
      retval->next = NULL;
      *head = retval;
   } else {
      var = *head;
      last = *head;
      while (var) {
         if (length < (var->length - var->used)) {
            retval = var;
            break;
         } else if (var->used == 0) { 
            /* If it's empty and too small, replace. */
            sfree(var->state_array);
            var->state_array = (Uint1*) malloc(chunksize*sizeof(Uint1));
            var->length = chunksize;
            retval = var;
            break;
         }
         last = var;
         var = var->next;
      }
      
      if (var == NULL)
      {
         retval = (GapStateArrayStruct*) malloc(sizeof(GapStateArrayStruct));
         retval->state_array = (Uint1*) malloc(chunksize*sizeof(Uint1));
         retval->length = chunksize;
         retval->used = 0;
         retval->next = NULL;
         last->next = retval;
      }
   }

#ifdef ERR_POST_EX_DEFINED
   if (retval->state_array == NULL)
      ErrPostEx(SEV_ERROR, 0, 0, "state array is NULL");
#endif
		
   return retval;

}

/** Remove a state from a state structure */
static Boolean
GapPurgeState(GapStateArrayStruct* state_struct)
{
   while (state_struct)
   {
      /*
	memset(state_struct->state_array, 0, state_struct->used);
      */
      state_struct->used = 0;
      state_struct = state_struct->next;
   }
   
   return TRUE;
}

/** Greatest common divisor */
static Int4 gcd(Int4 a, Int4 b)
{
    Int4 c;
    if (a < b) {
        c = a;
        a = b; b = c;
    }
    while ((c = a % b) != 0) {
        a = b; b = c;
    }

    return b;
}

/** Divide 3 numbers by their greatest common divisor
 * @return The GCD
 */
static Int4 gdb3(Int4* a, Int4* b, Int4* c)
{
    Int4 g;
    if (*b == 0) g = gcd(*a, *c);
    else g = gcd(*a, gcd(*b, *c));
    if (g > 1) {