blast_kappa.c

ncbi源码

    self->worstEvalue = self->heapArray[1].bestEvalue;
    KAPPA_ASSERT(SWheapIsValid(self->heapArray, 1, self->n));
  }
  /* end else "self" is currently a heap. */

  /* The matchRecord->alignments pointer is no longer valid */
  matchRecord->alignments = NULL;
}


/** Return true if only matches with evalue <= self->ecutoff 
 *  may be inserted. 
 * @param self heap containing data [in]
 */
static Boolean
SWheapWillAcceptOnlyBelowCutoff(SWheap * self)
{
  return self->n >= self->heapThreshold && self->worstEvalue <= self->ecutoff;
}


/** Initialize a new SWheap; parameters to this function correspond
 * directly to fields in the SWheap 
 * @param self the object to be filled [in|out]
 * @param capacity size of heap [in]
 * @param heapThreshold  items always inserted if fewer than this number in heap [in]
 * @param ecutoff items with a expect value less than this will always be inserted into heap [in]
 */
static void
SWheapInitialize(SWheap * self,
                 Int4 capacity,
                 Int4 heapThreshold,
                 double ecutoff)
{
  self->n             = 0;
  self->heapThreshold = heapThreshold;
  self->ecutoff       = ecutoff;
  self->heapArray     = NULL;
  self->capacity      = 0;
  self->worstEvalue   = 0;
  /* Begin life as a list */
  self->array =
    (SWheapRecord *) calloc(1, (capacity + 1) * sizeof(SWheapRecord));
  self->capacity      = capacity;
}


/** Release the storage associated with the fields of a SWheap. Don't
 * delete the SWheap structure itself.
 * @param self record to be cleared [in][out]
 */
static void
SWheapRelease(SWheap * self)
{
  if(self->heapArray) free(self->heapArray);
  if(self->array) free(self->array);

  self->n = self->capacity = self->heapThreshold = 0;
  self->heapArray = NULL;
}


/** Remove and return the element in the SWheap with largest (worst) evalue
 * @param self heap that contains record to be removed [in]
 * @return record that was removed
 */
static SWResults *
SWheapPop(SWheap * self)
{
  SWResults *results = NULL;

  ConvertToHeap(self);
  if(self->n > 0) { /* The heap is not empty */
    SWheapRecord *first, *last; /* The first and last elements of the
                                 * array that represents the heap. */
    first = &self->heapArray[1];
    last  = &self->heapArray[self->n];

    results = first->theseAlignments;
    
    first->theseAlignments = last->theseAlignments;
    first->bestEvalue      = last->bestEvalue;

    SWheapifyDown(self->heapArray, 1, --self->n);
  }
  
  KAPPA_ASSERT(SWheapIsValid(self->heapArray, 1, self->n));

  return results;
}


/** Convert a SWheap to a flat list of SWResults. Note that there
 * may be more than one alignment per match.  The list of all
 * alignments are sorted by the following keys:
 * - First by the evalue the best alignment between the query and a
 *   particular matching sequence;
 * - Second by the subject_index of the matching sequence; and
 * - Third by the evalue of each individual alignment.
 * @param self heap to be "flattened" [in]
 * @return "flattened" version of the input 
 */
static SWResults *
SWheapToFlatList(SWheap * self)
{
  SWResults *list = NULL;       /* the new list of SWResults */
  SWResults *result;            /* the next list of alignments to be
                                   prepended to "list" */

  while(NULL != (result = SWheapPop(self))) {
    SWResults *head, *remaining;     /* The head and remaining
                                        elements in a list of
                                        alignments to be prepended to
                                        "list" */
    remaining = result;
    while(NULL != (head = remaining)) {
      remaining   = head->next;
      head->next  = list;
      list        = head;
    }
  }

  return list;
}

/** keeps one row of the Smith-Waterman matrix
 */
typedef struct SWpairs {
  Int4 noGap;
  Int4 gapExists;
} SWpairs;


/** computes Smith-Waterman local alignment score and returns the
 * evalue
 *
 * @param matchSeq is a database sequence matched by this query [in]
 * @param matchSeqLength is the length of matchSeq in amino acids [in]
 * @param query is the input query sequence [in]
 * @param queryLength is the length of query [in]
 * @param matrix is the position-specific matrix associated with query [in]
 * @param gapOpen is the cost of opening a gap [in]
 * @param gapExtend is the cost of extending an existing gap by 1 position [in]
 * @param matchSeqEnd returns the final position in the matchSeq of an optimal
 *  local alignment [in]
 * @param queryEnd returns the final position in query of an optimal
 *  local alignment [in]
 * matchSeqEnd and queryEnd can be used to run the local alignment in reverse
 *  to find optimal starting positions [in]
 * @param score is used to pass back the optimal score [in]
 * @param kbp holds the Karlin-Altschul parameters [in]
 * @param effSearchSpace effective search space for calculation of expect value [in]
 * @param positionSpecific determines whether matrix is position specific or not [in]
 * @return the expect value of the alignment
*/

