~ubuntu-branches/ubuntu/precise/ncbi-tools6/precise

static char const rcsid[] = "$Id: blastkar.c,v 6.102 2004/09/28 16:04:19 papadopo Exp $";

 

 * $Revision: 6.102 $

 * Revision 6.102  2004/09/28 16:04:19  papadopo

 * From Michael Gertz:

 * 1. Pass the effective database size into BlastSmallGapSumE,

 *     BlastLargeGapSumE and BlastUnevenGapSumE.  The routines use this

 *     value in a simplified formula to compute the e-value of singleton sets.

 * 2. Caused all routines for calculating the significance of multiple

 *     distinct alignments (BlastSmallGapSumE, BlastLargeGapSumE and

 *     BlastUnevenGapSumE) to use

 *

 *        sum_{i in linked_set} (\lambda_i s_i - \ln K_i)

 *

 *     as the weighted sum score, where (\lambda_i, K_i) are taken from

 *     the appropriate query context.

 *

 * Revision 6.101  2004/06/07 20:03:23  coulouri

 * use floating point constants for comparisons with floating point variables

 *

 * Revision 6.100  2004/04/28 14:36:00  madden

 * Changes from Mike Gertz:

 * - I created the new routine BlastGapDecayDivisor that computes a

 *   divisor used to weight the evalue of a collection of distinct

 *   alignments.

 * - I removed  BlastGapDecay and BlastGapDecayInverse which had become

 *   redundant.

 * - I modified the BlastCutoffs routine so that it uses the value

 *   returned by BlastGapDecayDivisor to weight evalues.

 * - I modified BlastSmallGapSumE, BlastLargeGapSumE and

 *   BlastUnevenGapSumE no longer refer to the gap_prob parameter.

 *   Replaced the gap_decay_rate parameter of each of these routines with

 *   a weight_divisor parameter.  Added documentation.

 *

 * Revision 6.99  2004/04/23 13:49:43  madden

 * Cleaned up ifndef in BlastKarlinLHtoK

 *

 * Revision 6.98  2004/04/23 13:19:53  madden

 * Rewrote BlastKarlinLHtoK to do the following and more:

 * 1. fix a bug whereby the wrong formula was used when high score == 1

 *    and low score == -1;

 * 2. fix a methodological error of truncating the first sum

 *    and trying to make it converge quickly by adding terms

 *    of a geometric progression, even though the geometric progression

 *    estimate is not correct in all cases;

 *    the old adjustment code is left in for historical purposes but

 *    #ifdef'd out

 * 3. Eliminate the Boolean bi_modal_score variable.  The old test that

 *    set the value of bi_modal_score would frequently fail to choose the

 *    correct value due to rounding error.

 * 4. changed numerous local variable names to make them more meaningful;

 * 5. added substantial comments to explain what the procedure

 *    is doing and what each variable represents

 *

 * Revision 6.97  2004/04/01 13:43:08  lavr

 * Spell "occurred", "occurrence", and "occurring"

 *

 * Revision 6.96  2004/03/31 17:58:51  papadopo

 * Mike Gertz' changes for length adjustment calculations

 *

 * Revision 6.95  2003/12/12 16:00:34  madden

 * Add gap_decay_rate to BlastCutoffs, remove BlastCutoffs_simple, protection against overflow, removal of defunct _real variables (all from Mike Gertz)

 *

 * Revision 6.94  2003/11/30 03:36:38  camacho

 * Fix compilation error

 *

 * Revision 6.93  2003/11/28 22:39:40  camacho

 * +static keyword to BlastKarlinLtoH

 *

 * Revision 6.92  2003/11/28 15:16:38  camacho

 * Combine newkar.c's contents with blastkar.c

 *

 * Revision 6.91  2003/11/26 19:08:10  madden

 * simplified BlastKarlinLtoH and BlastKarlinLHtoK and provided better protection against overflow, new function NlmKarlinLambdaNR (all from Mike Gertz)

 *

 *

 * Revision 6.90  2003/06/30 20:01:32  dondosha

 * Correction in logic of finding matrices by BLASTMAT environment variable

 *

 * Revision 6.89  2003/05/30 17:25:36  coulouri

 * add rcsid

 *

 * Revision 6.88  2003/05/06 18:54:02  dondosha

 * Set all gap/residue scores in blastn matrix to INT4_MIN/2

 *

 * Revision 6.87  2003/03/07 22:33:25  bealer

 * - Fix UMR when alphabet_type is not set.

 *

 * Revision 6.86  2003/02/27 19:07:56  madden

 * Add functions PrintMatrixMessage and PrintAllowedValuesMessage

 *

 * Revision 6.85  2003/02/26 18:23:49  madden

 * Add functions BlastKarlinkGapBlkFill and BlastKarlinReportAllowedValues, call from BlastKarlinBlkGappedCalcEx

 *

 * Revision 6.84  2002/10/24 22:52:14  dondosha

 * When checking config file for matrices path, allow aa or nt subdirectories too

 *

 * Revision 6.83  2002/07/22 20:10:12  dondosha

 * Correction: previous change did not work for proteins

 *

 * Revision 6.82  2002/07/19 18:34:58  dondosha

 * Ignore bits higher than 4 when computing frequencies - needed for megablast

 *

 * Revision 6.81  2002/05/17 20:30:37  madden

 * Add comments on adding new matrix values

 *

 * Revision 6.80  2002/04/09 18:44:19  madden

 * Do not return if status of BlastScoreBlkMatCreate is non-zero

 *

 * Revision 6.79  2002/02/26 22:25:21  dondosha

 * Return error as soon as it is found that matrix name is not supported

 *

