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

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

生物技术

开发平台：

C/C++

phi_lookup.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: phi_lookup.c,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:08:29  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.14
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: phi_lookup.c,v 1000.2 2004/06/01 18:08:29 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 offical 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: Ilya Dondoshansky
  36.  *
  37.  */
  39. /** @file phi_lookup.c
  40.  * Functions for accessing the lookup table for PHI-BLAST
  41.  * @todo FIXME needs doxygen comments and lines shorter than 80 characters
  42.  */
  44. static char const rcsid[] = 
  45.     "$Id: phi_lookup.c,v 1000.2 2004/06/01 18:08:29 gouriano Exp $";
  47. #include <algo/blast/core/blast_def.h>
  48. #include <algo/blast/core/blast_util.h>
  49. #include <algo/blast/core/pattern.h>
  50. #include <algo/blast/core/phi_lookup.h>
  51. #include <algo/blast/core/blast_message.h>
  54. #define seedepsilon 0.00001
  55. #define allone  ((1 << ALPHABET_SIZE) - 1)
  57. typedef struct seedSearchItems {
  58.     double  charMultiple[ALPHABET_SIZE];
  59.     double  paramC; /*used in e-value computation*/
  60.     double  paramLambda; /*used in e-value computation*/
  61.     double  paramK; /*used in the bit score computation*/
  62.     Int4         cutoffScore; /*lower bound for what is a hit*/
  63.     double  standardProb[ALPHABET_SIZE]; /*probability of each letter*/
  64.     char         order[ASCII_SIZE];
  65.     char         pchars[ALPHABET_SIZE+1];
  66.     char         name_space[BUF_SIZE];  /*name of a pattern*/
  67.     char         pat_space[PATTERN_SPACE_SIZE];  /*string description
  68.                                                    of pattern*/
  69. } seedSearchItems;
  71. /*Initialize the order of letters in the alphabet, the score matrix,
  72. and the row sums of the score matrix. matrixToFill is the
  73. score matrix, program_flag says which variant of the program is
  74. used; is_dna tells whether the strings are DNA or protein*/
  75. static void init_order(Int4 **matrix, Int4 program_flag, Boolean is_dna, seedSearchItems *seedSearch)
  76. {
  77.     Uint1 i, j; /*loop indices for matrix*/ 
  78.     Int4 *matrixRow; /*row of matrixToFill*/ 
  79.     double rowSum; /*sum of scaled substitution probabilities on matrixRow*/
  80.     
  81.     if (is_dna) {
  82.       seedSearch->order['A'] = 0; 
  83.       seedSearch->order['C'] = 1;
  84.       seedSearch->order['G'] = 2; 
  85.       seedSearch->order['T'] = 3;
  86.     } 
  87.     else {
  88.       for (i = 0; i < ALPHABET_SIZE; i++) 
  89. seedSearch->order[(Uint1)seedSearch->pchars[i]] = i;
  90.     }
  91.     if (program_flag == SEED_FLAG) {
  92.       for (i = 0; i < ALPHABET_SIZE; i++) 
  93.         seedSearch->charMultiple[i] = 0;
  94.       for (i = 0; i < ALPHABET_SIZE; i++) {
  95. if (seedSearch->standardProb[i] > seedepsilon) {
  96.   matrixRow = matrix[i];
  97.   rowSum= 0;
  98.   for (j = 0; j < ALPHABET_SIZE; j++) {
  99.     if (seedSearch->standardProb[j] > seedepsilon) 
  100.       rowSum += seedSearch->standardProb[j]*exp(-(seedSearch->paramLambda)*matrixRow[j]);
  101.   }
  102.   seedSearch->charMultiple[i] = rowSum;
  103. }
  104.       }
  105.     }
  106. }
  108. /*Initialize occurrence probabilities for each amino acid*/
  109. static void initProbs(seedSearchItems * seedSearch)
  110. {
  111.    double totalCount;  /*for Robinson frequencies*/
  112.    seedSearch->pchars[0] = '-';
  113.    seedSearch->pchars[1] = 'A';
  114.    seedSearch->pchars[2] = 'B';
  115.    seedSearch->pchars[3] = 'C';
  116.    seedSearch->pchars[4] = 'D';
  117.    seedSearch->pchars[5] = 'E';
  118.    seedSearch->pchars[6] = 'F';
  119.    seedSearch->pchars[7] = 'G';
  120.    seedSearch->pchars[8] = 'H';
  121.    seedSearch->pchars[9] = 'I';
  122.    seedSearch->pchars[10] = 'K';
  123.    seedSearch->pchars[11] = 'L';
  124.    seedSearch->pchars[12] = 'M';
  125.    seedSearch->pchars[13] = 'N';
  126.    seedSearch->pchars[14] = 'P';
  127.    seedSearch->pchars[15] = 'Q';
  128.    seedSearch->pchars[16] = 'R';
  129.    seedSearch->pchars[17] = 'S';
  130.    seedSearch->pchars[18] = 'T';
  131.    seedSearch->pchars[19] = 'V';
  132.    seedSearch->pchars[20] = 'W';
  133.    seedSearch->pchars[21] = 'X';
  134.    seedSearch->pchars[22] = 'Y';
  135.    seedSearch->pchars[23] = 'Z';
  136.    seedSearch->pchars[24] = 'U';
  137.    seedSearch->pchars[25] = '*';
  138.    totalCount = 78.0 + 19.0 + 54.0 + 63.0 + 39.0 +
  139.     74.0 + 22.0 + 52.0 + 57.0 + 90.0 + 22.0 + 45.0 + 52.0 +
  140.      43.0 + 51.0 + 71.0 + 59.0 + 64.0 + 13.0 + 32.0;
  141.    seedSearch->standardProb[0] = 0.0;
  142.    seedSearch->standardProb[1] = 78.0/totalCount; /*A*/
  143.    seedSearch->standardProb[2] = 0.0;
  144.    seedSearch->standardProb[3] = 19.0/totalCount; /*C*/
  145.    seedSearch->standardProb[4] = 54.0/totalCount; /*D*/
  146.    seedSearch->standardProb[5] = 63.0/totalCount; /*E*/
  147.    seedSearch->standardProb[6] = 39.0/totalCount; /*F*/
  148.    seedSearch->standardProb[7] = 74.0/totalCount; /*G*/
  149.    seedSearch->standardProb[8] = 22.0/totalCount; /*H*/
  150.    seedSearch->standardProb[9] = 52.0/totalCount; /*I*/
  151.    seedSearch->standardProb[10] = 57.0/totalCount; /*K*/
  152.    seedSearch->standardProb[11] = 90.0/totalCount; /*L*/
  153.    seedSearch->standardProb[12] = 22.0/totalCount; /*M*/
  154.    seedSearch->standardProb[13] = 45.0/totalCount; /*N*/
  155.    seedSearch->standardProb[14] = 52.0/totalCount; /*P*/
  156.    seedSearch->standardProb[15] = 43.0/totalCount; /*Q*/
  157.    seedSearch->standardProb[16] = 51.0/totalCount; /*R*/
  158.    seedSearch->standardProb[17] = 71.0/totalCount; /*S*/
  159.    seedSearch->standardProb[18] = 59.0/totalCount; /*T*/
  160.    seedSearch->standardProb[19] = 64.0/totalCount; /*V*/
  161.    seedSearch->standardProb[20] = 13.0/totalCount; /*W*/
  162.    seedSearch->standardProb[21] = 0.0;   /*X*/
  163.    seedSearch->standardProb[22] = 32.0/totalCount;   /*Y*/
  164.    seedSearch->standardProb[23] = 0.0;   /*Z*/
  165.    seedSearch->standardProb[24] = 0.0;   /*U*/
  166. }
  168. /*set uo matches for words that encode 4 DNA characters; figure out
  169.   for each of 256 possible DNA 4-mers, where a prefix matches the pattern
  170.  and where a suffix matches the pattern; store in prefixPos and
  171.  suffixPos; mask has 1 bits for whatever lengths of string
  172. the pattern can match, mask2 has 4 1 bits corresponding to
  173. the last 4 positions of a match; they are used to
  174. do the prefixPos and suffixPos claculations with bit arithmetic*/
  175. static void setting_tt(Int4* S, Int4 mask, Int4 mask2, Uint4* prefixPos, 
  176.        Uint4* suffixPos)
  177. {
  178.   Int4 i; /*index over possible DNA encoded words, 4 bases per word*/
  179.   Int4 tmp; /*holds different mask combinations*/
  180.   Int4 maskLeftPlusOne; /*mask shifted left 1 plus 1; guarantees 1
  181.                            1 character match effectively */
  182.   Uint1 a1, a2, a3, a4;  /*four bases packed into an integer*/
  184.   maskLeftPlusOne = (mask << 1)+1;
  185.   for (i = 0; i < ASCII_SIZE; i++) {
  186.     /*find out the 4 bases packed in integer i*/
  187.     a1 = NCBI2NA_UNPACK_BASE(i, 3);
  188.     a2 = NCBI2NA_UNPACK_BASE(i, 2);
  189.     a3 = NCBI2NA_UNPACK_BASE(i, 1);
  190.     a4 = NCBI2NA_UNPACK_BASE(i, 0);
  191.     /*what positions match a prefix of a4 followed by a3*/
  192.     tmp = ((S[a4]>>1) | mask) & S[a3];
  193.     /*what positions match a prefix of a4 followed by a3 followed by a2*/
  194.     tmp = ((tmp >>1) | mask) & S[a2];
  195.     /*what positions match a prefix of a4, a3, a2,a1*/
  196.     prefixPos[i] = mask2 & ((tmp >>1) | mask) & S[a1];
  197.     
  198.     /*what positions match a suffix of a2, a1*/
  199.     tmp = ((S[a1]<<1) | maskLeftPlusOne) & S[a2];
  200.     /* what positions match a suffix of a3, a2, a1*/
  201.     tmp = ((tmp <<1) | maskLeftPlusOne) & S[a3];
  202.     /*what positions match a suffix of a4, a3, a2, a1*/
  203.     suffixPos[i] = ((((tmp <<1) | maskLeftPlusOne) & S[a4]) << 1) | maskLeftPlusOne;
  204.   }
  205. }
  207. /*initialize mask and other arrays for DNA patterns*/
  208. static void init_pattern_DNA(patternSearchItems *patternSearch)
  209. {
  210.   Int4 mask1; /*mask for one word in a set position*/
  211.   Int4 compositeMask; /*superimposed mask1 in 4 adjacent positions*/
  212.   Int4 wordIndex; /*index over words in pattern*/
  214.   if (patternSearch->flagPatternLength != ONE_WORD_PATTERN) {
  215.     for (wordIndex = 0; wordIndex < patternSearch->numWords; wordIndex++) {
  216.       mask1 = patternSearch->match_maskL[wordIndex];
  217.       compositeMask = mask1 + (mask1>>1)+(mask1>>2)+(mask1>>3);
  218.       setting_tt(patternSearch->SLL[wordIndex], 
  219.       patternSearch->match_maskL[wordIndex], 
  220.  compositeMask, patternSearch->DNAprefixSLL[wordIndex], patternSearch->DNAsuffixSLL[wordIndex]);
  221.     }
  222.   } 
  223.   else {
  224.     compositeMask = patternSearch->match_mask + 
  225.       (patternSearch->match_mask>>1) + 
  226.       (patternSearch->match_mask>>2) + (patternSearch->match_mask>>3); 
  227.     patternSearch->DNAwhichPrefixPosPtr = patternSearch->DNAwhichPrefixPositions; 
  228.     patternSearch->DNAwhichSuffixPosPtr = patternSearch->DNAwhichSuffixPositions;
  229.     setting_tt(patternSearch->whichPositionsByCharacter, 
  230.     patternSearch->match_mask, compositeMask, 
  231.     patternSearch->DNAwhichPrefixPositions, patternSearch->DNAwhichSuffixPositions);
  232.   }
  233. }
  235. /*length is the length of inputPattern, maxLength is a limit on
  236.    how long inputPattern can get;
  237.    return the final length of the pattern or -1 if too long*/
  238. static Int4 expanding(Int4 *inputPatternMasked, Uint1 *inputPattern, 
  239.       Int4 length, Int4 maxLength)
  240. {
  241.     Int4 i, j; /*pattern indices*/
  242.     Int4 numPos; /*number of positions index*/
  243.     Int4  k, t; /*loop indices*/
  244.     Int4 recReturnValue1, recReturnValue2; /*values returned from
  245.                                              recursive calls*/
  246.     Int4 thisPlaceMasked; /*value of one place in inputPatternMasked*/
  247.     Int4 tempPatternMask[MaxP]; /*used as a local representation of
  248.                                part of inputPatternMasked*/
  249.     Uint1 tempPattern[MaxP]; /*used as a local representation of part of
  250.                                inputPattern*/
  252.     for (i = 0; i < length; i++) {
  253.       thisPlaceMasked = -inputPatternMasked[i];
  254.       if (thisPlaceMasked > 0) {  /*represented variable wildcard*/
  255. inputPatternMasked[i] = allone;
  256. for (j = 0; j < length; j++) {
  257.   /*use this to keep track of pattern*/
  258.   tempPatternMask[j] = inputPatternMasked[j]; 
  259.   tempPattern[j] = inputPattern[j];
  260. }
  261. recReturnValue2 = recReturnValue1 = 
  262.   expanding(inputPatternMasked, inputPattern, length, maxLength);
  263. if (recReturnValue1 == -1)
  264.   return -1;
  265. for (numPos = 0; numPos <= thisPlaceMasked; numPos++) {
  266.   if (numPos == 1)
  267.     continue;
  268.   for (k = 0; k < length; k++) {
  269.     if (k == i) {
  270.       for (t = 0; t < numPos; t++) {
  271. inputPatternMasked[recReturnValue1++] = allone;
  272.                 if (recReturnValue1 >= maxLength)
  273.                   return(-1);
  274.       }
  275.     }
  276.     else {
  277.       inputPatternMasked[recReturnValue1] = tempPatternMask[k];
  278.       inputPattern[recReturnValue1++] = tempPattern[k];
  279.               if (recReturnValue1 >= maxLength)
  280.                   return(-1);
  281.     }
  282.     if (recReturnValue1 >= maxLength) 
  283.       return (-1);
  284.   }
  285.   recReturnValue1 = 
  286.     expanding(&inputPatternMasked[recReturnValue2], 
  287.       &inputPattern[recReturnValue2], 
  288.       length + numPos - 1, 
  289.       maxLength - recReturnValue2);
  290.   if (recReturnValue1 == -1) 
  291.     return -1;
  292.   recReturnValue2 += recReturnValue1; 
  293.   recReturnValue1 = recReturnValue2;
  294. }
  295. return recReturnValue1;
  296.       }
  297.     }
  298.     return length;
  299. }
  301. /*Pack the next length bytes of inputPattern into a bit vector
  302.   where the bit is 1 if and only if the byte is non-0 
  303.   Returns packed bit vector*/
  304. static Int4 packing(Uint1 *inputPattern, Int4 length)
  305. {
  306.     Int4 i; /*loop index*/
  307.     Int4 returnValue = 0; /*value to return*/
  308.     for (i = 0; i < length; i++) {
  309.       if (inputPattern[i])
  310. returnValue += (1 << i);
  311.     }
  312.     return returnValue;
  313. }
  315. /*Pack the bit representation of the inputPattern into
  316.    the array patternSearch->match_maskL
  317.    numPlaces is the number of positions in
  318.    inputPattern
  319.    Also packs patternSearch->bitPatternByLetter  */
  320. static void longpacking(Int4 numPlaces, Uint1 *inputPattern, patternSearchItems *patternSearch)
  321. {
  322.     Int4 charIndex; /*index over characters in alphabet*/
  323.     Int4 bitPattern; /*bit pattern for one word to pack*/
  324.     Int4 i;  /*loop index over places*/
  325.     Int4 wordIndex; /*loop counter over words to pack into*/
  326.     
  327.     patternSearch->numWords = (numPlaces-1) / BITS_PACKED_PER_WORD +1;
  329.     for (wordIndex = 0; wordIndex < patternSearch->numWords; wordIndex++) {
  330.       bitPattern = 0;
  331.       for (i = 0; i < BITS_PACKED_PER_WORD; i++) {
  332. if (inputPattern[wordIndex*BITS_PACKED_PER_WORD+i]) 
  333.   bitPattern += (1 << i);
  334.       }
  335.       patternSearch->match_maskL[wordIndex] = bitPattern;
  336.     }
  337.     for (charIndex = 0; charIndex < ALPHABET_SIZE; charIndex++) {
  338.       for (wordIndex = 0; wordIndex < patternSearch->numWords; wordIndex++) {
  339. bitPattern = 0;
  340. for (i = 0; i < BITS_PACKED_PER_WORD; i++) {
  341.   if ((1<< charIndex) & patternSearch->inputPatternMasked[wordIndex*BITS_PACKED_PER_WORD + i]) 
  342.     bitPattern = bitPattern | (1 << i);
  343. }
  344. patternSearch->bitPatternByLetter[charIndex][wordIndex] = 
  345.   bitPattern;
  346. }
  347.     }
  348. }
  350. /*Return the number of 1 bits in the base 2 representation of a*/
  351. static Int4 num_of_one(Int4 a)
  352. {
  353.   Int4 returnValue;
  354.   returnValue = 0;
  355.   while (a > 0) {
  356.     if (a % 2 == 1) 
  357.       returnValue++;
  358.     a = (a >> 1);
  359.   }
  360.   return returnValue;
  361. }
  363. /*Sets up field in patternSearch when pattern is very long*/
  364. static void longpacking2(Int4 *inputPatternMasked, Int4 numPlacesInPattern, patternSearchItems *patternSearch)
  365. {
  366.     Int4 placeIndex; /*index over places in pattern rep.*/
  367.     Int4 wordIndex; /*index over words*/
  368.     Int4 placeInWord, placeInWord2;  /*index for places in a single word*/
  369.     Int4 charIndex; /*index over characters in alphabet*/
  370.     Int4 oneWordMask; /*mask of matching characters for one word in
  371.                         pattern representation*/
  372.     double patternWordProbability;
  373.     double  most_specific; /*lowest probability of a word in the pattern*/
  374.     Int4 *oneWordSLL; /*holds patternSearch->SLL for one word*/
  376.     most_specific = 1.0; 
  377.     patternSearch->whichMostSpecific = 0; 
  378.     patternWordProbability = 1.0;
  379.     for (placeIndex = 0, wordIndex = 0, placeInWord=0; 
  380.  placeIndex <= numPlacesInPattern;   placeIndex++, placeInWord++) {
  381.       if (placeIndex==numPlacesInPattern || inputPatternMasked[placeIndex] < 0 
  382.   || placeInWord == BITS_PACKED_PER_WORD ) {
  383. patternSearch->match_maskL[wordIndex] = 1 << (placeInWord-1);
  384. oneWordSLL = patternSearch->SLL[wordIndex];
  385. for (charIndex = 0; charIndex < ALPHABET_SIZE; charIndex++) {
  386.   oneWordMask = 0;
  387.   for (placeInWord2 = 0; placeInWord2 < placeInWord; placeInWord2++) {
  388.     if ((1<< charIndex) & 
  389. inputPatternMasked[placeIndex-placeInWord+placeInWord2]) 
  390.       oneWordMask |= (1 << placeInWord2);
  391.   }
  392.   oneWordSLL[charIndex] = oneWordMask;
  393. }
  394. patternSearch->numPlacesInWord[wordIndex] = placeInWord;
  395. if (patternWordProbability < most_specific) {
  396.   most_specific = patternWordProbability;
  397.   patternSearch->whichMostSpecific = wordIndex;
  398. }
  399. if (placeIndex == numPlacesInPattern) 
  400.   patternSearch->spacing[wordIndex++] = 0; 
  401. else 
  402.   if (inputPatternMasked[placeIndex] < 0) { 
  403.     patternSearch->spacing[wordIndex++] = -inputPatternMasked[placeIndex];
  404.   }
  405.   else { 
  406.     placeIndex--; 
  407.     patternSearch->spacing[wordIndex++] = 0;
  408.   }
  409. placeInWord = -1; 
  410. patternWordProbability = 1.0;
  411.       }
  412.       else {
  413. patternWordProbability *= (double) 
  414.   num_of_one(inputPatternMasked[placeIndex])/ (double) ALPHABET_SIZE;
  415. }
  416.     }
  417.     patternSearch->numWords = wordIndex;
  418. }
  420. /*pattern is a string describing the pattern to search for;
  421.   is_dna is a boolean describing the strings are DNA or protein*/
  422. static Int4 
  423. init_pattern(Uint1 *pattern, Boolean is_dna, BlastScoreBlk* sbp, 
  424.              patternSearchItems* *pattern_info,
  425.              Blast_Message* *error_msg)
  426. {
  427.     Uint4 i; /*index over string describing the pattern*/
  428.     Uint4 j; /*index for position in pattern*/
  429.     Int4 charIndex; /*index over characters in alphabet*/
  430.     Int4 secondIndex; /*second index into pattern*/
  431.     Int4 numIdentical; /*number of consec. positions with identical specification*/
  432.     Uint4 charSetMask;  /*index over masks for specific characters*/
  433.     Int4 currentSetMask, prevSetMask ; /*mask for current and previous character positions*/    
  434.     Int4 thisMask;    /*integer representing a bit pattern for a 
  435.                         set of characters*/
  436.     Int4 minWildcard, maxWildcard; /*used for variable number of wildcard
  437.                                      positions*/
  438.     Uint4  tj=0; /*temporary copy of j*/
  439.     Int4 tempInputPatternMasked[MaxP]; /*local copy of parts
  440.             of inputPatternMasked*/
  441.     Uint1 c;  /*character occurring in pattern*/
  442.     Uint1 localPattern[MaxP]; /*local variable to hold
  443.                                for each position whether it is
  444.                                last in pattern (1) or not (0) */
  445.     double positionProbability; /*probability of a set of characters
  446.                                     allowed in one position*/
  447.     Int4 currentWildcardProduct; /*product of wildcard lengths for
  448.                                    consecutive character positions that
  449.                                    overlap*/
  450.     seedSearchItems *seedSearch;
  451.     patternSearchItems* patternSearch;
  453.     seedSearch = (seedSearchItems*) calloc(1, sizeof(seedSearchItems));
  454.     patternSearch = *pattern_info = 
  455.        (patternSearchItems*) calloc(1, sizeof(patternSearchItems));
  457.     initProbs(seedSearch);
  458.     init_order(sbp->matrix, PAT_SEED_FLAG, FALSE, seedSearch);
  460.     patternSearch->flagPatternLength = ONE_WORD_PATTERN; 
  461.     patternSearch->patternProbability = 1.0;
  462.     patternSearch->minPatternMatchLength = 0;
  463.     patternSearch->wildcardProduct = 1;
  464.     currentWildcardProduct = 1;
  465.     prevSetMask = 0;
  466.     currentSetMask = 0;
  468.     for (i = 0 ; i < MaxP; i++) {
  469.       patternSearch->inputPatternMasked[i] = 0; 
  470.       localPattern[i] = 0;
  471.     }
  472.     for (i = 0, j = 0; i < strlen((Char *) pattern); i++) {
  473.       if ((c=pattern[i]) == '-' || c == 'n' || c == '.' || c =='>' || c ==' ' 
  474. || c == '<')  /*spacers that mean nothing*/
  475. continue;
  476.       if ( c != '[' && c != '{') { /*not the start of a set of characters*/
  477. if (c == 'x' || c== 'X') {  /*wild-card character matches anything*/
  478.           /*next line checks to see if wild card is for multiple positions*/
  479.   if (pattern[i+1] == '(') {
  480.     i++;
  481.     secondIndex = i;
  482.             /*find end of description of how many positions are wildcarded
  483.                will look like x(2) or x(2,5) */
  484.     while (pattern[secondIndex] != ',' && pattern[secondIndex] != ')')
  485.       secondIndex++;
  486.     if (pattern[secondIndex] == ')') {  /*fixed number of positions wildcarded*/
  487.       i -= 1; 
  488.               /*wildcard, so all characters are allowed*/
  489.       charSetMask=allone; 
  490.       positionProbability = 1;
  491.     }
  492.     else { /*variable number of positions wildcarded*/   
  493.       sscanf((Char*) &pattern[++i], "%d,%d", &minWildcard, &maxWildcard);
  494.       maxWildcard = maxWildcard - minWildcard;
  495.          currentWildcardProduct *= (maxWildcard + 1);
  496.          if (currentWildcardProduct > patternSearch->wildcardProduct)
  497.             patternSearch->wildcardProduct = currentWildcardProduct;
  498.          patternSearch->minPatternMatchLength += minWildcard;
  499.       while (minWildcard-- > 0) { 
  500.             /*use one position each for the minimum number of
  501.               wildcard spaces required */
  502.             patternSearch->inputPatternMasked[j++] = allone; 
  503.             if (j >= MaxP) {
  504.                Blast_MessageWrite(error_msg, BLAST_SEV_WARNING, 2, 1, 
  505.                                   "pattern too long");
  506.                return(-1);
  507.             }
  508.       }
  509.       if (maxWildcard != 0) {
  510.             /*negative masking used to indicate variability
  511.               in number of wildcard spaces; e.g., if pattern looks
  512.               like x(3,5) then variability is 2 and there will
  513.               be three wildcard positions with mask allone followed
  514.               by a single position with mask -2*/
  515.             patternSearch->inputPatternMasked[j++] = -maxWildcard;
  516.             patternSearch->patternProbability *= maxWildcard;
  517.       }
  518.          /*now skip over wildcard description with the i index*/
  519.       while (pattern[++i] != ')') ; 
  520.       continue;
  521.     }
  522.   }
  523.   else {  /*wild card is for one position only*/
  524.     charSetMask=allone; 
  525.     positionProbability =1;
  526.   }
  527. } 
  528. else {
  529.   if (c == 'U') {   /*look for special U character*/
  530.     charSetMask = allone*2+1;
  531.     positionProbability = 1; 
  532.   }
  533.   else { 
  534.         /*exactly one character matches*/
  535.         prevSetMask = currentSetMask;
  536.         currentSetMask =  charSetMask = (1 << seedSearch->order[c]);
  537.         if (!(prevSetMask & currentSetMask)) /*character sets don't overlap*/
  538.            currentWildcardProduct = 1;
  539.         positionProbability = 
  540.            seedSearch->standardProb[(Uint1)seedSearch->order[c]];
  541.   }
  542. }
  543.       }
  544.       else {
  545. if (c == '[') {  /*start of a set of characters allowed*/
  546.   charSetMask = 0;
  547.   positionProbability = 0;
  548.   /*For each character in the set add it to the mask and
  549.             add its probability to positionProbability*/
  550.   while ((c=pattern[++i]) != ']') { /*end of set*/
  551.         if ((c < 'A') || (c > 'Z') || (c == '')) {
  552.            Blast_MessageWrite(error_msg, BLAST_SEV_WARNING, 2, 1, 
  553.                               "pattern description has a non-alphabetic"
  554.                               "character inside a bracket");
  555.            
  556.            return(-1);
  557.         }
  558.         charSetMask = charSetMask | (1 << seedSearch->order[c]);
  559.         positionProbability += seedSearch->standardProb[(Uint1)seedSearch->order[c]];
  560.   }
  561.      prevSetMask = currentSetMask;
  562.      currentSetMask = charSetMask;
  563.   if (!(prevSetMask & currentSetMask)) /*character sets don't overlap*/
  564.       currentWildcardProduct = 1;
  565.   } 
  566. else {   /*start of a set of characters forbidden*/
  567.   /*For each character forbidden remove it to the mask and
  568.             subtract its probability from positionProbability*/
  569.   charSetMask = allone; 
  570.   positionProbability = 1;
  571.   while ((c=pattern[++i]) != '}') { /*end of set*/
  572.     charSetMask = charSetMask -  (charSetMask & (1 << seedSearch->order[c]));
  573.     positionProbability -= seedSearch->standardProb[(Uint1)seedSearch->order[c]];
  574.   }
  575.           prevSetMask = currentSetMask;
  576.           currentSetMask = charSetMask;
  577.   if (!(prevSetMask & currentSetMask)) /*character sets don't overlap*/
  578.       currentWildcardProduct = 1;
  579. }
  580.       }
  581.       /*handle a number of positions that are the same */
  582.       if (pattern[i+1] == '(') {  /*read opening paren*/
  583. i++;
  584. numIdentical = atoi((Char *) &pattern[++i]);  /*get number of positions*/
  585.         patternSearch->minPatternMatchLength += numIdentical;
  586. while (pattern[++i] != ')') ;  /*skip over piece in pattern*/
  587. while ((numIdentical--) > 0) {
  588.   /*set up mask for these positions*/
  589.   patternSearch->inputPatternMasked[j++] = charSetMask;
  590.   patternSearch->patternProbability *= positionProbability; 
  591. }
  592.       } 
  593.       else {   /*specification is for one posiion only*/
  594.          patternSearch->inputPatternMasked[j++] = charSetMask;
  595.          patternSearch->minPatternMatchLength++;
  596.          patternSearch->patternProbability *= positionProbability;
  597.       }
  598.       if (j >= MaxP) {
  599.          Blast_MessageWrite(error_msg, BLAST_SEV_WARNING, 2, 1, 
  600.                             "pattern too long");
  601.       }
  602.     }
  603.     localPattern[j-1] = 1;
  604.     if (patternSearch->patternProbability > 1.0)
  605.       patternSearch->patternProbability = 1.0;
  607.     for (i = 0; i < j; i++) {
  608.       tempInputPatternMasked[i] = patternSearch->inputPatternMasked[i]; 
  609.       tj = j;
  610.     }
  611.     j = expanding(patternSearch->inputPatternMasked, localPattern, j, MaxP);
  612.     if ((j== -1) || ((j > BITS_PACKED_PER_WORD) && is_dna)) {
  613.       patternSearch->flagPatternLength = PATTERN_TOO_LONG;
  614.       longpacking2(tempInputPatternMasked, tj, patternSearch);
  615.       for (i = 0; i < tj; i++) 
  616.          patternSearch->inputPatternMasked[i] = tempInputPatternMasked[i];
  617.       patternSearch->highestPlace = tj;
  618.       if (is_dna) 
  619.          init_pattern_DNA(patternSearch);
  620.       return 1;
  621.     }
  622.     if (j > BITS_PACKED_PER_WORD) {
  623.       patternSearch->flagPatternLength = MULTI_WORD_PATTERN;
  624.       longpacking(j, localPattern, patternSearch);
  625.       return j;
  626.     } 
  627.     /*make a bit mask out of local pattern of length j*/
  628.     patternSearch->match_mask = packing(localPattern, j);
  629.     /*store for each character a bit mask of which positions
  630.       that character can occur in*/
  631.     for (charIndex = 0; charIndex < ALPHABET_SIZE; charIndex++) {
  632.       thisMask = 0;
  633.       for (charSetMask = 0; charSetMask < j; charSetMask++) {
  634. if ((1<< charIndex) & patternSearch->inputPatternMasked[charSetMask]) 
  635.   thisMask |= (1 << charSetMask);
  636.       }
  637.       patternSearch->whichPositionsByCharacter[charIndex] = thisMask;
  638.     }
  639.     patternSearch->whichPositionPtr = patternSearch->whichPositionsByCharacter;
  640.     if (is_dna) 
  641.       init_pattern_DNA(patternSearch);
  642.     return j; /*return number of places for pattern representation*/
  643. }
  645. Int2 PHILookupTableNew(const LookupTableOptions* opt, PHILookupTable* * lut,
  646.                        Boolean is_dna, BlastScoreBlk* sbp)
  647. {
  648.    PHILookupTable* lookup = *lut = 
  649.       (PHILookupTable*) malloc(sizeof(PHILookupTable));
  650.    Blast_Message* error_msg = NULL;
  652.    if (!lookup)
  653.       return -1;
  655.    lookup->is_dna = is_dna;
  656.    init_pattern((Uint1*)opt->phi_pattern, is_dna, sbp, &lookup->pattern_info,
  657.                 &error_msg); 
  658.    lookup->num_matches = 0;
  659.    lookup->allocated_size = MIN_PHI_LOOKUP_SIZE;
  660.    if ((lookup->start_offsets = 
  661.        (Int4*) malloc(MIN_PHI_LOOKUP_SIZE*sizeof(Int4))) == NULL)
  662.       return -1;
  663.    if ((lookup->lengths = (Int4*) malloc(MIN_PHI_LOOKUP_SIZE*sizeof(Int4)))
  664.         == NULL)
  665.       return -1;
  667.    return 0;
  668. }
  670. PHILookupTable* PHILookupTableDestruct(PHILookupTable* lut)
  671. {
  672.    sfree(lut->pattern_info);
  673.    sfree(lut->start_offsets);
  674.    sfree(lut->lengths);
  675.    sfree(lut);
  676.    return NULL;
  677. }
  679. static Int2 PHIBlastAddPatternHit(PHILookupTable* lookup, Int4 offset, 
  680.                                   Int4 length)
  681. {
  682.    if (lookup->num_matches >= lookup->allocated_size) {
  683.       lookup->start_offsets = (Int4*) realloc(lookup->start_offsets, 
  684.                                               2*lookup->allocated_size);
  685.       lookup->lengths = (Int4*) realloc(lookup->lengths, 
  686.                                               2*lookup->allocated_size);
  687.       if (!lookup->start_offsets || !lookup->lengths)
  688.          return -1;
  689.       lookup->allocated_size *= 2;
  690.    }      
  691.       
  692.    lookup->start_offsets[lookup->num_matches] = offset;
  693.    lookup->lengths[lookup->num_matches] = length;
  694.    ++lookup->num_matches;
  695.    return 0;
  696. }
  698. Int4 PHIBlastIndexQuery(PHILookupTable* lookup, 
  699.         BLAST_SequenceBlk* query, ListNode* location, Boolean is_dna)
  700. {
  701.    ListNode* loc;
  702.    Int4 from, to;
  703.    Int4 loc_length;
  704.    Uint1* sequence;
  705.    patternSearchItems* pattern_info = lookup->pattern_info;
  706.    Int4* hitArray;
  707.    Int4 i, twiceNumHits;
  708.    
  709.    hitArray = (Int4 *) calloc(2*query->length, sizeof(Int4));
  711.    for(loc=location; loc; loc=loc->next) {
  712.       from = ((SSeqRange*) loc->ptr)->left;
  713.       to = ((SSeqRange*) loc->ptr)->right;
  714.       loc_length = to - from + 1;
  715.       sequence = query->sequence + from;
  716.       
  717.       twiceNumHits = FindPatternHits(hitArray, sequence, loc_length, is_dna,
  718.                                      pattern_info);
  719.       
  720.       for (i = 0; i < twiceNumHits; i += 2) {
  721.          PHIBlastAddPatternHit(lookup, hitArray[i+1]+from, 
  722.                                hitArray[i]-hitArray[i+1]+1);
  723.       }
  724.    }
  726.    return lookup->num_matches;
  727. }
  729. static Boolean 
  730. PHIBlastMatchPatterns(Uint1* subject, Uint1* query, Int4 length)
  731. {
  732.    Int4 index;
  734.    for (index = 0; index < length; ++index) {
  735.       if (subject[index] != query[index])
  736.          break;
  737.    }
  738.    return (index == length);
  739. }
  741. Int4 PHIBlastScanSubject(const LookupTableWrap* lookup_wrap,
  742.         const BLAST_SequenceBlk *query_blk, 
  743.         const BLAST_SequenceBlk *subject_blk, 
  744.         Int4* offset_ptr, Uint4 * query_offsets, Uint4 * subject_offsets, 
  745.         Int4 array_size)
  746. {
  747.    Uint1* subject, *query;
  748.    PHILookupTable* lookup = (PHILookupTable*) lookup_wrap->lut;
  749.    Int4 index, count = 0, twiceNumHits, i;
  750.    Int4 *start_offsets = lookup->start_offsets;
  751.    Int4 *pat_lengths = lookup->lengths;
  752.    Int4 offset, length;
  753.    Int4 hitArray[MAX_HIT];
  755.    query = query_blk->sequence;
  756.    subject = subject_blk->sequence;
  757.    /* It must be guaranteed that all pattern matches for a given 
  758.       subject sequence are processed in one call to this function.
  759.    */
  760.    *offset_ptr = subject_blk->length;
  762.    twiceNumHits = FindPatternHits(hitArray, subject, subject_blk->length, 
  763.                                   lookup->is_dna, lookup->pattern_info);
  766.    for (i = 0; i < twiceNumHits; i += 2) {
  767. #if 0
  768.       if (count > array_size - lookup->num_matches)
  769.             break;
  770. #endif
  771.       length = hitArray[i] - hitArray[i+1] + 1;
  772.       offset = hitArray[i+1];
  774.       for (index = 0; index < lookup->num_matches; ++index) {
  775.          /* Match pattern lengths in subject in query first; then 
  776.             check for identical match of pattern */
  777.          if (length == pat_lengths[index] &&
  778.              PHIBlastMatchPatterns(subject+offset, query+start_offsets[index], 
  779.                                    length))
  780.          {
  781.             /* Pattern has matched completely. Save index into the array
  782.                of pattern start offsets in query (so pattern length will
  783.                be accessible in the word finder later), and the subject
  784.                offset. */
  785.             query_offsets[count] = index;
  786.             subject_offsets[count] = offset;
  787.             ++count;
  788.          }
  789.       }
  790.    }
  791.    return count;
  792. }