Calcul matriciel, hommage à Margaret Dayhoff

1 Opérations R sur les matrices

1.1 Création à partir d’un vecteur

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

1.2 Création à partir d’un data.frame

class(iris)
## [1] "data.frame"
iris
B = as.matrix( iris[,-5] )
class(B)
## [1] "matrix" "array"
head( B )
##      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\) :

I3 = diag(3)
I3
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Accès à la diagonale :

A
##      [,1] [,2] [,3]
## [1,]   11   12   13
## [2,]   21   22   23
## [3,]   31   32   33
diag(A)
## [1] 11 22 33

Modification de la diagonale :

diag(A) = c(11,44,66)
A
##      [,1] [,2] [,3]
## [1,]   11   12   13
## [2,]   21   44   23
## [3,]   31   32   66

Création à partir de la diagonale :

diag(6:9)
##      [,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

V = c(0.1, 10, 1000)
class(V)
## [1] "numeric"

Obtention de la matrice vecteur colonne :

as.matrix(V)
##       [,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 :

t( as.matrix(V) )
##      [,1] [,2] [,3]
## [1,]  0.1   10 1000

1.5 Opérations sur les matrices

  • Addition +, soustraction -, multiplication *, …
  • 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)

Comparaison de la multiplication et du produit matriciel :

A * V
##         [,1]    [,2]    [,3]
## [1,]     1.1     1.2     1.3
## [2,]   210.0   440.0   230.0
## [3,] 31000.0 32000.0 66000.0
A %*% V
##         [,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.

Arbre de parenté avec séquences ancêtres Mutations acceptées

2.1.1 Matrice de comptage count.mat

Cette matrice de comptage correspond à la figure 80 :

Fig. 80
Fig. 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.

