Confused about the calculation of BLOSUM matrix

Question

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?

Answer 1

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).