PssmQuery.h

   1 /**
   2  * PSSM queries
   3  *
   4  * Input is a Position Frequency Matrix (PFM, see e.g. http://jaspar.genereg.net/)
   5  * and a threshold:
   6  *
   7  *   contains('PSSM 0.1 3 4 19 0 16 0 20 0 1 0 0 0 0')
   8  *
   9  * Where
  10  *     'PSSM' is the query keyword.
  11  *      '0.1' is the minimum threshold for log-odds score.
  12  *        '3' is the number of columsn in the given PFM.
  13  *
  14  * Rest of the numbers contain the PFM in row-wise order,
  15  * e.g. in the above example we have:
  16  *    A [ 4 19  0 ]
  17  *    C [16  0 20 ]
  18  *    G [ 0  1  0 ]
  19  *    T [ 0  0  0 ]
  20  */
  21
  22 #ifndef _PssmQuery_H_
  23 #define _PssmQuery_H_
  24
  25 #include "Query.h"
  26
  27 #include <string>
  28 #include <cstdlib>
  29 #include <cstring>
  30 #include <cmath>
  31
  32 namespace SXSI
  33 {
  34
  35 class PssmQuery : public Query
  36 {
  37 public:
  38     PssmQuery(TextCollection const *tc_, double thr)
  39         : Query(tc_), threshold(thr)
  40     { }
  41
  42     virtual ~PssmQuery()
  43     { }
  44
  45 protected:
  46     void firstStep(void)
  47     {
  48         // Parse input into PFM
  49         uchar const *tmp = pat;
  50         int l = parseNextInt(tmp);
  51         if (l <= 0)
  52         {
  53             std::cerr << "PssmQuery::firstStep(): invalid length given: length must be greater than 0." << std::endl;
  54             std::exit(1);
  55         }
  56
  57         P = std::vector<double>();
  58         P.reserve(4 * l);
  59         while (*tmp != 0)
  60             P.push_back(parseNextInt(tmp));
  61
  62         if ((int)P.size() != 4 * l)
  63         {
  64             std::cerr << "PssmQuery::firstStep(): invalid length or matrix given: length was " << l << " * 4 == " << l * 4 << "  but the matrix had " << P.size() << " values." << std::endl;
  65             std::exit(1);
  66         }
  67
  68         // And into log-odds:
  69         P = convertPfmToLogOdds(P, l);
  70         calculateMaxSums(P, l);
  71         pat = 0;
  72         patlen = l;
  73         PssmMatch(patlen-1, 0);
  74     }
  75
  76 private:
  77
  78    int parseNextInt(uchar const *&tmp)
  79    {
  80        int r = std::atoi((char const *)tmp);
  81        while (*tmp != ' ' && *tmp != 0) ++tmp;
  82        if (*tmp == ' ')
  83            ++tmp;
  84        return r;
  85    }
  86
  87    void PssmMatch(const unsigned pos, const double score)
  88    {
  89        // ALPHABET has been defined in Query.cpp
  90        for (int i = 0; i < 4; ++i)
  91        {
  92            const char c = ALPHABET[i];
  93            if (!pushChar(c)) continue;
  94
  95            const double nval = P[pos + i * patlen];
  96
  97            if (pos == 0)
  98            {
  99                if (score + nval >= threshold)
 100                    newMatch();
 101            }
 102            else
 103            {
 104                const double achievableScore = score + nval + maxSum[pos-1];
 105                if (achievableScore >= threshold)
 106                    PssmMatch(pos - 1, score + nval);
 107            }
 108            popChar();
 109        }
 110    }
 111
 112     inline void newMatch()
 113     {
 114         ulong min = smin.top();
 115         ulong max = smax.top();
 116         for (ulong i = min; i <= max; i++)
 117             result.insert(tc->getDoc(i));
 118     }
 119
 120
 121     /**
 122      * Convert Position Frequency Matrix to Log-Odds
 123      *
 124      * Calculates log-odds and pseudocount into a PSSM.
 125      *
 126      * @param P Position Frequency Matrix.
 127      * @param l width of the input matrix.
 128      * @param threshold a threshold score.
 129      * @return a PSSM.
 130      */
 131     std::vector<double> convertPfmToLogOdds(std::vector<double> &P, unsigned l)
 132     {
 133         // Frequencies of different characters in human genome (calculated from NCBI Build 36.3)
 134         const double freq[] = { 790648081, 549365549, 549611686, 791675782 };
 135         const double freqSum = freq[0] + freq[1] + freq[2] + freq[3];
 136         const double logSum = std::log(freqSum);
 137
 138         double logArr[4];
 139         logArr[0] = std::log(freq[0]) - logSum;  // equals log (freq['A']/sum)
 140         logArr[1] = std::log(freq[1]) - logSum;
 141         logArr[2] = std::log(freq[2]) - logSum;
 142         logArr[3] = std::log(freq[3]) - logSum;
 143
 144         std::vector<double> R = std::vector<double>(l*4);
 145         for(unsigned i=0; i < l; i++)
 146         {
 147             // Sum + pseudocount
 148             double sum = P[i] + P[i+l] + P[i + 2*l] + P[i + 3*l] + 1;
 149             for(int j =0; j<4; j++)
 150             {
 151                 // Propability p' (background probability, used as a pseudocount)
 152                 double backProb = freq[j] / freqSum;
 153                 // Log of probability p with pseudocount added
 154                 double logProb = std::log(P[i+j*l] + backProb) - log(sum);  // equals log((P + pseudo)/sum)
 155
 156                 R[i+j*l] = logProb - logArr[j]; // equals log(p/p')
 157
 158 //                cout << "i = " << i+j*l << ", R = " << R[i+j*l] << ", back = " << backProb << ", logp= " << logProb << ", logA = " << logArr[j] << endl;
 159             }
 160         }
 161         return R;
 162     }
 163
 164     void calculateMaxSums(std::vector<double> P, unsigned l)
 165     {
 166         maxSum = std::vector<double>(l);
 167         // Max sum so far
 168         double sum = 0;
 169         // Iterate through columns in P
 170         for (unsigned i = 0; i < l;  ++i)
 171         {
 172             double value = P[i];
 173             // Scan the rest of column i
 174             for (unsigned j = 1; j < 4; j ++)
 175                 if (P[i + j*l] > value)
 176                     value = P[i + j*l]; // new maximum value in column i
 177             sum += value;
 178             maxSum[i] = sum;
 179 //                cout << "maxsum[" << i << "] = " << sum << endl;
 180         }
 181     }
 182
 183
 184     std::vector<double> P;
 185     double threshold;
 186     std::vector<double> maxSum;
 187
 188     PssmQuery();
 189     // No copy constructor or assignment
 190     PssmQuery(PssmQuery const&);
 191     PssmQuery& operator = (PssmQuery const&);
 192 };
 193
 194 } // namespace
 195
 196 #endif // _PssmQuery_H_