silico.biotoul.fr
 

M1 MABS BBS Math TD Calcul Matriciel

From silico.biotoul.fr

Jump to: navigation, search

Contents

Création d'une matrice

A = matrix( c(11, 12, 13, 21, 22, 23, 31, 32, 33), ncol=3, byrow=TRUE )

A partir d'un fichier Media:gold.metadata.txt :

G=read.table("Gold.metadata.txt", sep="\t", header=TRUE, row.names=1)
class(G)

ou à partir de l'URL du fichier :

G=read.table("http://silico.biotoul.fr/site/images/2/21/Gold.metadata.txt", sep="\t", header=TRUE, row.names=1)
class(G)


G est un data frame ; la première colonne contient le nom de l'organisme. Pour ne garder que les colonnes numériques et convertir le data.frame en matrice, on fait :

G=as.matrix( G[ , 2:11] )

Cas de la matrice unité d'ordre n (notée In), exemple avec n = 3 :

I3 = diag(3)

Cas des vecteurs colonne ou ligne :

V = c(0.1, 10, 1000)
class(V)
# obtention de la matrice vecteur colonne :
as.matrix(V)
# pour obtenir la matrice vecteur ligne, on fait la transposée du résultat pécédent :
t( as.matrix(V) )

Opérations sur les matrices

  • addition +, soustraction -, multiplication par un nombre *,
  • produit matriciel :
    • entre 2 matrices A et I3 : A %*% I3,
    • entre une matrice A et un vecteur v : A %*% v
  • transposition : t(M)
  • trace : sum( diag(M) )
  • déterminant : det(M)
  • inversion d'une matrice : solve(M)
  • calcul des valeurs propres et des vecteurs propres : eigen(M)


Comparez ce que vous obtenez avec la multiplication et le produit matriciel :

A * V
A %*% V

Application aux matrices de substitution pour les alignements de séquences

Rappel : Quelques diapositives sur l'alignement de 2 séquences et les matrices de substitution

Calcul d'une matrice de substitution à partir d'alignements

Nous allons construire une matrice de type Dayhoff (PAM) pour l'alignement de séquences protéiques, ceci à partir d'alignements de séquences. La méthode en quelques mots consiste à

  • compter le nombre de substitutions observées pour chaque acide aminé (a.a.),
  • à transformer ce nombre d'occurrence en fréquence,
  • puis à calculer la log odds matrix : la matrice contenant le coût attribué à chaque substitution (match et mismatch).

Source : A Model of Evolutionary Change in Proteins, M.O Dayhoff, R.M. Schwartz and B.C. Orcutt, Atlas of Proteins Sequence and Structure, 1978

La première étape (déjà réalisée pour vous) consiste à réaliser les alignements et à dénombrer le nombre de passage d'un acide aminé à un autre.

Fichier contenant, le nombre de substitutions observées entre chaque paire d'acides aminés dans un lot de séquences préalablement alignées File:Count matrix.txt.

count.mat=as.matrix(read.table("count_matrix.txt", sep="\t", header=TRUE, row.names=1))

Matrice de mutation (mutation.mat)

Après avoir chargé cette matrice dans R, transformez-la en une matrice mutation.mat qui contient des fréquences (aidez-vous des fonctions apply et sum). Ces fréquences sont interprétées comme des probabilités de passer d'un acide aminé en un autre : mutation.mat[i,j] = p_{j \rightarrow i} correspond à la probabilité que l'a.a. j soit remplacé par l'a.a. i.

Entraînez-vous avec la fonction apply pour faire la somme/moyenne/écart-type des lignes ou des colonnes :

apply(count.mat, 2, sum)
apply(count.mat, 1, sum) # pareil car c'est symétrique
sums = apply(count.mat, 2, sum)
sums

La fréquence correspond au nombre d’occurrences / total de la colonne :

round(   count.mat[ 1 ,   ] / sums[1]    ,3) # ligne A / sum(A) → colonne A de mutation.mat (mais vecteur !)
round(   count.mat[ 2 ,   ] / sums[2]    ,3) # ligne R / sum(R) → colonne R de mutation.mat
round(   count.mat[ 1:2 , ] / sums[1:2]  ,3) # lignes A & R → colonnes A et R