library(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
count.mat = read_csv('dayhoff.fig80.counts.csv')
## 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 :

count.mat

Il va falloir la rendre symétrique en recopiant les valeurs.

Avec la fonction replace_na, remplacez les NA par des 0.

count.mat = count.mat %>% replace(is.na(.), 0)

Vous devriez avoir :

count.mat

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

count.mat = count.mat %>% 
  select(-aa) %>%
  as.matrix
rownames(count.mat) = colnames(count.mat)
count.mat = count.mat+t(count.mat)

Vous devriez avoir :

count.mat
##     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 :

sum(count.mat) / 10 / 2
## [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 :

Fig. 81
Fig. 81

Et, elle nous est donnée pour chaque aa dans la table 21 :

Tab. 21
Tab. 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.

mut.rel = read_csv('dayhoff.mutabilities.csv' )
## 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.
foo = mut.rel$mut
names(foo) = mut.rel$aa
aa.mut = foo[ colnames(count.mat)]

Vous devriez obtenir :

aa.mut
##   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 :

Tab. 22
Tab. 22

Chargez le fichier dayhoff.frequencies.csv correspondant et créez le vecteur aa.freq.

foo = read_csv('dayhoff.frequencies.csv')
## 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.
aa.freq = foo$freq
names(aa.freq) = foo$aa

Vous devriez obtenir :

aa.freq
##     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
aa.freq =  aa.freq * 1/sum(aa.freq)

Mettre à l’échelle les fréquences pour que leur somme soit 1.

sum(aa.freq)
## [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 :

Fig. 82
Fig. 82

Ainsi que les formules pour son calcul :

formules
formules

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

lambda = 0.01 / sum(aa.freq * aa.mut)

Vous devriez obtenir :

lambda
## [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\)

#mm.1 = lambda * aa.mut * count.mat / sum(count.mat)
mm.1 = lambda * 1000 * count.mat / sum(count.mat)
diag(mm.1) = 1 - lambda * aa.mut
mutation.mat = mm.1

Vous devriez obtenir :

round( mutation.mat * 10000)
##      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 :

colA = mutation.mat[,'A']
round(colA * 10000)
##    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 :

(1 - colA['A']) # à atteindre
##          A 
## 0.01330303
sum(colA) - colA['A'] # somme actuelle
##          A 
## 0.01542356
(1 - colA['A']) / (sum(colA) - colA['A']) # facteur d'échelle
##         A 
## 0.8625137

Calculez les facteur de mise à l’échelle pour l’ensemble des colonnes.

scales = (1-diag(mm.1)) / (colSums(mm.1)-diag(mm.1))

Vous devriez obtenir :

scales %>% round(2)
##    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.

mm.1 = t(mm.1 * scales)
diag(mm.1) = 1 - lambda * aa.mut
mutation.mat = mm.1

Vous devriez obtenir :

round(mutation.mat * 10000)
##      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.

mutation.rate = function(m, f) sum(f * (1 - diag(m)))

Utilisation :

mutation.rate(mutation.mat, aa.freq)
## [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 :

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
E = eigen(A)
E
## 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
Ev = E$vectors

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 :

P = Ev

ainsi \(P\) diagonalise \(A\) :

A
##      [,1] [,2] [,3]
## [1,]    1    2    0
## [2,]    0    3    0
## [3,]    2   -4    2
solve(P) %*% A %*% P
##      [,1] [,2] [,3]
## [1,]    3    0    0
## [2,]    0    2    0
## [3,]    0    0    1

la diagonale obtenue contient les valeurs propres.

ou encore :

D = diag( eigen(A)$values )
A
##      [,1] [,2] [,3]
## [1,]    1    2    0
## [2,]    0    3    0
## [3,]    2   -4    2
P %*% D %*% solve(P)
##      [,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 :

P = eigen(mutation.mat)$vectors
D = diag(eigen(mutation.mat)$values)
mutation.mat.250 = P %*% D^250 %*% solve(P)
colnames(mutation.mat.250) = colnames(mm.1)
rownames(mutation.mat.250) = colnames(mm.1)

Vous devriez obtenir :

round(mutation.mat.250 * 100)
##    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 :

Probabilité de mutation pour 250-PAM
Probabilité de mutation pour 250-PAM

A quel taux de mutations cette matrice correspond ?

mutation.rate(mutation.mat.250, aa.freq)
## [1] 0.8038141

Est-ce que la somme des colonnes vaut bien 1 ?

colSums(mutation.mat.250)
## 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\).

PAM250= log10 ( mutation.mat.250 / aa.freq )

Vous devriez obtenir la matrice de scores de substitutions suivante :

round(PAM250 * 10)
##    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 :

PAM250
PAM250
pam.order = strsplit('CSTPAGNDEQHRKMILVFYW','')[[1]]
pam.order
##  [1] "C" "S" "T" "P" "A" "G" "N" "D" "E" "Q" "H" "R" "K" "M" "I" "L" "V" "F" "Y" "W"
round(PAM250[pam.order, pam.order]*10)
##    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.

PAM250 = read_csv('PAM250.csv') %>% as.matrix
## 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.
rownames(PAM250) = colnames(PAM250)
PAM250
##      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.

mutation.mat = 10**((PAM250 - PAM250[,1])/10)

On obtient donc la matrice suivante :

round(   mutation.mat , 2)
##   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=t(mutation.mat)
b=rep(1,20)

A est elle inversible ?

det(A)
## [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\)

x=as.vector( solve(A) %*% b )
round( x ,4)
##  [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 :

apply(mutation.mat * x, MARGIN=2, sum)
## 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

mutation.mat=mutation.mat*x

Vous devriez obtenir :

round(mutation.mat,2)
##      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 :

colSums(mutation.mat)
## 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

PAM250.freq=mutation.mat[,1]/mutation.mat[1,]
names(PAM250.freq)=names(aa.freq)

On obtient :

round( PAM250.freq , 2)
##    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

PAM250.freq=PAM250.freq/sum(PAM250.freq)
sum(PAM250.freq)
## [1] 1
#round( aa.freq*100,2 )
#round( PAM250.freq*100,2 )

Vous devriez obtenir les fréquences suivantes :

round( PAM250.freq *100 , 2)
##    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 :

mutation.rate(mutation.mat, PAM250.freq)
## [1] 0.822772