← Back to branch summary

~ubuntu-branches/ubuntu/edgy/ncbi-tools6/edgy

« back to all changes in this revision

Viewing changes to tools/blastutl.c

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-19 23:28:07 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20060719232807-et3cdmcjgmnyleyx
Tags: 6.1.20060507-3ubuntu1
Re-merge with Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/blastutl.c
1
 
static char const rcsid[] = "$Id: blastutl.c,v 6.464 2005/12/01 15:10:23 madden Exp $";
 
1
static char const rcsid[] = "$Id: blastutl.c,v 6.465 2006/02/15 18:23:47 madden Exp $";
2
2
 
3
3
/* ===========================================================================
4
4
*
32
32
 
33
33
Contents: Utilities for BLAST
34
34
 
35
 
$Revision: 6.464 $
 
35
$Revision: 6.465 $
36
36
 
37
37
******************************************************************************/
38
38
/*
39
39
 *
40
40
* $Log: blastutl.c,v $
 
41
* Revision 6.465  2006/02/15 18:23:47  madden
 
42
* Made changes so that CheckStartForGappedAlignment by default
 
43
* checks ungapped alignments of length 11, rather than length 10.
 
44
* Made changes to the rules used when the starting point is close to
 
45
* the edge of the preliminary gapped alignment.
 
46
* (from Mike Gertz)
 
47
*
41
48
* Revision 6.464  2005/12/01 15:10:23  madden
42
49
* Gave BLASTCheckHSPInclusion external linkage (i.e. removed the static specifier).
43
50
*
7984
7991
}
7985
7992
 
7986
7993
/*
7987
 
        Function to check that the highest scoring region in an HSP still gives a positive
7988
 
        score.  This value was originally calcualted by GetStartForGappedAlignment but it
7989
 
        may have changed due to the introduction of ambiguity characters.  Such a change
7990
 
        can lead to 'strange' results from ALIGN. 
 
7994
   Check whether the starting point for gapped alignment lies in
 
7995
   region that has positive score.  This routine is called after a
 
7996
   preliminary gapped alignment has been computed, but before the
 
7997
   traceback is computed.  The score of the region containing the
 
7998
   starting point may have changed due to the introduction of
 
7999
   ambiguity characters, further filtering of the sequences or the
 
8000
   application of composition based statistics.
 
8001
 
 
8002
   Usually, we check an ungapped alignment of length 11 about the
 
8003
   starting point: 5 characters to the left and 5 to the right.
 
8004
   However, the actual region checked is occassionally shorter because
 
8005
   we don't check characters before the start, or after the end, of
 
8006
   the preliminarily aligned regions in the query or subject.
7991
8007
*/
7992
8008
Boolean
7993
 
CheckStartForGappedAlignment (BlastSearchBlkPtr search, BLAST_HSPPtr hsp, Uint1Ptr query, Uint1Ptr subject, Int4Ptr PNTR matrix)
 
8009
CheckStartForGappedAlignment (BlastSearchBlkPtr search, BLAST_HSPPtr hsp,
 
8010
                              Uint1Ptr query, Uint1Ptr subject,
 
8011
                              Int4Ptr PNTR matrix)
7994
8012
{
7995
 
        Int4 index1, score, start, end, width;
7996
 
        Uint1Ptr query_var, subject_var;
7997
 
        Boolean positionBased = 
7998
 
           (search->positionBased && search->sbp->posMatrix);
7999
 
 
8000
 
        width = MIN((hsp->query.gapped_start-hsp->query.offset), HSP_MAX_WINDOW/2);
8001
 
        start = hsp->query.gapped_start - width;
8002
 
        end = MIN(hsp->query.end, hsp->query.gapped_start + width);
8003
 
        /* Assures that the start of subject is above zero. */
8004
 
        if ((hsp->subject.gapped_start + start - hsp->query.gapped_start) < 0)
8005
 
                start -= hsp->subject.gapped_start + (start - hsp->query.gapped_start);
8006
 
        query_var = query + start;
8007
 
        subject_var = subject + hsp->subject.gapped_start + (start - hsp->query.gapped_start);
8008
 
        score=0;
8009
 
        for (index1=start; index1<end; index1++)
8010
 
        {
8011
 
          if (!positionBased)
8012
 
            score += matrix[*query_var][*subject_var];
8013
 
          else
8014
 
            score += search->sbp->posMatrix[index1][*subject_var];
8015
 
          query_var++; subject_var++;
8016
 
        }
8017
 
        if (score <= 0)
8018
 
                return FALSE;
8019
 
 
8020
 
        return TRUE;
 
8013
    Int4 left, right;       /* Number of aligned characters to the
 
8014
                               left and right of the starting point */
 
8015
    Int4 score;             /* Score of the word alignment */
 
8016
    Uint1Ptr subject_var;   /* Current character in the subject sequence */
 
8017
    Uint1Ptr subject_right; /* last character to be considered in the subject
 
8018
                               sequence */
 
8019
    Boolean positionBased =
 
8020
        (search->positionBased && search->sbp->posMatrix);
 
8021
 
 
8022
    /* Compute the number of characters to the left of the start
 
8023
       to include in the word */
 
8024
    left = -HSP_MAX_WINDOW/2;
 
8025
    if (left < hsp->query.offset - hsp->query.gapped_start) {
 
8026
        left = hsp->query.offset - hsp->query.gapped_start;
 
8027
    }
 
8028
    if (left < hsp->subject.offset - hsp->subject.gapped_start) {
 
8029
        left = hsp->subject.offset - hsp->subject.gapped_start;
 
8030
    }
 
8031
 
 
8032
    /* Compute the number of characters to right to include in the word,
 
8033
       including the starting point itself. */
 
8034
    right = HSP_MAX_WINDOW/2 + 1;
 
8035
    if (right > hsp->query.end - hsp->query.gapped_start) {
 
8036
        right = hsp->query.end - hsp->query.gapped_start;
 
8037
    }
 
8038
    if (right > hsp->subject.end - hsp->subject.gapped_start) {
 
8039
        right = hsp->subject.end - hsp->subject.gapped_start;
 
8040
    }
 
8041
 
 
8042
    /* Calculate the score of the word */
 
8043
    score = 0;
 
8044
    subject_var   = subject + hsp->subject.gapped_start + left;
 
8045
    subject_right = subject + hsp->subject.gapped_start + right;
 
8046
    if ( !positionBased ) {
 
8047
        Uint1Ptr query_var;     /* Current character in the query */
 
8048
        query_var = query + hsp->query.gapped_start + left;
 
8049
        for ( ; subject_var < subject_right; subject_var++, query_var++) {
 
8050
           score += matrix[*query_var][*subject_var];
 
8051
        }
 
8052
    } else {
 
8053
        Int4 query_index;       /* Current position in the query */
 
8054
        query_index = hsp->query.gapped_start + left;
 
8055
        for ( ;  subject_var < subject_right;  subject_var++, query_index++) {
 
8056
            score += search->sbp->posMatrix[query_index][*subject_var];
 
8057
        }
 
8058
    }
 
8059
    if (score <= 0) {
 
8060
        return FALSE;
 
8061
    } else {
 
8062
        return TRUE;
 
8063
    }
8021
8064
}
 
8065
 
 
8066
 
8022
8067
/*
8023
8068
        Gets the ratio used to change an evalue calculated with the subject
8024
8069
        sequence length to one with a db length.