Vous devriez obtenir ceci (10 premières lignes et 10 premières colonne, arrondi à 2 chiffres après la virgule et multiplié par 100 pour l'affichage des questions de lisibilité. Il ne faut pas stocker les valeurs arrondies mais les valeurs avec le plus de précisions possibles pour les calculs dans mutation.mat) :

 round(mutation.mat[1:10,1:10] ,2)*100
   A  R  N  D  C  Q  E  G  H  I
A 54  3  3  3  8  4  5  5  2  2
R  2 57  3  1  2  6  3  1  4  1
N  1  2 48  6  1  3  3  3  5  1
D  2  1  9 57  1  4 11  2  3  0
C  1  0  0  0 57  0  0  0  0  0
Q  2  4  3  2  1 42  5  1  4  1
E  4  4  5 13  1 11 52  2  3  1
G  5  2  5  3  2  2  2 77  2  0
H  0  1  2  1  1  2  1  0 54  0
I  2  1  1  0  3  1  1  0  1 50

Vérifiez que les sommes des colonnes valent bien 1. Vous observerez que les sommes des lignes ne valent pas 1 car la matrice n'est pas symétrique.

Fréquences d'apparition des aa et matrice de scores de substitution

Pour un alignement de l'acide aminé AAi avec l'acide aminé AAj, on utilise la formule suivante :

\frac{P(\mathrm{alignement. AA_i. avec. AA_j. est. du. a. l'evolution})}{P(\mathrm{alignement. AA_i. avec. AA_j. est. du. au. hasard})} = \frac{F_i.M_{j,i}}{F_i.F_j} = \frac{F_j.M_{i,j}}{F_i.F_j}

avec F la fréquence de chacun des acides aminés (ici, les fréquences d'apparition dans les alignements de départ), et, Mi,j la probabilité que l'a.a. j soit remplacé par l'a.a. i (donc M correspond à la transposée de la matrice mutation.mat précédemment calculée).



On en déduit la vraisemblance d'un alignement entre AAi et AAj : L_{i,j} = \frac{F_j.M_{i,j}}{F_i.F_j} = \frac{M_{i,j}}{F_i}


En pratique, on utilise les logarithmes car log( S_1 \times S_2 \times \dots \times S_n ) = log S_1 + log S_2 + \dots + log S_n.


La matrice de substitution se calcule donc à l'aide de la formule suivante : S_{i,j} = 10 \times log_{10}\frac{M_{i,j}}{F_i}


Nous avons donc besoin de calculer la fréquences de chaque acide aminé dans une séquence (ici, fréquences d'apparition dans nos alignements de départ). Utilisez donc la matrice count.mat de départ qui contient le nombre de substitutions observées pour calculer le vecteur aa.freq contenant la fréquence de chacun de a.a.

Vous devriez obtenir le vecteur suivant :

round(aa.freq,2)
    A    R    N    D    C    Q    E    G    H    I    L    K    M    F    P    S    T    W    Y    V 
 0.08 0.05 0.04 0.06 0.01 0.03 0.07 0.08 0.02 0.07 0.09 0.07 0.02 0.04 0.04 0.06 0.06 0.01 0.03 0.08


Il n'y a plus qu'à appliquer la formule S_{i,j} = 10 \times log_{10}\frac{M_{i,j}}{F_i}.

Vous devriez obtenir la matrice de score/substitution suivante :

round(score.mat[1:10,1:10],1)
     A    R    N     D     C    Q     E     G    H     I
A  8.1 -5.3 -4.4  -4.5  -0.2 -3.0  -2.3  -2.1 -6.2  -6.0
R -5.3 10.5 -2.7  -6.0  -4.8  0.5  -2.9  -6.4 -1.6  -8.4
N -4.4 -2.7 10.9   1.9  -4.5 -1.4  -1.8  -1.8  0.7  -8.3
D -4.5 -6.0  1.9  10.0  -9.3 -1.7   2.9  -4.5 -2.2 -13.1
C -0.2 -4.8 -4.5  -9.3  16.6 -8.4 -11.3  -6.7 -6.1  -4.3
Q -3.0  0.5 -1.4  -1.7  -8.4 10.8   1.8  -6.5  0.1  -8.1
E -2.3 -2.9 -1.8   2.9 -11.3  1.8   8.8  -5.7 -3.3 -10.5
G -2.1 -6.4 -1.8  -4.5  -6.7 -6.5  -5.7   9.9 -7.1 -14.3
H -6.2 -1.6  0.7  -2.2  -6.1  0.1  -3.3  -7.1 14.1  -9.8
I -6.0 -8.4 -8.3 -13.1  -4.3 -8.1 -10.5 -14.3 -9.8   8.6

Calcul de PAM1

La matrice PAM1 est bien adaptée pour aligner des séquences dont la distance évolutive correspond au temps nécessaire pour qu'une mutation sur 100 a.a. soit fixée dans la population (PAM = Point Percent Accepted Mutation). En d'autres termes, chaque a.a. a une probabilité de 0.01 de muter. Pour apprécier la distance évolutive entre les séquences que nous avions au départ, nous allons calculer le pourcentage moyen de mutation d'un acide aminé. Vous devriez obtenir quelque chose comme 44,53783%.


Pour obtenir PAM1, il faut donc ajuster la probabilité de muter par λ : M_{i,j} = \lambda \frac{m_j A_{i,j} }{\sum_{i} A_{i,j} }

avec mj la mutabilité de l'a.a. j et Ai,j le nombre d'occurrences où l'a.a. i a été muté en a.a. j et λ un paramètre d'ajustement pour que l'on observe une mutation acceptée (au cours de l'évolution) sur 100 a.a..

Calcul de λ

On veut une probabilité de mutation de 1% pour PAM1 calculée par p = \lambda \sum m_i F_i = 0.01, donc  \lambda = \frac{ 0.01 }{ \sum_i m_i F_i}

Vous devriez obtenir 0.02245282.


Calcul de mutation.mat.1 pour PAM1 avec M_{i,j} = \lambda \frac{m_j A_{i,j} }{\sum_i A_{i,j} }

Vous pouvez donc calculer la matrice mutation.mat.1 et obtenir

round(mutation.mat.1,5)[,1:10]
        A       R       N       D       C       Q       E       G       H       I
A 0.01223 0.00056 0.00070 0.00067 0.00181 0.00096 0.00112 0.00116 0.00046 0.00047
R 0.00034 0.01276 0.00062 0.00029 0.00038 0.00128 0.00059 0.00026 0.00079 0.00017
N 0.00032 0.00046 0.01068 0.00133 0.00030 0.00063 0.00057 0.00057 0.00102 0.00013
D 0.00045 0.00033 0.00199 0.01277 0.00015 0.00087 0.00252 0.00046 0.00078 0.00006
C 0.00027 0.00009 0.00010 0.00003 0.01284 0.00004 0.00002 0.00006 0.00007 0.00010
Q 0.00039 0.00087 0.00057 0.00052 0.00011 0.00933 0.00119 0.00017 0.00080 0.00012
E 0.00091 0.00079 0.00103 0.00302 0.00011 0.00236 0.01163 0.00041 0.00072 0.00014
G 0.00108 0.00041 0.00117 0.00064 0.00038 0.00040 0.00048 0.01735 0.00034 0.00007
H 0.00011 0.00032 0.00055 0.00028 0.00011 0.00048 0.00022 0.00009 0.01209 0.00005
I 0.00038 0.00022 0.00023 0.00007 0.00057 0.00024 0.00014 0.00006 0.00016 0.01120
L 0.00059 0.00039 0.00026 0.00015 0.00057 0.00068 0.00023 0.00011 0.00045 0.00299
K 0.00061 0.00315 0.00129 0.00073 0.00023 0.00212 0.00142 0.00041 0.00105 0.00025
M 0.00020 0.00018 0.00011 0.00005 0.00019 0.00035 0.00010 0.00004 0.00018 0.00069
F 0.00014 0.00012 0.00011 0.00006 0.00042 0.00012 0.00004 0.00005 0.00062 0.00050
P 0.00060 0.00021 0.00023 0.00028 0.00011 0.00033 0.00032 0.00014 0.00039 0.00008
S 0.00167 0.00060 0.00128 0.00084 0.00115 0.00087 0.00073 0.00069 0.00055 0.00019
T 0.00084 0.00051 0.00099 0.00050 0.00102 0.00070 0.00066 0.00023 0.00046 0.00054
W 0.00002 0.00005 0.00002 0.00002 0.00013 0.00004 0.00002 0.00002 0.00015 0.00004
Y 0.00011 0.00014 0.00025 0.00010 0.00030 0.00020 0.00010 0.00004 0.00110 0.00013
V 0.00120 0.00027 0.00028 0.00012 0.00157 0.00045 0.00036 0.00012 0.00031 0.00453

Cette mise à l'échelle avec λ ne permet pas de calculer la diagonale de la matrice (non mutation). On peut la calculer avec M_{i,i} = 1 - \sum_{j=1, i \ne j}^{20} M_{i,j} ou plus simplement Mi,i = 1 − λmi

Vous devriez obtenir PAM1

round(mutation.mat.1,5)[,1:10]
        A       R       N       D       C       Q       E       G       H       I       L       K       M
A 0.98978 0.00056 0.00070 0.00067 0.00181 0.00096 0.00112 0.00116 0.00046 0.00047 0.00056 0.00080 0.00074
R 0.00034 0.99031 0.00062 0.00029 0.00038 0.00128 0.00059 0.00026 0.00079 0.00017 0.00023 0.00248 0.00039
N 0.00032 0.00046 0.98822 0.00133 0.00030 0.00063 0.00057 0.00057 0.00102 0.00013 0.00011 0.00076 0.00017
D 0.00045 0.00033 0.00199 0.99032 0.00015 0.00087 0.00252 0.00046 0.00078 0.00006 0.00010 0.00064 0.00011
C 0.00027 0.00009 0.00010 0.00003 0.99039 0.00004 0.00002 0.00006 0.00007 0.00010 0.00008 0.00004 0.00010
Q 0.00039 0.00087 0.00057 0.00052 0.00011 0.98688 0.00119 0.00017 0.00080 0.00012 0.00027 0.00113 0.00053
E 0.00091 0.00079 0.00103 0.00302 0.00011 0.00236 0.98918 0.00041 0.00072 0.00014 0.00018 0.00150 0.00030
G 0.00108 0.00041 0.00117 0.00064 0.00038 0.00040 0.00048 0.99490 0.00034 0.00007 0.00010 0.00050 0.00014
H 0.00011 0.00032 0.00055 0.00028 0.00011 0.00048 0.00022 0.00009 0.98964 0.00005 0.00010 0.00034 0.00016
I 0.00038 0.00022 0.00023 0.00007 0.00057 0.00024 0.00014 0.00006 0.00016 0.98875 0.00228 0.00026 0.00202
L 0.00059 0.00039 0.00026 0.00015 0.00057 0.00068 0.00023 0.00011 0.00045 0.00299 0.99102 0.00045 0.00419
K 0.00061 0.00315 0.00129 0.00073 0.00023 0.00212 0.00142 0.00041 0.00105 0.00025 0.00033 0.98868 0.00049
M 0.00020 0.00018 0.00011 0.00005 0.00019 0.00035 0.00010 0.00004 0.00018 0.00069 0.00109 0.00018 0.98698
F 0.00014 0.00012 0.00011 0.00006 0.00042 0.00012 0.00004 0.00005 0.00062 0.00050 0.00088 0.00008 0.00084
P 0.00060 0.00021 0.00023 0.00028 0.00011 0.00033 0.00032 0.00014 0.00039 0.00008 0.00017 0.00031 0.00012
S 0.00167 0.00060 0.00128 0.00084 0.00115 0.00087 0.00073 0.00069 0.00055 0.00019 0.00028 0.00079 0.00039
T 0.00084 0.00051 0.00099 0.00050 0.00102 0.00070 0.00066 0.00023 0.00046 0.00054 0.00035 0.00063 0.00065
W 0.00002 0.00005 0.00002 0.00002 0.00013 0.00004 0.00002 0.00002 0.00015 0.00004 0.00005 0.00001 0.00005
Y 0.00011 0.00014 0.00025 0.00010 0.00030 0.00020 0.00010 0.00004 0.00110 0.00013 0.00024 0.00012 0.00025
V 0.00120 0.00027 0.00028 0.00012 0.00157 0.00045 0.00036 0.00012 0.00031 0.00453 0.00158 0.00031 0.00137

Vérifiez que les sommes des colonnes valent 1.

On obtient ensuite la matrice de score correspondante à PAM1 avec la formule S_{i,j} = 10 \times log_{10}\frac{M_{i,j}}{F_i} :

round(PAM1,1)[,1:10]
       A     R     N     D     C     Q     E     G     H     I
 A  10.7 -21.8 -20.8 -21.0 -16.7 -19.5 -18.8 -18.6 -22.7 -22.5
 R -21.8  12.9 -19.2 -22.5 -21.3 -16.0 -19.4 -22.9 -18.1 -24.8
 N -20.8 -19.2  14.1 -14.6 -21.0 -17.9 -18.3 -18.3 -15.8 -24.7
 D -21.0 -22.5 -14.6  12.4 -25.8 -18.2 -13.6 -20.9 -18.7 -29.6
 C -16.7 -21.3 -21.0 -25.8  19.0 -24.9 -27.8 -23.2 -22.6 -20.8
 Q -19.5 -16.0 -17.9 -18.2 -24.9  14.5 -14.6 -23.0 -16.4 -24.6
 E -18.8 -19.4 -18.3 -13.6 -27.8 -14.6  11.6 -22.2 -19.8 -26.9
 G -18.6 -22.9 -18.3 -20.9 -23.2 -23.0 -22.2  11.0 -23.6 -30.8
 H -22.7 -18.1 -15.8 -18.7 -22.6 -16.4 -19.8 -23.6  16.8 -26.3
 I -22.5 -24.8 -24.7 -29.6 -20.8 -24.6 -26.9 -30.8 -26.3  11.6
 L -21.8 -23.5 -25.4 -27.8 -22.0 -21.2 -25.8 -29.2 -23.0 -14.7
 K -20.3 -13.1 -17.0 -19.5 -24.6 -14.9 -16.6 -22.0 -17.9 -24.1
 M -20.6 -21.2 -23.4 -27.1 -20.9 -18.2 -23.7 -27.4 -21.0 -15.3
 F -24.4 -25.0 -25.4 -28.2 -19.6 -24.9 -29.3 -28.8 -17.9 -18.8
 P -18.3 -22.8 -22.4 -21.6 -25.5 -20.9 -21.0 -24.7 -20.2 -26.8
 S -15.3 -19.8 -16.5 -18.3 -16.9 -18.1 -18.9 -19.2 -20.1 -24.7
 T -18.3 -20.5 -17.6 -20.6 -17.5 -19.1 -19.4 -23.9 -21.0 -20.2
 W -27.2 -22.7 -25.5 -27.2 -18.2 -23.3 -26.3 -25.6 -17.7 -23.2
 Y -24.1 -23.0 -20.5 -24.5 -19.6 -21.3 -24.6 -28.9 -14.1 -23.3
 V -18.1 -24.5 -24.5 -28.3 -17.0 -22.4 -23.4 -28.1 -24.0 -12.4

Autres matrices de la famille PAM : PAM120, ...

Pour obtenir les matrices de score correspondant à des séquences plus divergentes, il faut multiplier la matrice par elle-même (produit matriciel).

Élévation d'une matrice à une puissance : M^r = \underbrace{M \times M \times \dots \times M}_{r} (produits matriciels)

Si M est diagonalisable, il existe une matrice P et une matrice D vérifiant : M = P \times D \times P^{-1}

L'élévation à une puissance d'une matrice diagonale est facile :


\begin{pmatrix}
a & 0 & 0 \\
0 & b & 0 \\
0 & 0 & c \end{pmatrix} 
\times 
\begin{pmatrix}
d & 0 & 0 \\
0 & e & 0 \\
0 & 0 & f \end{pmatrix}
=
\begin{pmatrix}
a \times d & 0 & 0 \\
0 & b \times e & 0 \\
0 & 0 & c\times f \end{pmatrix}
et donc 
\begin{pmatrix}
a & 0 & 0 \\
0 & b & 0 \\
0 & 0 & c \end{pmatrix}^r
=
\begin{pmatrix}
a^r & 0   & 0 \\
0   & b^r & 0 \\
0   & 0   & c^r \end{pmatrix}


Et donc par exemple pour M^3 = P \times D \times P^{-1} \times P \times D \times P^{-1} \times P \times D \times P^{-1} = P \times D \times I \times D \times I \times D \times P =   P \times D^3 \times P^{-1}


Diagonalisation d'une matrice : P correspond à la matrice composée des vecteurs propres (en colonne) et D correspond à une matrice diagonale composée des valeurs propres.

Exemple avec une petite matrice :

A = matrix( c(1,0,2, 2,3,-4, 0,0,2), ncol=3)
A
     [,1] [,2] [,3]
[1,]    1    2    0
[2,]    0    3    0
[3,]    2   -4    2
 
eigen(A)
$values
[1] 3 2 1
 
$vectors
           [,1] [,2]       [,3]
[1,]  0.4082483    0  0.4472136
[2,]  0.4082483    0  0.0000000
[3,] -0.8164966    1 -0.8944272
 
P = eigen(A)$vectors
D = diag( eigen(A)$values )
 
P %*% D %*% solve(P)
     [,1] [,2] [,3]
[1,]    1    2    0
[2,]    0    3    0
[3,]    2   -4    2

Ainsi mutation.mat.120 sera

round(mutation.mat.120, 4)[,1:10]
       A      R      N      D      C      Q      E      G      H      I
A 0.3284 0.0501 0.0579 0.0563 0.0954 0.0635 0.0688 0.0748 0.0443 0.0472
R 0.0302 0.3414 0.0435 0.0323 0.0282 0.0660 0.0441 0.0254 0.0489 0.0184
N 0.0262 0.0327 0.2635 0.0605 0.0228 0.0372 0.0388 0.0359 0.0480 0.0131
D 0.0382 0.0363 0.0906 0.3472 0.0214 0.0580 0.1098 0.0387 0.0508 0.0128
C 0.0141 0.0069 0.0074 0.0047 0.3192 0.0056 0.0046 0.0059 0.0063 0.0093
Q 0.0260 0.0449 0.0337 0.0351 0.0154 0.2274 0.0508 0.0175 0.0391 0.0143
E 0.0559 0.0595 0.0696 0.1317 0.0255 0.1007 0.3145 0.0392 0.0537 0.0208
G 0.0697 0.0394 0.0739 0.0534 0.0370 0.0398 0.0450 0.5569 0.0352 0.0144
H 0.0109 0.0199 0.0260 0.0184 0.0105 0.0234 0.0162 0.0092 0.2967 0.0073
I 0.0380 0.0246 0.0232 0.0152 0.0503 0.0280 0.0205 0.0124 0.0239 0.3219
L 0.0485 0.0372 0.0304 0.0227 0.0541 0.0474 0.0290 0.0180 0.0416 0.1516
K 0.0449 0.1313 0.0690 0.0560 0.0288 0.0932 0.0751 0.0362 0.0626 0.0248
M 0.0142 0.0127 0.0098 0.0071 0.0153 0.0172 0.0096 0.0058 0.0133 0.0353
F 0.0149 0.0134 0.0133 0.0090 0.0288 0.0146 0.0093 0.0073 0.0396 0.0358
P 0.0373 0.0209 0.0224 0.0244 0.0162 0.0258 0.0262 0.0164 0.0273 0.0122
S 0.0726 0.0427 0.0623 0.0518 0.0595 0.0503 0.0488 0.0464 0.0395 0.0234
T 0.0521 0.0392 0.0542 0.0401 0.0584 0.0455 0.0445 0.0266 0.0360 0.0397
W 0.0026 0.0038 0.0030 0.0022 0.0083 0.0036 0.0024 0.0024 0.0097 0.0044
Y 0.0109 0.0126 0.0162 0.0102 0.0207 0.0149 0.0101 0.0060 0.0521 0.0151
V 0.0645 0.0303 0.0300 0.0218 0.0844 0.0377 0.0316 0.0192 0.0312 0.1783

Et les sommes des colonnes valent toujours 1.

Et la matrice de scores de substitutions correspondant à PAM120 :

round(PAM120, 1)[1:10,1:10]
      A    R    N    D    C    Q    E    G    H    I
 A  5.9 -2.3 -1.7 -1.8  0.5 -1.3 -0.9 -0.5 -2.8 -2.5
 R -2.3  8.3 -0.7 -2.0 -2.6  1.1 -0.6 -3.0 -0.2 -4.4
 N -1.7 -0.7  8.4  2.0 -2.3 -0.1  0.0 -0.3  1.0 -4.7
 D -1.8 -2.0  2.0  7.8 -4.3  0.0  2.8 -1.7 -0.5 -6.5
 C  0.5 -2.6 -2.3 -4.3 14.1 -3.5 -4.3 -3.3 -3.0 -1.3
 Q -1.3  1.1 -0.1  0.0 -3.5  8.2  1.7 -3.0  0.5 -3.9
 E -0.9 -0.6  0.0  2.8 -4.3  1.7  6.6 -2.4 -1.1 -5.2
 G -0.5 -3.0 -0.3 -1.7 -3.3 -3.0 -2.4  8.5 -3.5 -7.4
 H -2.8 -0.2  1.0 -0.5 -3.0  0.5 -1.1 -3.5 11.5 -4.5
 I -2.5 -4.4 -4.7 -6.5 -1.3 -3.9 -5.2 -7.4 -4.5  6.7

Calcul de la matrice des fréquences de mutation et calcul des fréquences d'apparition des acides aminés à partir d'une matrice PAM

Nous allons partir de la matrice PAM250 similaire à la matrice de substitution que nous avons calculée précédemment.


Comme vu dans la section précédente, nous savons que S_{i,j} = 10 \times \log_{10}\frac{M_{i,j}}{F_i}.

D'où la différence entre 2 valeurs d'une même ligne de la matrice (PAM250) :

S_{i,j} - S_{i,k} = 10 \times \log\frac{M_{i,j}}{F_i} - 10 \times \log\frac{M_{i,k}}{F_i}

S_{i,j} - S_{i,k} = 10 \times (\log\frac{M_{i,j}}{F_i} - \log\frac{M_{i,k}}{F_i})

avec \log a - \log b = \log \frac{a}{b}

S_{i,j} - S_{i,k} = 10 \times \log \frac{ \frac{M_{i,j}}{F_i} }{ \frac{M_{i,k}}{F_i} }= 10 \times (\log (\frac{M_{i,j}}{F_i} \times \frac{F_i}{M_{i,k}}))

les Fi se simplifient et on obtient :

S_{i,j} - S_{i,k} = 10 \times (\log \frac{M_{i,j}}{M_{i,k}})

d'où

\frac{S_{i,j} - S_{i,k}}{10} = \log \frac{M_{i,j}}{M_{i,k}}

d'où

10^{\frac{S_{i,j} - S_{i,k}}{10}} = \frac{M_{i,j}}{M_{i,k}}


En supposant que Mi,1 = 1 (première colonne de mutation.mat.250 à 1), on peut en déduire les autres valeurs de la matrice pour la ligne i de mutation.mat.250. En utilisant la formule ci-dessus avec k = 1 :

\frac{M_{i,j}}{M_{i,1}} = \frac{M_{i,j}}{1} = 10^{\frac{S_{i,j} - S_{i,1}}{10}}


On obtient donc la matrice suivante :

round(mutation.mat.250[1:10,1:10],2)
  A    R    N    D    C    Q    E    G    H    I
A 1 0.44 0.49 0.49 0.68 0.51 0.54 0.59 0.40 0.42
R 1 3.89 1.35 1.12 0.89 1.78 1.38 0.87 1.45 0.65
N 1 1.20 2.63 1.82 0.87 1.29 1.38 1.26 1.48 0.54
D 1 1.00 1.82 3.31 0.63 1.38 2.14 1.00 1.20 0.39
C 1 0.58 0.63 0.46 7.59 0.54 0.49 0.52 0.56 0.81
Q 1 1.51 1.23 1.32 0.71 2.09 1.55 0.74 1.29 0.63
E 1 1.12 1.26 1.95 0.62 1.48 2.19 0.79 1.02 0.46
G 1 0.65 1.05 0.83 0.60 0.65 0.72 4.57 0.59 0.26
H 1 1.58 1.82 1.48 0.95 1.66 1.38 0.87 6.92 0.72
I 1 0.68 0.63 0.46 1.32 0.78 0.59 0.37 0.69 3.31


Il s'agit maintenant d'ajuster ces fréquences pour que leur somme sur une colonne soit égale à 1.

S'il on transpose cette matrice (A = MT), cela correspond à un système linéaire de 20 équations à 20 inconnues : Ax = b avec bi = 1 pour i \in [1, 20] :

\left\{\begin{matrix}  a_{1,1}x_1+a_{1,2}x_2+...+a_{1,20}x_{20} = b_1 = 1 \\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,20}x_{20} = b_2 = 1 \\ \vdots \\ \vdots \\ a_{20,1}x_{1}+a_{20,2}x_{2}+...+a_{20,20}x_{20} = b_{20} = 1\end{matrix}\right.


que l'on représente sous forme matricielle : Ax = b

avec

A=M^t=\begin{pmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,20} \\
a_{2,1} & a_{2,2} & \cdots & a_{2,20} \\
\vdots & \vdots & \ddots & \vdots \\
a_{20,1} & a_{20,2} & \cdots & a_{20,20} \end{pmatrix} \qquad x=\begin{pmatrix} x_1 \\ x_2\\ \vdots \\ x_{20} \end{pmatrix}\quad\text{et}\quad b=\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_{20} \end{pmatrix}


Si A est une matrice carrée inversible, on a un système de Cramer que l'on peut résoudre de la manière suivante :

Ax = b

En multipliant des deux côtés par A − 1, on obtient :

A − 1Ax = A − 1b

or A − 1A = I d'où A − 1Ax = Ix = x = A − 1b


Et nous avons la solution x du système linéaire pour que la somme des lignes de A soit égale à 1 (donc la somme des colonnes de MM soit égale à 1).

Comme A correspond au système suivant :


\left\{\begin{matrix}  a_{1,1}x_1+a_{1,2}x_2+...+a_{1,20}x_{20} = b_1 = 1 \\ a_{2,1}x_1+a_{2,2}x_2+...+a_{2,20}x_{20} = b_2 = 1 \\ \vdots \\ \vdots \\ a_{20,1}x_{1}+a_{20,2}x_{2}+...+a_{20,20}x_{20} = b_{20} = 1\end{matrix}\right.


On remarque que cela correspond à multiplier la 1ère colonne de A par x1, la 2ème colonne de A par x2, ...

Et comme A = MT, cela correspond à multiplier la 1ère ligne de mutation.mat.250 par x1, la 2ème ligne de mutation.mat.250 par x2, ...


Après avoir vérifié que A est inversible, déterminez la matrice mutation.mat.250. Vous devriez obtenir :

round(mutation.mat.250[1:10,1:10],2)
     A    R    N    D    C    Q    E    G    H    I
A 0.15 0.07 0.08 0.08 0.10 0.08 0.08 0.09 0.06 0.06
R 0.04 0.17 0.06 0.05 0.04 0.08 0.06 0.04 0.06 0.03
N 0.04 0.05 0.10 0.07 0.03 0.05 0.05 0.05 0.06 0.02
D 0.05 0.05 0.09 0.16 0.03 0.07 0.11 0.05 0.06 0.02
C 0.02 0.01 0.01 0.01 0.12 0.01 0.01 0.01 0.01 0.01
Q 0.03 0.05 0.04 0.04 0.02 0.07 0.05 0.02 0.04 0.02
E 0.07 0.08 0.09 0.14 0.04 0.11 0.16 0.06 0.07 0.03
G 0.09 0.06 0.09 0.07 0.05 0.06 0.06 0.39 0.05 0.02
H 0.02 0.02 0.03 0.02 0.01 0.03 0.02 0.01 0.10 0.01
I 0.05 0.03 0.03 0.02 0.07 0.04 0.03 0.02 0.04 0.17


Calculons à présent la fréquence d'apparition des a.a. dans une séquence. On sait que

F_i \times M_{j,i} = F_j \times M_{i,j}

d'où

{F_i} = F_j \frac{M_{i,j}}{M_{j,i}}

Comme précédemment, on suppose que pour j = 1,Fj = F1 = 1, ce qui nous permet de calculer les autres fréquences : F_i = \frac{M_{i,1}}{M_{1,i}}

Sachant que la somme des fréquences des a.a. doit être égale à 1, vous pouvez maintenant mettre à l'échelle ces valeurs que l'on vient d'obtenir pour que leur somme soit 1.

Vous devriez obtenir les fréquences suivantes :

round(PAM250.freq,2)
   A    R    N    D    C    Q    E    G    H    I    L    K    M    F    P    S    T    W    Y    V 
0.08 0.05 0.04 0.05 0.01 0.03 0.07 0.08 0.02 0.07 0.09 0.06 0.03 0.04 0.04 0.05 0.06 0.01 0.03 0.08

Ce qui devrait correspondre aux fréquences d'a.a. calculées dans la première partie à partir de la matrice count.mat.