Calcul matriciel, hommage à Margaret Dayhoff
1 Opérations R sur les matrices
1.1 Création à partir d’un vecteur
## [,1] [,2] [,3]
## [1,] 11 12 13
## [2,] 21 22 23
## [3,] 31 32 33
1.2 Création à partir
d’un data.frame
## [1] "data.frame"
## [1] "matrix" "array"
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 5.1 3.5 1.4 0.2
## [2,] 4.9 3.0 1.4 0.2
## [3,] 4.7 3.2 1.3 0.2
## [4,] 4.6 3.1 1.5 0.2
## [5,] 5.0 3.6 1.4 0.2
## [6,] 5.4 3.9 1.7 0.4
1.3 Matrice identité et diagonale d’une matrice
Matrice identité d’ordre \(n\), notée \(I_n\) ; pour \(n = 3\) :
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
Accès à la diagonale :
## [,1] [,2] [,3]
## [1,] 11 12 13
## [2,] 21 22 23
## [3,] 31 32 33
## [1] 11 22 33
Modification de la diagonale :
## [,1] [,2] [,3]
## [1,] 11 12 13
## [2,] 21 44 23
## [3,] 31 32 66
Création à partir de la diagonale :
## [,1] [,2] [,3] [,4]
## [1,] 6 0 0 0
## [2,] 0 7 0 0
## [3,] 0 0 8 0
## [4,] 0 0 0 9
1.4 Vecteurs colonnes ou vecteurs lignes
Un vecteur
## [1] "numeric"
Obtention de la matrice vecteur colonne :
## [,1]
## [1,] 1e-01
## [2,] 1e+01
## [3,] 1e+03
Pour obtenir la matrice vecteur ligne, on fait la transposée du résultat précédent :
## [,1] [,2] [,3]
## [1,] 0.1 10 1000
1.5 Opérations sur les matrices
- Addition
+
, soustraction-
, multiplication*
, … - Produit matriciel
- entre 2 matrices
A
etI3
:A %*% I3
, - entre une matrice A et un vecteur v :
A %*% v
- entre 2 matrices
- 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)
Comparaison de la multiplication et du produit matriciel :
## [,1] [,2] [,3]
## [1,] 1.1 1.2 1.3
## [2,] 210.0 440.0 230.0
## [3,] 31000.0 32000.0 66000.0
## [,1]
## [1,] 13121.1
## [2,] 23442.1
## [3,] 66323.1
2 Matrices de scores de substitution : PAM
Rappel : Quelques diapositives sur l’alignement de 2 séquences et les matrices de substitution
2.1 Calcul d’une matrice de substitution PAM250 à partir d’alignements
Nous allons construire une matrice de type Dayhoff (la PAM250) pour l’alignement de séquences protéiques, ceci à partir d’alignements de séquences. La méthode utilise
- une matrice comptage du nombre de substitutions observées pour chaque acide aminé (aa),
- la mutabilité de chaque aa,, c’est-à-dire la probabilité qu’un aa a de muter,
- la fréquence d’apparition de chaque aa dans une séquence quelconque
A partir de ces informations, la matrice de probabilité de mutations d’un aa en un autre peut-être calculée pour une distance évolutive correspondant au temps qu’il faut pour qu’une mutation sur 100 aa soit fixée dans la population.
A partir de cette matrice, il est ensuite possible de calculer les matrices correspondant à des distances évolutives plus grandes. Et à partir de ces dernières matrices, on pourra calculer les matrices de scrores de substitutions utilisables pour faire l’alignement de séquences.
Nous allons donc suivre le cheminement décrit dans 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 globaux par familles de séquences, inférer les arbres de parenté correspondant afin de dénombrer le nombre de passage d’un acide aminé à un autre.
2.1.1 Matrice de comptage
count.mat
Cette matrice de comptage correspond à la figure 80 :
Cette matrice est symétrique.
Le fichier correspondant à la figure a été reconstitué : dayhoff.fig80.counts.csv.
Chargez la matrice count.mat
avec la fonction
read_csv
de tidyverse
.
## ── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Rows: 20 Columns: 21
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): aa
## dbl (19): A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y
## lgl (1): V
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Vous devriez avoir :
Il va falloir la rendre symétrique en recopiant les valeurs.
Avec la fonction replace_na
, remplacez les NA par des
0.
Vous devriez avoir :
A présent, vous pouvez ajouter cette matrice à sa transposée. Il vous faudra d’abord transformer ce tibble en matrice (et donc se débarrasser de la 1ère colonne).
Vous devriez avoir :
## A R N D C Q E G H I L K M F P S T W Y V
## A 0 30 109 154 33 93 266 579 21 66 95 57 29 20 345 772 590 0 20 365
## R 30 0 17 0 10 120 0 10 103 30 17 477 17 7 67 137 20 27 3 20
## N 109 17 0 532 0 50 94 156 226 36 37 322 0 7 27 432 169 3 36 13
## D 154 0 532 0 0 76 831 162 43 13 0 85 0 0 10 98 57 0 0 17
## C 33 10 0 0 0 0 0 10 10 17 0 0 0 0 10 117 10 0 30 33
## Q 93 120 50 76 0 0 422 30 243 8 75 147 20 0 93 47 37 0 0 27
## E 266 0 94 831 0 422 0 112 23 35 15 104 7 0 40 86 31 0 10 37
## G 579 10 156 162 10 30 112 0 10 0 17 60 7 17 49 450 50 0 0 97
## H 21 103 226 43 10 243 23 10 0 3 40 23 0 20 50 26 14 3 40 30
## I 66 30 36 13 17 8 35 0 3 0 253 43 57 90 7 20 129 0 13 661
## L 95 17 37 0 0 75 15 17 40 253 0 39 207 167 43 32 52 13 23 303
## K 57 477 322 85 0 147 104 60 23 43 39 0 90 0 43 168 200 0 10 17
## M 29 17 0 0 0 20 7 7 0 57 207 90 0 17 4 20 28 0 0 77
## F 20 7 7 0 0 0 0 17 20 90 167 0 17 0 7 40 10 10 260 10
## P 345 67 27 10 10 93 40 49 50 7 43 43 4 7 0 269 73 0 0 50
## S 772 137 432 98 117 47 86 450 26 20 32 168 20 40 269 0 696 17 22 43
## T 590 20 169 57 10 37 31 50 14 129 52 200 28 10 73 696 0 0 23 186
## W 0 27 3 0 0 0 0 0 3 0 13 0 0 10 0 17 0 0 6 0
## Y 20 3 36 0 30 0 10 0 40 13 23 10 0 260 0 22 23 6 0 17
## V 365 20 13 17 33 27 37 97 30 661 303 17 77 10 50 43 186 0 17 0
Remarque : Pensez à attribuer le nom des lignes avec
rownames(count.mat) = colnames(count.mat)
On peut vérifier que le nombre mutations correspond à celui présenté dans l’article :
## [1] 1571.5
Ce qui correspond à 1572 dans l’article (ici, à diviser par 10 car la table a été multipliée par 10 dans la figure ; et aussi par 2 pour ne pas compter 2 fois les mutations).
2.1.2 Mutabilité relative
aa.mut
La mutabilité est la probabilité qu’a un aa de muter. Son calcul est illustré dans la figure 81 :
Et, elle nous est donnée pour chaque aa dans la table 21 :
Le fichier correspondant a été reconstitué : dayhoff.mutabilities.csv
Chargez ce fichier avec la fonction read_csv
et créez le
vecteur correspondant aa.mut
. Il faudra gardez le même
ordre que dans la matrice count.mat
.
## Rows: 20 Columns: 2
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): aa
## dbl (1): mut
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## 100 65 134 106 20 93 102 49 66 96 40 56 94 41 56 120 97 18 41 74
2.1.3 Fréquences
d’apparition des aa aa.freq
Les fréquences d’apparition des aa dans les séquences sont données dans la table 22 :
Chargez le fichier dayhoff.frequencies.csv correspondant
et créez le vecteur aa.freq
.
## Rows: 20 Columns: 2
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): aa
## dbl (1): freq
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## 0.087 0.041 0.040 0.047 0.033 0.038 0.050 0.089 0.034 0.037 0.085 0.081 0.015 0.040 0.051 0.070 0.058 0.010 0.030 0.065
Vérifiez que la somme des fréquences vaut bien 1
On obtient :
## [1] 1.001
Mettre à l’échelle les fréquences pour que leur somme soit 1.
## [1] 1
2.1.4 Matrice de
probabilité de mutation mutation.mat
Il s’agit maintenant de calculer la matrice de probabilité de mutation correspondant à une distance évolutive de 1-PAM (temps qu’il faut pour qu’une mutation sur 100aa soit fixée dans la population). Elle nous est donnée en figure 82 :
Ainsi que les formules pour son calcul :
Les \(A_{ij}\) correspondent à notre
matrice count.mat
. Les \(m_j\) correspondent à
aa.mut
.
Il nous faut donc déterminer le facteur \(\lambda\) permettant de calculer
mutation.mat
correspondant à 1-PAM.
Nous allons utiliser les fréquences aa.freq
pour
considérer une séquence de 100aa. Elle sera donc composée de (arrondi à
1 chiffre après la virgule) :
## A R N D C Q E G H I L K M F P S T W Y V
## 8.7 4.1 4.0 4.7 3.3 3.8 5.0 8.9 3.4 3.7 8.5 8.1 1.5 4.0 5.1 7.0 5.8 1.0 3.0 6.5
Chacun de ces aa sera donc muté en fonction de sa mutabilité relative
aa.mut
et on souhaite mettre le total de mutations à
l’échelle (avec \(\lambda\)) pour avoir
1 seule mutation sur 100aa.
Le taux de mutations se calcule avec la formule suivante :
\(\mathrm{rate} = \lambda \sum f_i m_i\)
avec \(f_i\) la fréquence
d’apparition des aa (aa.freq
) et \(m_i\) la mutabilité de chaque aa
(aa.mut
).
On peut donc en déduire \(\lambda\) pour un taux de 1% : \(\lambda = \frac{0.01}{\sum f_i m_i}\)
Utilisez la formule précédente pour calculer \(\lambda\)
Vous devriez obtenir :
## [1] 0.0001330303
A présent, il est possible de remplir la matrice de probabilité de mutation avec les formules
en dehors de la diagonale : \(M_{ij} = \frac{\lambda \cancel{m_j} A_{ij}}{\sum_{i \ne j} A_{ij}}\) (il y a vraisemblablement une erreur dans l’article : il ne faut pas utiliser \(m_j\) ici. On pourra utiliser une constante à la place → 1000)
et
pour la diagonale : \(M_{ii} = 1 - \lambda m_i\)
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## A 9867 1 5 7 1 4 11 25 1 3 4 2 1 1 15 33 25 0 1 15
## R 1 9914 1 0 0 5 0 0 4 1 1 20 1 0 3 6 1 1 0 1
## N 5 1 9822 23 0 2 4 7 10 2 2 14 0 0 1 18 7 0 2 1
## D 7 0 23 9859 0 3 35 7 2 1 0 4 0 0 0 4 2 0 0 1
## C 1 0 0 0 9973 0 0 0 0 1 0 0 0 0 0 5 0 0 1 1
## Q 4 5 2 3 0 9876 18 1 10 0 3 6 1 0 4 2 2 0 0 1
## E 11 0 4 35 0 18 9864 5 1 1 1 4 0 0 2 4 1 0 0 2
## G 25 0 7 7 0 1 5 9935 0 0 1 3 0 1 2 19 2 0 0 4
## H 1 4 10 2 0 10 1 0 9912 0 2 1 0 1 2 1 1 0 2 1
## I 3 1 2 1 1 0 1 0 0 9872 11 2 2 4 0 1 5 0 1 28
## L 4 1 2 0 0 3 1 1 2 11 9947 2 9 7 2 1 2 1 1 13
## K 2 20 14 4 0 6 4 3 1 2 2 9926 4 0 2 7 8 0 0 1
## M 1 1 0 0 0 1 0 0 0 2 9 4 9875 1 0 1 1 0 0 3
## F 1 0 0 0 0 0 0 1 1 4 7 0 1 9945 0 2 0 0 11 0
## P 15 3 1 0 0 4 2 2 2 0 2 2 0 0 9926 11 3 0 0 2
## S 33 6 18 4 5 2 4 19 1 1 1 7 1 2 11 9840 29 1 1 2
## T 25 1 7 2 0 2 1 2 1 5 2 8 1 0 3 29 9871 0 1 8
## W 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 9976 0 0
## Y 1 0 2 0 1 0 0 0 2 1 1 0 0 11 0 1 1 0 9945 1
## V 15 1 1 1 1 1 2 4 1 28 13 1 3 0 2 2 8 0 1 9902
Vous noterez que cela diffère de la figure 82, tout particulièrement en dehors de la diagonale. Ceci est certainement dû aux arrondis effectués dans les tableaux de l’article.
Cette matrice représente la probabilité qu’un aa à de muter. Les sommes des colonnes devraient donc être de 1. Vérifiez que ce nest en fait pas le cas.
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P
## 1.0021205 0.9960597 0.9917650 0.9946941 0.9985245 0.9939263 0.9953744 1.0011679 0.9951478 0.9934976 1.0007229 1.0005287 0.9899501 0.9974324 0.9975744
## S T W Y V
## 0.9988166 0.9971485 0.9979398 0.9967171 0.9986336
Afin de ne pas modifier la diagonale (qui assure que le taux de mutations soit de 0.01), nous allons mettre à l’échelle les autres valeurs pour que la somme des colonnes s’approchent le plus possible de 1.
Si l’on décompose pour la 1ère colonne :
## A R N D C Q E G H I L K M F P S T W Y V
## 9867 1 5 7 1 4 11 25 1 3 4 2 1 1 15 33 25 0 1 15
La valeur à atteindre pour la somme des autres que A est donc 133. Et le facteur de mise à l’échelle est donc de 133 divisé par la somme des autres que A :
## A
## 0.01330303
## A
## 0.01542356
## A
## 0.8625137
Calculez les facteur de mise à l’échelle pour l’ensemble des colonnes.
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## 0.86 1.84 1.86 1.60 2.25 1.96 1.52 0.85 2.24 2.04 0.88 0.93 5.09 1.89 1.48 1.08 1.28 7.16 2.51 1.16
Appliquez ces mises à l’échelle pour obtenir la matrice
mutation.mat
.
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## A 9867 2 9 10 3 8 17 21 2 6 4 2 6 2 22 35 32 0 2 18
## R 1 9914 1 0 1 10 0 0 10 3 1 19 4 1 4 6 1 8 0 1
## N 4 1 9822 36 0 4 6 6 21 3 1 13 0 1 2 20 9 1 4 1
## D 6 0 42 9859 0 6 53 6 4 1 0 3 0 0 1 4 3 0 0 1
## C 1 1 0 0 9973 0 0 0 1 1 0 0 0 0 1 5 1 0 3 2
## Q 3 9 4 5 0 9876 27 1 23 1 3 6 4 0 6 2 2 0 0 1
## E 10 0 7 56 0 35 9864 4 2 3 1 4 2 0 3 4 2 0 1 2
## G 21 1 12 11 1 2 7 9935 1 0 1 2 2 1 3 21 3 0 0 5
## H 1 8 18 3 1 20 1 0 9912 0 1 1 0 2 3 1 1 1 4 1
## I 2 2 3 1 2 1 2 0 0 9872 9 2 12 7 0 1 7 0 1 32
## L 3 1 3 0 0 6 1 1 4 22 9947 2 45 13 3 1 3 4 2 15
## K 2 37 25 6 0 12 7 2 2 4 1 9926 19 0 3 8 11 0 1 1
## M 1 1 0 0 0 2 0 0 0 5 8 4 9875 1 0 1 2 0 0 4
## F 1 1 1 0 0 0 0 1 2 8 6 0 4 9945 0 2 1 3 28 0
## P 13 5 2 1 1 8 3 2 5 1 2 2 1 1 9926 12 4 0 0 2
## S 28 11 34 7 11 4 6 16 2 2 1 7 4 3 17 9840 38 5 2 2
## T 22 2 13 4 1 3 2 2 1 11 2 8 6 1 5 32 9871 0 2 9
## W 0 2 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 9976 1 0
## Y 1 0 3 0 3 0 1 0 4 1 1 0 0 21 0 1 1 2 9945 1
## V 13 2 1 1 3 2 2 3 3 57 11 1 17 1 3 2 10 0 2 9902
On notera que cela correspond davantage à la figure 82.
Testez à nouveau si les sommes des colonnes valent 1.
## A R N D C Q E G H I L K M F P S T W Y V
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Pour finir sur cette partie, écrivez une fonction
mutation.rate
qui prend en paramètre la matrice de
probabilité de mutation et les fréquences des aa et qui renvoie le taux
de mutations.
Utilisation :
## [1] 0.01
La matrice mutation.mat
contient les
probabilités de mutation d’un aa en un autre pour une distance évolutive
correspondant à 1-PAM, c’est-à-dire le temps qu’il faut pour qu’une
mutation sur 100aa soit fixée dans la population.
2.1.5 Matrice de mutation pour PAM250
A partir de la 1ère matrice mutation.mat
calculée
précédemment, on peut maintenant dériver d’autres matrices pour des
distances évolutives plus élevées. Pour 2-PAM, il faut multiplier
(produit matriciel) la matrice par elle-même :
mutation.mat %*% mutation.mat
. Pour PAM250, il faut donc le
faire 250 fois.
Élévation d’une matrice à une puissance : \(M^r = \underbrace{M \times M \times \dots \times M}_{r}\)
Élévation d’une matrice diagonalisable à une puissance https://fr.wikiversity.org/wiki/Initiation_aux_matrices/Puissance_d%27une_matrice
Si \(M\) est diagonalisable, il existe une matrice inversible \(P\) et une matrice diagonale \(D\) vérifiant : \(M = P \times D \times P^{−1}\)
L’élévation 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} \]
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 \cancel{P^{-1}} \times \cancel{P} \times D \times \cancel{P^{-1}} \times \cancel{P} \times D \times P^{-1} = P \times D \times \cancel{I} \times D \times \cancel{I} \times D \times P = P \times D^3 \times P^{-1}\)
Cette relation permet de calculer la matrice \(M^r\) : \(M^r = P \times D^r \times P^{−1}\) plutôt que d’enchaîner les produits matriciels : \(M^r = \underbrace{M \times M \times \dots \times M}_{r}\)
Diagonalisation d’une matrice https://fr.wikipedia.org/wiki/Diagonalisation
Exemple :
## [,1] [,2] [,3]
## [1,] 1 2 0
## [2,] 0 3 0
## [3,] 2 -4 2
## eigen() decomposition
## $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
Pour coller avec l’exemple wikipedia :
Ev[,1] = Ev[,1] * 1/Ev[1,1] # mise à l échelle de la 1ère colonne
Ev[,3] = Ev[,3] * 1/Ev[1,3] # mise à l échelle de la 3ème colonne
Ev
## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 1 0 0
## [3,] -2 1 -2
Matrice \(P\) correspond à la matrice des vecteurs propres en vecteur colonne :
ainsi \(P\) diagonalise \(A\) :
## [,1] [,2] [,3]
## [1,] 1 2 0
## [2,] 0 3 0
## [3,] 2 -4 2
## [,1] [,2] [,3]
## [1,] 3 0 0
## [2,] 0 2 0
## [3,] 0 0 1
la diagonale obtenue contient les valeurs propres.
ou encore :
## [,1] [,2] [,3]
## [1,] 1 2 0
## [2,] 0 3 0
## [3,] 2 -4 2
## [,1] [,2] [,3]
## [1,] 1 2 0
## [2,] 0 3 0
## [3,] 2 -4 2
Utilisez la méthode ci-dessus pour calculer la matrice de probabilités de mutation correspondant à 250-PAM.
Donc pour calculer MM250 et PAM250 :
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## A 13 6 9 9 5 8 9 12 6 8 6 7 7 4 11 11 11 2 4 9
## R 3 17 4 3 2 5 3 2 6 3 2 9 4 1 4 4 3 7 2 2
## N 4 4 6 7 2 5 6 4 6 3 2 5 3 2 4 5 4 2 3 3
## D 5 3 8 11 1 7 10 5 5 3 2 5 3 1 4 5 5 1 2 3
## C 2 1 1 1 52 1 1 2 2 2 1 1 1 1 2 3 2 1 4 2
## Q 3 5 5 6 1 10 7 3 8 2 3 5 3 1 4 3 3 1 2 2
## E 5 4 7 11 1 9 12 5 6 3 2 5 3 1 4 5 5 1 2 3
## G 12 5 10 10 4 7 9 27 5 5 3 6 5 3 8 11 9 2 3 7
## H 2 5 5 4 2 7 4 2 15 2 2 3 2 2 3 3 2 2 3 2
## I 3 2 2 2 2 2 2 2 2 10 6 2 6 5 2 3 4 1 3 9
## L 6 4 4 3 2 6 4 3 5 15 34 4 20 13 5 4 6 6 7 13
## K 6 18 10 8 2 10 8 5 8 5 4 24 9 2 6 8 8 4 3 5
## M 1 1 1 1 0 1 1 1 1 2 3 2 7 2 1 1 1 1 1 2
## F 2 1 2 1 1 1 1 1 3 5 6 1 4 32 1 2 2 4 20 3
## P 7 5 5 4 3 5 4 5 5 3 3 4 3 2 19 6 5 1 2 4
## S 9 6 8 7 7 6 7 9 6 5 4 7 5 3 9 10 9 4 4 6
## T 8 5 6 6 4 5 5 6 4 6 4 6 5 3 6 8 11 2 3 6
## W 0 2 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 55 1 0
## Y 1 1 2 1 3 1 1 1 3 2 2 1 2 15 1 2 2 3 31 2
## V 7 4 4 4 4 4 4 5 4 15 10 4 10 5 5 5 7 2 4 17
On pourra chercher les différences avec les valeurs de la figure 83 :
A quel taux de mutations cette matrice correspond ?
## [1] 0.8038141
Est-ce que la somme des colonnes vaut bien 1 ?
## A R N D C Q E G H I L K M F P S T W Y V
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Il est maintenant possible de calculer la PAM250.
2.1.6 Calcul de PAM250
Pour un alignement de l’acide aminé \(AA_i\) avec l’acide aminé \(AA_j\), 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_{ji}}{f_i.f_j} = \frac{f_j.M_{ij}}{f_i.f_j}\)
avec
- \(M\) la matrice
mutation.mat
précédemment calculée (mutation matrix) et - \(f_i\) la fréquence de chacun des
acides aminés
aa.freq
.
On en déduit la vraisemblance d’un alignement entre \(AA_i\) et \(AA_j\) : \(L_{i,j} = \frac{ \cancel{f_j}.M_{i,j}}{f_i.\cancel{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} = \log_{10}\frac{M_{i,j}}{f_i}\)
Il n’y a plus qu’à appliquer la formule \(S_{i,j} = \log_{10}\frac{M_{i,j}}{f_i}\) avec \(f_i\) la fréquence d’apparition des a.a. et \(M_{i,j}\) la probabilité que l’a.a. \(j\) soit remplacé par l’a.a. \(i\).
Vous devriez obtenir la matrice de scores de substitutions suivante :
## A R N D C Q E G H I L K M F P S T W Y V
## A 2 -2 0 0 -2 0 0 1 -1 -1 -2 -1 -1 -3 1 1 1 -6 -3 0
## R -2 6 0 -1 -4 1 -1 -3 2 -2 -3 3 0 -4 0 0 -1 2 -4 -3
## N 0 0 2 2 -4 1 1 0 2 -2 -3 1 -2 -3 0 1 0 -4 -2 -2
## D 0 -1 2 4 -5 2 3 1 1 -2 -4 0 -3 -6 -1 0 0 -7 -4 -2
## C -2 -4 -4 -5 12 -5 -5 -3 -3 -2 -6 -5 -5 -4 -3 0 -2 -8 0 -2
## Q 0 1 1 2 -5 4 3 -1 3 -2 -2 1 -1 -5 0 0 -1 -5 -4 -2
## E 0 -1 1 3 -5 2 4 0 1 -2 -3 0 -2 -5 -1 0 0 -7 -4 -2
## G 1 -3 0 1 -3 -1 0 5 -2 -3 -4 -2 -3 -5 -1 1 0 -7 -5 -1
## H -1 2 2 1 -3 3 1 -2 6 -2 -2 0 -2 -2 0 -1 -1 -3 0 -2
## I -1 -2 -2 -2 -2 -2 -2 -3 -2 5 2 -2 2 1 -2 -1 0 -5 -1 4
## L -2 -3 -3 -4 -6 -2 -3 -4 -2 2 6 -3 4 2 -3 -3 -2 -2 -1 2
## K -1 3 1 0 -5 1 0 -2 0 -2 -3 5 0 -5 -1 0 0 -3 -4 -2
## M -1 -1 -2 -3 -5 -1 -2 -3 -2 2 4 0 6 0 -2 -2 -1 -4 -3 2
## F -4 -4 -4 -6 -4 -5 -5 -5 -2 1 2 -5 0 9 -5 -3 -3 0 7 -1
## P 1 0 -1 -1 -3 0 -1 -1 0 -2 -3 -1 -2 -5 6 1 0 -6 -5 -1
## S 1 0 1 0 0 -1 0 1 -1 -1 -3 0 -2 -3 1 2 1 -2 -3 -1
## T 1 -1 0 0 -2 -1 0 0 -1 0 -2 0 -1 -3 0 1 3 -5 -3 0
## W -6 2 -4 -7 -8 -5 -7 -7 -2 -5 -2 -3 -4 1 -5 -2 -5 17 0 -6
## Y -3 -4 -2 -4 0 -4 -4 -5 0 -1 -1 -4 -2 7 -5 -3 -3 0 10 -2
## V 0 -3 -2 -2 -2 -2 -2 -1 -2 4 2 -2 2 -1 -1 -1 0 -6 -2 4
Pour pouvoir mieux la comparer à la figure 84, on réordonne les lignes et colonnes dans le même ordre que la figure :
## [1] "C" "S" "T" "P" "A" "G" "N" "D" "E" "Q" "H" "R" "K" "M" "I" "L" "V" "F" "Y" "W"
## C S T P A G N D E Q H R K M I L V F Y W
## C 12 0 -2 -3 -2 -3 -4 -5 -5 -5 -3 -4 -5 -5 -2 -6 -2 -4 0 -8
## S 0 2 1 1 1 1 1 0 0 -1 -1 0 0 -2 -1 -3 -1 -3 -3 -2
## T -2 1 3 0 1 0 0 0 0 -1 -1 -1 0 -1 0 -2 0 -3 -3 -5
## P -3 1 0 6 1 -1 -1 -1 -1 0 0 0 -1 -2 -2 -3 -1 -5 -5 -6
## A -2 1 1 1 2 1 0 0 0 0 -1 -2 -1 -1 -1 -2 0 -3 -3 -6
## G -3 1 0 -1 1 5 0 1 0 -1 -2 -3 -2 -3 -3 -4 -1 -5 -5 -7
## N -4 1 0 0 0 0 2 2 1 1 2 0 1 -2 -2 -3 -2 -3 -2 -4
## D -5 0 0 -1 0 1 2 4 3 2 1 -1 0 -3 -2 -4 -2 -6 -4 -7
## E -5 0 0 -1 0 0 1 3 4 2 1 -1 0 -2 -2 -3 -2 -5 -4 -7
## Q -5 0 -1 0 0 -1 1 2 3 4 3 1 1 -1 -2 -2 -2 -5 -4 -5
## H -3 -1 -1 0 -1 -2 2 1 1 3 6 2 0 -2 -2 -2 -2 -2 0 -3
## R -4 0 -1 0 -2 -3 0 -1 -1 1 2 6 3 0 -2 -3 -3 -4 -4 2
## K -5 0 0 -1 -1 -2 1 0 0 1 0 3 5 0 -2 -3 -2 -5 -4 -3
## M -5 -2 -1 -2 -1 -3 -2 -3 -2 -1 -2 -1 0 6 2 4 2 0 -3 -4
## I -2 -1 0 -2 -1 -3 -2 -2 -2 -2 -2 -2 -2 2 5 2 4 1 -1 -5
## L -6 -3 -2 -3 -2 -4 -3 -4 -3 -2 -2 -3 -3 4 2 6 2 2 -1 -2
## V -2 -1 0 -1 0 -1 -2 -2 -2 -2 -2 -3 -2 2 4 2 4 -1 -2 -6
## F -4 -3 -3 -5 -4 -5 -4 -6 -5 -5 -2 -4 -5 0 1 2 -1 9 7 0
## Y 0 -3 -3 -5 -3 -5 -2 -4 -4 -4 0 -4 -4 -2 -1 -1 -2 7 10 0
## W -8 -2 -5 -5 -6 -7 -4 -7 -7 -5 -2 2 -3 -4 -5 -2 -6 1 0 17
Le résultat est proche de la figure 84.
2.1.7 Pour aller plus loin
Il est possible de faire la même chose en python. Ce serait
l’occasion de s’entraîner, notamment avec le module numpy
ou encore l’environnement jupyter
:
Notebook jupyter : PAM.jupyter.ipynb et sa sortie HTML : PAM.jupyter.html
3 Calcul du taux de mutations à partir d’une matrice PAM
Nous allons partir de la matrice PAM250 la plus récente similaire à la matrice de substitution que nous avons calculée précédemment.
Chargez le fichier PAM250.csv avec la
fonction read_csv
.
## Rows: 20 Columns: 20
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (20): A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## A R N D C Q E G H I L K M F P S T W Y V
## A 2.7 -0.9 -0.4 -0.4 1.0 -0.2 0.0 0.4 -1.3 -1.1 -1.3 -0.5 -1.0 -2.5 0.6 1.2 0.5 -3.4 -2.6 -0.2
## R -0.9 5.0 0.4 -0.4 -1.4 1.6 0.5 -1.5 0.7 -2.8 -2.3 3.0 -1.5 -2.9 -1.4 -0.1 -0.4 -2.4 -2.0 -2.5
## N -0.4 0.4 3.8 2.2 -1.0 0.7 1.0 0.6 1.3 -3.1 -3.0 0.9 -2.3 -2.9 -1.1 0.8 0.4 -3.1 -1.4 -2.6
## D -0.4 -0.4 2.2 4.8 -2.4 1.0 2.9 -0.4 0.4 -4.5 -4.1 0.6 -3.3 -4.4 -0.8 0.5 -0.3 -4.3 -2.9 -3.6
## C 1.0 -1.4 -1.0 -2.4 9.8 -1.7 -2.1 -1.8 -1.5 0.1 -0.7 -1.8 -0.5 -0.2 -2.2 0.7 0.8 0.7 -0.2 1.0
## Q -0.2 1.6 0.7 1.0 -1.7 3.0 1.7 -1.5 0.9 -2.2 -1.7 1.7 -0.8 -2.5 -0.7 0.3 0.0 -2.6 -1.6 -1.8
## E 0.0 0.5 1.0 2.9 -2.1 1.7 3.4 -1.0 0.1 -3.4 -3.2 1.2 -2.4 -4.1 -0.6 0.4 -0.1 -3.9 -2.9 -2.5
## G 0.4 -1.5 0.6 -0.4 -1.8 -1.5 -1.0 7.0 -1.9 -5.4 -5.1 -1.1 -4.2 -5.4 -2.2 0.2 -1.7 -4.3 -4.8 -4.2
## H -1.3 0.7 1.3 0.4 -1.5 0.9 0.1 -1.9 7.1 -2.7 -1.8 0.6 -1.3 1.1 -0.6 -0.5 -0.7 1.4 3.0 -2.4
## I -1.1 -2.8 -3.1 -4.5 0.1 -2.2 -3.4 -5.4 -2.7 4.1 2.8 -2.7 2.3 0.8 -3.4 -2.2 -0.5 -1.4 -1.0 3.4
## L -1.3 -2.3 -3.0 -4.1 -0.7 -1.7 -3.2 -5.1 -1.8 2.8 4.0 -2.3 3.0 2.0 -2.7 -2.1 -1.1 -0.7 0.2 2.0
## K -0.5 3.0 0.9 0.6 -1.8 1.7 1.2 -1.1 0.6 -2.7 -2.3 3.2 -1.5 -3.3 -0.8 0.2 -0.1 -3.5 -2.3 -2.3
## M -1.0 -1.5 -2.3 -3.3 -0.5 -0.8 -2.4 -4.2 -1.3 2.3 3.0 -1.5 4.4 1.6 -2.7 -1.5 -0.6 -0.8 0.1 1.6
## F -2.5 -2.9 -2.9 -4.4 -0.2 -2.5 -4.1 -5.4 1.1 0.8 2.0 -3.3 1.6 7.1 -4.0 -2.7 -2.1 4.8 6.0 0.0
## P 0.6 -1.4 -1.1 -0.8 -2.2 -0.7 -0.6 -2.2 -0.6 -3.4 -2.7 -0.8 -2.7 -4.0 8.6 0.5 -0.8 -5.1 -3.7 -2.4
## S 1.2 -0.1 0.8 0.5 0.7 0.3 0.4 0.2 -0.5 -2.2 -2.1 0.2 -1.5 -2.7 0.5 2.3 1.5 -2.4 -2.3 -1.4
## T 0.5 -0.4 0.4 -0.3 0.8 0.0 -0.1 -1.7 -0.7 -0.5 -1.1 -0.1 -0.6 -2.1 -0.8 1.5 3.3 -3.2 -2.1 0.1
## W -3.4 -2.4 -3.1 -4.3 0.7 -2.6 -3.9 -4.3 1.4 -1.4 -0.7 -3.5 -0.8 4.8 -5.1 -2.4 -3.2 14.8 5.4 -1.8
## Y -2.6 -2.0 -1.4 -2.9 -0.2 -1.6 -2.9 -4.8 3.0 -1.0 0.2 -2.3 0.1 6.0 -3.7 -2.3 -2.1 5.4 8.4 -1.5
## V -0.2 -2.5 -2.6 -3.6 1.0 -1.8 -2.5 -4.2 -2.4 3.4 2.0 -2.3 1.6 0.0 -2.4 -1.4 0.1 -1.8 -1.5 3.5
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{M_{i,j}}{\cancel{f_i}} \times \frac{\cancel{f_i}}{M_{i,k}}))\)
les \(f_i\) 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 \(M_{i,1} = 1\)
(première colonne de mutation.mat
à 1), on peut en déduire
les autres valeurs de la matrice pour la ligne \(i\) de mutation.mat
. 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}}\)
Utilisez la formule précédente pour obtenir la matrice de probabilités.
On obtient donc la matrice suivante :
## A R N D C Q E G H I L K M F P S T W Y V
## A 1 0.44 0.49 0.49 0.68 0.51 0.54 0.59 0.40 0.42 0.40 0.48 0.43 0.30 0.62 0.71 0.60 0.25 0.30 0.51
## R 1 3.89 1.35 1.12 0.89 1.78 1.38 0.87 1.45 0.65 0.72 2.45 0.87 0.63 0.89 1.20 1.12 0.71 0.78 0.69
## N 1 1.20 2.63 1.82 0.87 1.29 1.38 1.26 1.48 0.54 0.55 1.35 0.65 0.56 0.85 1.32 1.20 0.54 0.79 0.60
## D 1 1.00 1.82 3.31 0.63 1.38 2.14 1.00 1.20 0.39 0.43 1.26 0.51 0.40 0.91 1.23 1.02 0.41 0.56 0.48
## C 1 0.58 0.63 0.46 7.59 0.54 0.49 0.52 0.56 0.81 0.68 0.52 0.71 0.76 0.48 0.93 0.95 0.93 0.76 1.00
## Q 1 1.51 1.23 1.32 0.71 2.09 1.55 0.74 1.29 0.63 0.71 1.55 0.87 0.59 0.89 1.12 1.05 0.58 0.72 0.69
## E 1 1.12 1.26 1.95 0.62 1.48 2.19 0.79 1.02 0.46 0.48 1.32 0.58 0.39 0.87 1.10 0.98 0.41 0.51 0.56
## G 1 0.65 1.05 0.83 0.60 0.65 0.72 4.57 0.59 0.26 0.28 0.71 0.35 0.26 0.55 0.95 0.62 0.34 0.30 0.35
## H 1 1.58 1.82 1.48 0.95 1.66 1.38 0.87 6.92 0.72 0.89 1.55 1.00 1.74 1.17 1.20 1.15 1.86 2.69 0.78
## I 1 0.68 0.63 0.46 1.32 0.78 0.59 0.37 0.69 3.31 2.45 0.69 2.19 1.55 0.59 0.78 1.15 0.93 1.02 2.82
## L 1 0.79 0.68 0.52 1.15 0.91 0.65 0.42 0.89 2.57 3.39 0.79 2.69 2.14 0.72 0.83 1.05 1.15 1.41 2.14
## K 1 2.24 1.38 1.29 0.74 1.66 1.48 0.87 1.29 0.60 0.66 2.34 0.79 0.52 0.93 1.17 1.10 0.50 0.66 0.66
## M 1 0.89 0.74 0.59 1.12 1.05 0.72 0.48 0.93 2.14 2.51 0.89 3.47 1.82 0.68 0.89 1.10 1.05 1.29 1.82
## F 1 0.91 0.91 0.65 1.70 1.00 0.69 0.51 2.29 2.14 2.82 0.83 2.57 9.12 0.71 0.95 1.10 5.37 7.08 1.78
## P 1 0.63 0.68 0.72 0.52 0.74 0.76 0.52 0.76 0.40 0.47 0.72 0.47 0.35 6.31 0.98 0.72 0.27 0.37 0.50
## S 1 0.74 0.91 0.85 0.89 0.81 0.83 0.79 0.68 0.46 0.47 0.79 0.54 0.41 0.85 1.29 1.07 0.44 0.45 0.55
## T 1 0.81 0.98 0.83 1.07 0.89 0.87 0.60 0.76 0.79 0.69 0.87 0.78 0.55 0.74 1.26 1.91 0.43 0.55 0.91
## W 1 1.26 1.07 0.81 2.57 1.20 0.89 0.81 3.02 1.58 1.86 0.98 1.82 6.61 0.68 1.26 1.05 66.07 7.59 1.45
## Y 1 1.15 1.32 0.93 1.74 1.26 0.93 0.60 3.63 1.45 1.91 1.07 1.86 7.24 0.78 1.07 1.12 6.31 12.59 1.29
## V 1 0.59 0.58 0.46 1.32 0.69 0.59 0.40 0.60 2.29 1.66 0.62 1.51 1.05 0.60 0.76 1.07 0.69 0.74 2.34
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 = M^T\)), cela correspond à un système linéaire de 20 équations à 20 inconnues : \(Ax=b\) avec \(b_i = 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=\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^{-1} A x = A^{-1} b\)
or \(A^{-1}A = I\) d’où \(A^{-1} A x = Ix = x = A^{-1}b\)
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 mutation.mat
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 \(x_1\), la 2ème colonne de \(A\) par \(x_2\), …
Et comme \(A=M^T\), cela correspond
à multiplier la 1ère ligne de mutation.mat
par
\(x_1\), la 2ème ligne de
mutation.mat
par \(x_2\),
…
Après avoir vérifié que \(A\) est
inversible, déterminez la matrice mutation.mat
.
mutation.mat
suppose première colonne à 1, maintenant il
faut que chaque colonne ait comme somme 1. on transpose et on a un
système linéaire de 20 équations à 20 inconnues : Système de Cramer
A est elle inversible ?
## [1] 6128934
déterminant non nul (6128934) donc là, c’est ok
on peut faire \(A^{-1} A x = A^{-1} b\) d’où \(I x = A^{-1} b\) d’où \(x = A^{-1} b\)
## [1] 0.1545 0.0426 0.0384 0.0498 0.0159 0.0330 0.0727 0.0858 0.0151 0.0517 0.0649 0.0555 0.0207 0.0210 0.0470 0.0701 0.0652 0.0040 0.0155 0.0767
on vérifie :
## A R N D C Q E G H I L K M F P S T W Y V
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
on est bon on peut mettre à jour mutation.mat
Vous devriez obtenir :
## A R N D C Q E G H I L K M F P S T W Y V
## A 0.15 0.07 0.08 0.08 0.10 0.08 0.08 0.09 0.06 0.06 0.06 0.07 0.07 0.05 0.10 0.11 0.09 0.04 0.05 0.08
## R 0.04 0.17 0.06 0.05 0.04 0.08 0.06 0.04 0.06 0.03 0.03 0.10 0.04 0.03 0.04 0.05 0.05 0.03 0.03 0.03
## N 0.04 0.05 0.10 0.07 0.03 0.05 0.05 0.05 0.06 0.02 0.02 0.05 0.02 0.02 0.03 0.05 0.05 0.02 0.03 0.02
## D 0.05 0.05 0.09 0.16 0.03 0.07 0.11 0.05 0.06 0.02 0.02 0.06 0.03 0.02 0.05 0.06 0.05 0.02 0.03 0.02
## C 0.02 0.01 0.01 0.01 0.12 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.01 0.01 0.02
## Q 0.03 0.05 0.04 0.04 0.02 0.07 0.05 0.02 0.04 0.02 0.02 0.05 0.03 0.02 0.03 0.04 0.03 0.02 0.02 0.02
## E 0.07 0.08 0.09 0.14 0.04 0.11 0.16 0.06 0.07 0.03 0.03 0.10 0.04 0.03 0.06 0.08 0.07 0.03 0.04 0.04
## G 0.09 0.06 0.09 0.07 0.05 0.06 0.06 0.39 0.05 0.02 0.02 0.06 0.03 0.02 0.05 0.08 0.05 0.03 0.03 0.03
## H 0.02 0.02 0.03 0.02 0.01 0.03 0.02 0.01 0.10 0.01 0.01 0.02 0.02 0.03 0.02 0.02 0.02 0.03 0.04 0.01
## I 0.05 0.03 0.03 0.02 0.07 0.04 0.03 0.02 0.04 0.17 0.13 0.04 0.11 0.08 0.03 0.04 0.06 0.05 0.05 0.15
## L 0.06 0.05 0.04 0.03 0.07 0.06 0.04 0.03 0.06 0.17 0.22 0.05 0.17 0.14 0.05 0.05 0.07 0.07 0.09 0.14
## K 0.06 0.12 0.08 0.07 0.04 0.09 0.08 0.05 0.07 0.03 0.04 0.13 0.04 0.03 0.05 0.07 0.06 0.03 0.04 0.04
## M 0.02 0.02 0.02 0.01 0.02 0.02 0.01 0.01 0.02 0.04 0.05 0.02 0.07 0.04 0.01 0.02 0.02 0.02 0.03 0.04
## F 0.02 0.02 0.02 0.01 0.04 0.02 0.01 0.01 0.05 0.04 0.06 0.02 0.05 0.19 0.01 0.02 0.02 0.11 0.15 0.04
## P 0.05 0.03 0.03 0.03 0.02 0.03 0.04 0.02 0.04 0.02 0.02 0.03 0.02 0.02 0.30 0.05 0.03 0.01 0.02 0.02
## S 0.07 0.05 0.06 0.06 0.06 0.06 0.06 0.06 0.05 0.03 0.03 0.06 0.04 0.03 0.06 0.09 0.08 0.03 0.03 0.04
## T 0.07 0.05 0.06 0.05 0.07 0.06 0.06 0.04 0.05 0.05 0.05 0.06 0.05 0.04 0.05 0.08 0.12 0.03 0.04 0.06
## W 0.00 0.01 0.00 0.00 0.01 0.00 0.00 0.00 0.01 0.01 0.01 0.00 0.01 0.03 0.00 0.01 0.00 0.26 0.03 0.01
## Y 0.02 0.02 0.02 0.01 0.03 0.02 0.01 0.01 0.06 0.02 0.03 0.02 0.03 0.11 0.01 0.02 0.02 0.10 0.19 0.02
## V 0.08 0.05 0.04 0.04 0.10 0.05 0.05 0.03 0.05 0.18 0.13 0.05 0.12 0.08 0.05 0.06 0.08 0.05 0.06 0.18
Vérification que les sommes des colonnes valent 1 :
## A R N D C Q E G H I L K M F P S T W Y V
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
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, f_j = f_1 = 1\), ce qui nous permet de calculer les autres fréquences : \(f_i = \frac{M_{i,j}}{M_{j,i}}\), exemple \(f_2 = \frac{M_{2,1}}{M_{1,2}}\)
calcul des fréquences. on suppose \(f_1\) (\(f_A\)) = 1
On obtient :
## A R N D C Q E G H I L K M F P S T W Y V
## 1.00 0.63 0.51 0.66 0.15 0.42 0.88 0.94 0.25 0.80 1.06 0.75 0.31 0.45 0.49 0.64 0.70 0.11 0.34 0.97
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.
on s’arrange pour que la somme des fréquences soit de 1
## [1] 1
Vous devriez obtenir les fréquences suivantes :
## A R N D C Q E G H I L K M F P S T W Y V
## 8.30 5.24 4.21 5.46 1.26 3.46 7.27 7.83 2.04 6.67 8.76 6.23 2.60 3.73 4.09 5.32 5.81 0.87 2.82 8.04
Avec la fonction mutation.rate
définie précédemment
:
## [1] 0.822772
4 Liens, références et sources utiles à la préparation de ce TP
Merci à http://www.biorecipes.com/Dayhoff/code.html
Fichiers
- https://github.com/telliott99/Dayhoff/blob/master/README.md
- http://www.deduveinstitute.be/~opperd/private/matrices.html
autres