blast_kappa.c

ncbi源码

/** computes Smith-Waterman local alignment score and returns the
 *  evalue assuming some positions are forbidden
 *  matchSeqEnd and query can be used to run the local alignment in reverse
 *  to find optimal starting positions
 * @param matchSeq is the matchSeq sequence [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 either 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 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]
 * @param score is used to pass back the optimal score [out]
 * @param kbp holds the Karlin-Altschul parameters  [in]
 * @param effSearchSpace effective search space [in]
 * @param numForbidden number of forbidden ranges [in]
 * @param forbiddenRanges lists areas that should not be aligned [in]
 * @param positionSpecific determines whether matrix is position specific or not [in]
*/


static double BLspecialSmithWatermanScoreOnly(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,
   Int4 *numForbidden, Int4 ** forbiddenRanges, 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 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*/
   double returnEvalue; /*e-value to return*/
   Boolean forbidden; /*is this position forbidden?*/
   Int4 f; /*index over forbidden positions*/


   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*/
       forbidden = FALSE;
       for(f = 0; f < numForbidden[queryPos]; f++) {
         if ((matchSeqPos >= forbiddenRanges[queryPos][2 * f]) &&
	     (matchSeqPos <= forbiddenRanges[queryPos][2*f + 1])) {
	   forbidden = TRUE;
	   break;
	 }
       }
       if (forbidden)
         newScore = BLAST_SCORE_MIN;
       else
	 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.   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 the matchSeq sequence [in]
 * @param matchSeqLength is the length of matchSeq in amino acids [in]
 * @param query is the sequence corresponding to some matrix profile [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 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 is passed in to check when it has
 *  been reached going backwards [in]
 * @param matchSeqStart optimal starting point [in]
 * @param queryStart optimal starting point [in]
 * @param numForbidden array of regions not to be aligned. [in]
 * @param numForbidden array of regions not to be aligned. [in]
 * @param forbiddenRanges regions not to be aligned. [in]
 * @param positionSpecific determines whether matrix is position specific or not
 * @return the score found
*/

static Int4 BLspecialSmithWatermanFindStart(Uint1 * matchSeq, 
   Int4 matchSeqLength, Uint1 *query, Int4 **matrix, 
   Int4 gapOpen, Int4 gapExtend,  Int4 matchSeqEnd, Int4 queryEnd, Int4 score,
   Int4 *matchSeqStart, Int4 *queryStart, Int4 *numForbidden, 
   Int4 ** forbiddenRanges, 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*/
   Boolean forbidden; /*is this position forbidden?*/
   Int4 f; /*index over forbidden positions*/

   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*/
       forbidden = FALSE;
       for(f = 0; f < numForbidden[queryPos]; f++) {
         if ((matchSeqPos >= forbiddenRanges[queryPos][2 * f]) &&
	     (matchSeqPos <= forbiddenRanges[queryPos][2*f + 1])) {
	   forbidden = TRUE;
	   break;
	 }
       }
       if (forbidden)
         newScore = BLAST_SCORE_MIN;
       else
	 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);
}


/** converts the list of Smith-Waterman alignments to a corresponding list
 * of HSP's. kbp stores parameters for computing the score
 * Code is adapted from procedure output_hits of pseed3.c
 * @param SWAligns List of Smith-Waterman alignments [in]
 * @param hitlist BlastHitList that is filled in [in|out]
 */
static Int2 newConvertSWalignsUpdateHitList(SWResults * SWAligns, BlastHitList* hitList)
{
    BlastHSPList* hspList=NULL;
    SWResults* curSW;
    
    if (SWAligns == NULL)
       return 0;

    curSW = SWAligns;
    while (curSW != NULL) {
        if (hspList == NULL)
        {
             hspList = Blast_HSPListNew(0); 
             hspList->oid = curSW->subject_index;
        }

        Blast_HSPListSaveHSP(hspList, curSW->hsp);
        curSW->hsp = NULL; /* Saved on the hitlist, will be deleted there. */

        /* Changing OID being worked on. */
        if (curSW->next == NULL ||
              curSW->subject_index != curSW->next->subject_index)
        {
             Blast_HitListUpdate(hitList, hspList);
             hspList = NULL;
        }

        curSW = curSW->next;
    }

    return 0;
}


/** allocates  a score matrix with numPositions positions and initializes some 
 * positions on the side
 * @param numPositions length of matrix (or query) [in]
 * @return matrix (Int4**)  
 */
static Int4 **allocateScaledMatrix(Int4 numPositions)
{
  Int4 **returnMatrix; /*allocated matrix to return*/
  Int4 c; /*loop index over characters*/

  returnMatrix = (Int4**) _PSIAllocateMatrix(numPositions+1, BLASTAA_SIZE, sizeof(Int4));
  for(c = 0; c < BLASTAA_SIZE; c++)
    returnMatrix[numPositions][c] = BLAST_SCORE_MIN;
  return(returnMatrix);
}

/** allocate a frequency ratio matrix with numPositions positions and initialize
 *  some positions.
 * @param numPositions the length of matrix or query [in]
 * @return frequency matrix (double**)
 */
static double **allocateStartFreqs(Int4 numPositions)
{
  double **returnMatrix; /*allocated matrix to return*/
  Int4 c; /*loop index over characters*/

  returnMatrix = (double**) _PSIAllocateMatrix(numPositions+1, BLASTAA_SIZE, sizeof(double));
  for(c = 0; c < BLASTAA_SIZE; c++)
    returnMatrix[numPositions][c] = BLAST_SCORE_MIN;
  return(returnMatrix);
}

#if 0 
FIXME delte if not needed
/*deallocate a frequency ratio matrix*/
static void freeStartFreqs(double **matrix, Int4 numPositions)
{
  int row; /*loop index*/

  for(row = 0; row <= numPositions; row++)
    sfree(matrix[row]);
  sfree(matrix);
}
#endif

/*matrix is a position-specific score matrix with matrixLength positions
  queryProbArray is an array containing the probability of occurrence
  of each residue in the query
  scoreArray is an array of probabilities for each score that is
    to be used as a field in return_sfp
  return_sfp is a the structure to be filled in and returned
  range is the size of scoreArray and is an upper bound on the
   difference between maximum score and minimum score in the matrix
  the routine posfillSfp computes the probability of each score weighted
   by the probability of each query residue and fills those probabilities
   into scoreArray and puts scoreArray as a field in