2
$\begingroup$

I am trying to understand how the BLOSUM matrix is calculated.

I found a very nice tutorial here: http://www.cs.columbia.edu/4761/assignments/assignment1/reference1.pdf

Also i read through the wikipedia article, and a book "Understanding bioinformatics" from Marketa Zvelebim and Jeremy O. Baum

After i thought i understand the method i tried myself with this "alignment" (i will refer to the variable as they are at the link i provided)

Q
Q
S

So it has only one column, and 2 type of amino acids.

The Count pair frequencies (C) for Q -> Q substitution is 1 (2*1/2) and for Q -> S it should be 2 (2*1).

To normalise it T = 3 as 3*2/2

The q values then seem to be q(Q->Q) = 1/3, q(Q->S) = 2/3

The p(Q) is then 1/3 + 2/3/2 = 4/6

p(S) is 2/3+1/3/2 = 2/6

The e(Q->Q) = (4/6)^2 = 4/9

e(Q->S) = 4/6*2/6*2 = 4/9

Thus

S(Q->Q) = round(log2((1/3)/(4/9)) = -1

S(Q->S) = round(log2((2/3)/(4/9)) = 1

So from this alignment it is more likely to substitute a Q to an S rather then Q to Q? This makes no sense for me.

What am i missing?

$\endgroup$
2
+50
$\begingroup$

Ok I'll try to go through it myself and see what happens.

First, we get the pair counts: $$ c_{ij}=\begin{cases} \displaystyle\binom{n_i}{2} & i=j\\ n_i n_j & i\ne j \end{cases} \;\;\;\;\;\implies\;\;\;\;\; c_{11}=1, c_{12}=2, c_{22}=0 $$ Clearly, the normalizer is just $T=\sum_{i\geq j} c_{ij}=3$ and $c_{ij}\equiv c_{ij}^{(1)}$ since there's only one column. Then our pair distribution is: $$ q_{ij} = \frac{c_{ij}}{T} \;\;\;\;\;\implies\;\;\;\;\; q_{11}=\frac{1}{3},\;q_{12}=\frac{2}{3},\; q_{22}= 0 $$ Hence, the occurrence probability per residue: $$ p_i = q_{ii} + \sum_{i\ne j} \frac{q_{ij}}{2} \;\;\;\;\;\implies\;\;\;\;\; p_1 = \frac{2}{3},\;p_2=\frac{1}{3} $$ and the expected pair frequency: $$ e_{ij} = \begin{cases} p_i^2 & i=j\\ 2p_ip_j & i\ne j \end{cases}\;\;\;\;\;\implies\;\;\;\;\; e_{11}=e_{12}=\frac{4}{9},\; e_{22}=\frac{1}{9} $$ Now we can finally compute the odds ratios: $$ r_{ij} = \frac{q_{ij}}{e_{ij}} \;\;\;\;\;\implies\;\;\;\;\; r_{11}=\frac{3}{4},\; r_{12}=\frac{3}{2}, r_{22}=0 $$ Computing the log-odds $s_{ij} = \log_2( r_{ij} )$ gives the same values as you get.

So why is $r_{11} > r_{12}$? Well, notice they have the same normalizer, i.e. $e_{11}=e_{12}$. So the difference comes from the fact that $ q_{11} > q_{12} $, the empirical pair distribution. Essentially, because of the way we computed $q_{ij}$, we consider there to be only one pair QQ but two pairs of QS. I would also note that sometimes results are a little unintuitive for very small samples (data sparsity throws off the statistics essentially).

$\endgroup$
2
  • 1
    $\begingroup$ Well I don't know whether this is right or wrong (I suspect right) but considering the irrelevant nonsense that gets up voted for other questions this surely deserves a +1. $\endgroup$
    – David
    Sep 16 '17 at 18:39
  • 1
    $\begingroup$ This is what i thought, thabks for the confirmation! $\endgroup$ Sep 17 '17 at 7:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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