static double BLbasicSmithWatermanScoreOnly(Uint1 * matchSeq, 
   Int4 matchSeqLength, Uint1 *query, Int4 queryLength, Int4 **matrix, 
   Int4 gapOpen, Int4 gapExtend,  Int4 *matchSeqEnd, Int4 *queryEnd, Int4 *score,
   Blast_KarlinBlk* kbp, Int8 effSearchSpace, Boolean positionSpecific)
{

   Int4 bestScore; /*best score seen so far*/
   Int4 newScore;  /* score of next entry*/
   Int4 bestMatchSeqPos, bestQueryPos; /*position ending best score in
                           matchSeq and query sequences*/
   SWpairs *scoreVector; /*keeps one row of the Smith-Waterman matrix
                           overwrite old row with new row*/
   Int4 *matrixRow; /*one row of score matrix*/
   Int4 newGapCost; /*cost to have a gap of one character*/
   Int4 prevScoreNoGapMatchSeq; /*score one row and column up
                               with no gaps*/
   Int4 prevScoreGapMatchSeq;   /*score if a gap already started in matchSeq*/
   Int4 continueGapScore; /*score for continuing a gap in matchSeq*/
   Int4 matchSeqPos, queryPos; /*positions in matchSeq and query*/
   double returnEvalue; /*e-value to return*/


   scoreVector = (SWpairs *) calloc(1, matchSeqLength * sizeof(SWpairs));
   bestMatchSeqPos = 0;
   bestQueryPos = 0;
   bestScore = 0;
   newGapCost = gapOpen + gapExtend;
   for (matchSeqPos = 0; matchSeqPos < matchSeqLength; matchSeqPos++) {
     scoreVector[matchSeqPos].noGap = 0;
     scoreVector[matchSeqPos].gapExists = -(gapOpen);
   }
   for(queryPos = 0; queryPos < queryLength; queryPos++) {  
     if (positionSpecific)
       matrixRow = matrix[queryPos];
     else
       matrixRow = matrix[query[queryPos]];
     newScore = 0;
     prevScoreNoGapMatchSeq = 0;
     prevScoreGapMatchSeq = -(gapOpen);
     for(matchSeqPos = 0; matchSeqPos < matchSeqLength; matchSeqPos++) {
       /*testing scores with a gap in matchSeq, either starting a new
         gap or extending an existing gap*/
       if ((newScore = newScore - newGapCost) > 
	   (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
         prevScoreGapMatchSeq = newScore;
       /*testing scores with a gap in query, either starting a new
         gap or extending an existing gap*/
       if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
           (continueGapScore = scoreVector[matchSeqPos].gapExists - gapExtend))
         continueGapScore = newScore;
       /*compute new score extending one position in matchSeq and query*/
       newScore = prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
       if (newScore < 0)
       newScore = 0; /*Smith-Waterman locality condition*/
       /*test two alternatives*/
       if (newScore < prevScoreGapMatchSeq)
         newScore = prevScoreGapMatchSeq;
       if (newScore < continueGapScore)
         newScore = continueGapScore;
       prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap; 
       scoreVector[matchSeqPos].noGap = newScore;
       scoreVector[matchSeqPos].gapExists = continueGapScore;
       if (newScore > bestScore) {
         bestScore = newScore;
         bestQueryPos = queryPos;
         bestMatchSeqPos = matchSeqPos;
       }
     }
   }
   sfree(scoreVector);
   if (bestScore < 0)
     bestScore = 0;
   *matchSeqEnd = bestMatchSeqPos;
   *queryEnd = bestQueryPos;
   *score = bestScore;
   returnEvalue = BLAST_KarlinStoE_simple(bestScore,kbp, effSearchSpace);
   return(returnEvalue);
}

/** computes where optimal Smith-Waterman local alignment starts given the
 *  ending positions and score
 *  matchSeqEnd and queryEnd can be used to run the local alignment in reverse
 *  to find optimal starting positions
 *  these are passed back in matchSeqStart and queryStart
 *  the optimal score is passed in to check when it has
 *  been reached going backwards
 * the score is also returned
 * @param matchSeq is a database sequence matched by this query [in]
 * @param matchSeqLength is the length of matchSeq in amino acids [in]
 * @param query is the input query sequence  [in]
 * @param matrix is the position-specific matrix associated with query
 *      or the standard matrix [in]
 * @param gapOpen is the cost of opening a gap [in]
 * @param gapExtend is the cost of extending an existing gap by 1 position [in]
 * @param matchSeqEnd is the final position in the matchSeq of an optimal
 *      local alignment [in]
 * @param queryEnd is the final position in query of an optimal
 *  local alignment [in]
 * @param score optimal score to be obtained [in]
 * @param matchSeqStart starting point of optimal alignment [out]
 * @param queryStart starting point of optimal alignment [out]
 * @param positionSpecific determines whether matrix is position specific or not
*/
  
static Int4 BLSmithWatermanFindStart(Uint1 * matchSeq, 
   Int4 matchSeqLength, Uint1 *query, Int4 **matrix, 
   Int4 gapOpen, Int4 gapExtend,  Int4 matchSeqEnd, Int4 queryEnd, Int4 score,
   Int4 *matchSeqStart, Int4 *queryStart, Boolean positionSpecific)
{

   Int4 bestScore; /*best score seen so far*/
   Int4 newScore;  /* score of next entry*/
   Int4 bestMatchSeqPos, bestQueryPos; /*position starting best score in
                           matchSeq and database sequences*/
   SWpairs *scoreVector; /*keeps one row of the Smith-Waterman matrix
                           overwrite old row with new row*/
   Int4 *matrixRow; /*one row of score matrix*/
   Int4 newGapCost; /*cost to have a gap of one character*/
   Int4 prevScoreNoGapMatchSeq; /*score one row and column up
                               with no gaps*/
   Int4 prevScoreGapMatchSeq;   /*score if a gap already started in matchSeq*/
   Int4 continueGapScore; /*score for continuing a gap in query*/
   Int4 matchSeqPos, queryPos; /*positions in matchSeq and query*/

   scoreVector = (SWpairs *) calloc(1, matchSeqLength * sizeof(SWpairs));
   bestMatchSeqPos = 0;
   bestQueryPos = 0;
   bestScore = 0;
   newGapCost = gapOpen + gapExtend;
   for (matchSeqPos = 0; matchSeqPos < matchSeqLength; matchSeqPos++) {
     scoreVector[matchSeqPos].noGap = 0;
     scoreVector[matchSeqPos].gapExists = -(gapOpen);
   }
   for(queryPos = queryEnd; queryPos >= 0; queryPos--) {  
     if (positionSpecific)
       matrixRow = matrix[queryPos];
     else
       matrixRow = matrix[query[queryPos]];
     newScore = 0;
     prevScoreNoGapMatchSeq = 0;
     prevScoreGapMatchSeq = -(gapOpen);
     for(matchSeqPos = matchSeqEnd; matchSeqPos >= 0; matchSeqPos--) {
       /*testing scores with a gap in matchSeq, either starting a new
         gap or extending an existing gap*/
       if ((newScore = newScore - newGapCost) > 
	   (prevScoreGapMatchSeq = prevScoreGapMatchSeq - gapExtend))
         prevScoreGapMatchSeq = newScore;
       /*testing scores with a gap in query, either starting a new
         gap or extending an existing gap*/
       if ((newScore = scoreVector[matchSeqPos].noGap - newGapCost) >
           (continueGapScore = scoreVector[matchSeqPos].gapExists - gapExtend))
         continueGapScore = newScore;
       /*compute new score extending one position in matchSeq and query*/
       newScore = prevScoreNoGapMatchSeq + matrixRow[matchSeq[matchSeqPos]];
       if (newScore < 0)
       newScore = 0; /*Smith-Waterman locality condition*/
       /*test two alternatives*/
       if (newScore < prevScoreGapMatchSeq)
         newScore = prevScoreGapMatchSeq;
       if (newScore < continueGapScore)
         newScore = continueGapScore;
       prevScoreNoGapMatchSeq = scoreVector[matchSeqPos].noGap; 
       scoreVector[matchSeqPos].noGap = newScore;
       scoreVector[matchSeqPos].gapExists = continueGapScore;
       if (newScore > bestScore) {
         bestScore = newScore;
         bestQueryPos = queryPos;
         bestMatchSeqPos = matchSeqPos;
       }
       if (bestScore >= score)
         break;
     }
     if (bestScore >= score)
       break;
   }
   sfree(scoreVector);
   if (bestScore < 0)
     bestScore = 0;
   *matchSeqStart = bestMatchSeqPos;
   *queryStart = bestQueryPos;
   return(bestScore);
}