blast_stat.h
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:14k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_stat.h,v $
  4.  * PRODUCTION Revision 1000.6  2004/06/01 18:03:52  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.45
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: blast_stat.h,v 1000.6 2004/06/01 18:03:52 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's official duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Author:  Tom Madden
  35.  *
  36.  */
  37. /** @file blast_stat.h
  38.  * Definitions and prototypes used by blast_stat.c to calculate BLAST
  39.  * statistics. @todo FIXME: needs doxygen comments
  40.  */
  41. #ifndef __BLAST_STAT__
  42. #define __BLAST_STAT__
  43. #include <algo/blast/core/blast_def.h>
  44. #include <algo/blast/core/blast_message.h>
  45. #ifdef __cplusplus
  46. extern "C" {
  47. #endif
  48. /*
  49. Defines for the matrix 'preferences' (as specified by S. Altschul).
  50. */
  51. #define BLAST_MATRIX_NOMINAL 0
  52. #define BLAST_MATRIX_PREFERRED 1
  53. #define BLAST_MATRIX_BEST 2
  54. /****************************************************************************
  55. For more accuracy in the calculation of K, set K_SUMLIMIT to 0.00001.
  56. For high speed in the calculation of K, use a K_SUMLIMIT of 0.001
  57. Note:  statistical significance is often not greatly affected by the value
  58. of K, so high accuracy is generally unwarranted.
  59. *****************************************************************************/
  60. /* K_SUMLIMIT_DEFAULT == sumlimit used in BlastKarlinLHtoK() */
  61. #define BLAST_KARLIN_K_SUMLIMIT_DEFAULT 0.0001
  62. /* LAMBDA_ACCURACY_DEFAULT == accuracy to which Lambda should be calc'd */
  63. #define BLAST_KARLIN_LAMBDA_ACCURACY_DEFAULT    (1.e-5)
  64. /* LAMBDA_ITER_DEFAULT == no. of iterations in LambdaBis = ln(accuracy)/ln(2)*/
  65. #define BLAST_KARLIN_LAMBDA_ITER_DEFAULT        17
  66. /* Initial guess for the value of Lambda in BlastKarlinLambdaNR */
  67. #define BLAST_KARLIN_LAMBDA0_DEFAULT    0.5
  68. #define BLAST_KARLIN_K_ITER_MAX 100
  69. #define BLAST_SUMP_EPSILON_DEFAULT 0.002 /* accuracy for SumP calculations */
  70. /* 
  71. Where are the BLAST matrices located?
  72. */
  73. #define BLASTMAT_DIR "/usr/ncbi/blast/matrix"
  74. /*************************************************************************
  75. Structure to the Karlin-Blk parameters.
  76. This structure was (more or less) copied from the old
  77. karlin.h.
  78. **************************************************************************/
  79. typedef struct Blast_KarlinBlk {
  80. double Lambda; /* Lambda value used in statistics */
  81. double K, logK; /* K value used in statistics */
  82. double H; /* H value used in statistics */
  83. double paramC; /* for use in seed. */
  84. } Blast_KarlinBlk;
  85. /********************************************************************
  86. Structures relating to scoring or the BlastScoreBlk
  87. ********************************************************************/
  88. /*
  89. SCORE_MIN is (-2**31 + 1)/2 because it has been observed more than once that
  90. a compiler did not properly calculate the quantity (-2**31)/2.  The list
  91. of compilers which failed this operation have at least at some time included:
  92. NeXT and a version of AIX/370's MetaWare High C R2.1r.
  93. For this reason, SCORE_MIN is not simply defined to be LONG_MIN/2.
  94. */
  95. #define BLAST_SCORE_MIN INT2_MIN
  96. #define BLAST_SCORE_MAX INT2_MAX
  97. #if defined(OS_DOS) || defined(OS_MAC)
  98. #define BLAST_SCORE_1MIN (-100)
  99. #define BLAST_SCORE_1MAX ( 100)
  100. #else
  101. #define BLAST_SCORE_1MIN (-10000)
  102. #define BLAST_SCORE_1MAX ( 10000)
  103. #endif
  104. #define BLAST_SCORE_RANGE_MAX (BLAST_SCORE_1MAX - BLAST_SCORE_1MIN)
  105. typedef struct Blast_ScoreFreq {
  106.     Int4 score_min, score_max;
  107.     Int4 obs_min, obs_max;
  108.     double score_avg;
  109.     double* sprob0,* sprob;
  110. } Blast_ScoreFreq;
  111. #define BLAST_MATRIX_SIZE 32
  112. /* Remove me */
  113. typedef struct SBLASTMatrixStructure {
  114.     Int4 *matrix[BLAST_MATRIX_SIZE];
  115.     Int4 long_matrix[BLAST_MATRIX_SIZE*BLAST_MATRIX_SIZE]; /* not used */
  116. } SBLASTMatrixStructure;
  117. typedef struct BlastScoreBlk {
  118. Boolean protein_alphabet; /* TRUE if alphabet_code is for a 
  119. protein alphabet (e.g., ncbistdaa etc.), FALSE for nt. alphabets. */
  120. Uint1 alphabet_code; /* NCBI alphabet code. */
  121. Int2  alphabet_size;  /* size of alphabet. */
  122. Int2  alphabet_start;  /* numerical value of 1st letter. */
  123. SBLASTMatrixStructure* matrix_struct; /* Holds info about matrix. */
  124. Int4 **matrix;  /* Substitution matrix */
  125. Int4 **posMatrix;  /* Sub matrix for position depend BLAST. */
  126.    double karlinK; /* Karlin-Altschul parameter associated with posMatrix */
  127. Int2 mat_dim1, mat_dim2; /* dimensions of matrix. */
  128. Int4 *maxscore; /* Max. score for each letter */
  129. Int4 loscore, hiscore; /* Min. & max. substitution scores */
  130. Int4 penalty, reward; /* penalty and reward for blastn. */
  131.         double  scale_factor; /* multiplier for all cutoff and dropoff scores */
  132. Boolean read_in_matrix; /* If TRUE, matrix is read in, otherwise
  133. produce one from penalty and reward above. */
  134. Blast_ScoreFreq** sfp; /* score frequencies. */
  135. double **posFreqs; /*matrix of position specific frequencies*/
  136. /* kbp & kbp_gap are ptrs that should be set to kbp_std, kbp_psi, etc. */
  137. Blast_KarlinBlk** kbp;  /* Karlin-Altschul parameters. */
  138. Blast_KarlinBlk** kbp_gap; /* K-A parameters for gapped alignments. */
  139. /* Below are the Karlin-Altschul parameters for non-position based ('std')
  140. and position based ('psi') searches. */
  141. Blast_KarlinBlk **kbp_std,
  142.                     **kbp_psi,
  143.                     **kbp_gap_std,
  144.                     **kbp_gap_psi;
  145. Blast_KarlinBlk*  kbp_ideal; /* Ideal values (for query with average database composition). */
  146. Int4 number_of_contexts; /* Used by sfp and kbp, how large are these*/
  147. char*  name; /* name of matrix. */
  148. Uint1*  ambiguous_res; /* Array of ambiguous res. (e.g, 'X', 'N')*/
  149. Int2 ambig_size, /* size of array above. */
  150. ambig_occupy; /* How many occupied? */
  151. ListNode* comments; /* Comments about matrix. */
  152. Int4     query_length;   /* the length of the query. */
  153. Int8 effective_search_sp; /* product of above two */
  154. } BlastScoreBlk;
  155. typedef struct Blast_ResFreq {
  156.     Uint1 alphabet_code;
  157.     double* prob; /* probs, (possible) non-zero offset. */
  158.     double* prob0; /* probs, zero offset. */
  159. } Blast_ResFreq;
  160. BlastScoreBlk* BlastScoreBlkNew (Uint1 alphabet, Int4 number_of_contexts);
  161. Int2 BlastScoreBlkMatrixLoad(BlastScoreBlk* sbp);
  162. BlastScoreBlk* BlastScoreBlkFree (BlastScoreBlk* sbp);
  163. Int2 BLAST_ScoreSetAmbigRes (BlastScoreBlk* sbp, char ambiguous_res);
  164. Int2 BLAST_ScoreBlkFill (BlastScoreBlk* sbp, char* string, Int4 length, Int4 context_number);
  165.  
  166. /** This function fills in the BlastScoreBlk structure.  
  167.  * Tasks are:
  168.  * -read in the matrix
  169.  * -set maxscore
  170.  * @param sbp Scoring block [in] [out]
  171.  * @param matrix Full path to the matrix in the directory structure [in]
  172. */
  173. Int2 BLAST_ScoreBlkMatFill (BlastScoreBlk* sbp, char* matrix);
  174.  
  175. /*
  176. Functions taken from the OLD karlin.c
  177. */
  178. Blast_KarlinBlk* Blast_KarlinBlkCreate (void);
  179. /** Deallocates the KarlinBlk
  180.  * @param kbp KarlinBlk to be deallocated [in]
  181. */
  182. Blast_KarlinBlk* Blast_KarlinBlkDestruct(Blast_KarlinBlk* kbp);
  183. Int2 Blast_KarlinBlkGappedCalc (Blast_KarlinBlk* kbp, Int4 gap_open, 
  184.         Int4 gap_extend, Int4 decline_align, char* matrix_name, 
  185.         Blast_Message** error_return);
  186. /** Calculates the standard Karlin parameters.  This is used
  187.  *       if the query is translated and the calculated (real) Karlin
  188.  *       parameters are bad, as they're calculated for non-coding regions.
  189.  * @param sbp ScoreBlk used to calculate "ideal" values. [in]
  190. */
  191. Blast_KarlinBlk* Blast_KarlinBlkIdealCalc(BlastScoreBlk* sbp);
  192. Int2 Blast_KarlinBlkStandardCalc(BlastScoreBlk* sbp, Int4 context_start, 
  193.                                  Int4 context_end);
  194. /*
  195.         Attempts to fill KarlinBlk for given gap opening, extensions etc.
  196.         Will return non-zero status if that fails.
  197.         return values:  -1 if matrix_name is NULL;
  198.                         1 if matrix not found
  199.                         2 if matrix found, but open, extend etc. values not supported.
  200. */
  201. Int2 Blast_KarlinkGapBlkFill(Blast_KarlinBlk* kbp, Int4 gap_open, Int4 gap_extend, Int4 decline_align, char* matrix_name);
  202. /* Prints a messages about the allowed matrices, BlastKarlinkGapBlkFill should return 1 before this is called. */
  203. char* BLAST_PrintMatrixMessage(const char *matrix);
  204. /* Prints a messages about the allowed open etc values for the given matrix, 
  205. BlastKarlinkGapBlkFill should return 2 before this is called. */
  206. char* BLAST_PrintAllowedValues(const char *matrix, Int4 gap_open, Int4 gap_extend, Int4 decline_align);
  207. /** Calculates the parameter Lambda given an initial guess for its value */
  208. double
  209. Blast_KarlinLambdaNR(Blast_ScoreFreq* sfp, double initialLambdaGuess);
  210. /** Calculates the Expect value based upon the search space and some Karlin-Altschul 
  211.  * parameters.  It is "simple" as it does not use sum-statistics.
  212.  * @param S the score of the alignment. [in]
  213.  * @param kbp the Karlin-Altschul parameters. [in]
  214.  * @param searchsp total search space to be used [in]
  215.  * @return the expect value
  216.  */
  217. double BLAST_KarlinStoE_simple (Int4 S, Blast_KarlinBlk* kbp, Int8  searchsp);
  218. double BLAST_GapDecayDivisor(double decayrate, unsigned nsegs );
  219. /** Calculate the cutoff score from the expected number of HSPs or vice versa.
  220.  * @param S The calculated score [in] [out]
  221.  * @param E The calculated e-value [in] [out]
  222.  * @param kbp The Karlin-Altschul statistical parameters [in]
  223.  * @param searchsp The effective search space [in]
  224.  * @param dodecay Use gap decay feature? [in]
  225.  * @param gap_decay_rate Gap decay rate to use, if dodecay is set [in]
  226.  */
  227. Int2 BLAST_Cutoffs (Int4 *S, double* E, Blast_KarlinBlk* kbp, 
  228.                     Int8 searchsp, Boolean dodecay, double gap_decay_rate);
  229. /* Functions to calculate SumE (for large and small gaps). */
  230. double BLAST_SmallGapSumE (Blast_KarlinBlk* kbp, Int4 gap, Int2 num,  double xsum, Int4 query_length, Int4 subject_length, double weight_divisor);
  231. double BLAST_UnevenGapSumE (Blast_KarlinBlk* kbp, Int4 p_gap, Int4 n_gap, Int2 num,  double xsum, Int4 query_length, Int4 subject_length, double weight_divisor);
  232. double BLAST_LargeGapSumE (Blast_KarlinBlk* kbp, Int2 num,  double xsum, Int4 query_length, Int4 subject_length, double weight_divisor );
  233. /*
  234. Obtains arrays of the allowed opening and extension penalties for gapped BLAST for
  235. the given matrix.  Also obtains arrays of Lambda, K, and H.  The pref_flags field is
  236. used for display purposes, with the defines: BLAST_MATRIX_NOMINAL, BLAST_MATRIX_PREFERRED, and
  237. BLAST_MATRIX_BEST.
  238. Any of these fields that
  239. are not required should be set to NULL.  The Int2 return value is the length of the
  240. arrays.
  241. */
  242. void BLAST_GetAlphaBeta (char* matrixName, double *alpha,
  243. double *beta, Boolean gapped, Int4 gap_open, Int4 gap_extend);
  244. Int4 ** RPSCalculatePSSM(double scalingFactor, Int4 rps_query_length, 
  245.                    Uint1 * rps_query_seq, Int4 db_seq_length, Int4 **posMatrix);
  246. Int4
  247. BLAST_ComputeLengthAdjustment(double K,
  248.                               double logK,
  249.                               double alpha_d_lambda,
  250.                               double beta,
  251.                               Int4 query_length,
  252.                               Int8 db_length,
  253.                               Int4 db_num_seqs,
  254.                               Int4 * length_adjustment);
  255. /** Allocates a new Blast_ResFreq structure and fills in the prob element
  256.     based upon the contents of sbp.
  257.  * @param sbp The BlastScoreBlk* used to init prob [in]
  258. */
  259. Blast_ResFreq* Blast_ResFreqNew(const BlastScoreBlk* sbp);
  260. /** Deallocates Blast_ResFreq and prob0 element.
  261.  * @param rfp the Blast_ResFreq to be deallocated.
  262. */
  263. Blast_ResFreq* Blast_ResFreqDestruct(Blast_ResFreq* rfp);
  264. /** Calculates residues frequencies given a standard distribution.
  265.  * @param sbp the BlastScoreBlk provides information on alphabet.
  266.  * @param rfp the prob element on this Blast_ResFreq is used.
  267. */
  268. Int2 Blast_ResFreqStdComp(const BlastScoreBlk* sbp, Blast_ResFreq* rfp);
  269. /** Creates a new structure to keep track of score frequencies for a scoring
  270.  * system.
  271.  * @param score_min Minimum score [in]
  272.  * @param score_max Maximum score [in]
  273.  */
  274. Blast_ScoreFreq*
  275. Blast_ScoreFreqNew(Int4 score_min, Int4 score_max);
  276. /** Deallocates the score frequencies structure 
  277.  * @param sfp the structure to deallocate [in]
  278.  * @return NULL
  279.  */
  280. Blast_ScoreFreq*
  281. Blast_ScoreFreqDestruct(Blast_ScoreFreq* sfp);
  282. /** Fills a buffer with the 'standard' alphabet 
  283.  * (given by STD_AMINO_ACID_FREQS[index].ch).
  284.  *
  285.  * @return Number of residues in alphabet or negative returns upon error.
  286.  */
  287. Int2
  288. Blast_GetStdAlphabet(Uint1 alphabet_code, Uint1* residues, 
  289.                      Uint4 residues_size);
  290. /* Please see comment on blast_stat.c */
  291. Int2
  292. Blast_KarlinBlkCalc(Blast_KarlinBlk* kbp, Blast_ScoreFreq* sfp);
  293. #ifdef __cplusplus
  294. }
  295. #endif
  296. #endif /* !__BLAST_STAT__ */