开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
blast_stat.c
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：113k
源码类别：

生物技术

开发平台：

C/C++

blast_stat.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_stat.c,v $
  4.  * PRODUCTION Revision 1000.6  2004/06/01 18:07:50  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.77
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: blast_stat.c,v 1000.6 2004/06/01 18:07:50 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's official duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Author: Tom Madden
  36.  *
  37.  */
  39. /** @file blast_stat.c
  40.  * Functions to calculate BLAST probabilities etc.
  41.  * Detailed Contents:
  42.  *
  43.  * - allocate and deallocate structures used by BLAST to calculate
  44.  * probabilities etc.
  45.  *
  46.  * - calculate residue frequencies for query and "average" database.
  47.  *
  48.  * - read in matrix or load it from memory.
  49.  *
  50.  *  - calculate sum-p from a collection of HSP's, for both the case
  51.  *   of a "small" gap and a "large" gap, when give a total score and the
  52.  *   number of HSP's.
  53.  *
  54.  * - calculate expect values for p-values.
  55.  *
  56.  * - calculate pseuod-scores from p-values.
  57.  *
  58.  * @todo FIXME needs doxygen comments
  59.  */
  61. static char const rcsid[] = 
  62.     "$Id: blast_stat.c,v 1000.6 2004/06/01 18:07:50 gouriano Exp $";
  64. #include <algo/blast/core/blast_stat.h>
  65. #include <algo/blast/core/blast_util.h>
  66. #include <util/tables/raw_scoremat.h>
  67. #include <algo/blast/core/blast_encoding.h>
  68. #include "blast_psi_priv.h"
  70. /* OSF1 apparently doesn't like this. */
  71. #if defined(HUGE_VAL) && !defined(OS_UNIX_OSF1)
  72. #define BLASTKAR_HUGE_VAL HUGE_VAL
  73. #else
  74. #define BLASTKAR_HUGE_VAL 1.e30
  75. #endif
  78. /* Allocates and Deallocates the two-dimensional matrix. */
  79. static SBLASTMatrixStructure* BlastMatrixAllocate (Int2 alphabet_size);
  81. /* performs sump calculation, used by BlastSumPStd */
  82. static double BlastSumPCalc (int r, double s);
  84. #define COMMENT_CHR '#'
  85. #define TOKSTR " tnr"
  86. #define BLAST_MAX_ALPHABET 40 /* ncbistdaa is only 26, this should be enough */
  88. /*
  89. How many of the first bases are not ambiguous
  90. (it's four, of course).
  91. */
  92. #define NUMBER_NON_AMBIG_BP 4
  94. /* Used in BlastKarlinBlkGappedCalc */
  95. typedef double array_of_8[8];
  97. /* Used to temporarily store matrix values for retrieval. */
  99. typedef struct MatrixInfo {
  100. char* name; /* name of matrix (e.g., BLOSUM90). */
  101. array_of_8  *values; /* The values (below). */
  102. Int4 *prefs; /* Preferences for display. */
  103. Int4 max_number_values; /* number of values (e.g., BLOSUM90_VALUES_MAX). */
  104. } MatrixInfo;
  107. /**************************************************************************************
  109. How the statistical parameters for the matrices are stored:
  110. -----------------------------------------------------------
  111. They parameters are stored in a two-dimensional array double (i.e., 
  112. doubles, which has as it's first dimensions the number of different 
  113. gap existence and extension combinations and as it's second dimension 8.
  114. The eight different columns specify:
  116. 1.) gap existence penalty (INT2_MAX denotes infinite).
  117. 2.) gap extension penalty (INT2_MAX denotes infinite).
  118. 3.) decline to align penalty (INT2_MAX denotes infinite).
  119. 4.) Lambda
  120. 5.) K
  121. 6.) H
  122. 7.) alpha
  123. 8.) beta
  125. (Items 4-8 are explained in:
  126. Altschul SF, Bundschuh R, Olsen R, Hwa T.
  127. The estimation of statistical parameters for local alignment score distributions.
  128. Nucleic Acids Res. 2001 Jan 15;29(2):351-61.).
  130. Take BLOSUM45 (below) as an example.  Currently (5/17/02) there are
  131. 14 different allowed combinations (specified by "#define BLOSUM45_VALUES_MAX 14").
  132. The first row in the array "blosum45_values" has INT2_MAX (i.e., 32767) for gap 
  133. existence, extension, and decline-to-align penalties.  For all practical purposes
  134. this value is large enough to be infinite, so the alignments will be ungapped.
  135. BLAST may also use this value (INT2_MAX) as a signal to skip gapping, so a
  136. different value should not be used if the intent is to have gapless extensions.
  137. The next row is for the gap existence penalty 13 and the extension penalty 3.
  138. The decline-to-align penalty is only supported in a few cases, so it is normally
  139. set to INT2_MAX.
  142. How to add a new matrix to blast_stat.c:
  143. --------------------------------------
  144. To add a new matrix to blast_stat.c it is necessary to complete 
  145. four steps.  As an example consider adding the matrix
  146. called TESTMATRIX
  148. 1.) add a define specifying how many different existence and extensions
  149. penalties are allowed, so it would be necessary to add the line:
  151. #define TESTMATRIX_VALUES_MAX 14
  153. if 14 values were to be allowed.
  155. 2.) add a two-dimensional array to contain the statistical parameters:
  157. static double testmatrix_values[TESTMATRIX_VALUES_MAX][8] ={ ...
  159. 3.) add a "prefs" array that should hint about the "optimal" 
  160. gap existence and extension penalties:
  162. static Int4 testmatrix_prefs[TESTMATRIX_VALUES_MAX] = {
  163. BLAST_MATRIX_NOMINAL,
  164. ...
  165. };
  167. 4.) Go to the function BlastLoadMatrixValues (in this file) and
  168. add two lines before the return at the end of the function: 
  170.         matrix_info = MatrixInfoNew("TESTMATRIX", testmatrix_values, testmatrix_prefs, TESTMATRIX_VALUES_MAX);
  171.         ListNodeAddPointer(&retval, 0, matrix_info);
  175. ***************************************************************************************/
  181. #define BLOSUM45_VALUES_MAX 14
  182. static double  blosum45_values[BLOSUM45_VALUES_MAX][8] = {
  183.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.2291, 0.0924, 0.2514, 0.9113, -5.7},
  184.     {13, 3, (double) INT2_MAX, 0.207, 0.049, 0.14, 1.5, -22},
  185.     {12, 3, (double) INT2_MAX, 0.199, 0.039, 0.11, 1.8, -34},
  186.     {11, 3, (double) INT2_MAX, 0.190, 0.031, 0.095, 2.0, -38},
  187.     {10, 3, (double) INT2_MAX, 0.179, 0.023, 0.075, 2.4, -51},
  188.     {16, 2, (double) INT2_MAX, 0.210, 0.051, 0.14, 1.5, -24},
  189.     {15, 2, (double) INT2_MAX, 0.203, 0.041, 0.12, 1.7, -31},
  190.     {14, 2, (double) INT2_MAX, 0.195, 0.032, 0.10, 1.9, -36},
  191.     {13, 2, (double) INT2_MAX, 0.185, 0.024, 0.084, 2.2, -45},
  192.     {12, 2, (double) INT2_MAX, 0.171, 0.016, 0.061, 2.8, -65},
  193.     {19, 1, (double) INT2_MAX, 0.205, 0.040, 0.11, 1.9, -43},
  194.     {18, 1, (double) INT2_MAX, 0.198, 0.032, 0.10, 2.0, -43},
  195.     {17, 1, (double) INT2_MAX, 0.189, 0.024, 0.079, 2.4, -57},
  196.     {16, 1, (double) INT2_MAX, 0.176, 0.016, 0.063, 2.8, -67},
  197. };
  199. static Int4 blosum45_prefs[BLOSUM45_VALUES_MAX] = {
  200. BLAST_MATRIX_NOMINAL,
  201. BLAST_MATRIX_NOMINAL,
  202. BLAST_MATRIX_NOMINAL,
  203. BLAST_MATRIX_NOMINAL,
  204. BLAST_MATRIX_NOMINAL,
  205. BLAST_MATRIX_NOMINAL,
  206. BLAST_MATRIX_NOMINAL,
  207. BLAST_MATRIX_BEST,
  208. BLAST_MATRIX_NOMINAL,
  209. BLAST_MATRIX_NOMINAL,
  210. BLAST_MATRIX_NOMINAL,
  211. BLAST_MATRIX_NOMINAL,
  212. BLAST_MATRIX_NOMINAL,
  213. BLAST_MATRIX_NOMINAL
  214. };
  217. #define BLOSUM50_VALUES_MAX 16
  218. static double  blosum50_values[BLOSUM50_VALUES_MAX][8] = {
  219.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.2318, 0.112, 0.3362, 0.6895, -4.0},
  220.     {13, 3, (double) INT2_MAX, 0.212, 0.063, 0.19, 1.1, -16},
  221.     {12, 3, (double) INT2_MAX, 0.206, 0.055, 0.17, 1.2, -18},
  222.     {11, 3, (double) INT2_MAX, 0.197, 0.042, 0.14, 1.4, -25},
  223.     {10, 3, (double) INT2_MAX, 0.186, 0.031, 0.11, 1.7, -34},
  224.     {9, 3, (double) INT2_MAX, 0.172, 0.022, 0.082, 2.1, -48},
  225.     {16, 2, (double) INT2_MAX, 0.215, 0.066, 0.20, 1.05, -15},
  226.     {15, 2, (double) INT2_MAX, 0.210, 0.058, 0.17, 1.2, -20},
  227.     {14, 2, (double) INT2_MAX, 0.202, 0.045, 0.14, 1.4, -27},
  228.     {13, 2, (double) INT2_MAX, 0.193, 0.035, 0.12, 1.6, -32},
  229.     {12, 2, (double) INT2_MAX, 0.181, 0.025, 0.095, 1.9, -41},
  230.     {19, 1, (double) INT2_MAX, 0.212, 0.057, 0.18, 1.2, -21},
  231.     {18, 1, (double) INT2_MAX, 0.207, 0.050, 0.15, 1.4, -28},
  232.     {17, 1, (double) INT2_MAX, 0.198, 0.037, 0.12, 1.6, -33},
  233.     {16, 1, (double) INT2_MAX, 0.186, 0.025, 0.10, 1.9, -42},
  234.     {15, 1, (double) INT2_MAX, 0.171, 0.015, 0.063, 2.7, -76},
  235. };
  237. static Int4 blosum50_prefs[BLOSUM50_VALUES_MAX] = {
  238. BLAST_MATRIX_NOMINAL,
  239. BLAST_MATRIX_NOMINAL,
  240. BLAST_MATRIX_NOMINAL,
  241. BLAST_MATRIX_NOMINAL,
  242. BLAST_MATRIX_NOMINAL,
  243. BLAST_MATRIX_NOMINAL,
  244. BLAST_MATRIX_NOMINAL,
  245. BLAST_MATRIX_NOMINAL,
  246. BLAST_MATRIX_NOMINAL,
  247. BLAST_MATRIX_BEST,
  248. BLAST_MATRIX_NOMINAL,
  249. BLAST_MATRIX_NOMINAL,
  250. BLAST_MATRIX_NOMINAL,
  251. BLAST_MATRIX_NOMINAL,
  252. BLAST_MATRIX_NOMINAL,
  253. BLAST_MATRIX_NOMINAL
  254. };
  256. #define BLOSUM62_VALUES_MAX 12
  257. static double  blosum62_values[BLOSUM62_VALUES_MAX][8] = {
  258.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.3176, 0.134, 0.4012, 0.7916, -3.2},
  259.     {11, 2, (double) INT2_MAX, 0.297, 0.082, 0.27, 1.1, -10},
  260.     {10, 2, (double) INT2_MAX, 0.291, 0.075, 0.23, 1.3, -15},
  261.     {9, 2, (double) INT2_MAX, 0.279, 0.058, 0.19, 1.5, -19},
  262.     {8, 2, (double) INT2_MAX, 0.264, 0.045, 0.15, 1.8, -26},
  263.     {7, 2, (double) INT2_MAX, 0.239, 0.027, 0.10, 2.5, -46},
  264.     {6, 2, (double) INT2_MAX, 0.201, 0.012, 0.061, 3.3, -58},
  265.     {13, 1, (double) INT2_MAX, 0.292, 0.071, 0.23, 1.2, -11},
  266.     {12, 1, (double) INT2_MAX, 0.283, 0.059, 0.19, 1.5, -19},
  267.     {11, 1, (double) INT2_MAX, 0.267, 0.041, 0.14, 1.9, -30},
  268.     {10, 1, (double) INT2_MAX, 0.243, 0.024, 0.10, 2.5, -44},
  269.     {9, 1, (double) INT2_MAX, 0.206, 0.010, 0.052, 4.0, -87},
  270. };
  272. static Int4 blosum62_prefs[BLOSUM62_VALUES_MAX] = {
  273.     BLAST_MATRIX_NOMINAL,
  274.     BLAST_MATRIX_NOMINAL,
  275.     BLAST_MATRIX_NOMINAL,
  276.     BLAST_MATRIX_NOMINAL,
  277.     BLAST_MATRIX_NOMINAL,
  278.     BLAST_MATRIX_NOMINAL,
  279.     BLAST_MATRIX_NOMINAL,
  280.     BLAST_MATRIX_NOMINAL,
  281.     BLAST_MATRIX_NOMINAL,
  282.     BLAST_MATRIX_BEST,
  283.     BLAST_MATRIX_NOMINAL,
  284.     BLAST_MATRIX_NOMINAL,
  285. };
  288. #define BLOSUM80_VALUES_MAX 10
  289. static double  blosum80_values[BLOSUM80_VALUES_MAX][8] = {
  290.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.3430, 0.177, 0.6568, 0.5222, -1.6},
  291.     {25, 2, (double) INT2_MAX, 0.342, 0.17, 0.66, 0.52, -1.6},
  292.     {13, 2, (double) INT2_MAX, 0.336, 0.15, 0.57, 0.59, -3},
  293.     {9, 2, (double) INT2_MAX, 0.319, 0.11, 0.42, 0.76, -6},
  294.     {8, 2, (double) INT2_MAX, 0.308, 0.090, 0.35, 0.89, -9},
  295.     {7, 2, (double) INT2_MAX, 0.293, 0.070, 0.27, 1.1, -14},
  296.     {6, 2, (double) INT2_MAX, 0.268, 0.045, 0.19, 1.4, -19},
  297.     {11, 1, (double) INT2_MAX, 0.314, 0.095, 0.35, 0.90, -9},
  298.     {10, 1, (double) INT2_MAX, 0.299, 0.071, 0.27, 1.1, -14},
  299.     {9, 1, (double) INT2_MAX, 0.279, 0.048, 0.20, 1.4, -19},
  300. };
  302. static Int4 blosum80_prefs[BLOSUM80_VALUES_MAX] = {
  303.     BLAST_MATRIX_NOMINAL,
  304.     BLAST_MATRIX_NOMINAL,
  305.     BLAST_MATRIX_NOMINAL,
  306.     BLAST_MATRIX_NOMINAL,
  307.     BLAST_MATRIX_NOMINAL,
  308.     BLAST_MATRIX_NOMINAL,
  309.     BLAST_MATRIX_NOMINAL,
  310.     BLAST_MATRIX_BEST,
  311.     BLAST_MATRIX_NOMINAL
  312. };
  314. #define BLOSUM90_VALUES_MAX 8
  315. static double  blosum90_values[BLOSUM90_VALUES_MAX][8] = {
  316.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.3346, 0.190, 0.7547, 0.4434, -1.4},
  317.     {9, 2, (double) INT2_MAX, 0.310, 0.12, 0.46, 0.67, -6},
  318.     {8, 2, (double) INT2_MAX, 0.300, 0.099, 0.39, 0.76, -7},
  319.     {7, 2, (double) INT2_MAX, 0.283, 0.072, 0.30, 0.93, -11},
  320.     {6, 2, (double) INT2_MAX, 0.259, 0.048, 0.22, 1.2, -16},
  321.     {11, 1, (double) INT2_MAX, 0.302, 0.093, 0.39, 0.78, -8},
  322.     {10, 1, (double) INT2_MAX, 0.290, 0.075, 0.28, 1.04, -15},
  323.     {9, 1, (double) INT2_MAX, 0.265, 0.044, 0.20, 1.3, -19},
  324. };
  326. static Int4 blosum90_prefs[BLOSUM90_VALUES_MAX] = {
  327. BLAST_MATRIX_NOMINAL,
  328. BLAST_MATRIX_NOMINAL,
  329. BLAST_MATRIX_NOMINAL,
  330. BLAST_MATRIX_NOMINAL,
  331. BLAST_MATRIX_NOMINAL,
  332. BLAST_MATRIX_NOMINAL,
  333. BLAST_MATRIX_BEST,
  334. BLAST_MATRIX_NOMINAL
  335. };
  337. #define PAM250_VALUES_MAX 16
  338. static double  pam250_values[PAM250_VALUES_MAX][8] = {
  339.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.2252, 0.0868, 0.2223, 0.98, -5.0},
  340.     {15, 3, (double) INT2_MAX, 0.205, 0.049, 0.13, 1.6, -23},
  341.     {14, 3, (double) INT2_MAX, 0.200, 0.043, 0.12, 1.7, -26},
  342.     {13, 3, (double) INT2_MAX, 0.194, 0.036, 0.10, 1.9, -31},
  343.     {12, 3, (double) INT2_MAX, 0.186, 0.029, 0.085, 2.2, -41},
  344.     {11, 3, (double) INT2_MAX, 0.174, 0.020, 0.070, 2.5, -48},
  345.     {17, 2, (double) INT2_MAX, 0.204, 0.047, 0.12, 1.7, -28},
  346.     {16, 2, (double) INT2_MAX, 0.198, 0.038, 0.11, 1.8, -29},
  347.     {15, 2, (double) INT2_MAX, 0.191, 0.031, 0.087, 2.2, -44},
  348.     {14, 2, (double) INT2_MAX, 0.182, 0.024, 0.073, 2.5, -53},
  349.     {13, 2, (double) INT2_MAX, 0.171, 0.017, 0.059, 2.9, -64},
  350.     {21, 1, (double) INT2_MAX, 0.205, 0.045, 0.11, 1.8, -34},
  351.     {20, 1, (double) INT2_MAX, 0.199, 0.037, 0.10, 1.9, -35},
  352.     {19, 1, (double) INT2_MAX, 0.192, 0.029, 0.083, 2.3, -52},
  353.     {18, 1, (double) INT2_MAX, 0.183, 0.021, 0.070, 2.6, -60},
  354.     {17, 1, (double) INT2_MAX, 0.171, 0.014, 0.052, 3.3, -86},
  355. };
  357. static Int4 pam250_prefs[PAM250_VALUES_MAX] = {
  358. BLAST_MATRIX_NOMINAL,
  359. BLAST_MATRIX_NOMINAL,
  360. BLAST_MATRIX_NOMINAL,
  361. BLAST_MATRIX_NOMINAL,
  362. BLAST_MATRIX_NOMINAL,
  363. BLAST_MATRIX_NOMINAL,
  364. BLAST_MATRIX_NOMINAL,
  365. BLAST_MATRIX_NOMINAL,
  366. BLAST_MATRIX_BEST,
  367. BLAST_MATRIX_NOMINAL,
  368. BLAST_MATRIX_NOMINAL,
  369. BLAST_MATRIX_NOMINAL,
  370. BLAST_MATRIX_NOMINAL,
  371. BLAST_MATRIX_NOMINAL,
  372. BLAST_MATRIX_NOMINAL,
  373. BLAST_MATRIX_NOMINAL
  374. };
  376. #define PAM30_VALUES_MAX 7
  377. static double  pam30_values[PAM30_VALUES_MAX][8] = {
  378.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.3400, 0.283, 1.754, 0.1938, -0.3},
  379.     {7, 2, (double) INT2_MAX, 0.305, 0.15, 0.87, 0.35, -3},
  380.     {6, 2, (double) INT2_MAX, 0.287, 0.11, 0.68, 0.42, -4},
  381.     {5, 2, (double) INT2_MAX, 0.264, 0.079, 0.45, 0.59, -7},
  382.     {10, 1, (double) INT2_MAX, 0.309, 0.15, 0.88, 0.35, -3},
  383.     {9, 1, (double) INT2_MAX, 0.294, 0.11, 0.61, 0.48, -6},
  384.     {8, 1, (double) INT2_MAX, 0.270, 0.072, 0.40, 0.68, -10},
  385. };
  387. static Int4 pam30_prefs[PAM30_VALUES_MAX] = {
  388. BLAST_MATRIX_NOMINAL,
  389. BLAST_MATRIX_NOMINAL,
  390. BLAST_MATRIX_NOMINAL,
  391. BLAST_MATRIX_NOMINAL,
  392. BLAST_MATRIX_NOMINAL,
  393. BLAST_MATRIX_BEST,
  394. BLAST_MATRIX_NOMINAL,
  395. };
  398. #define PAM70_VALUES_MAX 7
  399. static double  pam70_values[PAM70_VALUES_MAX][8] = {
  400.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.3345, 0.229, 1.029, 0.3250,   -0.7},
  401.     {8, 2, (double) INT2_MAX, 0.301, 0.12, 0.54, 0.56, -5},
  402.     {7, 2, (double) INT2_MAX, 0.286, 0.093, 0.43, 0.67, -7},
  403.     {6, 2, (double) INT2_MAX, 0.264, 0.064, 0.29, 0.90, -12},
  404.     {11, 1, (double) INT2_MAX, 0.305, 0.12, 0.52, 0.59, -6},
  405.     {10, 1, (double) INT2_MAX, 0.291, 0.091, 0.41, 0.71, -9},
  406.     {9, 1, (double) INT2_MAX, 0.270, 0.060, 0.28, 0.97, -14},
  407. };
  409. static Int4 pam70_prefs[PAM70_VALUES_MAX] = {
  410. BLAST_MATRIX_NOMINAL,
  411. BLAST_MATRIX_NOMINAL,
  412. BLAST_MATRIX_NOMINAL,
  413. BLAST_MATRIX_NOMINAL,
  414. BLAST_MATRIX_NOMINAL,
  415. BLAST_MATRIX_BEST,
  416. BLAST_MATRIX_NOMINAL
  417. };
  421. #define BLOSUM62_20_VALUES_MAX 65
  422. static double  blosum62_20_values[BLOSUM62_20_VALUES_MAX][8] = {
  423.     {(double) INT2_MAX, (double) INT2_MAX, (double) INT2_MAX, 0.03391, 0.125, 0.4544, 0.07462, -3.2},
  424.     {100, 12, (double) INT2_MAX, 0.0300, 0.056, 0.21, 0.14, -15},
  425.     {95, 12, (double) INT2_MAX, 0.0291, 0.047, 0.18, 0.16, -20},
  426.     {90, 12, (double) INT2_MAX, 0.0280, 0.038, 0.15, 0.19, -28},
  427.     {85, 12, (double) INT2_MAX, 0.0267, 0.030, 0.13, 0.21, -31},
  428.     {80, 12, (double) INT2_MAX, 0.0250, 0.021, 0.10, 0.25, -39},
  429.     {105, 11, (double) INT2_MAX, 0.0301, 0.056, 0.22, 0.14, -16},
  430.     {100, 11, (double) INT2_MAX, 0.0294, 0.049, 0.20, 0.15, -17},
  431.     {95, 11, (double) INT2_MAX, 0.0285, 0.042, 0.16, 0.18, -25},
  432.     {90, 11, (double) INT2_MAX, 0.0271, 0.031, 0.14, 0.20, -28},
  433.     {85, 11, (double) INT2_MAX, 0.0256, 0.023, 0.10, 0.26, -46},
  434.     {115, 10, (double) INT2_MAX, 0.0308, 0.062, 0.22, 0.14, -20},
  435.     {110, 10, (double) INT2_MAX, 0.0302, 0.056, 0.19, 0.16, -26},
  436.     {105, 10, (double) INT2_MAX, 0.0296, 0.050, 0.17, 0.17, -27},
  437.     {100, 10, (double) INT2_MAX, 0.0286, 0.041, 0.15, 0.19, -32},
  438.     {95, 10, (double) INT2_MAX, 0.0272, 0.030, 0.13, 0.21, -35},
  439.     {90, 10, (double) INT2_MAX, 0.0257, 0.022, 0.11, 0.24, -40},
  440.     {85, 10, (double) INT2_MAX, 0.0242, 0.017, 0.083, 0.29, -51},
  441.     {115, 9, (double) INT2_MAX, 0.0306, 0.061, 0.24, 0.13, -14},
  442.     {110, 9, (double) INT2_MAX, 0.0299, 0.053, 0.19, 0.16, -23},
  443.     {105, 9, (double) INT2_MAX, 0.0289, 0.043, 0.17, 0.17, -23},
  444.     {100, 9, (double) INT2_MAX, 0.0279, 0.036, 0.14, 0.20, -31},
  445.     {95, 9, (double) INT2_MAX, 0.0266, 0.028, 0.12, 0.23, -37},
  446.     {120, 8, (double) INT2_MAX, 0.0307, 0.062, 0.22, 0.14, -18},
  447.     {115, 8, (double) INT2_MAX, 0.0300, 0.053, 0.20, 0.15, -19},
  448.     {110, 8, (double) INT2_MAX, 0.0292, 0.046, 0.17, 0.17, -23},
  449.     {105, 8, (double) INT2_MAX, 0.0280, 0.035, 0.14, 0.20, -31},
  450.     {100, 8, (double) INT2_MAX, 0.0266, 0.026, 0.12, 0.23, -37},
  451.     {125, 7, (double) INT2_MAX, 0.0306, 0.058, 0.22, 0.14, -18},
  452.     {120, 7, (double) INT2_MAX, 0.0300, 0.052, 0.19, 0.16, -23},
  453.     {115, 7, (double) INT2_MAX, 0.0292, 0.044, 0.17, 0.17, -24},
  454.     {110, 7, (double) INT2_MAX, 0.0279, 0.032, 0.14, 0.20, -31},
  455.     {105, 7, (double) INT2_MAX, 0.0267, 0.026, 0.11, 0.24, -41},
  456.     {120,10,5, 0.0298, 0.049, 0.19, 0.16, -21},
  457.     {115,10,5, 0.0290, 0.042, 0.16, 0.18, -25},
  458.     {110,10,5, 0.0279, 0.033, 0.13, 0.21, -32},
  459.     {105,10,5, 0.0264, 0.024, 0.10, 0.26, -46},
  460.     {100,10,5, 0.0250, 0.018, 0.081, 0.31, -56},
  461.     {125,10,4, 0.0301, 0.053, 0.18, 0.17, -25},
  462.     {120,10,4, 0.0292, 0.043, 0.15, 0.20, -33},
  463.     {115,10,4, 0.0282, 0.035, 0.13, 0.22, -36},
  464.     {110,10,4, 0.0270, 0.027, 0.11, 0.25, -41},
  465.     {105,10,4, 0.0254, 0.020, 0.079, 0.32, -60},
  466.     {130,10,3, 0.0300, 0.051, 0.17, 0.18, -27},
  467.     {125,10,3, 0.0290, 0.040, 0.13, 0.22, -38},
  468.     {120,10,3, 0.0278, 0.030, 0.11, 0.25, -44},
  469.     {115,10,3, 0.0267, 0.025, 0.092, 0.29, -52},
  470.     {110,10,3, 0.0252, 0.018, 0.070, 0.36, -70},
  471.     {135,10,2, 0.0292, 0.040, 0.13, 0.22, -35},
  472.     {130,10,2, 0.0283, 0.034, 0.10, 0.28, -51},
  473.     {125,10,2, 0.0269, 0.024, 0.077, 0.35, -71},
  474.     {120,10,2, 0.0253, 0.017, 0.059, 0.43, -90},
  475.     {115,10,2, 0.0234, 0.011, 0.043, 0.55, -121},
  476.     {100,14,3, 0.0258, 0.023, 0.087, 0.33, -59},
  477.     {105,13,3, 0.0263, 0.024, 0.085, 0.31, -57},
  478.     {110,12,3, 0.0271, 0.028, 0.093, 0.29, -54},
  479.     {115,11,3, 0.0275, 0.030, 0.10, 0.27, -49},
  480.     {125,9,3, 0.0283, 0.034, 0.12, 0.23, -38},
  481.     {130,8,3, 0.0287, 0.037, 0.12, 0.23, -40},
  482.     {125,7,3, 0.0287, 0.036, 0.12, 0.24, -44},
  483.     {140,6,3, 0.0285, 0.033, 0.12, 0.23, -40},
  484.     {105,14,3, 0.0270, 0.028, 0.10, 0.27, -46},
  485.     {110,13,3, 0.0279, 0.034, 0.10, 0.27, -50},
  486.     {115,12,3, 0.0282, 0.035, 0.12, 0.24, -42},
  487.     {120,11,3, 0.0286, 0.037, 0.12, 0.24, -44},
  488. };
  490. static Int4 blosum62_20_prefs[BLOSUM62_20_VALUES_MAX] = {
  491. BLAST_MATRIX_NOMINAL,
  492. BLAST_MATRIX_NOMINAL,
  493. BLAST_MATRIX_NOMINAL,
  494. BLAST_MATRIX_NOMINAL,
  495. BLAST_MATRIX_NOMINAL,
  496. BLAST_MATRIX_NOMINAL,
  497. BLAST_MATRIX_NOMINAL,
  498. BLAST_MATRIX_NOMINAL,
  499. BLAST_MATRIX_NOMINAL,
  500. BLAST_MATRIX_NOMINAL,
  501. BLAST_MATRIX_NOMINAL,
  502. BLAST_MATRIX_NOMINAL,
  503. BLAST_MATRIX_NOMINAL,
  504. BLAST_MATRIX_NOMINAL,
  505. BLAST_MATRIX_NOMINAL,
  506. BLAST_MATRIX_NOMINAL,
  507. BLAST_MATRIX_NOMINAL,
  508. BLAST_MATRIX_NOMINAL,
  509. BLAST_MATRIX_NOMINAL,
  510. BLAST_MATRIX_NOMINAL,
  511. BLAST_MATRIX_NOMINAL,
  512. BLAST_MATRIX_NOMINAL,
  513. BLAST_MATRIX_NOMINAL,
  514. BLAST_MATRIX_NOMINAL,
  515. BLAST_MATRIX_NOMINAL,
  516. BLAST_MATRIX_NOMINAL,
  517. BLAST_MATRIX_NOMINAL,
  518. BLAST_MATRIX_NOMINAL,
  519. BLAST_MATRIX_NOMINAL,
  520. BLAST_MATRIX_NOMINAL,
  521. BLAST_MATRIX_NOMINAL,
  522. BLAST_MATRIX_NOMINAL,
  523. BLAST_MATRIX_NOMINAL,
  524. BLAST_MATRIX_NOMINAL,
  525. BLAST_MATRIX_NOMINAL,
  526. BLAST_MATRIX_NOMINAL,
  527. BLAST_MATRIX_NOMINAL,
  528. BLAST_MATRIX_NOMINAL,
  529. BLAST_MATRIX_NOMINAL,
  530. BLAST_MATRIX_NOMINAL,
  531. BLAST_MATRIX_NOMINAL,
  532. BLAST_MATRIX_NOMINAL,
  533. BLAST_MATRIX_NOMINAL,
  534. BLAST_MATRIX_NOMINAL,
  535. BLAST_MATRIX_NOMINAL,
  536. BLAST_MATRIX_BEST,
  537. BLAST_MATRIX_NOMINAL,
  538. BLAST_MATRIX_NOMINAL,
  539. BLAST_MATRIX_NOMINAL,
  540. BLAST_MATRIX_NOMINAL,
  541. BLAST_MATRIX_NOMINAL,
  542. BLAST_MATRIX_NOMINAL,
  543. BLAST_MATRIX_NOMINAL,
  544. BLAST_MATRIX_NOMINAL,
  545. BLAST_MATRIX_NOMINAL,
  546. BLAST_MATRIX_NOMINAL,
  547. BLAST_MATRIX_NOMINAL,
  548. BLAST_MATRIX_NOMINAL,
  549. BLAST_MATRIX_NOMINAL,
  550. BLAST_MATRIX_NOMINAL,
  551. BLAST_MATRIX_NOMINAL,
  552. BLAST_MATRIX_NOMINAL,
  553. BLAST_MATRIX_NOMINAL,
  554. BLAST_MATRIX_NOMINAL,
  555. BLAST_MATRIX_NOMINAL
  556. };
  558. /*
  559. Allocates memory for the BlastScoreBlk*.
  560. */
  562. BlastScoreBlk*
  563. BlastScoreBlkNew(Uint1 alphabet, Int4 number_of_contexts)
  565. {
  566. BlastScoreBlk* sbp;
  568. sbp = (BlastScoreBlk*) calloc(1, sizeof(BlastScoreBlk));
  570. if (sbp != NULL)
  571. {
  572. sbp->alphabet_code = alphabet;
  573. if (alphabet != BLASTNA_SEQ_CODE)
  574. sbp->alphabet_size = BLASTAA_SIZE;
  575. else
  576. sbp->alphabet_size = BLASTNA_SIZE;
  578. /* set the alphabet type (protein or not). */
  579. switch (alphabet) {
  580. case BLASTAA_SEQ_CODE:
  581. sbp->protein_alphabet = TRUE;
  582. break;
  583. case BLASTNA_SEQ_CODE:
  584. sbp->protein_alphabet = FALSE;
  585. break;
  586. default:
  587. break;
  588. }
  590. sbp->matrix_struct = BlastMatrixAllocate(sbp->alphabet_size);
  591. if (sbp->matrix_struct == NULL)
  592. {
  593. sbp = BlastScoreBlkFree(sbp);
  594. return sbp;
  595. }
  596. sbp->matrix = sbp->matrix_struct->matrix;
  597. sbp->maxscore = (Int4 *) calloc(BLAST_MATRIX_SIZE, sizeof(Int4));
  598.         sbp->scale_factor = 1.0;
  599. sbp->number_of_contexts = number_of_contexts;
  600. sbp->sfp = (Blast_ScoreFreq**) 
  601.          calloc(sbp->number_of_contexts, sizeof(Blast_ScoreFreq*));
  602. sbp->kbp_std = (Blast_KarlinBlk**)
  603.          calloc(sbp->number_of_contexts, sizeof(Blast_KarlinBlk*));
  604. sbp->kbp_gap_std = (Blast_KarlinBlk**)
  605.          calloc(sbp->number_of_contexts, sizeof(Blast_KarlinBlk*));
  606. sbp->kbp_psi = (Blast_KarlinBlk**)
  607.          calloc(sbp->number_of_contexts, sizeof(Blast_KarlinBlk*));
  608. sbp->kbp_gap_psi = (Blast_KarlinBlk**)
  609.          calloc(sbp->number_of_contexts, sizeof(Blast_KarlinBlk*));
  610. }
  612. return sbp;
  613. }
  615. Blast_ScoreFreq*
  616. Blast_ScoreFreqDestruct(Blast_ScoreFreq* sfp)
  617. {
  618. if (sfp == NULL)
  619. return NULL;
  621. if (sfp->sprob0 != NULL)
  622. sfree(sfp->sprob0);
  623. sfree(sfp);
  624. return sfp;
  625. }
  627. /*
  628. Deallocates the Karlin Block.
  629. */
  630. Blast_KarlinBlk*
  631. Blast_KarlinBlkDestruct(Blast_KarlinBlk* kbp)
  633. {
  634. sfree(kbp);
  636. return kbp;
  637. }
  639. static SBLASTMatrixStructure*
  640. BlastMatrixDestruct(SBLASTMatrixStructure* matrix_struct)
  642. {
  644. if (matrix_struct == NULL)
  645. return NULL;
  647. sfree(matrix_struct);
  649. return matrix_struct;
  650. }
  652. BlastScoreBlk*
  653. BlastScoreBlkFree(BlastScoreBlk* sbp)
  655. {
  656.     Int4 index, rows;
  657.     if (sbp == NULL)
  658.         return NULL;
  659.     
  660.     for (index=0; index<sbp->number_of_contexts; index++) {
  661.         if (sbp->sfp)
  662.             sbp->sfp[index] = Blast_ScoreFreqDestruct(sbp->sfp[index]);
  663.         if (sbp->kbp_std)
  664.             sbp->kbp_std[index] = Blast_KarlinBlkDestruct(sbp->kbp_std[index]);
  665.         if (sbp->kbp_gap_std)
  666.             sbp->kbp_gap_std[index] = Blast_KarlinBlkDestruct(sbp->kbp_gap_std[index]);
  667.         if (sbp->kbp_psi)
  668.             sbp->kbp_psi[index] = Blast_KarlinBlkDestruct(sbp->kbp_psi[index]);
  669.         if (sbp->kbp_gap_psi)
  670.             sbp->kbp_gap_psi[index] = Blast_KarlinBlkDestruct(sbp->kbp_gap_psi[index]);
  671.     }
  672.     if (sbp->kbp_ideal)
  673.         sbp->kbp_ideal = Blast_KarlinBlkDestruct(sbp->kbp_ideal);
  674.     sfree(sbp->sfp);
  675.     sfree(sbp->kbp_std);
  676.     sfree(sbp->kbp_psi);
  677.     sfree(sbp->kbp_gap_std);
  678.     sfree(sbp->kbp_gap_psi);
  679.     sbp->matrix_struct = BlastMatrixDestruct(sbp->matrix_struct);
  680.     sfree(sbp->maxscore);
  681.     sbp->comments = ListNodeFreeData(sbp->comments);
  682.     sfree(sbp->name);
  683.     sfree(sbp->ambiguous_res);
  684.     
  685.     /* Removing posMatrix and posFreqs if any */
  686.     rows = sbp->query_length + 1;
  687.     if(sbp->posMatrix != NULL) {
  688.         for (index=0; index < rows; index++) {
  689.             sfree(sbp->posMatrix[index]);
  690.         }
  691.         sfree(sbp->posMatrix);
  692.     }
  693.     
  694.     if(sbp->posFreqs != NULL) {
  695.         for (index = 0; index < rows; index++) {
  696.             sfree(sbp->posFreqs[index]);
  697.         }
  698.         sfree(sbp->posFreqs);
  699.     }
  700.     
  701.     sfree(sbp);
  702.     
  703.     return NULL;
  704. }
  706. /* 
  707. Set the ambiguous residue (e.g, 'N', 'X') in the BlastScoreBlk*.
  708. Convert from ncbieaa to sbp->alphabet_code (i.e., ncbistdaa) first.
  709. */
  710. Int2
  711. BLAST_ScoreSetAmbigRes(BlastScoreBlk* sbp, char ambiguous_res)
  713. {
  714. Int2 index;
  715. Uint1* ambig_buffer;
  717. if (sbp == NULL)
  718. return 1;
  720. if (sbp->ambig_occupy >= sbp->ambig_size)
  721. {
  722. sbp->ambig_size += 5;
  723. ambig_buffer = (Uint1 *) calloc(sbp->ambig_size, sizeof(Uint1));
  724. for (index=0; index<sbp->ambig_occupy; index++)
  725. {
  726. ambig_buffer[index] = sbp->ambiguous_res[index];
  727. }
  728. sfree(sbp->ambiguous_res);
  729. sbp->ambiguous_res = ambig_buffer;
  730. }
  732. if (sbp->alphabet_code == BLASTAA_SEQ_CODE)
  733. {
  734. sbp->ambiguous_res[sbp->ambig_occupy] = 
  735.          AMINOACID_TO_NCBISTDAA[toupper(ambiguous_res)];
  736. }
  737. else {
  738.       if (sbp->alphabet_code == BLASTNA_SEQ_CODE)
  739.          sbp->ambiguous_res[sbp->ambig_occupy] = 
  740.             IUPACNA_TO_BLASTNA[toupper(ambiguous_res)];
  741.       else if (sbp->alphabet_code == NCBI4NA_SEQ_CODE)
  742.          sbp->ambiguous_res[sbp->ambig_occupy] = 
  743.             IUPACNA_TO_NCBI4NA[toupper(ambiguous_res)];
  744.    }
  745. (sbp->ambig_occupy)++;
  748. return 0;
  749. }
  751. /* 
  752. Fill in the matrix for blastn using the penaly and rewards
  754. The query sequence alphabet is blastna, the subject sequence
  755. is ncbi2na.  The alphabet blastna is defined in blastkar.h
  756. and the first four elements of blastna are identical to ncbi2na.
  758. The query is in the first index, the subject is the second.
  759.         if matrix==NULL, it is allocated and returned.
  760. */
  762. static Int4 **BlastScoreBlkMatCreateEx(Int4 **matrix,Int4 penalty, 
  763.                                        Int4 reward)
  764. {
  766. Int2 index1, index2, degen;
  767. Int2 degeneracy[BLASTNA_SIZE+1];
  769.         if(!matrix) {
  770.             SBLASTMatrixStructure* matrix_struct;
  771.             matrix_struct =BlastMatrixAllocate((Int2) BLASTNA_SIZE);
  772.             matrix = matrix_struct->matrix;
  773.         }
  774. for (index1 = 0; index1<BLASTNA_SIZE; index1++) /* blastna */
  775. for (index2 = 0; index2<BLASTNA_SIZE; index2++) /* blastna */
  776. matrix[index1][index2] = 0;
  778. /* In blastna the 1st four bases are A, C, G, and T, exactly as it is ncbi2na. */
  779. /* ncbi4na gives them the value 1, 2, 4, and 8.  */
  781. /* Set the first four bases to degen. one */
  782. for (index1=0; index1<NUMBER_NON_AMBIG_BP; index1++)
  783. degeneracy[index1] = 1;
  785. for (index1=NUMBER_NON_AMBIG_BP; index1<BLASTNA_SIZE; index1++) /* blastna */
  786. {
  787. degen=0;
  788. for (index2=0; index2<NUMBER_NON_AMBIG_BP; index2++) /* ncbi2na */
  789. {
  790. if (BLASTNA_TO_NCBI4NA[index1] & BLASTNA_TO_NCBI4NA[index2])
  791. degen++;
  792. }
  793. degeneracy[index1] = degen;
  794. }
  797. for (index1=0; index1<BLASTNA_SIZE; index1++) /* blastna */
  798. {
  799. for (index2=index1; index2<BLASTNA_SIZE; index2++) /* blastna */
  800. {
  801. if (BLASTNA_TO_NCBI4NA[index1] & BLASTNA_TO_NCBI4NA[index2])
  802. { /* round up for positive scores, down for negatives. */
  803. matrix[index1][index2] = BLAST_Nint( (double) ((degeneracy[index2]-1)*penalty + reward))/degeneracy[index2];
  804. if (index1 != index2)
  805. {
  806.       matrix[index2][index1] = matrix[index1][index2];
  807. }
  808. }
  809. else
  810. {
  811. matrix[index1][index2] = penalty;
  812. matrix[index2][index1] = penalty;
  813. }
  814. }
  815. }
  817.         /* The value of 15 is a gap, which is a sentinel between strands in 
  818.            the ungapped extension algorithm */
  819. for (index1=0; index1<BLASTNA_SIZE; index1++) /* blastna */
  820.            matrix[BLASTNA_SIZE-1][index1] = INT4_MIN / 2;
  821. for (index1=0; index1<BLASTNA_SIZE; index1++) /* blastna */
  822.            matrix[index1][BLASTNA_SIZE-1] = INT4_MIN / 2;
  824. return matrix;
  825. }
  826. /* 
  827. Fill in the matrix for blastn using the penaly and rewards on
  828. the BlastScoreBlk*.
  830. The query sequence alphabet is blastna, the subject sequence
  831. is ncbi2na.  The alphabet blastna is defined in blastkar.h
  832. and the first four elements of blastna are identical to ncbi2na.
  834. The query is in the first index, the subject is the second.
  835. */
  836. static Int2 BlastScoreBlkMatCreate(BlastScoreBlk* sbp)
  837. {
  838.    sbp->matrix = BlastScoreBlkMatCreateEx(sbp->matrix,sbp->penalty, sbp->reward);
  839. sbp->mat_dim1 = BLASTNA_SIZE;
  840. sbp->mat_dim2 = BLASTNA_SIZE;
  842. return 0;
  843. }
  845. /* 
  846. Read in the matrix from the FILE *fp.
  848. This function ASSUMES that the matrices are in the ncbistdaa
  849. format.  BLAST should be able to use any alphabet, though it
  850. is expected that ncbistdaa will be used.
  851. */
  853. static Int2
  854. BlastScoreBlkMatRead(BlastScoreBlk* sbp, FILE *fp)
  855. {
  856.     char buf[512+3];
  857.     char temp[512];
  858.     char* cp,* lp;
  859.     char ch;
  860.     Int4 ** matrix;
  861.     Int4 * m;
  862.     Int4 score;
  863.     Uint4 a1cnt = 0, a2cnt = 0;
  864.     char    a1chars[BLAST_MAX_ALPHABET], a2chars[BLAST_MAX_ALPHABET];
  865.     long lineno = 0;
  866.     double xscore;
  867.     register int index1, index2;
  868.     Int2 status;
  869.     
  870.     matrix = sbp->matrix;
  871.     
  872.     if (sbp->alphabet_code != BLASTNA_SEQ_CODE) {
  873.         for (index1 = 0; index1 < sbp->alphabet_size; index1++)
  874.             for (index2 = 0; index2 < sbp->alphabet_size; index2++)
  875.                 matrix[index1][index2] = BLAST_SCORE_MIN;
  876.     } else {
  877.         /* Fill-in all the defaults ambiguity and normal codes */
  878.         status=BlastScoreBlkMatCreate(sbp); 
  879.         if(status != 0)
  880. {
  881.          return status;
  882. }
  883.     }
  884.     
  885.     /* Read the residue names for the second alphabet */
  886.     while (fgets(buf, sizeof(buf), fp) != NULL) {
  887.         ++lineno;
  888.         if (strchr(buf, 'n') == NULL) {
  889.             return 2;
  890.         }
  892.         if (buf[0] == COMMENT_CHR) {
  893.             /* save the comment line in a linked list */
  894.             *strchr(buf, 'n') = NULLB;
  895.             ListNodeCopyStr(&sbp->comments, 0, buf+1);
  896.             continue;
  897.         }
  898.         if ((cp = strchr(buf, COMMENT_CHR)) != NULL)
  899.             *cp = NULLB;
  900.         lp = (char*)strtok(buf, TOKSTR);
  901.         if (lp == NULL) /* skip blank lines */
  902.             continue;
  903.         while (lp != NULL) {
  904.            if (sbp->alphabet_code == BLASTAA_SEQ_CODE)
  905.               ch = AMINOACID_TO_NCBISTDAA[toupper(*lp)];
  906.            else if (sbp->alphabet_code == BLASTNA_SEQ_CODE) {
  907.               ch = IUPACNA_TO_BLASTNA[toupper(*lp)];
  908.            } else {
  909.               ch = *lp;
  910.            }
  911.             a2chars[a2cnt++] = ch;
  912.             lp = (char*)strtok(NULL, TOKSTR);
  913.         }
  914.         
  915.         break; /* Exit loop after reading one line. */
  916.     }
  917.     
  918.     if (a2cnt <= 1) { 
  919.         return 2;
  920.     }
  922.     if (sbp->alphabet_code != BLASTNA_SEQ_CODE) {
  923.         sbp->mat_dim2 = a2cnt;
  924.     }
  925.     while (fgets(buf, sizeof(buf), fp) != NULL)  {
  926.         ++lineno;
  927.         if ((cp = strchr(buf, 'n')) == NULL) {
  928.             return 2;
  929.         }
  930.         if ((cp = strchr(buf, COMMENT_CHR)) != NULL)
  931.             *cp = NULLB;
  932.         if ((lp = (char*)strtok(buf, TOKSTR)) == NULL)
  933.             continue;
  934.         ch = *lp;
  935.         cp = (char*) lp;
  936.         if ((cp = strtok(NULL, TOKSTR)) == NULL) {
  937.             return 2;
  938.         }
  939.         if (a1cnt >= DIM(a1chars)) {
  940.             return 2;
  941.         }
  943.         if (sbp->alphabet_code == BLASTAA_SEQ_CODE) {
  944.            ch = AMINOACID_TO_NCBISTDAA[toupper(ch)];
  945.         } else {
  946.             if (sbp->alphabet_code == BLASTNA_SEQ_CODE) {
  947.                 ch = IUPACNA_TO_BLASTNA[toupper(ch)];
  948.             }
  949.         }
  950.         a1chars[a1cnt++] = ch;
  951.         m = &matrix[(int)ch][0];
  952.         index2 = 0;
  953.         while (cp != NULL) {
  954.             if (index2 >= (int) a2cnt) {
  955.                 return 2;
  956.             }
  957.             strcpy(temp, cp);
  958.             
  959.             if (strcasecmp(temp, "na") == 0)  {
  960.                 score = BLAST_SCORE_1MIN;
  961.             } else  {
  962.                 if (sscanf(temp, "%lg", &xscore) != 1) {
  963.                     return 2;
  964.                 }
  965. /*xscore = MAX(xscore, BLAST_SCORE_1MIN);*/
  966.                 if (xscore > BLAST_SCORE_1MAX || xscore < BLAST_SCORE_1MIN) {
  967.                     return 2;
  968.                 }
  969.                 xscore += (xscore >= 0. ? 0.5 : -0.5);
  970.                 score = (Int4)xscore;
  971.             }
  972.             
  973.             m[(int)a2chars[index2++]] = score;
  974.             
  975.             cp = strtok(NULL, TOKSTR);
  976.         }
  977.     }
  978.     
  979.     if (a1cnt <= 1) {
  980.         return 2;
  981.     }
  982.     
  983.     if (sbp->alphabet_code != BLASTNA_SEQ_CODE) {
  984.         sbp->mat_dim1 = a1cnt;
  985.     }
  986.     
  987.     return 0;
  988. }
  990. static Int2
  991. BlastScoreBlkMaxScoreSet(BlastScoreBlk* sbp)
  992. {
  993. Int4 score, maxscore;
  994. Int4 ** matrix; 
  995. Int2 index1, index2;
  997. sbp->loscore = BLAST_SCORE_1MAX;
  998.         sbp->hiscore = BLAST_SCORE_1MIN;
  999. matrix = sbp->matrix;
  1000. for (index1=0; index1<sbp->alphabet_size; index1++)
  1001. {
  1002. maxscore=BLAST_SCORE_MIN;
  1003. for (index2=0; index2<sbp->alphabet_size; index2++)
  1004. {
  1005. score = matrix[index1][index2];
  1006. if (score <= BLAST_SCORE_MIN || score >= BLAST_SCORE_MAX)
  1007. continue;
  1008. if (score > maxscore)
  1009. {
  1010. maxscore = score;
  1011. }
  1012. if (sbp->loscore > score)
  1013. sbp->loscore = score;
  1014. if (sbp->hiscore < score)
  1015. sbp->hiscore = score;
  1016. }
  1017. sbp->maxscore[index1] = maxscore;
  1018. }
  1019. /* If the lo/hi-scores are BLAST_SCORE_MIN/BLAST_SCORE_MAX, (i.e., for
  1020. gaps), then use other scores. */
  1022. if (sbp->loscore < BLAST_SCORE_1MIN)
  1023. sbp->loscore = BLAST_SCORE_1MIN;
  1024. if (sbp->hiscore > BLAST_SCORE_1MAX)
  1025. sbp->hiscore = BLAST_SCORE_1MAX;
  1027. return 0;
  1028. }
  1030. Int2
  1031. BlastScoreBlkMatrixLoad(BlastScoreBlk* sbp)
  1032. {
  1033.     Int2 status = 0;
  1034.     SNCBIPackedScoreMatrix* psm;
  1035.     Int4** matrix = sbp->matrix;
  1036.     int i, j;   /* loop indices */
  1038.     ASSERT(sbp);
  1040.     if (strcasecmp(sbp->name, "BLOSUM62") == 0) {
  1041.         psm = (SNCBIPackedScoreMatrix*) &NCBISM_Blosum62;
  1042.     } else if (strcasecmp(sbp->name, "BLOSUM45") == 0) {
  1043.         psm = (SNCBIPackedScoreMatrix*) &NCBISM_Blosum45;
  1044.     } else if (strcasecmp(sbp->name, "BLOSUM80") == 0) {
  1045.         psm = (SNCBIPackedScoreMatrix*) &NCBISM_Blosum80;
  1046.     } else if (strcasecmp(sbp->name, "PAM30") == 0) {
  1047.         psm = (SNCBIPackedScoreMatrix*) &NCBISM_Pam30;
  1048.     } else if (strcasecmp(sbp->name, "PAM70") == 0) {
  1049.         psm = (SNCBIPackedScoreMatrix*) &NCBISM_Pam70;
  1050.     } else {
  1051.         return -1;
  1052.     }
  1054.     /* Initialize with BLAST_SCORE_MIN */
  1055.     for (i = 0; i < sbp->alphabet_size; i++) {
  1056.         for (j = 0; j < sbp->alphabet_size; j++) {
  1057.             matrix[i][j] = BLAST_SCORE_MIN;
  1058.         }
  1059.     }
  1061.     for (i = 0; i < sbp->alphabet_size; i++) {
  1062.         for (j = 0; j < sbp->alphabet_size; j++) {
  1063.             /* skip selenocysteine and gap */
  1064.             if (i == AMINOACID_TO_NCBISTDAA['U'] || 
  1065.                 i == AMINOACID_TO_NCBISTDAA['-'] ||
  1066.                 j == AMINOACID_TO_NCBISTDAA['U'] || 
  1067.                 j == AMINOACID_TO_NCBISTDAA['-']) {
  1068.                 continue;
  1069.             }
  1070.             matrix[i][j] = NCBISM_GetScore((const SNCBIPackedScoreMatrix*) psm,
  1071.                                            i, j);
  1072.         }
  1073.     }
  1074.     /* Sets dimensions of matrix. */
  1075.     sbp->mat_dim1 = sbp->mat_dim2 = sbp->alphabet_size;
  1077.     return status;
  1078. }
  1080. Int2
  1081. BLAST_ScoreBlkMatFill(BlastScoreBlk* sbp, char* matrix_path)
  1082. {
  1083.     Int2 status = 0;
  1084.     
  1085.     if (sbp->read_in_matrix) {
  1086.         
  1087.         if (matrix_path && *matrix_path != NULLB) {
  1089.             FILE *fp = NULL;
  1090.             char* full_matrix_path = NULL;
  1091.             int path_len = strlen(matrix_path);
  1092.             int buflen = path_len + strlen(sbp->name);
  1094.             full_matrix_path = (char*) malloc((buflen + 1) * sizeof(char));
  1095.             if (!full_matrix_path) {
  1096.                 return -1;
  1097.             }
  1098.             strncpy(full_matrix_path, matrix_path, buflen);
  1099.             strncat(full_matrix_path, sbp->name, buflen - path_len);
  1101.             if ( (fp=fopen(full_matrix_path, "r")) == NULL) {
  1102.                return -1;
  1103.             }
  1104.             sfree(full_matrix_path);
  1106.             if ( (status=BlastScoreBlkMatRead(sbp, fp)) != 0) {
  1107.                fclose(fp);
  1108.                return status;
  1109.             }
  1110.             fclose(fp);
  1112.         } else {
  1113.             if ( (status = BlastScoreBlkMatrixLoad(sbp)) !=0) {
  1114.                 return status;
  1115.             }
  1116.         }
  1117.     } else {
  1118.        /* Nucleotide BLAST only! */
  1119.        if ( (status=BlastScoreBlkMatCreate(sbp)) != 0)
  1120.           return status;
  1121.     }
  1122.     
  1123.     if ( (status=BlastScoreBlkMaxScoreSet(sbp)) != 0)
  1124.        return status;
  1126.     return status;
  1127. }
  1129. Blast_ResFreq*
  1130. Blast_ResFreqDestruct(Blast_ResFreq* rfp)
  1131. {
  1132. if (rfp == NULL)
  1133. return NULL;
  1135. if (rfp->prob0 != NULL)
  1136. sfree(rfp->prob0);
  1138. sfree(rfp);
  1140. return rfp;
  1141. }
  1143. /*
  1144. Allocates the Blast_ResFreq* and then fills in the frequencies
  1145. in the probabilities.
  1146. */ 
  1148. Blast_ResFreq*
  1149. Blast_ResFreqNew(const BlastScoreBlk* sbp)
  1150. {
  1151. Blast_ResFreq* rfp;
  1153. if (sbp == NULL)
  1154. {
  1155. return NULL;
  1156. }
  1158. rfp = (Blast_ResFreq*) calloc(1, sizeof(Blast_ResFreq));
  1159. if (rfp == NULL)
  1160. return NULL;
  1162. rfp->alphabet_code = sbp->alphabet_code;
  1164. rfp->prob0 = (double*) calloc(sbp->alphabet_size, sizeof(double));
  1165. if (rfp->prob0 == NULL) 
  1166. {
  1167. rfp = Blast_ResFreqDestruct(rfp);
  1168. return rfp;
  1169. }
  1170. rfp->prob = rfp->prob0 - sbp->alphabet_start;
  1172. return rfp;
  1173. }
  1175. typedef struct BLAST_LetterProb {
  1176. char ch;
  1177. double p;
  1178. } BLAST_LetterProb;
  1180. #if 0
  1182. /* Unused for right now, but do not remove */
  1184. /*  M. O. Dayhoff amino acid background frequencies   */
  1185. static BLAST_LetterProb Dayhoff_prob[] = {
  1186. { 'A', 87.13 },
  1187. { 'C', 33.47 },
  1188. { 'D', 46.87 },
  1189. { 'E', 49.53 },
  1190. { 'F', 39.77 },
  1191. { 'G', 88.61 },
  1192. { 'H', 33.62 },
  1193. { 'I', 36.89 },
  1194. { 'K', 80.48 },
  1195. { 'L', 85.36 },
  1196. { 'M', 14.75 },
  1197. { 'N', 40.43 },
  1198. { 'P', 50.68 },
  1199. { 'Q', 38.26 },
  1200. { 'R', 40.90 },
  1201. { 'S', 69.58 },
  1202. { 'T', 58.54 },
  1203. { 'V', 64.72 },
  1204. { 'W', 10.49 },
  1205. { 'Y', 29.92 }
  1206. };
  1208. /* Stephen Altschul amino acid background frequencies */
  1209. static BLAST_LetterProb Altschul_prob[] = {
  1210. { 'A', 81.00 },
  1211. { 'C', 15.00 },
  1212. { 'D', 54.00 },
  1213. { 'E', 61.00 },
  1214. { 'F', 40.00 },
  1215. { 'G', 68.00 },
  1216. { 'H', 22.00 },
  1217. { 'I', 57.00 },
  1218. { 'K', 56.00 },
  1219. { 'L', 93.00 },
  1220. { 'M', 25.00 },
  1221. { 'N', 45.00 },
  1222. { 'P', 49.00 },
  1223. { 'Q', 39.00 },
  1224. { 'R', 57.00 },
  1225. { 'S', 68.00 },
  1226. { 'T', 58.00 },
  1227. { 'V', 67.00 },
  1228. { 'W', 13.00 },
  1229. { 'Y', 32.00 }
  1230. };
  1231. #endif
  1233. /* amino acid background frequencies from Robinson and Robinson */
  1234. static BLAST_LetterProb Robinson_prob[] = {
  1235. { 'A', 78.05 },
  1236. { 'C', 19.25 },
  1237. { 'D', 53.64 },
  1238. { 'E', 62.95 },
  1239. { 'F', 38.56 },
  1240. { 'G', 73.77 },
  1241. { 'H', 21.99 },
  1242. { 'I', 51.42 },
  1243. { 'K', 57.44 },
  1244. { 'L', 90.19 },
  1245. { 'M', 22.43 },
  1246. { 'N', 44.87 },
  1247. { 'P', 52.03 },
  1248. { 'Q', 42.64 },
  1249. { 'R', 51.29 },
  1250. { 'S', 71.20 },
  1251. { 'T', 58.41 },
  1252. { 'V', 64.41 },
  1253. { 'W', 13.30 },
  1254. { 'Y', 32.16 }
  1255. };
  1257. #define STD_AMINO_ACID_FREQS Robinson_prob
  1259. static BLAST_LetterProb nt_prob[] = {
  1260. { 'A', 25.00 },
  1261. { 'C', 25.00 },
  1262. { 'G', 25.00 },
  1263. { 'T', 25.00 }
  1264. };
  1266. /*
  1267. Normalize the frequencies to "norm".
  1268. */
  1269. static Int2
  1270. Blast_ResFreqNormalize(const BlastScoreBlk* sbp, Blast_ResFreq* rfp, double norm)
  1271. {
  1272. Int2 alphabet_stop, index;
  1273. double sum = 0., p;
  1275. if (norm == 0.)
  1276. return 1;
  1278. alphabet_stop = sbp->alphabet_start + sbp->alphabet_size;
  1279. for (index=sbp->alphabet_start; index<alphabet_stop; index++)
  1280. {
  1281. p = rfp->prob[index];
  1282. if (p < 0.)
  1283. return 1;
  1284. sum += p;
  1285. }
  1286. if (sum <= 0.)
  1287. return 0;
  1289. for (index=sbp->alphabet_start; index<alphabet_stop; index++)
  1290. {
  1291. rfp->prob[index] /= sum;
  1292. rfp->prob[index] *= norm;
  1293. }
  1294. return 0;
  1295. }
  1297. Int2
  1298. Blast_GetStdAlphabet(Uint1 alphabet_code, Uint1* residues, Uint4 residues_size)
  1299. {
  1300. Int2 index;
  1302. if (residues_size < DIM(STD_AMINO_ACID_FREQS))
  1303. return -2;
  1305. for (index=0; index<(int)DIM(STD_AMINO_ACID_FREQS); index++) 
  1306. {
  1307. if (alphabet_code == BLASTAA_SEQ_CODE)
  1308. {
  1309.   residues[index] = 
  1310.             AMINOACID_TO_NCBISTDAA[toupper(STD_AMINO_ACID_FREQS[index].ch)];
  1311. }
  1312. else
  1313. {
  1314. residues[index] = STD_AMINO_ACID_FREQS[index].ch;
  1315. }
  1316. }
  1318. return index;
  1319. }
  1321. Int2
  1322. Blast_ResFreqStdComp(const BlastScoreBlk* sbp, Blast_ResFreq* rfp)
  1323. {
  1324. Int2 retval;
  1325.         Uint4 index;
  1326. Uint1* residues;
  1328. if (sbp->protein_alphabet == TRUE)
  1329. {
  1330. residues = (Uint1*) calloc(DIM(STD_AMINO_ACID_FREQS), sizeof(Uint1));
  1331. retval = Blast_GetStdAlphabet(sbp->alphabet_code, residues, DIM(STD_AMINO_ACID_FREQS));
  1332. if (retval < 1)
  1333. return retval;
  1335. for (index=0; index<DIM(STD_AMINO_ACID_FREQS); index++) 
  1336. {
  1337. rfp->prob[residues[index]] = STD_AMINO_ACID_FREQS[index].p;
  1338. }
  1339. sfree(residues);
  1340. }
  1341. else
  1342. { /* beginning of blastna and ncbi2na are the same. */
  1343. /* Only blastna used  for nucleotides. */
  1344. for (index=0; index<DIM(nt_prob); index++) 
  1345. {
  1346. rfp->prob[index] = nt_prob[index].p;
  1347. }
  1348. }
  1350. Blast_ResFreqNormalize(sbp, rfp, 1.0);
  1352. return 0;
  1353. }
  1355. typedef struct Blast_ResComp {
  1356.     Uint1 alphabet_code;
  1357.     Int4* comp;  /* composition of alphabet, array starts at beginning of alphabet. */
  1358.     Int4*   comp0; /* Same array as above, starts at zero. */
  1359. } Blast_ResComp;
  1361. static Blast_ResComp*
  1362. BlastResCompDestruct(Blast_ResComp* rcp)
  1363. {
  1364. if (rcp == NULL)
  1365. return NULL;
  1367. if (rcp->comp0 != NULL)
  1368. sfree(rcp->comp0);
  1370. sfree(rcp);
  1371. return NULL;
  1372. }
  1374. /* 
  1375. Allocated the Blast_ResComp* for a given alphabet.  Only the
  1376. alphabets ncbistdaa and ncbi4na should be used by BLAST.
  1377. */
  1378. static Blast_ResComp*
  1379. BlastResCompNew(BlastScoreBlk* sbp)
  1380. {
  1381. Blast_ResComp* rcp;
  1383. rcp = (Blast_ResComp*) calloc(1, sizeof(Blast_ResComp));
  1384. if (rcp == NULL)
  1385. return NULL;
  1387. rcp->alphabet_code = sbp->alphabet_code;
  1389. /* comp0 has zero offset, comp starts at 0, only one 
  1390. array is allocated.  */
  1391. rcp->comp0 = (Int4*) calloc(BLAST_MATRIX_SIZE, sizeof(Int4));
  1392. if (rcp->comp0 == NULL) 
  1393. {
  1394. rcp = BlastResCompDestruct(rcp);
  1395. return rcp;
  1396. }
  1398. rcp->comp = rcp->comp0 - sbp->alphabet_start;
  1400. return rcp;
  1401. }
  1403. /*
  1404. Store the composition of a (query) string.  
  1405. */
  1406. static Int2
  1407. BlastResCompStr(BlastScoreBlk* sbp, Blast_ResComp* rcp, char* str, Int4 length)
  1408. {
  1409. char* lp,* lpmax;
  1410. Int2 index;
  1411.         Uint1 mask;
  1413. if (sbp == NULL || rcp == NULL || str == NULL)
  1414. return 1;
  1416. if (rcp->alphabet_code != sbp->alphabet_code)
  1417. return 1;
  1419.         /* For megablast, check only the first 4 bits of the sequence values */
  1420.         mask = (sbp->protein_alphabet ? 0xff : 0x0f);
  1422. /* comp0 starts at zero and extends for "num", comp is the same array, but 
  1423. "start_at" offset. */
  1424. for (index=0; index<(sbp->alphabet_size); index++)
  1425. rcp->comp0[index] = 0;
  1427. for (lp = str, lpmax = lp+length; lp < lpmax; lp++)
  1428. {
  1429. ++rcp->comp[(int)(*lp & mask)];
  1430. }
  1432. /* Don't count ambig. residues. */
  1433. for (index=0; index<sbp->ambig_occupy; index++)
  1434. {
  1435. rcp->comp[sbp->ambiguous_res[index]] = 0;
  1436. }
  1438. return 0;
  1439. }
  1441. static Int2
  1442. Blast_ResFreqClr(const BlastScoreBlk* sbp, Blast_ResFreq* rfp)
  1443. {
  1444. Int2 alphabet_max, index;
  1445.  
  1446. alphabet_max = sbp->alphabet_start + sbp->alphabet_size;
  1447. for (index=sbp->alphabet_start; index<alphabet_max; index++)
  1448.                 rfp->prob[index] = 0.0;
  1449.  
  1450.         return 0;
  1451. }
  1453. /*
  1454. Calculate the residue frequencies associated with the provided ResComp
  1455. */
  1456. static Int2
  1457. Blast_ResFreqResComp(BlastScoreBlk* sbp, Blast_ResFreq* rfp, 
  1458.                      Blast_ResComp* rcp)
  1459. {
  1460. Int2 alphabet_max, index;
  1461. double sum = 0.;
  1463. if (rfp == NULL || rcp == NULL)
  1464. return 1;
  1466. if (rfp->alphabet_code != rcp->alphabet_code)
  1467. return 1;
  1469. alphabet_max = sbp->alphabet_start + sbp->alphabet_size;
  1470. for (index=sbp->alphabet_start; index<alphabet_max; index++)
  1471. sum += rcp->comp[index];
  1473. if (sum == 0.) {
  1474. Blast_ResFreqClr(sbp, rfp);
  1475. return 0;
  1476. }
  1478. for (index=sbp->alphabet_start; index<alphabet_max; index++)
  1479. rfp->prob[index] = rcp->comp[index] / sum;
  1481. return 0;
  1482. }
  1484. static Int2
  1485. Blast_ResFreqString(BlastScoreBlk* sbp, Blast_ResFreq* rfp, char* string, Int4 length)
  1486. {
  1487. Blast_ResComp* rcp;
  1489. rcp = BlastResCompNew(sbp);
  1491. BlastResCompStr(sbp, rcp, string, length);
  1493. Blast_ResFreqResComp(sbp, rfp, rcp);
  1495. rcp = BlastResCompDestruct(rcp);
  1497. return 0;
  1498. }
  1500. static Int2
  1501. BlastScoreChk(Int4 lo, Int4 hi)
  1502. {
  1503. if (lo >= 0 || hi <= 0 ||
  1504. lo < BLAST_SCORE_1MIN || hi > BLAST_SCORE_1MAX)
  1505. return 1;
  1507. if (hi - lo > BLAST_SCORE_RANGE_MAX)
  1508. return 1;
  1510. return 0;
  1511. }
  1513. Blast_ScoreFreq*
  1514. Blast_ScoreFreqNew(Int4 score_min, Int4 score_max)
  1515. {
  1516. Blast_ScoreFreq* sfp;
  1517. Int4 range;
  1519. if (BlastScoreChk(score_min, score_max) != 0)
  1520. return NULL;
  1522. sfp = (Blast_ScoreFreq*) calloc(1, sizeof(Blast_ScoreFreq));
  1523. if (sfp == NULL)
  1524. return NULL;
  1526. range = score_max - score_min + 1;
  1527. sfp->sprob = (double*) calloc(range, sizeof(double));
  1528. if (sfp->sprob == NULL) 
  1529. {
  1530. Blast_ScoreFreqDestruct(sfp);
  1531. return NULL;
  1532. }
  1534. sfp->sprob0 = sfp->sprob;
  1535. sfp->sprob -= score_min;        /* center around 0 */
  1536. sfp->score_min = score_min;
  1537. sfp->score_max = score_max;
  1538. sfp->obs_min = sfp->obs_max = 0;
  1539. sfp->score_avg = 0.0;
  1540. return sfp;
  1541. }
  1543. static Int2
  1544. BlastScoreFreqCalc(BlastScoreBlk* sbp, Blast_ScoreFreq* sfp, Blast_ResFreq* rfp1, Blast_ResFreq* rfp2)
  1545. {
  1546. Int4 ** matrix;
  1547. Int4 score, obs_min, obs_max;
  1548. double score_sum, score_avg;
  1549. Int2  alphabet_start, alphabet_end, index1, index2;
  1551. if (sbp == NULL || sfp == NULL)
  1552. return 1;
  1554. if (sbp->loscore < sfp->score_min || sbp->hiscore > sfp->score_max)
  1555. return 1;
  1557. for (score = sfp->score_min; score <= sfp->score_max; score++)
  1558. sfp->sprob[score] = 0.0;
  1560. matrix = sbp->matrix;
  1562. alphabet_start = sbp->alphabet_start;
  1563. alphabet_end = alphabet_start + sbp->alphabet_size;
  1564. for (index1=alphabet_start; index1<alphabet_end; index1++)
  1565. {
  1566. for (index2=alphabet_start; index2<alphabet_end; index2++)
  1567. {
  1568. score = matrix[index1][index2];
  1569. if (score >= sbp->loscore) 
  1570. {
  1571. sfp->sprob[score] += rfp1->prob[index1] * rfp2->prob[index2];
  1572. }
  1573. }
  1574. }
  1576. score_sum = 0.;
  1577. obs_min = obs_max = BLAST_SCORE_MIN;
  1578. for (score = sfp->score_min; score <= sfp->score_max; score++)
  1579. {
  1580. if (sfp->sprob[score] > 0.) 
  1581. {
  1582. score_sum += sfp->sprob[score];
  1583. obs_max = score;
  1584. if (obs_min == BLAST_SCORE_MIN)
  1585. obs_min = score;
  1586. }
  1587. }
  1588. sfp->obs_min = obs_min;
  1589. sfp->obs_max = obs_max;
  1591. score_avg = 0.0;
  1592. if (score_sum > 0.0001 || score_sum < -0.0001)
  1593. {
  1594. for (score = obs_min; score <= obs_max; score++) 
  1595. {
  1596. sfp->sprob[score] /= score_sum;
  1597. score_avg += score * sfp->sprob[score];
  1598. }
  1599. }
  1600. sfp->score_avg = score_avg;
  1602. return 0;
  1603. }
  1606. #define DIMOFP0 (iterlimit*range + 1)
  1607. #define DIMOFP0_MAX (BLAST_KARLIN_K_ITER_MAX*BLAST_SCORE_RANGE_MAX+1)
  1610. #define SMALL_LAMBDA_THRESHOLD 20 /*defines special case in K computation*/
  1611.                                 /*threshold is on exp(-Lambda)*/
  1613. /*The following procedure computes K. The input includes Lambda, H,
  1614.  *  and an array of probabilities for each score.
  1615.  *  There are distinct closed form for three cases:
  1616.  *  1. high score is 1 low score is -1
  1617.  *  2. high score is 1 low score is not -1
  1618.  *  3. low score is -1, high score is not 1
  1619.  *
  1620.  * Otherwise, in most cases the value is computed as:
  1621.  * -exp(-2.0*outerSum) / ((H/lambda)*(exp(-lambda) - 1)
  1622.  * The last term (exp(-lambda) - 1) can be computed in two different
  1623.  * ways depending on whether lambda is small or not.
  1624.  * outerSum is a sum of the terms
  1625.  * innerSum/j, where j is denoted by iterCounter in the code.
  1626.  * The sum is truncated when the new term innersum/j i sufficiently small.
  1627.  * innerSum is a weighted sum of the probabilities of
  1628.  * of achieving a total score i in a gapless alignment,
  1629.  * which we denote by P(i,j).
  1630.  * of exactly j characters. innerSum(j) has two parts
  1631.  * Sum over i < 0  P(i,j)exp(-i * lambda) +
  1632.  * Sum over i >=0  P(i,j)
  1633.  * The terms P(i,j) are computed by dynamic programming.
  1634.  * An earlier version was flawed in that ignored the special case 1
  1635.  * and tried to replace the tail of the computation of outerSum
  1636.  * by a geometric series, but the base of the geometric series
  1637.  * was not accurately estimated in some cases.
  1638.  */
  1640. static double
  1641. BlastKarlinLHtoK(Blast_ScoreFreq* sfp, double lambda, double H)
  1642. {
  1643.     /*The next array stores the probabilities of getting each possible
  1644.       score in an alignment of fixed length; the array is shifted
  1645.       during part of the computation, so that
  1646.       entry 0 is for score 0.  */
  1647.     double         *alignmentScoreProbabilities = NULL;
  1648.     Int4            low;    /* Lowest score (must be negative) */
  1649.     Int4            high;   /* Highest score (must be positive) */
  1650.     Int4            range;  /* range of scores, computed as high - low*/
  1651.     double          K;      /* local copy of K  to return*/
  1652.     int             i;   /*loop index*/
  1653.     int             iterCounter; /*counter on iterations*/
  1654.     Int4            divisor; /*candidate divisor of all scores with
  1655.                                non-zero probabilities*/
  1656.     /*highest and lowest possible alignment scores for current length*/
  1657.     Int4            lowAlignmentScore, highAlignmentScore;
  1658.     Int4            first, last; /*loop indices for dynamic program*/
  1659.     register double innerSum;
  1660.     double          oldsum, oldsum2;  /* values of innerSum on previous
  1661.                                          iterations*/
  1662.     double          outerSum;        /* holds sum over j of (innerSum
  1663.                                         for iteration j/j)*/
  1665.     double          score_avg; /*average score*/
  1666.     /*first term to use in the closed form for the case where
  1667.       high == 1 or low == -1, but not both*/
  1668.     double          firstTermClosedForm;  /*usually store H/lambda*/
  1669.     int             iterlimit; /*upper limit on iterations*/
  1670.     double          sumlimit; /*lower limit on contributions
  1671.                                 to sum over scores*/
  1673.     /*array of score probabilities reindexed so that low is at index 0*/
  1674.     double         *probArrayStartLow;
  1676.     /*pointers used in dynamic program*/
  1677.     double         *ptrP, *ptr1, *ptr2, *ptr1e;
  1678.     double          expMinusLambda; /*e^^(-Lambda) */
  1680.     if (lambda <= 0. || H <= 0.) {
  1681.         /* Theory dictates that H and lambda must be positive, so
  1682.          * return -1 to indicate an error */
  1683.         return -1.;
  1684.     }
  1686.     /*Karlin-Altschul theory works only if the expected score
  1687.       is negative*/
  1688.     if (sfp->score_avg >= 0.0) {
  1689.         return -1.;
  1690.     }
  1692.     low   = sfp->obs_min;
  1693.     high  = sfp->obs_max;
  1694.     range = high - low;
  1696.     probArrayStartLow = &sfp->sprob[low];
  1697.     /* Look for the greatest common divisor ("delta" in Appendix of PNAS 87 of
  1698.        Karlin&Altschul (1990) */
  1699.     for (i = 1, divisor = -low; i <= range && divisor > 1; ++i) {
  1700.         if (probArrayStartLow[i])
  1701.             divisor = BLAST_Gcd(divisor, i);
  1702.     }
  1704.     high   /= divisor;
  1705.     low    /= divisor;
  1706.     lambda *= divisor;
  1708.     range = high - low;
  1710.     firstTermClosedForm = H/lambda;
  1711.     expMinusLambda      = exp((double) -lambda);
  1713.     if (low == -1 && high == 1) {
  1714.         K = (sfp->sprob[low] - sfp->sprob[high]) *
  1715.             (sfp->sprob[low] - sfp->sprob[high]) / sfp->sprob[low];
  1716.         return(K);
  1717.     }
  1719.     if (low == -1 || high == 1) {
  1720.         if (high == 1)
  1721.             ;
  1722.         else {
  1723.             score_avg = sfp->score_avg / divisor;
  1724.             firstTermClosedForm
  1725.                 = (score_avg * score_avg) / firstTermClosedForm;
  1726.         }
  1727.         return firstTermClosedForm * (1.0 - expMinusLambda);
  1728.     }
  1730.     sumlimit  = BLAST_KARLIN_K_SUMLIMIT_DEFAULT;
  1731.     iterlimit = BLAST_KARLIN_K_ITER_MAX;
  1733.     if (DIMOFP0 > DIMOFP0_MAX) {
  1734.         return -1.;
  1735.     }
  1736.     alignmentScoreProbabilities =
  1737.         (double *)calloc(DIMOFP0, sizeof(*alignmentScoreProbabilities));
  1738.     if (alignmentScoreProbabilities == NULL)
  1739.         return -1.;
  1741.     outerSum = 0.;
  1742.     lowAlignmentScore = highAlignmentScore = 0;
  1743.     alignmentScoreProbabilities[0] = innerSum = oldsum = oldsum2 = 1.;
  1745.     for (iterCounter = 0;
  1746.          ((iterCounter < iterlimit) && (innerSum > sumlimit));
  1747.          outerSum += innerSum /= ++iterCounter) {
  1748.         first = last = range;
  1749.         lowAlignmentScore  += low;
  1750.         highAlignmentScore += high;
  1751.         /*dynamic program to compute P(i,j)*/
  1752.         for (ptrP = alignmentScoreProbabilities +
  1753.                  (highAlignmentScore-lowAlignmentScore);
  1754.              ptrP >= alignmentScoreProbabilities;
  1755.              *ptrP-- =innerSum) {
  1756.             ptr1  = ptrP - first;
  1757.             ptr1e = ptrP - last;
  1758.             ptr2  = probArrayStartLow + first;
  1759.             for (innerSum = 0.; ptr1 >= ptr1e; )
  1760.                 innerSum += *ptr1--  *  *ptr2++;
  1761.             if (first)
  1762.                 --first;
  1763.             if (ptrP - alignmentScoreProbabilities <= range)
  1764.                 --last;
  1765.         }
  1766.         /* Horner's rule */
  1767.         innerSum = *++ptrP;
  1768.         for( i = lowAlignmentScore + 1; i < 0; i++ ) {
  1769.             innerSum = *++ptrP + innerSum * expMinusLambda;
  1770.         }
  1771.         innerSum *= expMinusLambda;
  1773.         for (; i <= highAlignmentScore; ++i)
  1774.             innerSum += *++ptrP;
  1775.         oldsum2 = oldsum;
  1776.         oldsum  = innerSum;
  1777.     }
  1779. #ifdef ADD_GEOMETRIC_TERMS_TO_K
  1780.     /*old code assumed that the later terms in sum were
  1781.       asymptotically comparable to those of a geometric
  1782.       progression, and tried to speed up convergence by
  1783.       guessing the estimated ratio between sucessive terms
  1784.       and using the explicit terms of a geometric progression
  1785.       to speed up convergence. However, the assumption does not
  1786.       always hold, and convergenece of the above code is fast
  1787.       enough in practice*/
  1788.     /* Terms of geometric progression added for correction */
  1789.     {
  1790.         double     ratio;  /* fraction used to generate the
  1791.                                    geometric progression */
  1793.         ratio = oldsum / oldsum2;
  1794.         if (ratio >= (1.0 - sumlimit*0.001)) {
  1795.             K = -1.;
  1796.             if (alignmentScoreProbabilities != NULL)
  1797.                 sfree(alignmentScoreProbabilities);
  1798.             return K;
  1799.         }
  1800.         sumlimit *= 0.01;
  1801.         while (innerSum > sumlimit) {
  1802.             oldsum   *= ratio;
  1803.             outerSum += innerSum = oldsum / ++iterCounter;
  1804.         }
  1805.     }
  1806. #endif
  1808.     if (expMinusLambda <  SMALL_LAMBDA_THRESHOLD ) {
  1809.         K = -exp((double)-2.0*outerSum) /
  1810.             (firstTermClosedForm*(expMinusLambda - 1.0));
  1811.     } else {
  1812.         K = -exp((double)-2.0*outerSum) /
  1813.             (firstTermClosedForm*BLAST_Expm1(-(double)lambda));
  1814.     }
  1816.     if (alignmentScoreProbabilities != NULL)
  1817.         sfree(alignmentScoreProbabilities);
  1819.     return K;
  1820. }
  1823. /**
  1824.  * Find positive solution to sum_{i=low}^{high} exp(i lambda) = 1.
  1825.  * 
  1826.  * @param probs probabilities of a score occurring 
  1827.  * @param d the gcd of the possible scores. This equals 1 if the scores
  1828.  * are not a lattice
  1829.  * @param low the lowest possible score
  1830.  * @param high the highest possible score
  1831.  * @param lambda0 an initial value for lambda
  1832.  * @param tolx the tolerance to which lambda must be computed
  1833.  * @param itmax the maximum number of times the function may be
  1834.  * evaluated
  1835.  * @param maxNewton the maximum permissible number of Newton
  1836.  * iteration. After that the computation will proceed by bisection.
  1837.  * @param itn a pointer to an integer that will receive the actually
  1838.  * number of iterations performed.
  1839.  *
  1840.  * Let phi(lambda) =  sum_{i=low}^{high} exp(i lambda) - 1. Then phi(lambda)
  1841.  * may be written
  1842.  *
  1843.  *     phi(lamdba) = exp(u lambda) p( exp(-lambda) )
  1844.  *
  1845.  * where p(x) is a polynomial that has exactly two zeros, one at x = 1
  1846.  * and one at y = exp(-lamdba). It is simpler, more numerically
  1847.  * efficient and stable to apply Newton's method to p(x) than to
  1848.  * phi(lambda).
  1849.  *
  1850.  * We define a safeguarded Newton iteration as follows. Let the
  1851.  * initial interval of uncertainty be [0,1]. If p'(x) >= 0, we bisect
  1852.  * the interval. Otherwise we try a Newton step. If the Newton iterate
  1853.  * lies in the current interval of uncertainty and it reduces the
  1854.  * value of | p(x) | by at least 10%, we accept the new
  1855.  * point. Otherwise, we bisect the current interval of uncertainty.
  1856.  * It is clear that this method converges to a zero of p(x).  Since
  1857.  * p'(x) > 0 in an interval containing x = 1, the method cannot
  1858.  * converge to x = 1 and therefore converges to the only other zero,
  1859.  * y.
  1860.  */
  1862. static double 
  1863. NlmKarlinLambdaNR( double* probs, Int4 d, Int4 low, Int4 high, double lambda0, double tolx,
  1864.  Int4 itmax, Int4 maxNewton, Int4 * itn ) 
  1865. {
  1866.   Int4 k;
  1867.   double x0, x, a = 0, b = 1;
  1868.   double f = 4;  /* Larger than any possible value of the poly in [0,1] */
  1869.   Int4 isNewton = 0; /* we haven't yet taken a Newton step. */
  1871.   assert( d > 0 );
  1873. x0 = exp( -lambda0 );
  1874.   x = ( 0 < x0 && x0 < 1 ) ? x0 : .5;
  1875.   
  1876.   for( k = 0; k < itmax; k++ ) { /* all iteration indices k */
  1877.     Int4 i;
  1878.     double g, fold = f;
  1879.     Int4 wasNewton = isNewton; /* If true, then the previous step was a */
  1880.                               /* Newton step */
  1881.     isNewton  = 0;            /* Assume that this step is not */
  1882.     
  1883.     /* Horner's rule for evaluating a polynomial and its derivative */
  1884.     g = 0;
  1885.     f = probs[low];
  1886.     for( i = low + d; i < 0; i += d ) {
  1887.       g = x * g + f;
  1888.       f = f * x + probs[i];
  1889.     }
  1890.     g = x * g + f;
  1891.     f = f * x + probs[0] - 1;
  1892.     for( i = d; i <= high; i += d ) {
  1893.       g = x * g + f;
  1894.       f = f * x + probs[i];
  1895.     }
  1896.     /* End Horner's rule */
  1898.     if( f > 0 ) {
  1899.       a = x; /* move the left endpoint */
  1900.     } else if( f < 0 ) { 
  1901.       b = x; /* move the right endpoint */
  1902.     } else { /* f == 0 */
  1903.       break; /* x is an exact solution */
  1904.     }
  1905.     if( b - a < 2 * a * ( 1 - b ) * tolx ) {
  1906.       /* The midpoint of the interval converged */
  1907.       x = (a + b) / 2; break;
  1908.     }
  1910.     if( k >= maxNewton ||
  1911.         /* If convergence of Newton's method appears to be failing; or */
  1912. ( wasNewton && fabs( f ) > .9 * fabs(fold) ) ||  
  1913.         /* if the previous iteration was a Newton step but didn't decrease 
  1914.          * f sufficiently; or */
  1915.         g >= 0 
  1916.         /* if a Newton step will move us away from the desired solution */
  1917.         ) { /* then */
  1918.       /* bisect */
  1919.       x = (a + b)/2;
  1920.     } else {
  1921.       /* try a Newton step */
  1922.       double p = - f/g;
  1923.       double y = x + p;
  1924.       if( y <= a || y >= b ) { /* The proposed iterate is not in (a,b) */
  1925.         x = (a + b)/2;
  1926.       } else { /* The proposed iterate is in (a,b). Accept it. */
  1927.         isNewton = 1;
  1928.         x = y;
  1929.         if( fabs( p ) < tolx * x * (1-x) ) break; /* Converged */
  1930.       } /* else the proposed iterate is in (a,b) */
  1931.     } /* else try a Newton step. */ 
  1932.   } /* end for all iteration indices k */
  1933. *itn = k; 
  1934.   return -log(x)/d;
  1935. }
  1938. double
  1939. Blast_KarlinLambdaNR(Blast_ScoreFreq* sfp, double initialLambdaGuess)
  1940. {
  1941. Int4 low; /* Lowest score (must be negative)  */
  1942. Int4 high; /* Highest score (must be positive) */
  1943. Int4 itn;
  1944. Int4 i, d;
  1945. double*  sprob;
  1946. double returnValue;
  1948. low = sfp->obs_min;
  1949. high = sfp->obs_max;
  1950. if (sfp->score_avg >= 0.) { /* Expected score must be negative */
  1951. return -1.0;
  1952. }
  1953. if (BlastScoreChk(low, high) != 0) return -1.;
  1955. sprob = sfp->sprob;
  1956. /* Find greatest common divisor of all scores */
  1957. for (i = 1, d = -low; i <= high-low && d > 1; ++i) {
  1958. if (sprob[i+low] != 0) {
  1959. d = BLAST_Gcd(d, i);
  1960. }
  1961. }
  1962. returnValue =
  1963. NlmKarlinLambdaNR( sprob, d, low, high,
  1964.                            initialLambdaGuess,
  1965.                            BLAST_KARLIN_LAMBDA_ACCURACY_DEFAULT,
  1966.    20, 20 + BLAST_KARLIN_LAMBDA_ITER_DEFAULT, &itn );
  1969. return returnValue;
  1970. }
  1972. /*
  1973. BlastKarlinLtoH
  1975. Calculate H, the relative entropy of the p's and q's
  1976. */
  1977. static double 
  1978. BlastKarlinLtoH(Blast_ScoreFreq* sfp, double lambda)
  1979. {
  1980. Int4 score;
  1981. double H, etonlam, sum, scale;
  1983. double *probs = sfp->sprob;
  1984. Int4 low   = sfp->obs_min,  high  = sfp->obs_max;
  1986. if (lambda < 0.) {
  1987. return -1.;
  1988. }
  1989. if (BlastScoreChk(low, high) != 0) return -1.;
  1991. etonlam = exp( - lambda );
  1992.   sum = low * probs[low];
  1993.   for( score = low + 1; score <= high; score++ ) {
  1994.     sum = score * probs[score] + etonlam * sum;
  1995.   }
  1997.   scale = BLAST_Powi( etonlam, high );
  1998.   if( scale > 0 ) {
  1999.     H = lambda * sum/scale;
  2000.   } else { /* Underflow of exp( -lambda * high ) */
  2001.     H = lambda * exp( lambda * high + log(sum) );
  2002.   }
  2003. return H;
  2004. }
  2006. /**************** Statistical Significance Parameter Subroutine ****************
  2008.     Version 1.0     February 2, 1990
  2009.     Version 1.2     July 6,     1990
  2011.     Program by:     Stephen Altschul
  2013.     Address:        National Center for Biotechnology Information
  2014.                     National Library of Medicine
  2015.                     National Institutes of Health
  2016.                     Bethesda, MD  20894
  2018.     Internet:       altschul@ncbi.nlm.nih.gov
  2020. See:  Karlin, S. & Altschul, S.F. "Methods for Assessing the Statistical
  2021.     Significance of Molecular Sequence Features by Using General Scoring
  2022.     Schemes,"  Proc. Natl. Acad. Sci. USA 87 (1990), 2264-2268.
  2024.     Computes the parameters lambda and K for use in calculating the
  2025.     statistical significance of high-scoring segments or subalignments.
  2027.     The scoring scheme must be integer valued.  A positive score must be
  2028.     possible, but the expected (mean) score must be negative.
  2030.     A program that calls this routine must provide the value of the lowest
  2031.     possible score, the value of the greatest possible score, and a pointer
  2032.     to an array of probabilities for the occurrence of all scores between
  2033.     these two extreme scores.  For example, if score -2 occurs with
  2034.     probability 0.7, score 0 occurs with probability 0.1, and score 3
  2035.     occurs with probability 0.2, then the subroutine must be called with
  2036.     low = -2, high = 3, and pr pointing to the array of values
  2037.     { 0.7, 0.0, 0.1, 0.0, 0.0, 0.2 }.  The calling program must also provide
  2038.     pointers to lambda and K; the subroutine will then calculate the values
  2039.     of these two parameters.  In this example, lambda=0.330 and K=0.154.
  2041.     The parameters lambda and K can be used as follows.  Suppose we are
  2042.     given a length N random sequence of independent letters.  Associated
  2043.     with each letter is a score, and the probabilities of the letters
  2044.     determine the probability for each score.  Let S be the aggregate score
  2045.     of the highest scoring contiguous segment of this sequence.  Then if N
  2046.     is sufficiently large (greater than 100), the following bound on the
  2047.     probability that S is greater than or equal to x applies:
  2049.             P( S >= x )   <=   1 - exp [ - KN exp ( - lambda * x ) ].
  2051.     In other words, the p-value for this segment can be written as
  2052.     1-exp[-KN*exp(-lambda*S)].
  2054.     This formula can be applied to pairwise sequence comparison by assigning
  2055.     scores to pairs of letters (e.g. amino acids), and by replacing N in the
  2056.     formula with N*M, where N and M are the lengths of the two sequences
  2057.     being compared.
  2059.     In addition, letting y = KN*exp(-lambda*S), the p-value for finding m
  2060.     distinct segments all with score >= S is given by:
  2062.                            2             m-1           -y
  2063.             1 - [ 1 + y + y /2! + ... + y   /(m-1)! ] e
  2065.     Notice that for m=1 this formula reduces to 1-exp(-y), which is the same
  2066.     as the previous formula.
  2068. *******************************************************************************/
  2069. Int2
  2070. Blast_KarlinBlkCalc(Blast_KarlinBlk* kbp, Blast_ScoreFreq* sfp)
  2071. {
  2074. if (kbp == NULL || sfp == NULL)
  2075. return 1;
  2077. /* Calculate the parameter Lambda */
  2079. kbp->Lambda = Blast_KarlinLambdaNR(sfp, BLAST_KARLIN_LAMBDA0_DEFAULT);
  2080. if (kbp->Lambda < 0.)
  2081. goto ErrExit;
  2084. /* Calculate H */
  2086. kbp->H = BlastKarlinLtoH(sfp, kbp->Lambda);
  2087. if (kbp->H < 0.)
  2088. goto ErrExit;
  2091. /* Calculate K and log(K) */
  2093. kbp->K = BlastKarlinLHtoK(sfp, kbp->Lambda, kbp->H);
  2094. if (kbp->K < 0.)
  2095. goto ErrExit;
  2096. kbp->logK = log(kbp->K);
  2098. /* Normal return */
  2099. return 0;
  2101. ErrExit:
  2102. kbp->Lambda = kbp->H = kbp->K = -1.;
  2103. #ifdef BLASTKAR_HUGE_VAL
  2104. kbp->logK = BLASTKAR_HUGE_VAL;
  2105. #else
  2106. kbp->logK = 1.e30;
  2107. #endif
  2108. return 1;
  2109. }
  2111. /*
  2112.   Calculate the Karlin parameters.  This function should be called once
  2113.   for each context, or frame translated.
  2114.   
  2115.   The rfp and stdrfp are calculated for each context, this should be
  2116.   fixed. 
  2117. */
  2119. Int2
  2120. BLAST_ScoreBlkFill(BlastScoreBlk* sbp, char* query, Int4 query_length, Int4 context_number)
  2121. {
  2122. Blast_ResFreq* rfp,* stdrfp;
  2123. Int2 retval=0;
  2125. rfp = Blast_ResFreqNew(sbp);
  2126. stdrfp = Blast_ResFreqNew(sbp);
  2127. Blast_ResFreqStdComp(sbp, stdrfp);
  2128. Blast_ResFreqString(sbp, rfp, query, query_length);
  2129. sbp->sfp[context_number] = Blast_ScoreFreqNew(sbp->loscore, sbp->hiscore);
  2130. BlastScoreFreqCalc(sbp, sbp->sfp[context_number], rfp, stdrfp);
  2131. sbp->kbp_std[context_number] = Blast_KarlinBlkCreate();
  2132. retval = Blast_KarlinBlkCalc(sbp->kbp_std[context_number], sbp->sfp[context_number]);
  2133. if (retval)
  2134. {
  2135. rfp = Blast_ResFreqDestruct(rfp);
  2136. stdrfp = Blast_ResFreqDestruct(stdrfp);
  2137. return retval;
  2138. }
  2139. sbp->kbp_psi[context_number] = Blast_KarlinBlkCreate();
  2140. retval = Blast_KarlinBlkCalc(sbp->kbp_psi[context_number], sbp->sfp[context_number]);
  2141. rfp = Blast_ResFreqDestruct(rfp);
  2142. stdrfp = Blast_ResFreqDestruct(stdrfp);
  2144. return retval;
  2145. }
  2147. /*
  2148. Calculates the standard Karlin parameters.  This is used
  2149. if the query is translated and the calculated (real) Karlin
  2150. parameters are bad, as they're calculated for non-coding regions.
  2151. */
  2153. Blast_KarlinBlk*
  2154. Blast_KarlinBlkIdealCalc(BlastScoreBlk* sbp)
  2156. {
  2157. Blast_KarlinBlk* kbp_ideal;
  2158. Blast_ResFreq* stdrfp;
  2159. Blast_ScoreFreq* sfp;
  2161. stdrfp = Blast_ResFreqNew(sbp);
  2162. Blast_ResFreqStdComp(sbp, stdrfp);
  2163. sfp = Blast_ScoreFreqNew(sbp->loscore, sbp->hiscore);
  2164. BlastScoreFreqCalc(sbp, sfp, stdrfp, stdrfp);
  2165. kbp_ideal = Blast_KarlinBlkCreate();
  2166. Blast_KarlinBlkCalc(kbp_ideal, sfp);
  2167. stdrfp = Blast_ResFreqDestruct(stdrfp);
  2169. sfp = Blast_ScoreFreqDestruct(sfp);
  2171. return kbp_ideal;
  2172. }
  2174. Int2
  2175. Blast_KarlinBlkStandardCalc(BlastScoreBlk* sbp, Int4 context_start, Int4 context_end)
  2176. {
  2178. Blast_KarlinBlk* kbp_ideal,* kbp;
  2179. Int4 index;
  2181. kbp_ideal = Blast_KarlinBlkIdealCalc(sbp);
  2182. /* Replace the calculated values with ideal ones for blastx, tblastx. */
  2183. for (index=context_start; index<=context_end; index++)
  2184. {
  2185. kbp = sbp->kbp[index];
  2186.       if (!kbp)
  2187.          continue;
  2188. if (kbp->Lambda >= kbp_ideal->Lambda)
  2189. {
  2190. kbp->Lambda = kbp_ideal->Lambda;
  2191. kbp->K = kbp_ideal->K;
  2192. kbp->logK = kbp_ideal->logK;
  2193. kbp->H = kbp_ideal->H;
  2194. }
  2195. }
  2196. kbp_ideal = Blast_KarlinBlkDestruct(kbp_ideal);
  2198. return 0;
  2199. }
  2201. /*
  2202. Creates the Karlin Block.
  2203. */
  2205. Blast_KarlinBlk*
  2206. Blast_KarlinBlkCreate(void)
  2208. {
  2209. Blast_KarlinBlk* kbp;
  2211. kbp = (Blast_KarlinBlk*) calloc(1, sizeof(Blast_KarlinBlk));
  2213. return kbp;
  2214. }
  2216. static SBLASTMatrixStructure*
  2217. BlastMatrixAllocate(Int2 alphabet_size)
  2219. {
  2220. SBLASTMatrixStructure* matrix_struct;
  2221. Int2 index;
  2223. if (alphabet_size <= 0 || alphabet_size >= BLAST_MATRIX_SIZE)
  2224. return NULL;
  2226. matrix_struct = (SBLASTMatrixStructure*) calloc(1, sizeof(SBLASTMatrixStructure));
  2228. if (matrix_struct == NULL)
  2229. return NULL;
  2231. for (index=0; index<BLAST_MATRIX_SIZE-1; index++)
  2232. {
  2233. matrix_struct->matrix[index] = matrix_struct->long_matrix + index*BLAST_MATRIX_SIZE;
  2234. }
  2236. return matrix_struct;
  2237. }
  2239. /*
  2240. Deallocates MatrixInfo*
  2241. */
  2243. static MatrixInfo*
  2244. MatrixInfoDestruct(MatrixInfo* matrix_info)
  2246. {
  2247. if (matrix_info == NULL)
  2248. return NULL;
  2250. sfree(matrix_info->name);
  2251. sfree(matrix_info);
  2252. return NULL;
  2253. }
  2255. /*
  2256. Makes New MatrixInfo*
  2257. */
  2259. static MatrixInfo*
  2260. MatrixInfoNew(char* name, array_of_8 *values, Int4* prefs, Int4 max_number)
  2262. {
  2263. MatrixInfo* matrix_info;
  2265. matrix_info = (MatrixInfo*) calloc(1, sizeof(MatrixInfo));
  2266. matrix_info->name = strdup(name);
  2267. matrix_info->values = values;
  2268. matrix_info->prefs = prefs;
  2269. matrix_info->max_number_values = max_number;
  2271. return matrix_info;
  2272. }
  2274. static ListNode*
  2275. BlastMatrixValuesDestruct(ListNode* vnp)
  2277. {
  2278. ListNode* head;
  2280. head = vnp;
  2281. while (vnp)
  2282. {
  2283. MatrixInfoDestruct((MatrixInfo*) vnp->ptr);
  2284. vnp = vnp->next;
  2285. }
  2287. head = ListNodeFree(head);
  2289. return head;
  2290. }
  2292. /*
  2293. ListNode* BlastLoadMatrixValues (void)
  2295. Loads all the matrix values, returns a ListNode* chain that contains
  2296. MatrixInfo*'s.
  2298. */
  2299. static ListNode* 
  2300. BlastLoadMatrixValues (void)
  2302. {
  2303. MatrixInfo* matrix_info;
  2304. ListNode* retval=NULL;
  2306. matrix_info = MatrixInfoNew("BLOSUM80", blosum80_values, blosum80_prefs, BLOSUM80_VALUES_MAX);
  2307. ListNodeAddPointer(&retval, 0, matrix_info);
  2309. matrix_info = MatrixInfoNew("BLOSUM62", blosum62_values, blosum62_prefs, BLOSUM62_VALUES_MAX);
  2310. ListNodeAddPointer(&retval, 0, matrix_info);
  2312. matrix_info = MatrixInfoNew("BLOSUM50", blosum50_values, blosum50_prefs, BLOSUM50_VALUES_MAX);
  2313. ListNodeAddPointer(&retval, 0, matrix_info);
  2315. matrix_info = MatrixInfoNew("BLOSUM45", blosum45_values, blosum45_prefs, BLOSUM45_VALUES_MAX);
  2316. ListNodeAddPointer(&retval, 0, matrix_info);
  2318. matrix_info = MatrixInfoNew("PAM250", pam250_values, pam250_prefs, PAM250_VALUES_MAX);
  2319. ListNodeAddPointer(&retval, 0, matrix_info);
  2321. matrix_info = MatrixInfoNew("BLOSUM62_20", blosum62_20_values, blosum62_20_prefs, BLOSUM62_20_VALUES_MAX);
  2322. ListNodeAddPointer(&retval, 0, matrix_info);
  2324. matrix_info = MatrixInfoNew("BLOSUM90", blosum90_values, blosum90_prefs, BLOSUM90_VALUES_MAX);
  2325. ListNodeAddPointer(&retval, 0, matrix_info);
  2327. matrix_info = MatrixInfoNew("PAM30", pam30_values, pam30_prefs, PAM30_VALUES_MAX);
  2328. ListNodeAddPointer(&retval, 0, matrix_info);
  2330. matrix_info = MatrixInfoNew("PAM70", pam70_values, pam70_prefs, PAM70_VALUES_MAX);
  2331. ListNodeAddPointer(&retval, 0, matrix_info);
  2333. return retval;
  2334. }
  2336. /*
  2337. Int2
  2338. BlastKarlinGetMatrixValuesEx2(char* matrix, Int4* open, Int4* extension, Int4* decline_align, double* lambda, double* K, double* H)
  2340. Obtains arrays of the allowed opening and extension penalties for gapped BLAST for
  2341. the given matrix.  Also obtains arrays of Lambda, K, and H.  Any of these fields that
  2342. are not required should be set to NULL.  The Int2 return value is the length of the
  2343. arrays.
  2344. */
  2346. static Int2
  2347. BlastKarlinGetMatrixValuesEx2(char* matrix, Int4** open, Int4** extension, Int4** decline_align, double** lambda, double** K, double** H, double** alpha, double** beta, Int4** pref_flags)
  2349. {
  2350. array_of_8 *values = NULL;
  2351. Boolean found_matrix=FALSE;
  2352. Int4 index, max_number_values=0;
  2353. Int4* open_array=NULL,* extension_array=NULL,* decline_align_array=NULL,* pref_flags_array=NULL,* prefs=NULL;
  2354. double* lambda_array=NULL,* K_array=NULL,* H_array=NULL,* alpha_array=NULL,* beta_array=NULL;
  2355. MatrixInfo* matrix_info;
  2356. ListNode* vnp,* head;
  2358. if (matrix == NULL)
  2359. return 0;
  2361. vnp = head = BlastLoadMatrixValues();
  2363. while (vnp)
  2364. {
  2365. matrix_info = vnp->ptr;
  2366. if (strcasecmp(matrix_info->name, matrix) == 0)
  2367. {
  2368. values = matrix_info->values;
  2369. max_number_values = matrix_info->max_number_values;
  2370. prefs = matrix_info->prefs;
  2371. found_matrix = TRUE;
  2372. break;
  2373. }
  2374. vnp = vnp->next;
  2375. }
  2377. if (found_matrix)
  2378. {
  2379. if (open)
  2380. *open = open_array = (Int4 *) calloc(max_number_values, sizeof(Int4));
  2381. if (extension)
  2382. *extension = extension_array = 
  2383.             (Int4 *) calloc(max_number_values, sizeof(Int4));
  2384. if (decline_align)
  2385. *decline_align = decline_align_array = 
  2386.             (Int4 *) calloc(max_number_values, sizeof(Int4));
  2387. if (lambda)
  2388. *lambda = lambda_array = 
  2389.             (double*) calloc(max_number_values, sizeof(double));
  2390. if (K)
  2391. *K = K_array = (double*) calloc(max_number_values, sizeof(double));
  2392. if (H)
  2393. *H = H_array = (double*) calloc(max_number_values, sizeof(double));
  2394. if (alpha)
  2395. *alpha = alpha_array = (double*) calloc(max_number_values, sizeof(double));
  2396. if (beta)
  2397. *beta = beta_array = (double*) calloc(max_number_values, sizeof(double));
  2398. if (pref_flags)
  2399. *pref_flags = pref_flags_array = 
  2400.             (Int4 *) calloc(max_number_values, sizeof(Int4));
  2402. for (index=0; index<max_number_values; index++)
  2403. {
  2404. if (open)
  2405. open_array[index] = (Int4) values[index][0];
  2406. if (extension)
  2407. extension_array[index] = (Int4) values[index][1];
  2408. if (decline_align)
  2409. decline_align_array[index] = (Int4) values[index][2];
  2410. if (lambda)
  2411. lambda_array[index] = values[index][3];
  2412. if (K)
  2413. K_array[index] = values[index][4];
  2414. if (H)
  2415. H_array[index] = values[index][5];
  2416. if (alpha)
  2417. alpha_array[index] = values[index][6];
  2418. if (beta)
  2419. beta_array[index] = values[index][7];
  2420. if (pref_flags)
  2421. pref_flags_array[index] = prefs[index];
  2422. }
  2423. }
  2425. BlastMatrixValuesDestruct(head);
  2427. return max_number_values;
  2428. }
  2430. /*Extract the alpha and beta settings for this matrixName, and these
  2431.   gap open and gap extension costs*/
  2432. void BLAST_GetAlphaBeta(char* matrixName, double *alpha,
  2433. double *beta, Boolean gapped, Int4 gap_open, Int4 gap_extend)
  2434. {
  2435.    Int4* gapOpen_arr,* gapExtend_arr,* pref_flags;
  2436.    double* alpha_arr,* beta_arr;
  2437.    Int2 num_values;
  2438.    Int4 i; /*loop index*/
  2440.    num_values = BlastKarlinGetMatrixValuesEx2(matrixName, &gapOpen_arr, 
  2441.      &gapExtend_arr, NULL, NULL, NULL, NULL,  &alpha_arr, &beta_arr, 
  2442.      &pref_flags);
  2444.    if (gapped) {
  2445.      if ((0 == gap_open) && (0 == gap_extend)) {
  2446.        for(i = 1; i < num_values; i++) {
  2447.  if(pref_flags[i]==BLAST_MATRIX_BEST) {
  2448.    (*alpha) = alpha_arr[i];
  2449.    (*beta) = beta_arr[i];
  2450.    break;
  2451.  }
  2452.        }
  2453.      }
  2454.      else {
  2455.        for(i = 1; i < num_values; i++) {
  2456.  if ((gapOpen_arr[i] == gap_open) &&
  2457.      (gapExtend_arr[i] == gap_extend)) {
  2458.    (*alpha) = alpha_arr[i];
  2459.    (*beta) = beta_arr[i];
  2460.    break;
  2461.  }
  2462.        }
  2463.      }
  2464.    }
  2465.    else if (num_values > 0) {
  2466.      (*alpha) = alpha_arr[0];
  2467.      (*beta) = beta_arr[0];
  2468.    }
  2470.    sfree(gapOpen_arr);
  2471.    sfree(gapExtend_arr);
  2472.    sfree(pref_flags);
  2473.    sfree(alpha_arr);
  2474.    sfree(beta_arr);
  2475. }
  2477. static Int2
  2478. BlastKarlinReportAllowedValues(const char *matrix_name, 
  2479.    Blast_Message** error_return)
  2480. {
  2481. array_of_8 *values = NULL;
  2482. Boolean found_matrix=FALSE;
  2483. char buffer[256];
  2484. Int4 index, max_number_values=0;
  2485. MatrixInfo* matrix_info;
  2486. ListNode* vnp,* head;
  2488. vnp = head = BlastLoadMatrixValues();
  2489. while (vnp)
  2490. {
  2491. matrix_info = vnp->ptr;
  2492. if (strcasecmp(matrix_info->name, matrix_name) == 0)
  2493. {
  2494. values = matrix_info->values;
  2495. max_number_values = matrix_info->max_number_values;
  2496. found_matrix = TRUE;
  2497. break;
  2498. }
  2499. vnp = vnp->next;
  2500. }
  2502. if (found_matrix)
  2503. {
  2504. for (index=0; index<max_number_values; index++)
  2505. {
  2506. if (BLAST_Nint(values[index][2]) == INT2_MAX)
  2507. sprintf(buffer, "Gap existence and extension values of %ld and %ld are supported", (long) BLAST_Nint(values[index][0]), (long) BLAST_Nint(values[index][1]));
  2508. else
  2509. sprintf(buffer, "Gap existence, extension and decline-to-align values of %ld, %ld and %ld are supported", (long) BLAST_Nint(values[index][0]), (long) BLAST_Nint(values[index][1]), (long) BLAST_Nint(values[index][2]));
  2510. Blast_MessageWrite(error_return, BLAST_SEV_ERROR, 0, 0, buffer);
  2511. }
  2512. }
  2514. BlastMatrixValuesDestruct(head);
  2516. return 0;
  2517. }
  2519. /*
  2520. Supplies lambda, H, and K values, as calcualted by Stephen Altschul 
  2521. in Methods in Enzy. (vol 266, page 474).
  2523. if kbp is NULL, then a validation is perfomed.
  2524. */
  2525. Int2
  2526. Blast_KarlinBlkGappedCalc(Blast_KarlinBlk* kbp, Int4 gap_open, Int4 gap_extend, Int4 decline_align, char* matrix_name, Blast_Message** error_return)
  2528. {
  2529. char buffer[256];
  2530. Int2 status = Blast_KarlinkGapBlkFill(kbp, gap_open, gap_extend, decline_align, matrix_name);
  2532. if (status && error_return)
  2533. {
  2534. if (status == 1)
  2535. {
  2536. MatrixInfo* matrix_info;
  2537. ListNode* vnp,* head; 
  2539. vnp = head = BlastLoadMatrixValues();
  2541. sprintf(buffer, "%s is not a supported matrix", matrix_name);
  2542. Blast_MessageWrite(error_return, BLAST_SEV_ERROR, 0, 0, buffer);
  2544. while (vnp)
  2545. {
  2546. matrix_info = vnp->ptr;
  2547. sprintf(buffer, "%s is a supported matrix", matrix_info->name);
  2548.             Blast_MessageWrite(error_return, BLAST_SEV_ERROR, 0, 0, buffer);
  2549. vnp = vnp->next;
  2550. }
  2552. BlastMatrixValuesDestruct(head);
  2553. }
  2554. else if (status == 2)
  2555. {
  2556. if (decline_align == INT2_MAX)
  2557. sprintf(buffer, "Gap existence and extension values of %ld and %ld not supported for %s", (long) gap_open, (long) gap_extend, matrix_name);
  2558. else
  2559. 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);
  2560. Blast_MessageWrite(error_return, BLAST_SEV_ERROR, 0, 0, buffer);
  2561. BlastKarlinReportAllowedValues(matrix_name, error_return);
  2562. }
  2563. }
  2565. return status;
  2566. }
  2568. /*
  2569. Attempts to fill KarlinBlk for given gap opening, extensions etc.  
  2570. Will return non-zero status if that fails.
  2572. return values: -1 if matrix_name is NULL;
  2573. 1 if matrix not found
  2574. 2 if matrix found, but open, extend etc. values not supported.
  2575. */
  2576. Int2
  2577. Blast_KarlinkGapBlkFill(Blast_KarlinBlk* kbp, Int4 gap_open, Int4 gap_extend, Int4 decline_align, char* matrix_name)
  2578. {
  2579. Boolean found_matrix=FALSE, found_values=FALSE;
  2580. array_of_8 *values;
  2581. Int2 status=0;
  2582. Int4 index, max_number_values=0;
  2583. MatrixInfo* matrix_info;
  2584. ListNode* vnp,* head;
  2586. if (matrix_name == NULL)
  2587. return -1;
  2589. values = NULL;
  2591. vnp = head = BlastLoadMatrixValues();
  2592. while (vnp)
  2593. {
  2594. matrix_info = vnp->ptr;
  2595. if (strcasecmp(matrix_info->name, matrix_name) == 0)
  2596. {
  2597. values = matrix_info->values;
  2598. max_number_values = matrix_info->max_number_values;
  2599. found_matrix = TRUE;
  2600. break;
  2601. }
  2602. vnp = vnp->next;
  2603. }
  2606. if (found_matrix)
  2607. {
  2608. for (index=0; index<max_number_values; index++)
  2609. {
  2610. if (BLAST_Nint(values[index][0]) == gap_open &&
  2611. BLAST_Nint(values[index][1]) == gap_extend &&
  2612. (BLAST_Nint(values[index][2]) == INT2_MAX || BLAST_Nint(values[index][2]) == decline_align))
  2613. {
  2614. if (kbp)
  2615. {
  2616. kbp->Lambda = values[index][3];
  2617. kbp->K = values[index][4];
  2618. kbp->logK = log(kbp->K);
  2619. kbp->H = values[index][5];
  2620. }
  2621. found_values = TRUE;
  2622. break;
  2623. }
  2624. }
  2626. if (found_values == TRUE)
  2627. {
  2628. status = 0;
  2629. }
  2630. else
  2631. {
  2632. status = 2;
  2633. }
  2634. }
  2635. else
  2636. {
  2637. status = 1;
  2638. }
  2640. BlastMatrixValuesDestruct(head);
  2642. return status;
  2643. }
  2645. char*
  2646. BLAST_PrintMatrixMessage(const char *matrix_name)
  2647. {
  2648. char* buffer= (char *) calloc(1024, sizeof(char));
  2649. char* ptr;
  2650. MatrixInfo* matrix_info;
  2651.         ListNode* vnp,* head;
  2653. ptr = buffer;
  2654.         sprintf(ptr, "%s is not a supported matrix, supported matrices are:n", matrix_name);
  2656. ptr += strlen(ptr);
  2658.         vnp = head = BlastLoadMatrixValues();
  2660.         while (vnp)
  2661.         {
  2662.          matrix_info = vnp->ptr;
  2663.          sprintf(ptr, "%s n", matrix_info->name);
  2664. ptr += strlen(ptr);
  2665. vnp = vnp->next;
  2666.         }
  2667.         BlastMatrixValuesDestruct(head);
  2669. return buffer;
  2670. }
  2672. char*
  2673. BLAST_PrintAllowedValues(const char *matrix_name, Int4 gap_open, Int4 gap_extend, Int4 decline_align) 
  2674. {
  2675. array_of_8 *values = NULL;
  2676. Boolean found_matrix=FALSE;
  2677. char* buffer,* ptr;
  2678. Int4 index, max_number_values=0;
  2679. MatrixInfo* matrix_info;
  2680. ListNode* vnp,* head;
  2682. ptr = buffer = (char *) calloc(2048, sizeof(char));
  2684.         if (decline_align == INT2_MAX)
  2685.               sprintf(ptr, "Gap existence and extension values of %ld and %ld not supported for %snsupported values are:n", 
  2686. (long) gap_open, (long) gap_extend, matrix_name);
  2687.         else
  2688.                sprintf(ptr, "Gap existence, extension and decline-to-align values of %ld, %ld and %ld not supported for %sn supported values are:n", 
  2689. (long) gap_open, (long) gap_extend, (long) decline_align, matrix_name);
  2691. ptr += strlen(ptr);
  2693. vnp = head = BlastLoadMatrixValues();
  2694. while (vnp)
  2695. {
  2696. matrix_info = vnp->ptr;
  2697. if (strcasecmp(matrix_info->name, matrix_name) == 0)
  2698. {
  2699. values = matrix_info->values;
  2700. max_number_values = matrix_info->max_number_values;
  2701. found_matrix = TRUE;
  2702. break;
  2703. }
  2704. vnp = vnp->next;
  2705. }
  2707. if (found_matrix)
  2708. {
  2709. for (index=0; index<max_number_values; index++)
  2710. {
  2711. if (BLAST_Nint(values[index][2]) == INT2_MAX)
  2712. sprintf(ptr, "%ld, %ldn", (long) BLAST_Nint(values[index][0]), (long) BLAST_Nint(values[index][1]));
  2713. else
  2714. sprintf(ptr, "%ld, %ld, %ldn", (long) BLAST_Nint(values[index][0]), (long) BLAST_Nint(values[index][1]), (long) BLAST_Nint(values[index][2]));
  2715. ptr += strlen(ptr);
  2716. }
  2717. }
  2719. BlastMatrixValuesDestruct(head);
  2721. return buffer;
  2722. }
  2725. /* Smallest float that might not cause a floating point exception in
  2726. S = (Int4) (ceil( log((double)(K * searchsp / E)) / Lambda ));
  2727. below.
  2728. */
  2729. #define BLASTKAR_SMALL_FLOAT 1.0e-297
  2730. static Int4
  2731. BlastKarlinEtoS_simple(double E, /* Expect value */
  2732. Blast_KarlinBlk* kbp,
  2733. Int8 searchsp) /* size of search space */
  2734. {
  2736. double Lambda, K, H; /* parameters for Karlin statistics */
  2737. Int4 S;
  2739. Lambda = kbp->Lambda;
  2740. K = kbp->K;
  2741. H = kbp->H;
  2742. if (Lambda < 0. || K < 0. || H < 0.) 
  2743. {
  2744. return BLAST_SCORE_MIN;
  2745. }
  2747. E = MAX(E, BLASTKAR_SMALL_FLOAT);
  2749. S = (Int4) (ceil( log((double)(K * searchsp / E)) / Lambda ));
  2750. return S;
  2751. }
  2753. /* Compute a divisor used to weight the evalue of a collection of
  2754.  * "nsegs" distinct alignments.  These divisors are used to compensate
  2755.  * for the effect of choosing the best among multiple collections of
  2756.  * alignments.  See
  2757.  *
  2758.  * Stephen F. Altschul. Evaluating the statitical significance of
  2759.  * multiple distinct local alignments. In Suhai, editior, Theoretical
  2760.  * and Computational Methods in Genome Research, pages 1-14. Plenum
  2761.  * Press, New York, 1997.
  2762.  *
  2763.  * The "decayrate" parameter of this routine is a value in the
  2764.  * interval (0,1). Typical values of decayrate are .1 and .5. */
  2766. double
  2767. BLAST_GapDecayDivisor(double decayrate, unsigned nsegs )
  2768. {
  2769.     return (1. - decayrate) * BLAST_Powi(decayrate, nsegs - 1);
  2770. }
  2773. /*
  2774. BlastCutoffs
  2775. Calculate the cutoff score, S, and the highest expected score.
  2776. */
  2777. Int2
  2778. BLAST_Cutoffs(Int4 *S, /* cutoff score */
  2779. double* E, /* expected no. of HSPs scoring at or above S */
  2780. Blast_KarlinBlk* kbp,
  2781. Int8 searchsp, /* size of search space. */
  2782. Boolean dodecay,  /* TRUE ==> use gapdecay feature */
  2783.    double gap_decay_rate)
  2784. {
  2785. Int4 s = *S, es;
  2786. double e = *E, esave;
  2787. Boolean s_changed = FALSE;
  2789. if (kbp->Lambda == -1. || kbp->K == -1. || kbp->H == -1.)
  2790. return 1;
  2792. /*
  2793. Calculate a cutoff score, S, from the Expected
  2794. (or desired) number of reported HSPs, E.
  2795. */
  2796. es = 1;
  2797. esave = e;
  2798. if (e > 0.) 
  2799. {
  2800.         if (dodecay) {
  2801.             /* Invert the adjustment to the e-value that will be applied
  2802.              * to compensate for the effect of choosing the best among
  2803.              * multiple alignments */
  2804.             if( gap_decay_rate > 0 && gap_decay_rate < 1 ) {
  2805.                 e *= BLAST_GapDecayDivisor(gap_decay_rate, 1);
  2806.             }
  2807.         }
  2808.         es = BlastKarlinEtoS_simple(e, kbp, searchsp);
  2809. }
  2810. /*
  2811. Pick the larger cutoff score between the user's choice
  2812. and that calculated from the value of E.
  2813. */
  2814. if (es > s) {
  2815. s_changed = TRUE;
  2816. *S = s = es;
  2817. }
  2819. /*
  2820. Re-calculate E from the cutoff score, if E going in was too high
  2821. */
  2822. if (esave <= 0. || !s_changed) 
  2823. {
  2824. e = BLAST_KarlinStoE_simple(s, kbp, searchsp);
  2825. if (dodecay) {
  2826.             /* Weight the e-value to compensate for the effect of
  2827.              * choosing the best of more than one collection of
  2828.              * distinct alignments */
  2829.             if( gap_decay_rate > 0 && gap_decay_rate < 1 ) {
  2830.                 e /= BLAST_GapDecayDivisor(gap_decay_rate, 1);
  2831.             }
  2832.         }
  2833. *E = e;
  2834. }
  2836. return 0;
  2837. }
  2839. /*
  2840. BlastKarlinStoE() -- given a score, return the associated Expect value
  2842. Error return value is -1.
  2843. */
  2844. double
  2845. BLAST_KarlinStoE_simple(Int4 S,
  2846. Blast_KarlinBlk* kbp,
  2847. Int8 searchsp) /* size of search space. */
  2848. {
  2849. double Lambda, K, H; /* parameters for Karlin statistics */
  2851. Lambda = kbp->Lambda;
  2852. K = kbp->K;
  2853. H = kbp->H;
  2854. if (Lambda < 0. || K < 0. || H < 0.) {
  2855. return -1.;
  2856. }
  2858. return (double) searchsp * exp((double)(-Lambda * S) + kbp->logK);
  2859. }
  2861. /*
  2862. BlastKarlinPtoE -- convert a P-value to an Expect value
  2863. */
  2864. static double
  2865. BlastKarlinPtoE(double p)
  2866. {
  2867.         if (p < 0. || p > 1.0) 
  2868. {
  2869.                 return INT4_MIN;
  2870.         }
  2872. if (p == 1)
  2873. return INT4_MAX;
  2875.         return -BLAST_Log1p(-p);
  2876. }
  2878. static double tab2[] = { /* table for r == 2 */
  2879. 0.01669,  0.0249,   0.03683,  0.05390,  0.07794,  0.1111,   0.1559,   0.2146,   
  2880. 0.2890,   0.3794,   0.4836,   0.5965,   0.7092,   0.8114,   0.8931,   0.9490,   
  2881. 0.9806,   0.9944,   0.9989
  2882. };
  2884. static double tab3[] = { /* table for r == 3 */
  2885. 0.9806,   0.9944,   0.9989,   0.0001682,0.0002542,0.0003829,0.0005745,0.0008587,
  2886. 0.001278, 0.001893, 0.002789, 0.004088, 0.005958, 0.008627, 0.01240,  0.01770,  
  2887. 0.02505,  0.03514,  0.04880,  0.06704,  0.09103,  0.1220,   0.1612,   0.2097,   
  2888. 0.2682,   0.3368,   0.4145,   0.4994,   0.5881,   0.6765,   0.7596,   0.8326,   
  2889. 0.8922,   0.9367,   0.9667,   0.9846,   0.9939,   0.9980
  2890. };
  2892. static double tab4[] = { /* table for r == 4 */
  2893. 2.658e-07,4.064e-07,6.203e-07,9.450e-07,1.437e-06,2.181e-06,3.302e-06,4.990e-06,
  2894. 7.524e-06,1.132e-05,1.698e-05,2.541e-05,3.791e-05,5.641e-05,8.368e-05,0.0001237,
  2895. 0.0001823,0.0002677,0.0003915,0.0005704,0.0008275,0.001195, 0.001718, 0.002457,
  2896. 0.003494, 0.004942, 0.006948, 0.009702, 0.01346,  0.01853,  0.02532,  0.03431,
  2897. 0.04607,  0.06128,  0.08068,  0.1051,   0.1352,   0.1719,   0.2157,   0.2669,
  2898. 0.3254,   0.3906,   0.4612,   0.5355,   0.6110,   0.6849,   0.7544,   0.8168,
  2899. 0.8699,   0.9127,   0.9451,   0.9679,   0.9827,   0.9915,   0.9963
  2900. };
  2902. static double* table[] = { tab2, tab3, tab4 };
  2903. static short tabsize[] = { DIM(tab2)-1, DIM(tab3)-1, DIM(tab4)-1 };
  2905. static double f (double,void*);
  2906. static double g (double,void*);
  2909. /*
  2910.     Estimate the Sum P-value by calculation or interpolation, as appropriate.
  2911. Approx. 2-1/2 digits accuracy minimum throughout the range of r, s.
  2912. r = number of segments
  2913. s = total score (in nats), adjusted by -r*log(KN)
  2914. */
  2915. static double
  2916. BlastSumP(Int4 r, double s)
  2917. {
  2918. Int4 i, r1, r2;
  2919. double a;
  2921. if (r == 1)
  2922. return -BLAST_Expm1(-exp(-s));
  2924. if (r <= 4) {
  2925. if (r < 1)
  2926. return 0.;
  2927. r1 = r - 1;
  2928. if (s >= r*r+r1) {
  2929. a = BLAST_LnGammaInt(r+1);
  2930. return r * exp(r1*log(s)-s-a-a);
  2931. }
  2932. if (s > -2*r) {
  2933. /* interpolate */
  2934. i = (Int4) (a = s+s+(4*r));
  2935. a -= i;
  2936. i = tabsize[r2 = r - 2] - i;
  2937. return a*table[r2][i-1] + (1.-a)*table[r2][i];
  2938. }
  2939. return 1.;
  2940. }
  2942. return BlastSumPCalc(r, s);
  2943. }
  2945. /*
  2946.     BlastSumPCalc
  2948.     Evaluate the following double integral, where r = number of segments
  2949.     and s = the adjusted score in nats:
  2951.                     (r-2)         oo           oo
  2952.      Prob(r,s) =   r              -            -   (r-2)
  2953.                  -------------   |   exp(-y)  |   x   exp(-exp(x - y/r)) dx dy
  2954.                  (r-1)! (r-2)!  U            U
  2955.                                 s            0
  2956. */
  2957. static double
  2958. BlastSumPCalc(int r, double s)
  2959. {
  2960. int r1, itmin;
  2961. double t, d, epsilon;
  2962. double est_mean, mean, stddev, stddev4;
  2963. double xr, xr1, xr2, logr;
  2964. double args[6];
  2966. epsilon = BLAST_SUMP_EPSILON_DEFAULT; /* accuracy for SumP calcs. */
  2968. if (r == 1) {
  2969. if (s > 8.)
  2970. return exp(-s);
  2971. return -BLAST_Expm1(-exp(-s));
  2972. }
  2973. if (r < 1)
  2974. return 0.;
  2976. xr = r;
  2978. if (r < 8) {
  2979. if (s <= -2.3*xr)
  2980. return 1.;
  2981. }
  2982. else if (r < 15) {
  2983. if (s <= -2.5*xr)
  2984. return 1.;
  2985. }
  2986. else if (r < 27) {
  2987. if (s <= -3.0*xr)
  2988. return 1.;
  2989. }
  2990. else if (r < 51) {
  2991. if (s <= -3.4*xr)
  2992. return 1.;
  2993. }
  2994. else if (r < 101) {
  2995. if (s <= -4.0*xr)
  2996. return 1.;
  2997. }
  2999. /* stddev in the limit of infinite r, but quite good for even small r */
  3000. stddev = sqrt(xr);
  3001. stddev4 = 4.*stddev;
  3002. xr1 = r1 = r - 1;
  3004. if (r > 100) {
  3005. /* Calculate lower bound on the mean using inequality log(r) <= r */
  3006. est_mean = -xr * xr1;
  3007. if (s <= est_mean - stddev4)
  3008. return 1.;
  3009. }
  3011. /* mean is rather close to the mode, and the mean is readily calculated */
  3012. /* mean in the limit of infinite r, but quite good for even small r */
  3013. logr = log(xr);
  3014. mean = xr * (1. - logr) - 0.5;
  3015. if (s <= mean - stddev4)
  3016. return 1.;
  3018. if (s >= mean) {
  3019. t = s + 6.*stddev;
  3020. itmin = 1;
  3021. }
  3022. else {
  3023. t = mean + 6.*stddev;
  3024. itmin = 2;
  3025. }
  3027. #define ARG_R args[0]
  3028. #define ARG_R2 args[1]
  3029. #define ARG_ADJ1 args[2]
  3030. #define ARG_ADJ2 args[3]
  3031. #define ARG_SDIVR args[4]
  3032. #define ARG_EPS args[5]
  3034. ARG_R = xr;
  3035. ARG_R2 = xr2 = r - 2;
  3036. ARG_ADJ1 = xr2*logr - BLAST_LnGammaInt(r1) - BLAST_LnGammaInt(r);
  3037. ARG_EPS = epsilon;
  3039. do {
  3040. d = BLAST_RombergIntegrate(g, args, s, t, epsilon, 0, itmin);
  3041. #ifdef BLASTKAR_HUGE_VAL
  3042. if (d == BLASTKAR_HUGE_VAL)
  3043. return d;
  3044. #endif
  3045. } while (s < mean && d < 0.4 && itmin++ < 4);
  3047. return (d < 1. ? d : 1.);
  3048. }
  3050. static double
  3051. g(double s, void* vp)
  3052. {
  3053. register double* args = vp;
  3054. double mx;
  3056. ARG_ADJ2 = ARG_ADJ1 - s;
  3057. ARG_SDIVR = s / ARG_R; /* = s / r */
  3058. mx = (s > 0. ? ARG_SDIVR + 3. : 3.);
  3059. return BLAST_RombergIntegrate(f, vp, 0., mx, ARG_EPS, 0, 1);
  3060. }
  3062. static double
  3063. f(double x, void* vp)
  3064. {
  3065. register double* args = vp;
  3066. register double y;
  3068. y = exp(x - ARG_SDIVR);
  3069. #ifdef BLASTKAR_HUGE_VAL
  3070. if (y == BLASTKAR_HUGE_VAL)
  3071. return 0.;
  3072. #endif
  3073. if (ARG_R2 == 0.)
  3074. return exp(ARG_ADJ2 - y);
  3075. if (x == 0.)
  3076. return 0.;
  3077. return exp(ARG_R2*log(x) + ARG_ADJ2 - y);
  3078. }
  3080. /*
  3081.     Calculates the e-value for alignments with "small" gaps (typically
  3082.     under fifty residues/basepairs) following ideas of Stephen Altschul's.
  3083. */
  3085. double
  3086. BLAST_SmallGapSumE(
  3087.     Blast_KarlinBlk * kbp,     /* statistical parameters */
  3088.     Int4 gap,                   /* maximum size of gaps between alignments */
  3089.     Int2 num,                   /* the number of distinct alignments in this
  3090.                                  * collection */
  3091.     double score_prime,         /* the sum of the scores of these alignments
  3092.                                  * each weighted by an appropriate value of
  3093.                                  * Lambda */
  3094.     Int4 query_length,          /* the effective len of the query seq */
  3095.     Int4 subject_length,        /* the effective len of the database seq */
  3096.     double weight_divisor)      /* a divisor used to weight the e-value
  3097.                                  * when multiple collections of alignments
  3098.                                  * are being considered by the calling
  3099.                                  * routine */
  3100. {
  3101.     double search_space;        /* The effective size of the search space */
  3102.     double sum_p;               /* The p-value of this set of alignments */
  3103.     double sum_e;               /* The e-value of this set of alignments */
  3105.     search_space = (double)subject_length*(double)query_length;
  3107.     score_prime -= kbp->logK +
  3108.         log(search_space) + (num-1)*(kbp->logK + 2*log((double)gap));
  3109.     score_prime -= BLAST_LnFactorial((double) num);
  3111.     sum_p = BlastSumP(num, score_prime);
  3113.     sum_e = BlastKarlinPtoE(sum_p);
  3115.     if( weight_divisor == 0 || (sum_e /= weight_divisor) > INT4_MAX ) {
  3116.         sum_e = INT4_MAX;
  3117.     }
  3119.     return sum_e;
  3120. }
  3122. /*
  3123.     Calculates the e-value of a collection multiple distinct
  3124.     alignments with asymmetric gaps between the alignments. The gaps
  3125.     in one (protein) sequence are typically small (like in
  3126.     BLAST_SmallGapSumE) gap an the gaps in the other (translated DNA)
  3127.     sequence are possibly large (up to 4000 bp.)  This routine is used
  3128.     for linking HSPs representing exons in the DNA sequence that are
  3129.     separated by introns.
  3130. */
  3132. double
  3133. BLAST_UnevenGapSumE(
  3134.     Blast_KarlinBlk * kbp,      /* statistical parameters */
  3135.     Int4 p_gap,                 /* maximum size of gaps between alignments,
  3136.                                  * in one sequence */
  3137.     Int4 n_gap,                 /* maximum size of gaps between alignments,
  3138.                                  * in the other sequence */
  3139.     Int2 num,                   /* the number of distinct alignments in this
  3140.                                  * collection */
  3141.     double score_prime,         /* the sum of the scores of these alignments
  3142.                                  * each weighted by an appropriate value of
  3143.                                  * Lambda */
  3144.     Int4 query_length,          /* the effective len of the query seq */
  3145.     Int4 subject_length,        /* the effective len of the database seq */
  3146.     double weight_divisor)      /* a divisor used to weight the e-value
  3147.                                  * when multiple collections of alignments
  3148.                                  * are being considered by the calling
  3149.                                  * routine */
  3150. {
  3151.     double search_space;   /* The effective size of the search space */
  3152.     double sum_p;          /* The p-value of this set of alignments */
  3153.     double sum_e;          /* The e-value of this set of alignments */
  3155.     search_space = (double)subject_length*(double)query_length;
  3157.     score_prime -=
  3158.         kbp->logK + log(search_space) +
  3159.         (num-1)*(kbp->logK + log((double)p_gap) + log((double)n_gap));
  3160.     score_prime -= BLAST_LnFactorial((double) num);
  3162.     sum_p = BlastSumP(num, score_prime);
  3164.     sum_e = BlastKarlinPtoE(sum_p);
  3166.     if( weight_divisor == 0 || (sum_e /= weight_divisor) > INT4_MAX ) {
  3167.         sum_e = INT4_MAX;
  3168.     }
  3170.     return sum_e;
  3171. }
  3174. /*
  3175.     Calculates the e-value if a collection of distinct alignments with
  3176.     arbitrarily large gaps between the alignments
  3177. */
  3179. double
  3180. BLAST_LargeGapSumE(
  3181.     Blast_KarlinBlk * kbp,      /* statistical parameters */
  3182.     Int2 num,                   /* the number of distinct alignments in this
  3183.                                  * collection */
  3184.     double      score_prime,    /* the sum of the scores of these alignments
  3185.                                  * each weighted by an appropriate value of
  3186.                                  * Lambda */
  3187.     Int4 query_length,          /* the effective len of the query seq */
  3188.     Int4 subject_length,        /* the effective len of the database seq */
  3189.     double weight_divisor)      /* a divisor used to weight the e-value
  3190.                                  * when multiple collections of alignments
  3191.                                  * are being considered by the calling
  3192.                                  * routine */
  3193. {
  3194.     double sum_p;               /* The p-value of this set of alignments */
  3195.     double sum_e;               /* The e-value of this set of alignments */
  3197. /* The next two variables are for compatability with Warren's code. */
  3198.     double lcl_subject_length;     /* query_length as a float */
  3199.     double lcl_query_length;       /* subject_length as a float */
  3201.     lcl_query_length = (double) query_length;
  3202.     lcl_subject_length = (double) subject_length;
  3204.     score_prime -= num*(kbp->logK + log(lcl_subject_length*lcl_query_length))
  3205.         - BLAST_LnFactorial((double) num);
  3207.     sum_p = BlastSumP(num, score_prime);
  3209.     sum_e = BlastKarlinPtoE(sum_p);
  3211.     if( weight_divisor == 0 || (sum_e /= weight_divisor) > INT4_MAX ) {
  3212.         sum_e = INT4_MAX;
  3213.     }
  3215.     return sum_e;
  3216. }
  3218. /*------------------- RPS BLAST functions --------------------*/
  3220. static double
  3221. RPSfindUngappedLambda(Char *matrixName)
  3222. {
  3223.     if (0 == strcmp(matrixName, "BLOSUM62"))
  3224.         return 0.3176;
  3225.     if (0 == strcmp(matrixName, "BLOSUM90"))
  3226.         return 0.3346;
  3227.     if (0 == strcmp(matrixName, "BLOSUM80"))
  3228.         return 0.3430;
  3229.     if (0 == strcmp(matrixName, "BLOSUM50"))
  3230.         return 0.232;
  3231.     if (0 == strcmp(matrixName, "BLOSUM45"))
  3232.         return 0.2291;
  3233.     if (0 == strcmp(matrixName, "PAM30"))
  3234.         return 0.340;
  3235.     if (0 == strcmp(matrixName, "PAM70"))
  3236.         return 0.3345;
  3237.     if (0 == strcmp(matrixName, "PAM250"))
  3238.         return 0.229;
  3239.     return(0);
  3240. }
  3242. /* matrix is a position-specific score matrix with matrixLength positions
  3243.    queryProbArray is an array containing the probability of occurrence
  3244.         of each residue in the query
  3245.    scoreArray is an array of probabilities for each score that is
  3246.         to be used as a field in return_sfp
  3247.    return_sfp is a the structure to be filled in and returned
  3248.    range is the size of scoreArray and is an upper bound on the
  3249.         difference between maximum score and minimum score in the matrix
  3250.    the routine fillSfp computes the probability of each score weighted
  3251.         by the probability of each query residue and fills those probabilities
  3252.         into scoreArray and puts scoreArray as a field in
  3253.         that in the structure that is returned
  3254.    for indexing convenience the field storing scoreArray points to the
  3255.         entry for score 0, so that referring to the -k index corresponds to
  3256.         score -k 
  3257. */
  3259. static void
  3260. RPSFillScores(Int4 **matrix, Int4 matrixLength, 
  3261.               double *queryProbArray, double *scoreArray,  
  3262.               Blast_ScoreFreq* return_sfp, Int4 range)
  3263. {
  3264.     Int4 minScore, maxScore;    /*observed minimum and maximum scores */
  3265.     Int4 i,j;                   /* indices */
  3266.     double recipLength;        /* 1/matrix length as a double */
  3268.     minScore = maxScore = 0;
  3270.     for (i = 0; i < matrixLength; i++) {
  3271.         for (j = 0 ; j < PSI_ALPHABET_SIZE; j++) {
  3272.             if (j == AMINOACID_TO_NCBISTDAA['X'])
  3273.                 continue;
  3274.             if ((matrix[i][j] > BLAST_SCORE_MIN) && 
  3275.                 (matrix[i][j] < minScore))
  3276.                 minScore = matrix[i][j];
  3277.             if (matrix[i][j] > maxScore)
  3278.                 maxScore = matrix[i][j];
  3279.         }
  3280.     }
  3282.     return_sfp->obs_min = minScore;
  3283.     return_sfp->obs_max = maxScore;
  3284.     for (i = 0; i < range; i++)
  3285.         scoreArray[i] = 0.0;
  3287.     return_sfp->sprob = &(scoreArray[-minScore]); /*center around 0*/
  3288.     recipLength = 1.0 / (double) matrixLength;
  3289.     for(i = 0; i < matrixLength; i++) {
  3290.         for (j = 0; j < PSI_ALPHABET_SIZE; j++) {
  3291.             if (j == AMINOACID_TO_NCBISTDAA['X'])
  3292.                 continue;
  3293.             if(matrix[i][j] >= minScore)
  3294.                 return_sfp->sprob[matrix[i][j]] += recipLength * 
  3295.                                                 queryProbArray[j];
  3296.         }
  3297.     }
  3299.     return_sfp->score_avg = 0;
  3300.     for(i = minScore; i <= maxScore; i++)
  3301.         return_sfp->score_avg += i * return_sfp->sprob[i];
  3302. }
  3304. /*  Given a sequence of 'length' amino acid residues, compute the
  3305.     probability of each residue and put that in the array resProb*/
  3307. static void 
  3308. RPSFillResidueProbability(Uint1 * sequence, Int4 length, double * resProb)
  3309. {
  3310.     Int4 frequency[PSI_ALPHABET_SIZE];  /*frequency of each letter*/
  3311.     Int4 i;                             /*index*/
  3312.     Int4 denominator;                   /*length not including X's*/
  3314.     denominator = length;
  3315.     for(i = 0; i < PSI_ALPHABET_SIZE; i++)
  3316.         frequency[i] = 0;
  3318.     for(i = 0; i < length; i++) {
  3319.         if (sequence[i] != AMINOACID_TO_NCBISTDAA['X'])
  3320.             frequency[sequence[i]]++;
  3321.         else
  3322.             denominator--;
  3323.     }
  3325.     for(i = 0; i < PSI_ALPHABET_SIZE; i++) {
  3326.         if (frequency[i] == 0)
  3327.             resProb[i] = 0.0;
  3328.         else
  3329.             resProb[i] = ((double) frequency[i]) /((double) denominator);
  3330.     }
  3331. }
  3333. /* Calculate a new PSSM, using composition-based statistics, for use
  3334.    with RPS BLAST. This function produces a PSSM for a single RPS DB
  3335.    sequence (of size db_seq_length) and incorporates information from 
  3336.    the RPS blast query. Each individual database sequence must call this
  3337.    function to retrieve its own PSSM. The matrix is returned (and must
  3338.    be freed elsewhere). posMatrix is the portion of the complete 
  3339.    concatenated PSSM that is specific to this DB sequence */
  3341. Int4 **
  3342. RPSCalculatePSSM(double scalingFactor, Int4 rps_query_length, 
  3343.                    Uint1 * rps_query_seq, Int4 db_seq_length, 
  3344.                    Int4 **posMatrix)
  3345. {
  3346.     double *scoreArray;         /*array of score probabilities*/
  3347.     double *resProb;            /*array of probabilities for each residue*/
  3348.     Blast_ScoreFreq * return_sfp;/*score frequency pointers to compute lambda*/
  3349.     Int4* * returnMatrix;        /*the PSSM to return */
  3350.     double initialUngappedLambda; 
  3351.     double scaledInitialUngappedLambda; 
  3352.     double correctUngappedLambda;
  3353.     double finalLambda;
  3354.     double temp;               /*intermediate variable for adjusting matrix*/
  3355.     Int4 index, inner_index; 
  3357.     resProb = (double *)malloc(PSI_ALPHABET_SIZE * sizeof(double));
  3358.     scoreArray = (double *)malloc(BLAST_SCORE_RANGE_MAX * sizeof(double));
  3359.     return_sfp = (Blast_ScoreFreq *)malloc(sizeof(Blast_ScoreFreq));
  3361.     RPSFillResidueProbability(rps_query_seq, rps_query_length, resProb);
  3363.     RPSFillScores(posMatrix, db_seq_length, resProb, scoreArray, 
  3364.                  return_sfp, BLAST_SCORE_RANGE_MAX);
  3366.     initialUngappedLambda = RPSfindUngappedLambda("BLOSUM62");
  3367.     scaledInitialUngappedLambda = initialUngappedLambda / scalingFactor;
  3368.     correctUngappedLambda = Blast_KarlinLambdaNR(return_sfp, 
  3369.                                               scaledInitialUngappedLambda);
  3370.     sfree(resProb);
  3371.     sfree(scoreArray);
  3372.     sfree(return_sfp);
  3373.     if(correctUngappedLambda == -1.0)
  3374.         return NULL;
  3376.     finalLambda = correctUngappedLambda/scaledInitialUngappedLambda;
  3378.     returnMatrix = (Int4 **)_PSIAllocateMatrix((db_seq_length+1),
  3379.                                                PSI_ALPHABET_SIZE,
  3380.                                                sizeof(Int4));
  3382.     for (index = 0; index < db_seq_length+1; index++) {
  3383.         for (inner_index = 0; inner_index < PSI_ALPHABET_SIZE; inner_index++) {
  3384.             if (posMatrix[index][inner_index] <= BLAST_SCORE_MIN || 
  3385.                 inner_index == AMINOACID_TO_NCBISTDAA['X']) {
  3386.                 returnMatrix[index][inner_index] = 
  3387.                     posMatrix[index][inner_index];
  3388.             }
  3389.             else {
  3390.                temp = ((double)(posMatrix[index][inner_index])) * finalLambda;
  3391.                returnMatrix[index][inner_index] = BLAST_Nint(temp);
  3392.            }
  3393.         }
  3394.     }
  3396.     return returnMatrix;
  3397. }
  3399. /** 
  3400.  * Computes the adjustment to the lengths of the query and database sequences
  3401.  * that is used to compensate for edge effects when computing evalues. 
  3402.  *
  3403.  * The length adjustment is an integer-valued approximation to the fixed
  3404.  * point of the function
  3405.  *
  3406.  *    f(ell) = beta + 
  3407.  *               (alpha/lambda) * (log K + log((m - ell)*(n - N ell)))
  3408.  *
  3409.  * where m is the query length n is the length of the database and N is the
  3410.  * number of sequences in the database. The values beta, alpha, lambda and
  3411.  * K are statistical, Karlin-Altschul parameters.
  3412.  * 
  3413.  * The value of the length adjustment computed by this routine, A, 
  3414.  * will always be an integer smaller than the fixed point of
  3415.  * f(ell). Usually, it will be the largest such integer.  However, the
  3416.  * computed length adjustment, A, will also be so small that 
  3417.  *
  3418.  *    K * (m - A) * (n - N * A) > MAX(m,n).
  3419.  *
  3420.  * Moreover, an iterative method is used to compute A, and under
  3421.  * unusual circumstances the iterative method may not converge. 
  3422.  *
  3423.  * @param K      the statistical parameter K
  3424.  * @param logK   the natural logarithm of K
  3425.  * @param alpha_d_lambda    the ratio of the statistical parameters 
  3426.  *                          alpha and lambda (for ungapped alignments, the
  3427.  *                          value 1/H should be used)
  3428.  * @param beta              the statistical parameter beta (for ungapped
  3429.  *                          alignments, beta == 0)
  3430.  * @param query_length      the length of the query sequence
  3431.  * @param db_length         the length of the database
  3432.  * @param db_num_seqs       the number of sequences in the database
  3433.  * @param length_adjustment the computed value of the length adjustment [out]
  3434.  *
  3435.  * @return   0 if length_adjustment is known to be the largest integer less
  3436.  *           than the fixed point of f(ell); 1 otherwise.
  3437.  */
  3438. Int4
  3439. BLAST_ComputeLengthAdjustment(double K,
  3440.                              double logK,
  3441.                              double alpha_d_lambda,
  3442.                              double beta,
  3443.                              Int4 query_length,
  3444.                              Int8 db_length,
  3445.                              Int4 db_num_seqs,
  3446.                              Int4 * length_adjustment)
  3447. {
  3448.     Int4 i;                     /* iteration index */
  3449.     const Int4 maxits = 20;     /* maximum allowed iterations */
  3450.     double m = (double) query_length;
  3451.     double n = (double) db_length;
  3452.     double N = (double) db_num_seqs;
  3454.     double ell;                 /* A float value of the length adjustment */
  3455.     double ss;                  /* effective size of the search space */
  3456.     double ell_min = 0, ell_max;   /* At each iteration i,
  3457.                                          * ell_min <= ell <= ell_max. */
  3458.     Boolean converged    = FALSE;       /* True if the iteration converged */
  3459.     double ell_next = 0;        /* Value the variable ell takes at iteration
  3460.                                  * i + 1 */
  3461.     /* Choose ell_max to be the largest nonnegative value that satisfies
  3462.      *
  3463.      *    K * (m - ell) * (n - N * ell) > MAX(m,n)
  3464.      *
  3465.      * Use quadratic formula: 2 c /( - b + sqrt( b*b - 4 * a * c )) */
  3466.     { /* scope of a, mb, and c, the coefficients in the quadratic formula
  3467.        * (the variable mb is -b) */
  3468.         double a  = N;
  3469.         double mb = m * N + n;
  3470.         double c  = n * m - MAX(m, n) / K;
  3472.         if(c < 0) {
  3473.             *length_adjustment = 0;
  3474.             return 1;
  3475.         } else {
  3476.             ell_max = 2 * c / (mb + sqrt(mb * mb - 4 * a * c));
  3477.         }
  3478.     } /* end scope of a, mb and c */
  3480.     for(i = 1; i <= maxits; i++) {      /* for all iteration indices */
  3481.         double ell_bar;         /* proposed next value of ell */
  3482.         ell      = ell_next;
  3483.         ss       = (m - ell) * (n - N * ell);
  3484.         ell_bar  = alpha_d_lambda * (logK + log(ss)) + beta;
  3485.         if(ell_bar >= ell) { /* ell is no bigger than the true fixed point */
  3486.             ell_min = ell;
  3487.             if(ell_bar - ell_min <= 1.0) {
  3488.                 converged = TRUE;
  3489.                 break;
  3490.             }
  3491.             if(ell_min == ell_max) { /* There are no more points to check */
  3492.                 break;
  3493.             }
  3494.         } else { /* else ell is greater than the true fixed point */
  3495.             ell_max = ell;
  3496.         }
  3497.         if(ell_min <= ell_bar && ell_bar <= ell_max) { 
  3498.           /* ell_bar is in range. Accept it */
  3499.             ell_next = ell_bar;
  3500.         } else { /* else ell_bar is not in range. Reject it */
  3501.             ell_next = (i == 1) ? ell_max : (ell_min + ell_max) / 2;
  3502.         }
  3503.     } /* end for all iteration indices */
  3504.     if(converged) { /* the iteration converged */
  3505.         /* If ell_fixed is the (unknown) true fixed point, then we
  3506.          * wish to set (*length_adjustment) to floor(ell_fixed).  We
  3507.          * assume that floor(ell_min) = floor(ell_fixed) */
  3508.         *length_adjustment = (Int4) ell_min;
  3509.         /* But verify that ceil(ell_min) != floor(ell_fixed) */
  3510.         ell = ceil(ell_min);
  3511.         if( ell <= ell_max ) {
  3512.           ss = (m - ell) * (n - N * ell);
  3513.           if(alpha_d_lambda * (logK + log(ss)) + beta >= ell) {
  3514.             /* ceil(ell_min) == floor(ell_fixed) */
  3515.             *length_adjustment = (Int4) ell;
  3516.           }
  3517.         }
  3518.     } else { /* else the iteration did not converge. */
  3519.         /* Use the best value seen so far */
  3520.         *length_adjustment = (Int4) ell_min;
  3521.     }
  3523.     return converged ? 0 : 1;
  3524. }
  3526. /*
  3527.  * ===========================================================================
  3528.  *
  3529.  * $Log: blast_stat.c,v $
  3530.  * Revision 1000.6  2004/06/01 18:07:50  gouriano
  3531.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.77
  3532.  *
  3533.  * Revision 1.77  2004/05/24 15:09:40  camacho
  3534.  * Fixed conflict
  3535.  *
  3536.  * Revision 1.76  2004/05/24 13:26:27  madden
  3537.  * Fix PC compiler warnings
  3538.  *
  3539.  * Revision 1.75  2004/05/20 16:29:30  madden
  3540.  * Make searchsp an Int8 consistent with rest of blast
  3541.  *
  3542.  * Revision 1.74  2004/05/19 15:34:38  dondosha
  3543.  * Moved Blast_ResComp definition from header file
  3544.  *
  3545.  * Revision 1.73  2004/05/19 14:52:03  camacho
  3546.  * 1. Added doxygen tags to enable doxygen processing of algo/blast/core
  3547.  * 2. Standardized copyright, CVS $Id string, $Log and rcsid formatting and i
  3548.  *    location
  3549.  * 3. Added use of @todo doxygen keyword
  3550.  *
  3551.  * Revision 1.72  2004/05/17 10:37:38  camacho
  3552.  * Rename BLAST_ScoreFreq, BLASTMatrixStructure and BLAST_ResComp to avoid conflicts with C toolkit
  3553.  *
  3554.  * Revision 1.71  2004/05/07 15:23:47  papadopo
  3555.  * add initialization of scale factor to ScoreBlkNew
  3556.  *
  3557.  * Revision 1.70  2004/05/06 15:59:29  camacho
  3558.  * Made Blast_KarlinBlkCalc non-static
  3559.  *
  3560.  * Revision 1.69  2004/05/06 15:05:13  camacho
  3561.  * Fix to previous commit
  3562.  *
  3563.  * Revision 1.68  2004/05/06 14:44:27  camacho
  3564.  * Made Blast_ScoreFreqDestruct non-static
  3565.  *
  3566.  * Revision 1.67  2004/05/05 21:16:24  camacho
  3567.  * Make Blast_GetStdAlphabet and Blast_ScoreFreqNew non-static
  3568.  *
  3569.  * Revision 1.66  2004/05/04 13:00:02  madden
  3570.  * Change BlastKarlinBlkStandardCalcEx to more descriptive Blast_KarlinBlkIdealCalc, make public
  3571.  *
  3572.  * Revision 1.65  2004/04/30 14:39:44  papadopo
  3573.  * 1. Remove unneeded #defines
  3574.  * 2. use BLAST_SCORE_RANGE_MAX during RPS PSSM creation instead of
  3575.  *  (possibly incompatible) RPS_SCORE_MAX
  3576.  * 3. return NULL instead of FALSE on an error
  3577.  *
  3578.  * Revision 1.64  2004/04/30 12:58:49  camacho
  3579.  * Replace RPSKarlinLambdaNR by Blast_KarlinLambdaNR
  3580.  *
  3581.  * Revision 1.63  2004/04/29 20:32:38  papadopo
  3582.  * remove RPS_SCORE_MIN, since it turned out to be a workaround for a bug that has since been fixed
  3583.  *
  3584.  * Revision 1.62  2004/04/29 19:58:03  camacho
  3585.  * Use generic matrix allocator/deallocator from blast_psi_priv.h
  3586.  *
  3587.  * Revision 1.61  2004/04/28 14:40:23  madden
  3588.  * Changes from Mike Gertz:
  3589.  * - I created the new routine BLAST_GapDecayDivisor that computes a
  3590.  *   divisor used to weight the evalue of a collection of distinct
  3591.  *   alignments.
  3592.  * - I removed  BLAST_GapDecay and BLAST_GapDecayInverse which had become
  3593.  *   redundant.
  3594.  * - I modified the BLAST_Cutoffs routine so that it uses the value
  3595.  *   returned by BLAST_GapDecayDivisor to weight evalues.
  3596.  * - I modified BLAST_SmallGapSumE, BLAST_LargeGapSumE and
  3597.  *   BLAST_UnevenGapSumE no longer refer to the gap_prob parameter.
  3598.  *   Replaced the gap_decay_rate parameter of each of these routines with
  3599.  *   a weight_divisor parameter.  Added documentation.
  3600.  *
  3601.  * Revision 1.60  2004/04/23 19:06:33  camacho
  3602.  * Do NOT use lowercase names for #defines
  3603.  *
  3604.  * Revision 1.59  2004/04/23 13:49:20  madden
  3605.  * Cleaned up ifndef in BlastKarlinLHtoK
  3606.  *
  3607.  * Revision 1.58  2004/04/23 13:21:25  madden
  3608.  * Rewrote BlastKarlinLHtoK to do the following and more:
  3609.  * 1. fix a bug whereby the wrong formula was used when high score == 1
  3610.  *    and low score == -1;
  3611.  * 2. fix a methodological error of truncating the first sum
  3612.  *    and trying to make it converge quickly by adding terms
  3613.  *    of a geometric progression, even though the geometric progression
  3614.  *    estimate is not correct in all cases;
  3615.  *    the old adjustment code is left in for historical purposes but
  3616.  *    #ifdef'd out
  3617.  * 3. Eliminate the Boolean bi_modal_score variable.  The old test that
  3618.  *    set the value of bi_modal_score would frequently fail to choose the
  3619.  *    correct value due to rounding error.
  3620.  * 4. changed numerous local variable names to make them more meaningful;
  3621.  * 5. added substantial comments to explain what the procedure
  3622.  *    is doing and what each variable represents
  3623.  *
  3624.  * Revision 1.57  2004/04/19 12:58:18  madden
  3625.  * Changed BLAST_KarlinBlk to Blast_KarlinBlk to avoid conflict with blastkar.h structure, renamed some functions to start with Blast_Karlin, made Blast_KarlinBlkDestruct public
  3626.  *
  3627.  * Revision 1.56  2004/04/12 18:57:31  madden
  3628.  * Rename BLAST_ResFreq to Blast_ResFreq, make Blast_ResFreqNew, Blast_ResFreqDestruct, and Blast_ResFreqStdComp non-static
  3629.  *
  3630.  * Revision 1.55  2004/04/08 13:53:10  papadopo
  3631.  * fix doxygen warning
  3632.  *
  3633.  * Revision 1.54  2004/04/07 03:06:16  camacho
  3634.  * Added blast_encoding.[hc], refactoring blast_stat.[hc]
  3635.  *
  3636. v * Revision 1.53  2004/04/05 18:53:35  madden
  3637.  * Set dimensions if matrix from memory
  3638.  *
  3639.  * Revision 1.52  2004/04/01 14:14:02  lavr
  3640.  * Spell "occurred", "occurrence", and "occurring"
  3641.  *
  3642.  * Revision 1.51  2004/03/31 17:50:09  papadopo
  3643.  * Mike Gertz' changes for length adjustment calculations
  3644.  *
  3645.  * Revision 1.50  2004/03/11 18:52:41  camacho
  3646.  * Remove THREADS_IMPLEMENTED
  3647.  *
  3648.  * Revision 1.49  2004/03/10 18:00:06  camacho
  3649.  * Remove outdated references to blastkar
  3650.  *
  3651.  * Revision 1.48  2004/03/05 17:52:33  papadopo
  3652.  * Allow 32-bit context numbers for queries
  3653.  *
  3654.  * Revision 1.47  2004/03/04 21:07:51  papadopo
  3655.  * add RPS BLAST functionality
  3656.  *
  3657.  * Revision 1.46  2004/02/19 21:16:48  dondosha
  3658.  * Use enum type for severity argument in Blast_MessageWrite
  3659.  *
  3660.  * Revision 1.45  2003/12/05 16:03:57  camacho
  3661.  * Remove compiler warnings
  3662.  *
  3663.  * Revision 1.44  2003/11/28 22:39:11  camacho
  3664.  * + static keyword to BlastKarlinLtoH
  3665.  *
  3666.  * Revision 1.43  2003/11/28 15:03:48  camacho
  3667.  * Added static keyword to BlastKarlinLtoH
  3668.  *
  3669.  * Revision 1.42  2003/11/26 19:12:13  madden
  3670.  * code to simplify some routines and use NlmKarlinLambdaNR in place of BlastKarlinLambdaBis (following Mike Gertzs changes to blastkar.c )
  3671.  *
  3672.  * Revision 1.41  2003/11/24 23:18:32  dondosha
  3673.  * Added gap_decay_rate argument to BLAST_Cutoffs; removed BLAST_Cutoffs_simple
  3674.  *
  3675.  * Revision 1.40  2003/11/19 15:17:42  dondosha
  3676.  * Removed unused members from Karlin block structure
  3677.  *
  3678.  * Revision 1.39  2003/10/16 15:55:22  coulouri
  3679.  * fix uninitialized variables
  3680.  *
  3681.  * Revision 1.38  2003/10/16 15:52:08  coulouri
  3682.  * fix uninitialized variables
  3683.  *
  3684.  * Revision 1.37  2003/10/15 16:59:43  coulouri
  3685.  * type correctness fixes
  3686.  *
  3687.  * Revision 1.36  2003/10/02 22:08:34  dondosha
  3688.  * Corrections for one-strand translated searches
  3689.  *
  3690.  * Revision 1.35  2003/09/26 19:01:59  madden
  3691.  * Prefix ncbimath functions with BLAST_
  3692.  *
  3693.  * Revision 1.34  2003/09/09 14:21:39  coulouri
  3694.  * change blastkar.h to blast_stat.h
  3695.  *
  3696.  * Revision 1.33  2003/09/02 21:12:07  camacho
  3697.  * Fix small memory leak
  3698.  *
  3699.  * Revision 1.32  2003/08/26 15:23:51  dondosha
  3700.  * Rolled back previous change as it is not necessary any more
  3701.  *
  3702.  * Revision 1.31  2003/08/25 22:29:07  dondosha
  3703.  * Default matrix loading is defined only in C++ toolkit
  3704.  *
  3705.  * Revision 1.30  2003/08/25 18:05:41  dondosha
  3706.  * Moved assert statement after variables declarations
  3707.  *
  3708.  * Revision 1.29  2003/08/25 16:23:33  camacho
  3709.  * +Loading protein scoring matrices from utils/tables
  3710.  *
  3711.  * Revision 1.28  2003/08/11 15:01:59  dondosha
  3712.  * Added algo/blast/core to all #included headers
  3713.  *
  3714.  * Revision 1.27  2003/08/01 17:27:04  dondosha
  3715.  * Renamed external functions to avoid collisions with ncbitool library; made other functions static
  3716.  *
  3717.  * Revision 1.26  2003/07/31 18:48:49  dondosha
  3718.  * Use Int4 instead of BLAST_Score
  3719.  *
  3720.  * Revision 1.25  2003/07/31 17:48:06  madden
  3721.  * Remove call to FileLength
  3722.  *
  3723.  * Revision 1.24  2003/07/31 14:31:41  camacho
  3724.  * Replaced Char for char
  3725.  *
  3726.  * Revision 1.23  2003/07/31 14:19:28  camacho
  3727.  * Replaced FloatHi for double
  3728.  *
  3729.  * Revision 1.22  2003/07/31 00:32:37  camacho
  3730.  * Eliminated Ptr notation
  3731.  *
  3732.  * Revision 1.21  2003/07/30 22:08:09  dondosha
  3733.  * Process of finding path to the matrix is moved out of the blast library
  3734.  *
  3735.  * Revision 1.20  2003/07/30 21:52:41  camacho
  3736.  * Follow conventional structure definition
  3737.  *
  3738.  * Revision 1.19  2003/07/30 19:39:14  camacho
  3739.  * Remove PNTRs
  3740.  *
  3741.  * Revision 1.18  2003/07/30 17:58:25  dondosha
  3742.  * Changed ValNode to ListNode
  3743.  *
  3744.  * Revision 1.17  2003/07/30 17:15:00  dondosha
  3745.  * Minor fixes for very strict compiler warnings
  3746.  *
  3747.  * Revision 1.16  2003/07/30 17:06:40  camacho
  3748.  * Removed old cvs log
  3749.  *
  3750.  * Revision 1.15  2003/07/30 16:32:02  madden
  3751.  * Use ansi functions when possible
  3752.  *
  3753.  * Revision 1.14  2003/07/30 15:29:37  madden
  3754.  * Removed MemSets
  3755.  *
  3756.  * Revision 1.13  2003/07/29 14:42:31  coulouri
  3757.  * use strdup() instead of StringSave()
  3758.  *
  3759.  * Revision 1.12  2003/07/28 19:04:15  camacho
  3760.  * Replaced all MemNews for calloc
  3761.  *
  3762.  * Revision 1.11  2003/07/28 03:41:49  camacho
  3763.  * Use f{open,close,gets} instead of File{Open,Close,Gets}
  3764.  *
  3765.  * Revision 1.10  2003/07/25 21:12:28  coulouri
  3766.  * remove constructions of the form "return sfree();" and "a=sfree(a);"
  3767.  *
  3768.  * Revision 1.9  2003/07/25 18:58:43  camacho
  3769.  * Avoid using StrUpper and StringHasNoText
  3770.  *
  3771.  * Revision 1.8  2003/07/25 17:25:43  coulouri
  3772.  * in progres:
  3773.  *  * use malloc/calloc/realloc instead of Malloc/Calloc/Realloc
  3774.  *  * add sfree() macro and __sfree() helper function to util.[ch]
  3775.  *  * use sfree() instead of MemFree()
  3776.  *
  3777.  * Revision 1.7  2003/07/24 22:37:33  dondosha
  3778.  * Removed some unused function parameters
  3779.  *
  3780.  * Revision 1.6  2003/07/24 22:01:44  camacho
  3781.  * Removed unused variables
  3782.  *
  3783.  * Revision 1.5  2003/07/24 21:31:06  dondosha
  3784.  * Changed to calls to BlastConstructErrorMessage to API from blast_message.h
  3785.  *
  3786.  * Revision 1.4  2003/07/24 20:38:30  dondosha
  3787.  * Removed LIBCALL etc. macros
  3788.  *
  3789.  * Revision 1.3  2003/07/24 17:37:46  dondosha
  3790.  * Removed MakeBlastScore function that is dependent on objalign.h
  3791.  *
  3792.  * Revision 1.2  2003/07/24 15:50:49  dondosha
  3793.  * Commented out mutex operations
  3794.  *
  3795.  * Revision 1.1  2003/07/24 15:18:09  dondosha
  3796.  * Copy of blastkar.h from ncbitools library, stripped of dependency on ncbiobj
  3797.  *
  3798.  * ===========================================================================
  3799.  */