Take the 2-minute tour ×
Biology Stack Exchange is a question and answer site for biology researchers, academics, and students. It's 100% free, no registration required.

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

share|improve this question
    
Maybe your sequence is a bad match? –  terdon Oct 11 '12 at 19:05
    
I edited my question with an example. See above ^ –  Omar Oct 11 '12 at 19:58
add comment

2 Answers 2

According to this page, 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

enter image description here

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

share|improve this answer
    
He is assuming by Msij he means the log(frequency/background.frequency). I dont have the frequency matrix. Just the probabilities. –  Omar Oct 11 '12 at 21:12
    
It shouldn't make a difference, the score for the sequence should still be the sum of the scores at each position. –  terdon Oct 11 '12 at 21:23
2  
@terdon actually it does make a difference, in probabilities it is a product but in log space it is mathematically equivalent to the sum... –  Bitwise Oct 12 '12 at 1:05
    
@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. –  terdon Oct 12 '12 at 9:21
    
@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. –  Bitwise Oct 12 '12 at 13:37
add comment

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.

share|improve this answer
    
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? –  Omar Oct 12 '12 at 5:03
    
Also, if I have real life background probabilities (i.e. determined from a large database) would I just take the product of everything? –  Omar Oct 12 '12 at 5:30
    
@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. –  Bitwise Oct 12 '12 at 13:43
    
@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. –  Bitwise Oct 12 '12 at 13:46
    
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 –  Omar Oct 12 '12 at 16:33
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.