How the statistical parameters for the matrices are stored:

-----------------------------------------------------------

They parameters are stored in a two-dimensional array FloatHi (i.e., 

doubles, which has as it's first dimensions the number of different 

gap existence and extension combinations and as it's second dimension 8.

The eight different columns specify:

 

1.) gap existence penalty (INT2_MAX denotes infinite).

2.) gap extension penalty (INT2_MAX denotes infinite).

3.) decline to align penalty (INT2_MAX denotes infinite).

4.) Lambda

5.) K

6.) H

7.) alpha

8.) beta

 

(Items 4-8 are explained in:

Altschul SF, Bundschuh R, Olsen R, Hwa T.

The estimation of statistical parameters for local alignment score distributions.

Nucleic Acids Res. 2001 Jan 15;29(2):351-61.).

 

Take BLOSUM45 (below) as an example.  Currently (5/17/02) there are

14 different allowed combinations (specified by "#define BLOSUM45_VALUES_MAX 14").

The first row in the array "blosum45_values" has INT2_MAX (i.e., 32767) for gap 

existence, extension, and decline-to-align penalties.  For all practical purposes

this value is large enough to be infinite, so the alignments will be ungapped.

BLAST may also use this value (INT2_MAX) as a signal to skip gapping, so a

different value should not be used if the intent is to have gapless extensions.

The next row is for the gap existence penalty 13 and the extension penalty 3.

The decline-to-align penalty is only supported in a few cases, so it is normally

set to INT2_MAX.

 

 

How to add a new matrix to blastkar.c:

--------------------------------------

To add a new matrix to blastkar.c it is necessary to complete 

four steps.  As an example consider adding the matrix

called TESTMATRIX

 

1.) add a define specifying how many different existence and extensions

penalties are allowed, so it would be necessary to add the line:

 

#define TESTMATRIX_VALUES_MAX 14

 

if 14 values were to be allowed.

 

2.) add a two-dimensional array to contain the statistical parameters:

 

static Nlm_FloatHi testmatrix_values[TESTMATRIX_VALUES_MAX][8] ={ ...

 

3.) add a "prefs" array that should hint about the "optimal" 

gap existence and extension penalties:

 

static Int4 testmatrix_prefs[TESTMATRIX_VALUES_MAX] = {

BLAST_MATRIX_NOMINAL,

...

};

 

4.) Go to the function BlastLoadMatrixValues (in this file) and

