11
$\begingroup$

I have a weight matrix of length 20 x 15 (amino acids x sequence positions). Each element of my weight matrix is a relative probability

If I have a sequence say "AAPGTGASMHSGLLW" how would I score it against the matrix? I tried taking the product of probabilities corresponding to the matrix, but I end up with a really small number

Any ideas?

Edit:

Consider the simple matrix:

    1    2   3     4
A 0.3 0.90 0.5 0.0001
B 0.2 0.05 0.4 0.2
C 0.5 0.05 0.1 0.8

The best match is, with a score of:

CAAC = 0.5 * 0.9 * 0.5 * 0.8 = 0.18

If you change the first letter to an B instead of C

you get a match, with a score of:

BAAC = 0.2 * 0.9 * 0.5 * 0.8 = 0.072

Which is a huge difference for such a small change... This is even worse with my larger matrix since the score is easily affected by small probabilities

Improve this question
$\endgroup$
2
  • $\begingroup$ Maybe your sequence is a bad match? $\endgroup$
    – terdon
    Oct 11, 2012 at 19:05
  • $\begingroup$ I edited my question with an example. See above ^ $\endgroup$
    – Omar Wagih
    Oct 11, 2012 at 19:58

2 Answers 2

Reset to default
12
$\begingroup$

The probabilities are correct. You must take the product (in log space this is equivalent to sum). The reason the probability looks small is just that you are perhaps thinking the score should be close to 1. However, this is not the case. To get a score of 1, you need the PWM to have 1/0/0/0 at all positions and get a perfect match.

So what should you compare to? What people usually do is compare this to a background distribution, the easiest being uniform, so the PWM is 0.25 everywhere. For your example, the score in this case will be 0.25^4 = ~0.004 and this is what you should expect by random.

This is why people usually look at the ratio between the score of the PWM relative to the score for the background model (and usually take the log2 of that), which in your case will be 0.18/0.004 = ~46 so the sequence you got is 46 times more than you would expect by random! And for your second example, 0.072/0.004 = ~18 times more than expected, so that is still high.

More conceptually, what you are doing is comparing two probabilistic models, your PWM and a background PWM, and comparing the probability to get your observed sequence according to each one of them. This is a common approach in general for comparing probabilistic models, even if they are more complicated.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you for that detailed answer! I am trying to get scored between 0 and 1 (i.e. a probability). I thought about taking the ratio of the score to the best match to the PWM. I'm guessing from what you just explained, I would have to divide by 0.004 for the sequence score and the best score. So using the score from BAAC : (0.072/0.004)/(0.18/0.004)? Would this work? $\endgroup$
    – Omar Wagih
    Oct 12, 2012 at 5:03
  • $\begingroup$ Also, if I have real life background probabilities (i.e. determined from a large database) would I just take the product of everything? $\endgroup$
    – Omar Wagih
    Oct 12, 2012 at 5:30
  • $\begingroup$ @Omar regarding normalizing the score between 0 and 1, you can normalize the score in the way you suggested (note that the background model probability cancels out so you don't need it). The interpretation then will be different, as it will not tell you anything about how likely it is to find your sequence by random. So I guess it depends what you want to use the scores for. Note that as long as you divide all your sequences by a constant score, be it background or the highest probability, the ratios between the scores of different sequences will be maintained. $\endgroup$
    – Bitwise
    Oct 12, 2012 at 13:43
  • $\begingroup$ @Omar regarding background models, if you are not using uniform background and you are making predictions on a given organism, it would make sense to calculate the A/C/G/T percent of the whole genome, and use that as the background model, i.e. PWM with the same length as your PWM of interest, just that all positions have the same genomic probabilities. $\endgroup$
    – Bitwise
    Oct 12, 2012 at 13:46
  • $\begingroup$ Thanks. I am using amino acid sequences so I guess the background model would be (1/20)^15 then. I run into the problem when dividing the product of my score by the product of the best score since they are both very small numbers: for example 10^-19 / 10^-12 = 10^-7, which is indeed between 0 and 1 but it does not give me a sense of how similar it is to the best. What I am looking to score is how well some peptide matches a PWM (ranging from 0 to 1). Does this make sense? Appreciate the help $\endgroup$
    – Omar Wagih
    Oct 12, 2012 at 16:33
3
$\begingroup$

According to [this page][1], you should take the sum and not the product:

Once a Profile has been derived from a set of functionally related sites, the Profile can be used to scan a query sequence for the presence of potential sites. Usually you run a window the length of the matrix along the sequence, and sum the coefficients from the matrix corresponding to each nucleotide in each position on the window sequence. Formally, the score of a matrix M for a site s of length l (s = s1, ... , sl, and sk being one of {A, C, G, T}) is computed as

$m_s=\sum\limits_{j=1}^lM_{s_{lj}}$

I highly recommend you read the rest of the page, the author, Roderic Guigó, is an authority on the subject.

Improve this answer
$\endgroup$
5
  • $\begingroup$ He is assuming by Msij he means the log(frequency/background.frequency). I dont have the frequency matrix. Just the probabilities. $\endgroup$
    – Omar Wagih
    Oct 11, 2012 at 21:12
  • $\begingroup$ It shouldn't make a difference, the score for the sequence should still be the sum of the scores at each position. $\endgroup$
    – terdon
    Oct 11, 2012 at 21:23
  • 3
    $\begingroup$ @terdon actually it does make a difference, in probabilities it is a product but in log space it is mathematically equivalent to the sum... $\endgroup$
    – Bitwise
    Oct 12, 2012 at 1:05
  • $\begingroup$ @Bitwise, all I am saying is that the summ is still sufficient to score the match. Said sum will not, of course, be a probability but it will be maximized at the region of the highest sequence similarity. In any case, your answer is much better and you are clearly more knowledgeable on the subject than I so I won't argue the point. $\endgroup$
    – terdon
    Oct 12, 2012 at 9:21
  • $\begingroup$ @terdon not arguing, it is just important to correct so the information will be clear. Consider a PWM that in first position has A with probability 1 and second position C with probability 1. The sequence AA has a probability of zero for binding this sequence (1*0), but the sum will give you 1, which is half the maximal score and thus clearly wrong. PWM is a probabilistic model of the binding affinity so it has to be manipulated as probabilities. $\endgroup$
    – Bitwise
    Oct 12, 2012 at 13:37

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .