silico.biotoul.fr
 

M1 MABS BBS Math TD Calcul Matriciel

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Application aux matrices de substitution pour les alignements de séquences)
m (Application aux matrices de substitution pour les alignements de séquences)
Line 43: Line 43:
= Application aux matrices de substitution pour les alignements de séquences =
= Application aux matrices de substitution pour les alignements de séquences =
-
== Calcul d'une matrice de substitution à partir des alignements ==
+
== 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 à  
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 à  
Line 130: Line 130:
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é... à vos claviers, vous devriez obtenir quelque chose comme 44,5%.
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é... à vos claviers, vous devriez obtenir quelque chose comme 44,5%.
 +
 +
 +
== 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 [[File:PAM250.txt|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 <math>S_{i,j} = 10 \times log_{10}\frac{M_{i,j}}{F_i}</math>.
 +
 +
D'où la différence entre 2 valeurs d'une même ligne de la matrice (PAM250) :
 +
 +
<math>S_{i,j} - S_{i,k} = 10 \times log\frac{M_{i,j}}{F_i} - 10 \times log\frac{M_{i,k}}{F_i}</math>
 +
 +
<math>S_{i,j} - S_{i,k} = 10 \times (log\frac{M_{i,j}}{F_i} - log\frac{M_{i,k}}{F_i})</math>
 +
 +
avec <math>log a - log b = log \frac{a}{b}</math>
 +
 +
<math>S_{i,j} - S_{i,k} = 10 \times (log (\frac{M_{i,j}}{F_i} \times \frac{F_i}{M_{i,k}}))</math>
 +
 +
les <math>F_i</math> se simplifient et on obtient :
 +
 +
<math>S_{i,j} - S_{i,k} = 10 \times (log \frac{M_{i,j}}{M_{i,k}})</math>
 +
 +
d'où
 +
 +
<math>\frac{S_{i,j} - S_{i,k}}{10} = log \frac{M_{i,j}}{M_{i,k}}</math>
 +
 +
d'où
 +
 +
<math>10^{\frac{S_{i,j} - S_{i,k}}{10}} = \frac{M_{i,j}}{M_{i,k}}</math>
 +
 +
 +
En supposant que <math>M_{i,1} = 1</math> (première colonne de <tt>MM</tt> à 1), on peut en déduire les autres valeurs de la matrice pour la ligne <math>i</math> de <tt>MM</tt>. En utilisant la formule ci-dessus avec <math>k=1</math> :
 +
 +
<math>\frac{M_{i,j}}{M_{i,1}} = \frac{M_{i,j}}{1} = 10^{\frac{S_{i,j} - S_{i,1}}{10}}</math>
 +
 +
 +
On obtient donc la matrice suivante :
 +
 +
<source lang="rsplus">
 +
round(MM[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
 +
</source>
 +
 +
 +
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 (<math>A = M^T</math>), cela correspond à un système linéaire de 20 équations à 20 inconnues : <math>A.X=B</math> avec <math>B_i =  1</math> pour <math>i \in [1, 20]</math>
 +
 +
 +
 +
Cela correspond donc à un système de Cramer de 20 équations à 20 inconnues :

Revision as of 12:07, 29 September 2011

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 File:Gold.metadata.txt :

G=read.table("gold.metadata.txt", sep="\t", header=TRUE)
class(G)

G est un data frame ; les deux premières colonnes contiennent l'identifiant et le nom de l'organisme. Pour extraire la matrice de données (colonnes numériques), on fait :

as.matrix( G[ , 3:12] )

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(1, 0, 5)
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 : entre 2 matrices A et I3 : A %*% I3, entre une matrice A et un vecteur V A %*% V
  • transposition : t(M)


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

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 les coûts attribuer à chaque substitution (match et mismatch).


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.

CM=read.table("count_matrix.txt", sep="\t", header=TRUE, row.names=1)


Après avoir chargé cette matrice dans R, transformez-la en une matrice MM (mutation matrix) 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.

Vous devriez obtenir ceci (10 premières lignes et 10 premières colonne) :

round(MM[1:10,1:10],2)
     A    R    N    D    C    Q    E    G    H    I
A 0.54 0.03 0.03 0.03 0.08 0.04 0.05 0.05 0.02 0.02
R 0.02 0.57 0.03 0.01 0.02 0.06 0.03 0.01 0.04 0.01
N 0.01 0.02 0.48 0.06 0.01 0.03 0.03 0.03 0.05 0.01
D 0.02 0.01 0.09 0.57 0.01 0.04 0.11 0.02 0.03 0.00
C 0.01 0.00 0.00 0.00 0.57 0.00 0.00 0.00 0.00 0.00
Q 0.02 0.04 0.03 0.02 0.01 0.42 0.05 0.01 0.04 0.01
E 0.04 0.04 0.05 0.13 0.01 0.11 0.52 0.02 0.03 0.01
G 0.05 0.02 0.05 0.03 0.02 0.02 0.02 0.77 0.02 0.00
H 0.00 0.01 0.02 0.01 0.01 0.02 0.01 0.00 0.54 0.00
I 0.02 0.01 0.01 0.00 0.03 0.01 0.01 0.00 0.01 0.50


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 M la matrice MM précédemment calculée (mutation matrix) et F la fréquence de chacun des acides aminés (ici, les fréquences d'apparition dans nos alignements de départ).


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 CM de départ qui contient le nombre de substitutions observées pour calculer le vecteur Freq contenant la fréquence de chacun de a.a.

Vous devriez obtenir le vecteur suivant :

round(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 substitution suivante :

round(S[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


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é... à vos claviers, vous devriez obtenir quelque chose comme 44,5%.


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 File:PAM250.txt 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{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 MM à 1), on peut en déduire les autres valeurs de la matrice pour la ligne i de MM. 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(MM[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 : A.X = B avec Bi = 1 pour i \in [1, 20]


Cela correspond donc à un système de Cramer de 20 équations à 20 inconnues :