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

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

生物技术

开发平台：

C/C++

blast_psi_priv.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_psi_priv.c,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:07:34  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. static char const rcsid[] =
  11.     "$Id: blast_psi_priv.c,v 1000.1 2004/06/01 18:07:34 gouriano Exp $";
  12. /* ===========================================================================
  13.  *
  14.  *                            PUBLIC DOMAIN NOTICE
  15.  *               National Center for Biotechnology Information
  16.  *
  17.  *  This software/database is a "United States Government Work" under the
  18.  *  terms of the United States Copyright Act.  It was written as part of
  19.  *  the author's official duties as a United States Government employee and
  20.  *  thus cannot be copyrighted.  This software/database is freely available
  21.  *  to the public for use. The National Library of Medicine and the U.S.
  22.  *  Government have not placed any restriction on its use or reproduction.
  23.  *
  24.  *  Although all reasonable efforts have been taken to ensure the accuracy
  25.  *  and reliability of the software and data, the NLM and the U.S.
  26.  *  Government do not and cannot warrant the performance or results that
  27.  *  may be obtained by using this software or data. The NLM and the U.S.
  28.  *  Government disclaim all warranties, express or implied, including
  29.  *  warranties of performance, merchantability or fitness for any particular
  30.  *  purpose.
  31.  *
  32.  *  Please cite the author in any work or product based on this material.
  33.  *
  34.  * ===========================================================================
  35.  *
  36.  * Author:  Alejandro Schaffer, ported by Christiam Camacho
  37.  *
  38.  */
  40. /** @file blast_psi_priv.c
  41.  * Defintions for functions in private interface for Position Iterated BLAST 
  42.  * API.
  43.  */
  45. #include "blast_psi_priv.h"
  46. #include "matrix_freq_ratios.h"
  48. /****************************************************************************/
  49. /* Constants */
  50. const unsigned int kQueryIndex = 0;
  51. const double kEpsilon = 0.0001;
  52. const double kDefaultEvalueForPosition = 1.0;
  53. const int kPsiScaleFactor = 200;
  55. /****************************************************************************/
  57. void**
  58. _PSIAllocateMatrix(unsigned int ncols, unsigned int nrows, 
  59.                    unsigned int data_type_sz)
  60. {
  61.     void** retval = NULL;
  62.     unsigned int i = 0;
  64.     if ( !(retval = (void**) malloc(sizeof(void*) * ncols)))
  65.         return NULL;
  67.     for (i = 0; i < ncols; i++) {
  68.         if ( !(retval[i] = (void*) calloc(nrows, data_type_sz))) {
  69.             retval = _PSIDeallocateMatrix(retval, i);
  70.             break;
  71.         }
  72.     }
  73.     return retval;
  74. }
  76. void**
  77. _PSIDeallocateMatrix(void** matrix, unsigned int ncols)
  78. {
  79.     unsigned int i = 0;
  81.     if (!matrix)
  82.         return NULL;
  84.     for (i = 0; i < ncols; i++) {
  85.         sfree(matrix[i]);
  86.     }
  87.     sfree(matrix);
  88.     return NULL;
  89. }
  91. void
  92. _PSICopyMatrix(double** dest, const double** src, 
  93.                unsigned int ncols, unsigned int nrows)
  94. {
  95.     unsigned int i = 0;
  96.     unsigned int j = 0;
  98.     for (i = 0; i < ncols; i++) {
  99.         for (j = 0; j < nrows; j++) {
  100.             dest[i][j] = src[i][j];
  101.         }
  102.     }
  103. }
  105. /****************************************************************************/
  107. PsiAlignedBlock*
  108. _PSIAlignedBlockNew(Uint4 num_positions)
  109. {
  110.     PsiAlignedBlock* retval = NULL;
  111.     Uint4 i = 0;
  113.     if ( !(retval = (PsiAlignedBlock*) calloc(1, sizeof(PsiAlignedBlock)))) {
  114.         return NULL;
  115.     }
  117.     retval->pos_extnt = (SSeqRange*) calloc(num_positions, sizeof(SSeqRange));
  118.     if ( !retval->pos_extnt ) {
  119.         return _PSIAlignedBlockFree(retval);
  120.     }
  122.     retval->size = (Uint4*) calloc(num_positions, sizeof(Uint4));
  123.     if ( !retval->size ) {
  124.         return _PSIAlignedBlockFree(retval);
  125.     }
  127.     for (i = 0; i < num_positions; i++) {
  128.         retval->pos_extnt[i].left = -1;
  129.         retval->pos_extnt[i].right = num_positions;
  130.     }
  131.     return retval;
  132. }
  134. PsiAlignedBlock*
  135. _PSIAlignedBlockFree(PsiAlignedBlock* aligned_blocks)
  136. {
  137.     if ( !aligned_blocks ) {
  138.         return NULL;
  139.     }
  141.     if (aligned_blocks->pos_extnt) {
  142.         sfree(aligned_blocks->pos_extnt);
  143.     }
  145.     if (aligned_blocks->size) {
  146.         sfree(aligned_blocks->size);
  147.     }
  149.     sfree(aligned_blocks);
  150.     return NULL;
  151. }
  153. PsiSequenceWeights*
  154. _PSISequenceWeightsNew(const PsiInfo* info, const BlastScoreBlk* sbp)
  155. {
  156.     PsiSequenceWeights* retval = NULL;
  158.     ASSERT(info);
  160.     retval = (PsiSequenceWeights*) calloc(1, sizeof(PsiSequenceWeights));
  161.     if ( !retval ) {
  162.         return NULL;
  163.     }
  165.     retval->row_sigma = (double*) calloc(info->num_seqs + 1, sizeof(double));
  166.     if ( !retval->row_sigma ) {
  167.         return _PSISequenceWeightsFree(retval);
  168.     }
  170.     retval->norm_seq_weights = (double*) calloc(info->num_seqs + 1, 
  171.                                                 sizeof(double));
  172.     if ( !retval->norm_seq_weights ) {
  173.         return _PSISequenceWeightsFree(retval);
  174.     }
  176.     retval->sigma = (double*) calloc(info->num_seqs + 1, sizeof(double));
  177.     if ( !retval->sigma ) {
  178.         return _PSISequenceWeightsFree(retval);
  179.     }
  181.     retval->match_weights = (double**) _PSIAllocateMatrix(info->query_sz + 1,
  182.                                                           PSI_ALPHABET_SIZE,
  183.                                                           sizeof(double));
  184.     retval->match_weights_size = info->query_sz + 1;
  185.     if ( !retval->match_weights ) {
  186.         return _PSISequenceWeightsFree(retval);
  187.     }
  189.     retval->std_prob = _PSIGetStandardProbabilities(sbp);
  190.     if ( !retval->std_prob ) {
  191.         return _PSISequenceWeightsFree(retval);
  192.     }
  194.     retval->info_content = (double*) calloc(info->query_sz, sizeof(double));
  195.     if ( !retval->info_content ) {
  196.         return _PSISequenceWeightsFree(retval);
  197.     }
  199.     return retval;
  200. }
  202. PsiSequenceWeights*
  203. _PSISequenceWeightsFree(PsiSequenceWeights* seq_weights)
  204. {
  205.     if ( !seq_weights ) {
  206.         return NULL;
  207.     }
  209.     if (seq_weights->row_sigma) {
  210.         sfree(seq_weights->row_sigma);
  211.     }
  213.     if (seq_weights->norm_seq_weights) {
  214.         sfree(seq_weights->norm_seq_weights);
  215.     }
  217.     if (seq_weights->sigma) {
  218.         sfree(seq_weights->sigma);
  219.     }
  221.     if (seq_weights->match_weights) {
  222.         _PSIDeallocateMatrix((void**) seq_weights->match_weights,
  223.                              seq_weights->match_weights_size);
  224.     }
  226.     if (seq_weights->std_prob) {
  227.         sfree(seq_weights->std_prob);
  228.     }
  230.     if (seq_weights->info_content) {
  231.         sfree(seq_weights->info_content);
  232.     }
  234.     sfree(seq_weights);
  236.     return NULL;
  237. }
  239. /****************************************************************************/
  240. /* Function prototypes */
  242. /* Purges any aligned segments which are identical to the query sequence */
  243. static void
  244. _PSIPurgeIdenticalAlignments(PsiAlignmentData* alignment);
  246. /* Keeps only one copy of any aligned sequences which are >PSI_NEAR_IDENTICAL%
  247.  * identical to one another */
  248. static void
  249. _PSIPurgeNearIdenticalAlignments(PsiAlignmentData* alignment);
  250. static void
  251. _PSIUpdatePositionCounts(PsiAlignmentData* alignment);
  253. /* FIXME: needs more descriptive name */
  254. static void
  255. _PSIPurgeSimilarAlignments(PsiAlignmentData* alignment,
  256.                            Uint4 seq_index1,
  257.                            Uint4 seq_index2,
  258.                            double max_percent_identity);
  259. /****************************************************************************/
  261. /**************** PurgeMatches stage of PSSM creation ***********************/
  262. int
  263. PSIPurgeBiasedSegments(PsiAlignmentData* alignment)
  264. {
  265.     if ( !alignment ) {
  266.         return PSIERR_BADPARAM;
  267.     }
  269.     _PSIPurgeIdenticalAlignments(alignment);
  270.     _PSIPurgeNearIdenticalAlignments(alignment);
  271.     _PSIUpdatePositionCounts(alignment);
  273.     return PSI_SUCCESS;
  274. }
  276. /** Remove those sequences which are identical to the query sequence 
  277.  * FIXME: Rename to _PSIPurgeSelfHits() ?
  278.  */
  279. static void
  280. _PSIPurgeIdenticalAlignments(PsiAlignmentData* alignment)
  281. {
  282.     Uint4 s = 0;        /* index on sequences */
  284.     ASSERT(alignment);
  286.     for (s = 0; s < alignment->dimensions->num_seqs + 1; s++) {
  287.         _PSIPurgeSimilarAlignments(alignment, kQueryIndex, s, PSI_IDENTICAL);
  288.     }
  289. }
  291. static void
  292. _PSIPurgeNearIdenticalAlignments(PsiAlignmentData* alignment)
  293. {
  294.     Uint4 i = 0;
  295.     Uint4 j = 0;
  297.     ASSERT(alignment);
  299.     for (i = 1; i < alignment->dimensions->num_seqs + 1; i++) { 
  300.         for (j = 1; (i + j) < alignment->dimensions->num_seqs + 1; j++) {
  301.             _PSIPurgeSimilarAlignments(alignment, i, (i + j), 
  302.                                        PSI_NEAR_IDENTICAL);
  303.         }
  304.     }
  305. }
  307. /** Counts the number of sequences matching the query per query position
  308.  * (columns of the multiple alignment) as well as the number of residues
  309.  * present in each position of the query.
  310.  * Should be called after multiple alignment data has been purged from biased
  311.  * sequences.
  312.  */
  313. static void
  314. _PSIUpdatePositionCounts(PsiAlignmentData* alignment)
  315. {
  316.     Uint4 s = 0;        /* index on aligned sequences */
  317.     Uint4 p = 0;        /* index on positions */
  319.     ASSERT(alignment);
  321.     for (s = kQueryIndex + 1; s < alignment->dimensions->num_seqs + 1; s++) {
  323.         if ( !alignment->use_sequences[s] )
  324.             continue;
  326.         for (p = 0; p < alignment->dimensions->query_sz; p++) {
  327.             if (alignment->desc_matrix[s][p].used) {
  328.                 const Uint1 res = alignment->desc_matrix[s][p].letter;
  329.                 if (res >= PSI_ALPHABET_SIZE) {
  330.                     continue;
  331.                 }
  332.                 alignment->res_counts[p][res]++;
  333.                 alignment->match_seqs[p]++;
  334.             }
  335.         }
  336.     }
  337. }
  339. /** This function compares the sequences in the alignment->desc_matrix
  340.  * structure indexed by sequence_index1 and seq_index2. If it finds aligned 
  341.  * regions that have a greater percent identity than max_percent_identity, 
  342.  * it removes the sequence identified by seq_index2.
  343.  */
  344. static void
  345. _PSIPurgeSimilarAlignments(PsiAlignmentData* alignment,
  346.                            Uint4 seq_index1,
  347.                            Uint4 seq_index2,
  348.                            double max_percent_identity)
  349. {
  350.     Uint4 i = 0;
  352.     /* Nothing to do if sequences are the same or not selected for further
  353.        processing */
  354.     if ( seq_index1 == seq_index2 ||
  355.          !alignment->use_sequences[seq_index1] ||
  356.          !alignment->use_sequences[seq_index2] ) {
  357.         return;
  358.     }
  360.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  362.         const PsiDesc* seq1 = alignment->desc_matrix[seq_index1];
  363.         const PsiDesc* seq2 = alignment->desc_matrix[seq_index2];
  366.         /* starting index of the aligned region */
  367.         Uint4 align_start = i;    
  368.         /* length of the aligned region */
  369.         Uint4 align_length = 0;   
  370.         /* # of identical residues in aligned region */
  371.         Uint4 nidentical = 0;     
  373.         /* both positions in the sequences must be used */
  374.         if ( !(seq1[i].used && seq2[i].used) ) {
  375.             continue;
  376.         }
  378.         /* examine the aligned region (FIXME: should we care about Xs?) */
  379.         while ( (i < alignment->dimensions->query_sz) && 
  380.                 (seq1[i].used && seq2[i].used)) {
  381.             if (seq1[i].letter == seq2[i].letter)
  382.                 nidentical++;
  383.             align_length++;
  384.             i++;
  385.         }
  386.         ASSERT(align_length != 0);
  388.         /* percentage of similarity of an aligned region between seq1 and 
  389.            seq2 */
  390.         {
  391.             double percent_identity = (double) nidentical / align_length;
  393.                 fprintf(stderr, "%f%% for seqs %d and %dn",
  394.                         percent_identity, seq_index1, seq_index2);
  395.             if (percent_identity >= max_percent_identity) {
  396.                 int rv = _PSIPurgeAlignedRegion(alignment, seq_index2, 
  397.                                                 align_start, 
  398.                                                 align_start+align_length);
  399.                 ASSERT(rv == PSI_SUCCESS);
  400.             }
  401.         }
  402.     }
  403. }
  405. /****************************************************************************/
  406. /* Function prototypes */
  407. static void
  408. _PSIGetLeftExtents(const PsiAlignmentData* alignment, Uint4 seq_index);
  409. static void
  410. _PSIGetRightExtents(const PsiAlignmentData* alignment, Uint4 seq_index);
  412. static void
  413. _PSIComputePositionExtents(const PsiAlignmentData* alignment, 
  414.                            Uint4 seq_index,
  415.                            PsiAlignedBlock* aligned_blocks);
  416. static void
  417. _PSIComputeAlignedRegionLengths(const PsiAlignmentData* alignment,
  418.                                 PsiAlignedBlock* aligned_blocks);
  420. /****************************************************************************/
  421. /******* Compute alignment extents stage of PSSM creation *******************/
  422. /* posComputeExtents in posit.c */
  423. int
  424. PSIComputeAlignmentBlocks(const PsiAlignmentData* alignment,        /* [in] */
  425.                           PsiAlignedBlock* aligned_blocks)          /* [out] */
  426. {
  427.     Uint4 s = 0;     /* index on aligned sequences */
  429.     if ( !alignment || !aligned_blocks ) {
  430.         return PSIERR_BADPARAM;
  431.     }
  433.     /* no need to compute extents for query sequence */
  434.     for (s = kQueryIndex + 1; s < alignment->dimensions->num_seqs + 1; s++) {
  435.         if ( !alignment->use_sequences[s] )
  436.             continue;
  438.         _PSIGetLeftExtents(alignment, s);
  439.         _PSIGetRightExtents(alignment, s);
  440.         _PSIComputePositionExtents(alignment, s, aligned_blocks);
  441.     }
  443.     _PSIComputeAlignedRegionLengths(alignment, aligned_blocks);
  445.     return PSI_SUCCESS;
  446. }
  448. static void
  449. _PSIGetLeftExtents(const PsiAlignmentData* alignment, Uint4 seq_index)
  450. {
  451.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  452.     PsiDesc* sequence_position = NULL;
  453.     Uint4 prev = 0;  /* index for the first and previous position */
  454.     Uint4 curr = 0;  /* index for the current position */
  456.     ASSERT(alignment);
  457.     ASSERT(seq_index < alignment->dimensions->num_seqs + 1);
  458.     ASSERT(alignment->use_sequences[seq_index]);
  460.     sequence_position = alignment->desc_matrix[seq_index];
  462.     if (sequence_position[prev].used && sequence_position[prev].letter != GAP) {
  463.         sequence_position[prev].extents.left = prev;
  464.     }
  466.     for (curr = prev + 1; curr < alignment->dimensions->query_sz; 
  467.          curr++, prev++) {
  469.         if ( !sequence_position[curr].used ) {
  470.             continue;
  471.         }
  473.         if (sequence_position[prev].used) {
  474.             sequence_position[curr].extents.left =
  475.                 sequence_position[prev].extents.left;
  476.         } else {
  477.             sequence_position[curr].extents.left = curr;
  478.         }
  479.     }
  480. }
  482. static void
  483. _PSIGetRightExtents(const PsiAlignmentData* alignment, Uint4 seq_index)
  484. {
  485.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  486.     PsiDesc* sequence_position = NULL;
  487.     Uint4 last = 0;      /* index for the last position */
  488.     Uint4 curr = 0;      /* index for the current position */
  490.     ASSERT(alignment);
  491.     ASSERT(seq_index < alignment->dimensions->num_seqs + 1);
  492.     ASSERT(alignment->use_sequences[seq_index]);
  494.     sequence_position = alignment->desc_matrix[seq_index];
  495.     last = alignment->dimensions->query_sz - 1;
  497.     if (sequence_position[last].used && sequence_position[last].letter != GAP) {
  498.         sequence_position[last].extents.right = last;
  499.     }
  501.     for (curr = last - 1; curr >= 0; curr--, last--) {
  503.         if ( !sequence_position[curr].used ) {
  504.             continue;
  505.         }
  507.         if (sequence_position[last].used) {
  508.             sequence_position[curr].extents.right =
  509.                 sequence_position[last].extents.right;
  510.         } else {
  511.             sequence_position[curr].extents.right = curr;
  512.         }
  513.     }
  514. }
  516. static void
  517. _PSIComputePositionExtents(const PsiAlignmentData* alignment, 
  518.                            Uint4 seq_index,
  519.                            PsiAlignedBlock* aligned_blocks)
  520. {
  521. #ifdef PSI_IGNORE_GAPS_IN_COLUMNS
  522.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  523. #endif
  524.     PsiDesc* sequence_position = NULL;
  525.     Uint4 i = 0;
  527.     ASSERT(aligned_blocks);
  528.     ASSERT(alignment);
  529.     ASSERT(seq_index < alignment->dimensions->num_seqs + 1);
  530.     ASSERT(alignment->use_sequences[seq_index]);
  532.     sequence_position = alignment->desc_matrix[seq_index];
  534.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  535. #ifdef PSI_IGNORE_GAPS_IN_COLUMNS
  536.         if (sequence_position[i].used && 
  537.             sequence_position[i].letter != GAP) {
  538. #else
  539.         if (sequence_position[i].used) {
  540. #endif
  541.             aligned_blocks->pos_extnt[i].left = 
  542.                 MAX(aligned_blocks->pos_extnt[i].left, 
  543.                     sequence_position[i].extents.left);
  544.             aligned_blocks->pos_extnt[i].right = 
  545.                 MIN(aligned_blocks->pos_extnt[i].right, 
  546.                     sequence_position[i].extents.right);
  547.         }
  548.     }
  549. }
  551. static void
  552. _PSIComputeAlignedRegionLengths(const PsiAlignmentData* alignment,
  553.                                 PsiAlignedBlock* aligned_blocks)
  554. {
  555.     const Uint1 X = AMINOACID_TO_NCBISTDAA['X'];
  556.     PsiDesc* query_seq = NULL;
  557.     Uint4 i = 0;
  558.     
  559.     ASSERT(alignment);
  560.     ASSERT(aligned_blocks);
  562.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  563.         aligned_blocks->size[i] = aligned_blocks->pos_extnt[i].right - 
  564.                                    aligned_blocks->pos_extnt[i].left;
  565.     }
  567.     query_seq = alignment->desc_matrix[kQueryIndex];
  569.     /* Do not include X's in aligned region lengths */
  570.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  572.         if (query_seq[i].letter == X) {
  574.             Uint4 idx = 0;
  575.             for (idx = 0; idx < i; idx++) {
  576.                 if ((Uint4)aligned_blocks->pos_extnt[idx].right >= i &&
  577.                     query_seq[idx].letter != X) {
  578.                     aligned_blocks->size[idx]--;
  579.                 }
  580.             }
  581.             for (idx = alignment->dimensions->query_sz - 1; idx > i; idx--) {
  582.                 if ((Uint4)aligned_blocks->pos_extnt[idx].left <= i &&
  583.                     query_seq[idx].letter != X) {
  584.                     aligned_blocks->size[idx]--;
  585.                 }
  586.             }
  587.         }
  589.     }
  590. }
  592. /****************************************************************************/
  593. static Uint4
  594. _PSIGetAlignedSequencesForPosition(
  595.     const PsiAlignmentData* alignment, 
  596.     Uint4 position,
  597.     Uint4* aligned_sequences);
  599. static void
  600. _PSICalculatePositionWeightsAndIntervalSigmas(
  601.     const PsiAlignmentData* alignment,     /* [in] */
  602.     const PsiAlignedBlock* aligned_blocks, /* [in] */
  603.     Uint4 position,                        /* [in] */
  604.     const Uint4* aligned_seqs,             /* [in] */
  605.     Uint4 num_aligned_seqs,                /* [in] */
  606.     PsiSequenceWeights* seq_weights,       /* [out] */
  607.     double* sigma,                         /* [out] */
  608.     double* interval_sigma);               /* [out] */
  610. static void
  611. _PSICalculateNormalizedSequenceWeights(
  612.     const PsiAlignedBlock* aligned_blocks, /* [in] */
  613.     Uint4 position,                        /* [in] */
  614.     const Uint4* aligned_seqs,             /* [in] */
  615.     Uint4 num_aligned_seqs,                /* [in] */
  616.     double sigma,                          /* [in] */
  617.     PsiSequenceWeights* seq_weights);      /* [out] */
  619. static void
  620. _PSICalculateMatchWeights(
  621.     const PsiAlignmentData* alignment,  /* [in] */
  622.     Uint4 position,                     /* [in] */
  623.     const Uint4* aligned_seqs,          /* [in] */
  624.     Uint4 num_aligned_seqs,             /* [in] */
  625.     PsiSequenceWeights* seq_weights);   /* [out] */
  627. static int
  628. _PSICheckSequenceWeights(
  629.     const PsiAlignmentData* alignment,         /* [in] */
  630.     PsiSequenceWeights* seq_weights);          /* [in] */
  632. /****************************************************************************/
  633. /******* Calculate sequence weights stage of PSSM creation ******************/
  634. /* Needs the PsiAlignedBlock structure calculated in previous stage as well
  635.  * as PsiAlignmentData structure */
  637. int
  638. PSIComputeSequenceWeights(const PsiAlignmentData* alignment,        /* [in] */
  639.                           const PsiAlignedBlock* aligned_blocks,    /* [in] */
  640.                           PsiSequenceWeights* seq_weights)          /* [out] */
  641. {
  642.     Uint4* aligned_seqs = NULL;     /* list of indices of sequences
  643.                                               which participate in an
  644.                                               aligned position */
  645.     Uint4 pos = 0;                  /* position index */
  647.     ASSERT(alignment);
  648.     ASSERT(aligned_blocks);
  649.     ASSERT(seq_weights);
  651.     aligned_seqs = (Uint4*) calloc(alignment->dimensions->num_seqs + 1,
  652.                                           sizeof(Uint4));
  653.     if ( !aligned_seqs ) {
  654.         return PSIERR_OUTOFMEM;
  655.     }
  657.     for (pos = 0; pos < alignment->dimensions->query_sz; pos++) {
  659.         Uint4 num_aligned_seqs = 0;
  660.         double sigma = 0.0; /**< number of different characters
  661.                                   occurring in matches within a multiple
  662.                                   alignment block, excluding identical
  663.                                   columns */
  664.         double interval_sigma = 0.0; /**< same as sigma but includes identical
  665.                                   columns */
  667.         /* ignore positions of no interest */
  668.         if (aligned_blocks->size[pos] == 0 || alignment->match_seqs[pos] <= 1) {
  669.             continue;
  670.         }
  672.         memset((void*)aligned_seqs, 0, 
  673.                sizeof(Uint4) * (alignment->dimensions->num_seqs + 1));
  675.         num_aligned_seqs = _PSIGetAlignedSequencesForPosition(alignment, pos,
  676.                                                               aligned_seqs);
  677.         if (num_aligned_seqs <= 1) {
  678.             continue;
  679.         }
  681.         /* Skipping optimization about redundant sets */
  683.         /* if (newSequenceSet) in posit.c */
  684.         memset((void*)seq_weights->norm_seq_weights, 0, 
  685.                sizeof(double)*(alignment->dimensions->num_seqs+1));
  686.         memset((void*)seq_weights->row_sigma, 0,
  687.                sizeof(double)*(alignment->dimensions->num_seqs+1));
  688.         _PSICalculatePositionWeightsAndIntervalSigmas(alignment, 
  689.                                      aligned_blocks, pos, aligned_seqs, 
  690.                                      num_aligned_seqs, seq_weights, 
  691.                                      &sigma, &interval_sigma);
  693.         seq_weights->sigma[pos] = interval_sigma;
  695.         /* Populates norm_seq_weights */
  696.         _PSICalculateNormalizedSequenceWeights(aligned_blocks, pos, 
  697.                                  aligned_seqs, num_aligned_seqs, sigma,
  698.                                  seq_weights);
  701.         /* Uses seq_weights->norm_seq_weights to populate match_weights */
  702.         _PSICalculateMatchWeights(alignment, pos, aligned_seqs,
  703.                                   num_aligned_seqs, seq_weights);
  704.     }
  706.     sfree(aligned_seqs);
  708.     /* Check that the sequence weights add up to 1 in each column */
  709.     _PSICheckSequenceWeights(alignment, seq_weights);
  711.     /* Return seq_weights->match_weigths, should free others? FIXME: need to
  712.      * keep sequence weights for diagnostics for structure group */
  713.     return PSI_SUCCESS;
  714. }
  716. /* Calculates the position based weights using a modified version of the
  717.  * Henikoff's algorithm presented in "Position-based sequence weights".
  718.  * More documentation pending */
  719. static void
  720. _PSICalculatePositionWeightsAndIntervalSigmas(
  721.     const PsiAlignmentData* alignment,     /* [in] */
  722.     const PsiAlignedBlock* aligned_blocks, /* [in] */
  723.     Uint4 position,                        /* [in] */
  724.     const Uint4* aligned_seqs,             /* [in] */
  725.     Uint4 num_aligned_seqs,                /* [in] */
  726.     PsiSequenceWeights* seq_weights,       /* [out] */
  727.     double* sigma,                         /* [out] */
  728.     double* interval_sigma)                /* [out] */
  729. {
  730.     /** keeps track of how many occurrences of each residue was found at a
  731.      * given position */
  732.     Uint4 residue_counts[PSI_ALPHABET_SIZE];
  733.     Uint4 num_distinct_residues = 0; /**< number of distinct 
  734.                                               residues found at a given 
  735.                                               position */
  737.     Uint4 i = 0;         /**< index into aligned block for requested
  738.                                    position */
  740.     ASSERT(seq_weights);
  741.     ASSERT(sigma);
  742.     ASSERT(interval_sigma);
  744.     *sigma = 0.0;
  745.     *interval_sigma = 0.0;
  747.     for (i = 0; i < sizeof(residue_counts); i++) {
  748.         residue_counts[i] = 0;
  749.     }
  751.     for (i = (Uint4) aligned_blocks->pos_extnt[position].left; 
  752.          i <= (Uint4) aligned_blocks->pos_extnt[position].right; i++) {
  754.         Uint4 asi = 0;   /**< index into array of aligned sequences */
  756.         /**** Count number of occurring residues at a position ****/
  757.         for (asi = 0; asi < num_aligned_seqs; asi++) {
  758.             const Uint4 seq_idx = aligned_seqs[asi];
  759.             const Uint1 residue = 
  760.                 alignment->desc_matrix[seq_idx][i].letter;
  762.             if (residue_counts[residue] == 0) {
  763.                 num_distinct_residues++;
  764.             }
  765.             residue_counts[residue]++;
  766.         }
  767.         /**** END: Count number of occurring residues at a position ****/
  769.         /* FIXME: see Alejandro's email about this */
  770.         (*interval_sigma) += num_distinct_residues;
  771.         if (num_distinct_residues > 1) {    /* if this is not 1 */
  772.             (*sigma) += num_distinct_residues;
  773.         }
  775.         /* Calculate row_sigma */
  776.         for (asi = 0; asi < num_aligned_seqs; asi++) {
  777.             const Uint4 seq_idx = aligned_seqs[asi];
  778.             const Uint1 residue =
  779.                 alignment->desc_matrix[seq_idx][i].letter;
  781.             seq_weights->row_sigma[seq_idx] += 
  782.                 (1.0 / (double) 
  783.                  (residue_counts[residue] * num_distinct_residues) );
  784.         }
  785.     }
  787.     return;
  788. }
  790. /** Calculates the normalized sequence weights for the requested position */
  791. static void
  792. _PSICalculateNormalizedSequenceWeights(
  793.     const PsiAlignedBlock* aligned_blocks, /* [in] */
  794.     Uint4 position,                        /* [in] */
  795.     const Uint4* aligned_seqs,             /* [in] */
  796.     Uint4 num_aligned_seqs,                /* [in] */
  797.     double sigma,                          /* [in] */
  798.     PsiSequenceWeights* seq_weights)       /* [out] */
  799. {
  800.     Uint4 asi = 0;   /**< index into array of aligned sequences */
  802.     if (sigma > 0) {
  803.         double weight_sum = 0.0;
  805.         for (asi = 0; asi < num_aligned_seqs; asi++) {
  806.             const Uint4 seq_idx = aligned_seqs[asi];
  807.             seq_weights->norm_seq_weights[seq_idx] = 
  808.                 seq_weights->row_sigma[seq_idx] / 
  809.                 (aligned_blocks->pos_extnt[position].right -
  810.                  aligned_blocks->pos_extnt[position].left + 1);
  811. #ifndef PSI_IGNORE_GAPS_IN_COLUMNS
  812.             weight_sum += seq_weights->norm_seq_weights[seq_idx];
  813. #endif
  814.         }
  816.         for (asi = 0; asi < num_aligned_seqs; asi++) {
  817.             const Uint4 seq_idx = aligned_seqs[asi];
  818.             seq_weights->norm_seq_weights[seq_idx] /= weight_sum;
  819.         }
  821.     } else {
  822.         for (asi = 0; asi < num_aligned_seqs; asi++) {
  823.             const Uint4 seq_idx = aligned_seqs[asi];
  824.             seq_weights->norm_seq_weights[seq_idx] = 
  825.                 (1.0/(double) num_aligned_seqs);
  826.         }
  827.     }
  829. }
  831. static void
  832. _PSICalculateMatchWeights(
  833.     const PsiAlignmentData* alignment,  /* [in] */
  834.     Uint4 position,                     /* [in] */
  835.     const Uint4* aligned_seqs,          /* [in] */
  836.     Uint4 num_aligned_seqs,             /* [in] */
  837.     PsiSequenceWeights* seq_weights)    /* [out] */
  838. {
  839.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  840.     Uint4 asi = 0;   /**< index into array of aligned sequences */
  842.     for (asi = 0; asi < num_aligned_seqs; asi++) {
  843.         const Uint4 seq_idx = aligned_seqs[asi];
  844.         const Uint1 residue =
  845.             alignment->desc_matrix[seq_idx][position].letter;
  847.         seq_weights->match_weights[position][residue] += 
  848.             seq_weights->norm_seq_weights[seq_idx];
  850.         /* FIXME: this field is populated but never used */
  851.         if (residue == GAP) {
  852.             /*seq_weights->gapless_column_weights[position] +=
  853.              * seq_weights->a[seq_idx]; */
  854.             ;
  855.         }
  856.     }
  857. }
  859. /** Finds the sequences aligned in a given position.
  860.  * @param alignment Multiple-alignment data [in]
  861.  * @param position position of interest [in]
  862.  * @param aligned_sequences array which will contain the indices of the
  863.  * sequences aligned at the requested position. This array must have size
  864.  * greater than or equal to the number of sequences + 1 in multiple alignment 
  865.  * data structure (alignment->dimensions->num_seqs + 1) [out]
  866.  * @return number of sequences aligned at the requested position
  867.  */
  868. static Uint4
  869. _PSIGetAlignedSequencesForPosition(const PsiAlignmentData* alignment, 
  870.                                    Uint4 position,
  871.                                    Uint4* aligned_sequences)
  872. {
  873. #ifdef PSI_IGNORE_GAPS_IN_COLUMNS
  874.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  875. #endif
  876.     Uint4 retval = 0;
  877.     Uint4 i = 0;
  879.     ASSERT(alignment);
  880.     ASSERT(position < alignment->dimensions->query_sz);
  881.     ASSERT(aligned_sequences);
  883.     for (i = 0; i < alignment->dimensions->num_seqs + 1; i++) {
  885.         if ( !alignment->use_sequences[i] ) {
  886.             continue;
  887.         }
  889. #ifdef PSI_IGNORE_GAPS_IN_COLUMNS
  890.         if (alignment->desc_matrix[i][position].used &&
  891.             alignment->desc_matrix[i][position].letter != GAP) {
  892. #else
  893.         if (alignment->desc_matrix[i][position].used) {
  894. #endif
  895.             aligned_sequences[retval++] = i;
  896.         }
  897.     }
  899.     return retval;
  900. }
  902. /* The second parameter is not really const, it's updated! */
  903. static int
  904. _PSICheckSequenceWeights(const PsiAlignmentData* alignment,
  905.                          PsiSequenceWeights* seq_weights)
  906. {
  907.     const Uint1 GAP = AMINOACID_TO_NCBISTDAA['-'];
  908.     const Uint1 X = AMINOACID_TO_NCBISTDAA['X'];
  909.     Uint4 pos = 0;   /* residue position (ie: column number) */
  910.     Uint4 res = 0;   /* residue */
  912.     ASSERT(alignment);
  913.     ASSERT(seq_weights);
  915.     for (pos = 0; pos < alignment->dimensions->query_sz; pos++) {
  917.         double running_total = 0.0;
  919.         if (alignment->match_seqs[pos] <= 1 ||
  920.             alignment->desc_matrix[kQueryIndex][pos].letter == X) {
  921.             continue;
  922.         }
  924.         for (res = 0; res < PSI_ALPHABET_SIZE; res++) {
  925.             running_total += seq_weights->match_weights[pos][res];
  926.         }
  927.         ASSERT(running_total < 0.99 || running_total > 1.01);
  929. #ifndef PSI_IGNORE_GAPS_IN_COLUMNS
  930.         /* Disperse method of spreading the gap weights */
  931.         for (res = 0; res < PSI_ALPHABET_SIZE; res++) {
  932.             if (seq_weights->std_prob[res] > kEpsilon) {
  933.                 seq_weights->match_weights[pos][res] += 
  934.                     (seq_weights->match_weights[pos][GAP] * 
  935.                      seq_weights->std_prob[res]);
  936.             }
  937.         }
  938. #endif
  939.         seq_weights->match_weights[pos][GAP] = 0.0;
  940.         running_total = 0.0;
  941.         for (res = 0; res < PSI_ALPHABET_SIZE; res++) {
  942.             running_total += seq_weights->match_weights[pos][res];
  943.         }
  945.         if (running_total < 0.99 || running_total > 1.01) {
  946.             return PSIERR_BADSEQWEIGHTS;
  947.         }
  948.     }
  950.     return PSI_SUCCESS;
  951. }
  953. /****************************************************************************/
  954. /******* Compute residue frequencies stage of PSSM creation *****************/
  955. /* port of posComputePseudoFreqs */
  956. int
  957. PSIComputeResidueFrequencies(const PsiAlignmentData* alignment,     /* [in] */
  958.                              const PsiSequenceWeights* seq_weights, /* [in] */
  959.                              const BlastScoreBlk* sbp,              /* [in] */
  960.                              const PsiAlignedBlock* aligned_blocks, /* [in] */
  961.                              const PSIBlastOptions* opts,           /* [in] */
  962.                              PsiMatrix* score_matrix)               /* [out] */
  963. {
  964.     const Uint1 X = AMINOACID_TO_NCBISTDAA['X'];
  965.     Uint4 i = 0;             /* loop index into query positions */
  966.     SFreqRatios* freq_ratios;       /* matrix-specific frequency ratios */
  968.     if ( !alignment || !seq_weights || !sbp || 
  969.          !aligned_blocks || !opts || !score_matrix ) {
  970.         return PSIERR_BADPARAM;
  971.     }
  973.     freq_ratios = _PSIMatrixFrequencyRatiosNew(sbp->name);
  975.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  977.         Uint4 j = 0;     /* loop index into alphabet */
  978.         double info_sum = 0.0;  /* for information content - FIXME calculate
  979.                                    separately */
  981.         /* If there is an 'X' in the query sequence at position i... */
  982.         if (alignment->desc_matrix[kQueryIndex][i].letter == X) {
  984.             for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  985.                 score_matrix->res_freqs[i][j] = 0.0;
  986.             }
  988.         } else {
  990.             for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  992.                 if (seq_weights->std_prob[j] > kEpsilon) {
  993.                     Uint4 interval_size = 0; /* length of a block */
  994.                     Uint4 idx = 0;  /* loop index */
  995.                     double sigma = 0.0;    /* number of chars in an interval */
  997.                     double pseudo = 0.0;            /* intermediate term */
  998.                     double numerator = 0.0;         /* intermediate term */
  999.                     double denominator = 0.0;       /* intermediate term */
  1000.                     double qOverPEstimate = 0.0;    /* intermediate term */
  1002.                     /* changed to matrix specific ratios here May 2000 */
  1003.                     for (idx = 0; idx < (Uint4) sbp->alphabet_size; idx++) {
  1004.                         if (sbp->matrix[j][idx] != BLAST_SCORE_MIN) {
  1005.                             pseudo += (seq_weights->match_weights[i][idx] *
  1006.                                        freq_ratios->data[j][idx]);
  1007.                         }
  1008.                     }
  1009.                     pseudo *= opts->pseudo_count;
  1011.                     /* FIXME: document where this formula is coming from
  1012.                      * (probably 2001 paper, p 2996) */
  1013.                     sigma = seq_weights->sigma[i];
  1014.                     interval_size = aligned_blocks->size[i];
  1016.                     numerator = pseudo + 
  1017.                         ((sigma/interval_size-1) * 
  1018.                          seq_weights->match_weights[i][j] / 
  1019.                          seq_weights->std_prob[j]);
  1021.                     denominator = (sigma/interval_size-1) +
  1022.                         opts->pseudo_count;
  1024.                     qOverPEstimate = numerator/denominator;
  1026.                     /* Note artificial multiplication by standard probability
  1027.                      * to normalize */
  1028.                     score_matrix->res_freqs[i][j] = qOverPEstimate *
  1029.                         seq_weights->std_prob[j];
  1031.                     if ( qOverPEstimate != 0.0 &&
  1032.                          (seq_weights->std_prob[j] > kEpsilon) ) {
  1033.                         info_sum += qOverPEstimate * seq_weights->std_prob[j] *
  1034.                             log(qOverPEstimate)/NCBIMATH_LN2;
  1035.                     }
  1037.                 } else {
  1038.                     score_matrix->res_freqs[i][j] = 0.0;
  1039.                 } /* END: if (seq_weights->std_prob[j] > kEpsilon) { */
  1040.             } /* END: for (j = 0; j < sbp->alphabet_size; j++) */
  1042.         }
  1043.         /* FIXME: Should move out the calculation of information content to its
  1044.          * own function (see posFreqsToInformation)! */
  1045.         seq_weights->info_content[i] = info_sum;
  1046.     }
  1048.     freq_ratios = _PSIMatrixFrequencyRatiosFree(freq_ratios);
  1050.     return PSI_SUCCESS;
  1051. }
  1053. /****************************************************************************/
  1054. /**************** Convert residue frequencies to PSSM stage *****************/
  1056. /* FIXME: Answer questions
  1057.    FIXME: need ideal_labmda, regular scoring matrix, length of query
  1058. */
  1059. int
  1060. PSIConvertResidueFreqsToPSSM(PsiMatrix* score_matrix,           /* [in|out] */
  1061.                              const Uint1* query,                /* [in] */
  1062.                              const BlastScoreBlk* sbp,          /* [in] */
  1063.                              const double* std_probs)           /* [in] */
  1064. {
  1065.     const Uint4 X = AMINOACID_TO_NCBISTDAA['X'];
  1066.     const Uint4 Star = AMINOACID_TO_NCBISTDAA['*'];
  1067.     Uint4 i = 0;
  1068.     Uint4 j = 0;
  1069.     SFreqRatios* std_freq_ratios = NULL;    /* only needed when there are not
  1070.                                                residue frequencies for a given 
  1071.                                                column */
  1072.     double ideal_lambda;
  1074.     if ( !score_matrix || !sbp || !std_probs )
  1075.         return PSIERR_BADPARAM;
  1077.     std_freq_ratios = _PSIMatrixFrequencyRatiosNew(sbp->name);
  1078.     ideal_lambda = sbp->kbp_ideal->Lambda;
  1080.     /* Each column is a position in the query */
  1081.     for (i = 0; i < score_matrix->ncols; i++) {
  1083.         /* True if all frequencies in column i are zero */
  1084.         Boolean is_unaligned_column = TRUE;
  1085.         const Uint4 query_res = query[i];
  1087.         for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  1089.             double qOverPEstimate = 0.0;
  1091.             /* Division compensates for multiplication in previous function */
  1092.             if (std_probs[j] > kEpsilon) {
  1093.                 qOverPEstimate = 
  1094.                     score_matrix->res_freqs[i][j] / std_probs[j];
  1095.             }
  1097.             if (is_unaligned_column && qOverPEstimate != 0.0) {
  1098.                 is_unaligned_column = FALSE;
  1099.             }
  1101.             /* Populate scaled matrix */
  1102.             if (qOverPEstimate == 0.0 || std_probs[j] < kEpsilon) {
  1103.                 score_matrix->scaled_pssm[i][j] = BLAST_SCORE_MIN;
  1104.             } else {
  1105.                 double tmp = log(qOverPEstimate)/ideal_lambda;
  1106.                 score_matrix->scaled_pssm[i][j] = (int)
  1107.                     BLAST_Nint(kPsiScaleFactor * tmp);
  1108.             }
  1110.             if ( (j == X || j == Star) &&
  1111.                  (sbp->matrix[query_res][X] != BLAST_SCORE_MIN) ) {
  1112.                 score_matrix->scaled_pssm[i][j] = 
  1113.                     sbp->matrix[query_res][j] * kPsiScaleFactor;
  1114.             }
  1115.         }
  1117.         if (is_unaligned_column) {
  1118.             for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  1120.                 score_matrix->pssm[i][j] = sbp->matrix[query_res][j];
  1122.                 if (sbp->matrix[query_res][j] != BLAST_SCORE_MIN) {
  1123.                     double tmp = 
  1124.                         kPsiScaleFactor * std_freq_ratios->bit_scale_factor *
  1125.                         log(std_freq_ratios->data[query_res][j])/NCBIMATH_LN2;
  1127.                     score_matrix->scaled_pssm[i][j] = BLAST_Nint(tmp);
  1128.                 } else {
  1129.                     score_matrix->scaled_pssm[i][j] = BLAST_SCORE_MIN;
  1130.                 }
  1131.             }
  1132.         }
  1133.     }
  1135.     std_freq_ratios = _PSIMatrixFrequencyRatiosFree(std_freq_ratios);
  1137.     /* Set the last column of the matrix to BLAST_SCORE_MIN (why?) */
  1138.     for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  1139.         score_matrix->scaled_pssm[score_matrix->ncols-1][j] = BLAST_SCORE_MIN;
  1140.     }
  1142.     return PSI_SUCCESS;
  1143. }
  1145. /****************************************************************************/
  1146. /************************* Scaling of PSSM stage ****************************/
  1148. /**
  1149.  * @param initial_lambda_guess value to be used when calculating lambda if this
  1150.  * is not null [in]
  1151.  * @param sbp Score block structure where the calculated lambda and K will be
  1152.  * returned
  1153.  */
  1154. void
  1155. _PSIUpdateLambdaK(const int** pssm,              /* [in] */
  1156.                   const Uint1* query,            /* [in] */
  1157.                   Uint4 query_length,            /* [in] */
  1158.                   const double* std_probs,       /* [in] */
  1159.                   double* initial_lambda_guess,  /* [in] */
  1160.                   BlastScoreBlk* sbp);           /* [in|out] */
  1162. /* FIXME: change so that only lambda is calculated inside the loop that scales
  1163.    the matrix and kappa is calculated before returning from this function.
  1164.    Scaling factor should be optional argument to accomodate kappa.c's needs?
  1165. */
  1166. int
  1167. PSIScaleMatrix(const Uint1* query,              /* [in] */
  1168.                Uint4 query_length,              /* [in] */
  1169.                const double* std_probs,         /* [in] */
  1170.                double* scaling_factor,          /* [in] */
  1171.                PsiMatrix* score_matrix,         /* [in|out] */
  1172.                BlastScoreBlk* sbp)              /* [in|out] */
  1173. {
  1174.     Boolean first_time = TRUE;
  1175.     Uint4 index = 0;     /* loop index */
  1176.     int** scaled_pssm = NULL;
  1177.     int** pssm = NULL;
  1178.     double factor;
  1179.     double factor_low;
  1180.     double factor_high;
  1181.     double new_lambda;      /* Karlin-Altschul parameter */
  1183.     const double kPositPercent = 0.05;
  1184.     const Uint4 kPositNumIterations = 10;
  1185.     Boolean too_high = TRUE;
  1186.     double ideal_lambda;
  1188.     if ( !score_matrix || !sbp || !query || !std_probs )
  1189.         return PSIERR_BADPARAM;
  1191.     ASSERT(sbp->kbp_psi[0]);
  1193.     scaled_pssm = score_matrix->scaled_pssm;
  1194.     pssm = score_matrix->pssm;
  1195.     ideal_lambda = sbp->kbp_ideal->Lambda;
  1197.     /* FIXME: need to take scaling_factor into account */
  1199.     factor = 1.0;
  1200.     while (1) {
  1201.         Uint4 i = 0;
  1202.         Uint4 j = 0;
  1204.         for (i = 0; i < score_matrix->ncols; i++) {
  1205.             for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  1206.                 if (scaled_pssm[i][j] != BLAST_SCORE_MIN) {
  1207.                     pssm[i][j] = 
  1208.                         BLAST_Nint(factor*scaled_pssm[i][j]/kPsiScaleFactor);
  1209.                 } else {
  1210.                     pssm[i][j] = BLAST_SCORE_MIN;
  1211.                 }
  1212.             }
  1213.         }
  1215.         if (scaling_factor) {
  1216.             double init_lambda_guess = 
  1217.                 sbp->kbp_psi[0]->Lambda / *scaling_factor;
  1218.             _PSIUpdateLambdaK((const int**)pssm, query, query_length, 
  1219.                               std_probs, &init_lambda_guess, sbp);
  1220.         } else {
  1221.             _PSIUpdateLambdaK((const int**)pssm, query, query_length, 
  1222.                               std_probs, NULL, sbp);
  1223.         }
  1225.         new_lambda = sbp->kbp_psi[0]->Lambda;
  1227.         if (new_lambda > ideal_lambda) {
  1228.             if (first_time) {
  1229.                 factor_high = 1.0 + kPositPercent;
  1230.                 factor = factor_high;
  1231.                 too_high = TRUE;
  1232.                 first_time = FALSE;
  1233.             } else {
  1234.                 if (too_high == FALSE) {
  1235.                     break;
  1236.                 }
  1237.                 factor_high += (factor_high - 1.0);
  1238.                 factor = factor_high;
  1239.             }
  1240.         } else if (new_lambda > 0) {
  1241.             if (first_time) {
  1242.                 factor_high = 1.0;
  1243.                 factor_low = 1.0 - kPositPercent;
  1244.                 factor = factor_low;
  1245.                 too_high = FALSE;
  1246.                 first_time = FALSE;
  1247.             } else {
  1248.                 if (too_high == TRUE) {
  1249.                     break;
  1250.                 }
  1251.                 factor_low += (factor_low - 1.0);
  1252.                 factor = factor_low;
  1253.             }
  1254.         } else {
  1255.             return PSIERR_POSITIVEAVGSCORE;
  1256.         }
  1257.     }
  1259.     /* Binary search for kPositNumIterations times */
  1260.     for (index = 0; index < kPositNumIterations; index++) {
  1261.         Uint4 i = 0;
  1262.         Uint4 j = 0;
  1264.         factor = (factor_high + factor_low)/2;
  1266.         for (i = 0; i < score_matrix->ncols; i++) {
  1267.             for (j = 0; j < (Uint4) sbp->alphabet_size; j++) {
  1268.                 if (scaled_pssm[i][j] != BLAST_SCORE_MIN) {
  1269.                     pssm[i][j] = 
  1270.                         BLAST_Nint(factor*scaled_pssm[i][j]/kPsiScaleFactor);
  1271.                 } else {
  1272.                     pssm[i][j] = BLAST_SCORE_MIN;
  1273.                 }
  1274.             }
  1275.         }
  1277.         if (scaling_factor) {
  1278.             double init_lambda_guess = 
  1279.                 sbp->kbp_psi[0]->Lambda / *scaling_factor;
  1280.             _PSIUpdateLambdaK((const int**)pssm, query, query_length, 
  1281.                               std_probs, &init_lambda_guess, sbp);
  1282.         } else {
  1283.             _PSIUpdateLambdaK((const int**)pssm, query, query_length, 
  1284.                               std_probs, NULL, sbp);
  1285.         }
  1287.         new_lambda = sbp->kbp_psi[0]->Lambda;
  1289.         if (new_lambda > ideal_lambda) {
  1290.             factor_low = factor;
  1291.         } else {
  1292.             factor_high = factor;
  1293.         }
  1294.     }
  1296.     /* FIXME: Why is this needed? */
  1297.     for (index = 0; index < (Uint4) sbp->alphabet_size; index++) {
  1298.         pssm[score_matrix->ncols-1][index] = BLAST_SCORE_MIN;
  1299.     }
  1301.     return PSI_SUCCESS;
  1302. }
  1304. Uint4
  1305. _PSISequenceLengthWithoutX(const Uint1* seq, Uint4 length)
  1306. {
  1307.     const Uint1 X = AMINOACID_TO_NCBISTDAA['X'];
  1308.     Uint4 retval = 0;       /* the return value */
  1309.     Uint4 i = 0;            /* loop index */
  1311.     ASSERT(seq);
  1313.     for(i = 0; i < length; i++) {
  1314.         if (seq[i] != X) {
  1315.             retval++;
  1316.         }
  1317.     }
  1319.     return retval;
  1320. }
  1322. Blast_ScoreFreq*
  1323. _PSIComputeScoreProbabilities(const int** pssm,                     /* [in] */
  1324.                               const Uint1* query,                   /* [in] */
  1325.                               Uint4 query_length,                   /* [in] */
  1326.                               const double* std_probs,              /* [in] */
  1327.                               const BlastScoreBlk* sbp)             /* [in] */
  1328. {
  1329.     const Uint1 X = AMINOACID_TO_NCBISTDAA['X'];
  1330.     Uint1 aa_alphabet[BLASTAA_SIZE];            /* ncbistdaa alphabet */
  1331.     Uint4 effective_length = 0;                 /* length of query w/o Xs */
  1332.     Uint4 p = 0;                                /* index on positions */
  1333.     Uint4 c = 0;                                /* index on characters */
  1334.     int s = 0;                                  /* index on scores */
  1335.     int min_score = 0;                          /* minimum score in pssm */
  1336.     int max_score = 0;                          /* maximum score in pssm */
  1337.     short rv = 0;                               /* temporary return value */
  1338.     Blast_ScoreFreq* score_freqs = NULL;        /* score frequencies */
  1340.     ASSERT(pssm);
  1341.     ASSERT(query);
  1342.     ASSERT(std_probs);
  1343.     ASSERT(sbp);
  1344.     ASSERT(sbp->alphabet_code == BLASTAA_SEQ_CODE);
  1346.     rv = Blast_GetStdAlphabet(sbp->alphabet_code, aa_alphabet, BLASTAA_SIZE);
  1347.     if (rv <= 0) {
  1348.         return NULL;
  1349.     }
  1350.     ASSERT(rv == sbp->alphabet_size);
  1352.     effective_length = _PSISequenceLengthWithoutX(query, query_length);
  1354.     /* Get the minimum and maximum scores */
  1355.     for (p = 0; p < query_length; p++) {
  1356.         for (c = 0; c < (Uint4) sbp->alphabet_size; c++) {
  1357.             const int kScore = pssm[p][aa_alphabet[c]];
  1359.             if (kScore <= BLAST_SCORE_MIN || kScore >= BLAST_SCORE_MAX) {
  1360.                 continue;
  1361.             }
  1362.             max_score = MAX(kScore, max_score);
  1363.             min_score = MIN(kScore, min_score);
  1364.         }
  1365.     }
  1367.     if ( !(score_freqs = Blast_ScoreFreqNew(min_score, max_score)))
  1368.         return NULL;
  1370.     score_freqs->obs_min = min_score;
  1371.     score_freqs->obs_max = max_score;
  1372.     for (p = 0; p < query_length; p++) {
  1373.         if (query[p] == X) {
  1374.             continue;
  1375.         }
  1377.         for (c = 0; c < (Uint4) sbp->alphabet_size; c++) {
  1378.             const Uint4 kScore = pssm[p][aa_alphabet[c]];
  1380.             if (kScore <= BLAST_SCORE_MIN || kScore >= BLAST_SCORE_MAX) {
  1381.                 continue;
  1382.             }
  1384.             /* Increment the weight for the score in position p, residue c */
  1385.             score_freqs->sprob[kScore] += 
  1386.                 (std_probs[aa_alphabet[c]]/effective_length);
  1387.         }
  1388.     }
  1390.     for (s = min_score; s < max_score; s++) {
  1391.         score_freqs->score_avg += (s * score_freqs->sprob[s]);
  1392.     }
  1394.     return score_freqs;
  1395. }
  1397. void
  1398. _PSIUpdateLambdaK(const int** pssm,              /* [in] */
  1399.                   const Uint1* query,            /* [in] */
  1400.                   Uint4 query_length,            /* [in] */
  1401.                   const double* std_probs,       /* [in] */
  1402.                   double* initial_lambda_guess,  /* [in] */
  1403.                   BlastScoreBlk* sbp)            /* [in|out] */
  1404. {
  1405.     Blast_ScoreFreq* score_freqs = 
  1406.         _PSIComputeScoreProbabilities(pssm, query, query_length, 
  1407.                                       std_probs, sbp);
  1409.     if (initial_lambda_guess) {
  1410.         sbp->kbp_psi[0]->Lambda = Blast_KarlinLambdaNR(score_freqs, 
  1411.                                                        *initial_lambda_guess);
  1413.     } else {
  1414.         /* Calculate lambda and K */
  1415.         Blast_KarlinBlkCalc(sbp->kbp_psi[0], score_freqs);
  1417.         /* Shouldn't this be in a function? */
  1418.         sbp->kbp_gap_psi[0]->K = 
  1419.             sbp->kbp_psi[0]->K * sbp->kbp_gap_std[0]->K / sbp->kbp_ideal->K;
  1420.         sbp->kbp_gap_psi[0]->logK = log(sbp->kbp_gap_psi[0]->K);
  1421.     }
  1423.     score_freqs = Blast_ScoreFreqDestruct(score_freqs);
  1424. }
  1427. /****************************************************************************/
  1428. /* Function definitions for auxiliary functions for the stages above */
  1429. int
  1430. _PSIPurgeAlignedRegion(PsiAlignmentData* alignment,
  1431.                        unsigned int seq_index,
  1432.                        unsigned int start,
  1433.                        unsigned int stop)
  1434. {
  1435.     PsiDesc* sequence_position = NULL;
  1436.     unsigned int i = 0;
  1438.     if (!alignment)
  1439.         return PSIERR_BADPARAM;
  1441.     /* Cannot remove the query sequence from multiple alignment data or
  1442.        bad index */
  1443.     if (seq_index == kQueryIndex || 
  1444.         seq_index > alignment->dimensions->num_seqs + 1 ||
  1445.         stop > alignment->dimensions->query_sz)
  1446.         return PSIERR_BADPARAM;
  1449.     sequence_position = alignment->desc_matrix[seq_index];
  1450.     fprintf(stderr, "NEW purge %d (%d-%d)n", seq_index, start, stop);
  1451.     for (i = start; i < stop; i++) {
  1452.         sequence_position[i].letter = (unsigned char) -1;
  1453.         sequence_position[i].used = FALSE;
  1454.         sequence_position[i].e_value = kDefaultEvalueForPosition;
  1455.         sequence_position[i].extents.left = (unsigned int) -1;
  1456.         sequence_position[i].extents.right = alignment->dimensions->query_sz;
  1457.     }
  1459.     _PSIDiscardIfUnused(alignment, seq_index);
  1461.     return PSI_SUCCESS;
  1462. }
  1464. /* Check if we still need this sequence */
  1465. void
  1466. _PSIDiscardIfUnused(PsiAlignmentData* alignment, unsigned int seq_index)
  1467. {
  1468.     Boolean contains_aligned_regions = FALSE;
  1469.     unsigned int i = 0;
  1471.     for (i = 0; i < alignment->dimensions->query_sz; i++) {
  1472.         if (alignment->desc_matrix[seq_index][i].used) {
  1473.             contains_aligned_regions = TRUE;
  1474.             break;
  1475.         }
  1476.     }
  1478.     if ( !contains_aligned_regions ) {
  1479.         alignment->use_sequences[seq_index] = FALSE;
  1480.         fprintf(stderr, "NEW Removing sequence %dn", seq_index);
  1481.     }
  1482. }
  1484. /****************************************************************************/
  1485. double*
  1486. _PSIGetStandardProbabilities(const BlastScoreBlk* sbp)
  1487. {
  1488.     Blast_ResFreq* standard_probabilities = NULL;
  1489.     Uint4 i = 0;
  1490.     double* retval = NULL;
  1492.     if ( !(retval = (double*) malloc(sbp->alphabet_size * sizeof(double))))
  1493.         return NULL;
  1495.     standard_probabilities = Blast_ResFreqNew(sbp);
  1496.     Blast_ResFreqStdComp(sbp, standard_probabilities);
  1498.     for (i = 0; i < (Uint4) sbp->alphabet_size; i++) {
  1499.         retval[i] = standard_probabilities->prob[i];
  1500.     }
  1502.     Blast_ResFreqDestruct(standard_probabilities);
  1503.     return retval;
  1504. }
  1506. PsiDiagnostics*
  1507. _PSISaveDiagnostics(const PsiAlignmentData* alignment,
  1508.                     const PsiAlignedBlock* aligned_block,
  1509.                     const PsiSequenceWeights* seq_weights)
  1510. {
  1511.     /* _PSICalculateInformationContent(seq_weights); */
  1512.     abort();
  1513.     return NULL;
  1514. }
  1517. /*
  1518.  * ===========================================================================
  1519.  * $Log: blast_psi_priv.c,v $
  1520.  * Revision 1000.1  2004/06/01 18:07:34  gouriano
  1521.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  1522.  *
  1523.  * Revision 1.6  2004/05/28 17:35:03  camacho
  1524.  * Fix msvc6 warnings
  1525.  *
  1526.  * Revision 1.5  2004/05/28 16:00:09  camacho
  1527.  * + first port of PSSM generation engine
  1528.  *
  1529.  * Revision 1.4  2004/05/19 14:52:02  camacho
  1530.  * 1. Added doxygen tags to enable doxygen processing of algo/blast/core
  1531.  * 2. Standardized copyright, CVS $Id string, $Log and rcsid formatting and i
  1532.  *    location
  1533.  * 3. Added use of @todo doxygen keyword
  1534.  *
  1535.  * Revision 1.3  2004/05/06 14:01:40  camacho
  1536.  * + _PSICopyMatrix
  1537.  *
  1538.  * Revision 1.2  2004/04/07 22:08:37  kans
  1539.  * needed to include blast_def.h for sfree prototype
  1540.  *
  1541.  * Revision 1.1  2004/04/07 19:11:17  camacho
  1542.  * Initial revision
  1543.  *
  1544.  * ===========================================================================
  1545.  */