add two lines before the return at the end of the function: 

 

        matrix_info = MatrixInfoNew("TESTMATRIX", testmatrix_values, testmatrix_prefs, TESTMATRIX_VALUES_MAX);

        ValNodeAddPointer(&retval, 0, matrix_info);

 

 

 

 

        for (index1=0; index1<BLASTNA_SIZE; index1++) /* blastna */

           matrix[index1][BLASTNA_SIZE-1] = INT4_MIN / 2;

    Char string[PATH_MAX] = "", alphabet_type[3] = "";

    CharPtr matrix_dir = NULL;

    Int2 status = 0;

    FILE *fp = NULL;

    

        /* 1. Try local directory */

        /* 2. Try configuration file */

           if (sbp->protein_alphabet)

              Nlm_StringNCpy(alphabet_type, "aa", 2);

           else

              Nlm_StringNCpy(alphabet_type, "nt", 2);

           alphabet_type[2] = NULLB;

 

           if(FindPath("ncbi", "ncbi", "data", string, PATH_MAX)) {

               matrix_dir = StringSave(string);

               sprintf(string, "%s%s", matrix_dir, matrix);

               if(FileLength(string) > 0) {

                  fp = FileOpen(string, "r");

               } else {

                  sprintf(string, "%s%s%s%s", matrix_dir, 

                          alphabet_type, DIRDELIMSTR, matrix);

                  if(FileLength(string) > 0)

                     fp = FileOpen(string, "r");

               }

               matrix_dir = MemFree(matrix_dir);

            sprintf(string, "data%s%s", DIRDELIMSTR, matrix);

        if(status != 0)

        {

                NlmMutexUnlock(read_matrix_mutex); 

                return status;

        }

        Uint1 mask;

        /* For megablast, check only the first 4 bits of the sequence values */

        mask = (sbp->protein_alphabet ? 0xff : 0x0f);

 

                ++rcp->comp[(int)(*lp & mask)];

   Int2 num_values;

   else if (num_values > 0) {

        Char buffer[256];

        Int2 status = BlastKarlinkGapBlkFill(kbp, gap_open, gap_extend, decline_align, matrix_name);

 

        if (status && error_return)

        {

                if (status == 1)

                {

                        MatrixInfoPtr matrix_info;

                        ValNodePtr vnp, head;           

 

                        vnp = head = BlastLoadMatrixValues();

 

                        sprintf(buffer, "%s is not a supported matrix", matrix_name);

                        BlastConstructErrorMessage("BlastKarlinBlkGappedCalc", buffer, 2, error_return);

 

                        while (vnp)

                        {

                                matrix_info = vnp->data.ptrvalue;

                                sprintf(buffer, "%s is a supported matrix", matrix_info->name);

                                BlastConstructErrorMessage("BlastKarlinBlkGappedCalc", buffer, 2, error_return);

                                vnp = vnp->next;

                        }

 

                        BlastMatrixValuesDestruct(head);

                }

                else if (status == 2)

                {

                        if (decline_align == INT2_MAX)

                                sprintf(buffer, "Gap existence and extension values of %ld and %ld not supported for %s", (long) gap_open, (long) gap_extend, matrix_name);

                        else

                                sprintf(buffer, "Gap existence, extension and decline-to-align values of %ld, %ld and %ld not supported for %s", (long) gap_open, (long) gap_extend, (long) decline_align, matrix_name);

                        BlastConstructErrorMessage("BlastKarlinBlkGappedCalc", buffer, 2, error_return);

                        BlastKarlinReportAllowedValues(matrix_name, error_return);

                }

        }

 

        return status;

}

 

/*

        Attempts to fill KarlinBlk for given gap opening, extensions etc.  

        Will return non-zero status if that fails.

 

        return values:  -1 if matrix_name is NULL;

                        1 if matrix not found

                        2 if matrix found, but open, extend etc. values not supported.

*/

Int2 LIBCALL

BlastKarlinkGapBlkFill(BLAST_KarlinBlkPtr kbp, Int4 gap_open, Int4 gap_extend, Int4 decline_align, CharPtr matrix_name)

{

        Boolean found_matrix=FALSE, found_values=FALSE;

        Int2 status=0;

                return -1;

                                        kbp->Lambda = values[index][3];

                                        kbp->K = values[index][4];

                                        kbp->logK = log(kbp->K);

                                        kbp->H = values[index][5];

                        status = 0;

                }

                else

                {

                        status = 2;

                status = 1;

        }

 

        BlastMatrixValuesDestruct(head);

 

        return status;

}

 

CharPtr

PrintMatrixMessage(const Char *matrix_name)

{

        CharPtr buffer=MemNew(1024*sizeof(Char));

        CharPtr ptr;

        MatrixInfoPtr matrix_info;

        ValNodePtr vnp, head;

 

        ptr = buffer;

        sprintf(ptr, "%s is not a supported matrix, supported matrices are:\n", matrix_name);

 

        ptr += StringLen(ptr);

 

        vnp = head = BlastLoadMatrixValues();

 

        while (vnp)

        {

                matrix_info = vnp->data.ptrvalue;

                sprintf(ptr, "%s \n", matrix_info->name);

                ptr += StringLen(ptr);

                vnp = vnp->next;

        }

        BlastMatrixValuesDestruct(head);

 

        return buffer;

}

 

CharPtr

PrintAllowedValuesMessage(const Char *matrix_name, Int4 gap_open, Int4 gap_extend, Int4 decline_align) 

{

        array_of_8 *values;