Devoir 5 - solutions
1. Parcimonie pondérée (15 points)
Donnez un algorithme qui minimise la valeur de parcimonie pondérée sur un arbre. Soit S une matrice de pondération de substitutions sur un alphabet A de taille m. (Le probléme de base considéré dans le cours utilise S(i,j)=0 si i=j et =1 sinon.) Soit T une phylogénie avec n feuilles où chaque noeud u est étiqueté par un caractère cu de A. La valeur de parcimonie sur cet arbre est égale au somme des valeurs S(cu, cv) sur les branches uv. L'algorithme doit trouver les meilleurs caractères cu ainsi que la valeur de parcimonie minimale sur l'arbre. Démontrer que l'algorithme est correct. Analysez son temps et espace de calcul en fonction de n et m.
Solution. Pour arriver à une solution
par programmation dynamique, on
a besoin d'une bonne recurrence
. Définissons
W(u,c) pour chaque noeud u
et caractère c comme
la valeur de parcimonie dans le sous-arbre raciné à u
quand u est étiqueté par u.
Les W(.,.) pour les noeuds internes peuvent être calculées par
récurrence:
W(u,c')=
min {W(v1,c)+S(c',c)}
+ min {W(v2,c)+S(c',c)}
+ ...
+ min {W(vk,c)+S(c',c)},
où v1, ..., vk
sont les enfants de u et
chaque min prend le minimum sur les choix de caractère pour c.
Pour les feuilles, on a W(u,c)
= 0 si u est étiqueté par c
et = infini sinon.
Pour démontrer que la récurrence est correcte, notons
que (1) le côté droit représente un étiquetage
possible [donc on a <= au moins] et
que (2) l'étiquetage optimal qui correspond à
W(u,c')
doit minimiser les
W(vi,c)+S(c',c)
[sinon, on pourrait en trouver un meilleur en
choisissant un meilleur c à l'enfant].
Il est aussi utile de définir
Cu(v,c')
comme le caractére pour lequel min atteint sa valeur
dans le calcul de W(u,c')
où u est le parent de v.
L'algorithme calcule les W ainsi que les C
à partir des feuilles
en montant vers la racine. La valeur de parcimonie pour l'arbre
avec racine u0
est min W(u0, c) où le min
est pris sur les choix possibles de c.
Le meilleur étiquetage est retrouvé
à l'aide des C
en descendant vers les feuilles.
Le temps de
calculer W(u,c) et Cu(.,c)
pour u et c fixés
est de O(mk) où k
est le nombre des enfants de u.
Calculer tous les W et C
prend O(nm2). L'espace de calcul (stockage
des W et C) est
de O(nm).
2. Analyse de données d'expression de gènes (20 points)
Dans cet exercice vous travaillez avec des données d'une
expérince classique
. Il s'agit d'une étude sur le cancer du sein,
décrite dans l'article de
Perou et al. Distinctive gene expression patterns in human mammary
epithelial cells and breast cancers
, Proceedings of the National Acadamy of Sciences
of the USA, vol. 96, pages 9212-9217, 1999. Vous pouvez
télécharger l'article
à partir du
site du journal.
Tout d'abord, lisez l'article. Les données sont accessibles au
site du
Stanford Breast Cancer Microarray Project.
a.
Les données pour chaque échantillon sont stockées dans un
fichier texte. Dans cet exercice vous travaillez avec
les fichiers pour les échantillons de
BC4 (copie locale),
BC14 (copie locale),
BC21 (copie locale),
-EGF (copie locale),
et +Matrigel (copie locale).
Chaque fichier et de taille 1.3M.
Le tableau de chaque fichier contient les données pour
une expérience: chaque rangée correspond à une sonde (spot
).
Familiarisez-vous avec le format
du tableau, en lisant la
déscription des colonnes.
On veut travailler avec la valeur x5/x3
où x5 est la valeur du signal Cy5 après
la sous-traction du bruit (background
), et
x3 et la même pour Cy3.
Quelle colonne donne cette valeur?
Soit v(s,p)
la valeur observée pour sonde s sur puce p
(p=BC4,BC14,BC21,-EGF,+Matrigel).
Compilez un tableau qui donne le logarithme (base 2) de v(s,p)
pour des sondes (s) avec nombres d'accession
R11480,
H46600,
N92646,
AA004872,
T63045,
H71488,
AA100116
(attention: il est possible qu'il y a plusieurs sondes
avec le même accession). Exemple (la colonne id
sert seulement à distinguer entre des sondes
avec accessions identiques) :
id accno BC4 BC14 BC21 -EGF +Matrigel 1 AA004872 2.4 5.7 6.0 0.15 12.1 2 T63045 1.1 2.6 10.2 18.2 1.15 ...
b. Trouvez les entrées de Genbank pour ces nombres d'accession: ce sont des ESTs (Expressed Sequence Tags), donnez les génes qui y correspondent (PUTATIVE ID). Trouvez les clusters auxquels ils appartiennent selon l'article (Figure 2, clusters A,B,C et D).
c.
Finalement, calculez un clustering hiérarchique agglomeratif sur les colonnes du tableau par
la méthode de Eisen (
Eisen et al. Cluster analysis and display of genome-wide expression patterns
,
PNAS vol. 95, pages 14863-14868, 1998).
Dans cette approche, la similarité de deux vecteurs (clusters ou feuilles)
est calculée par la correlation de Pearson. Donnez un tableau
de similarités entre les colonnes (tableau de taille 5*5).
Montrez chaque étape dans le calcul du clustering.
Est-ce que votre arbre correspond à celui de l'article?
Vous pouvez faire tous les calculs à l'aide d'un logiciel comme Excel. Pour extraire les rangées pour des sondes avec accession AA99999, vous pouvez taper awk -F '[\t]+' '$4 == "AA99999" {print ;}' < fichier.txt sous Unix.
Solution.
Les valeurs avec lesquelles on va travailler se trouvent dans la colonne
RAT2
. Voici le tableau des logarithmes des RAT2.
Nom (fichiers originaux) | Accno | BC4 | BC21 | -EGF | BC14 | +Matrigel |
---|---|---|---|---|---|---|
IMMUNOGLOBULIN MU | R11480 | 2.0 | 4.4 | 0.1 | 3.9 | -0.7 |
BIGLYCAN | H46600 | 2.2 | 2.6 | 0.2 | 2.5 | -0.8 |
BIGLYCAN | H46600 | 2.3 | 3.8 | 0.0 | 2.6 | -0.9 |
IMMUNOGLOBULIN CONSTANT HEAVY CHAIN REGION | N92646 | 2.4 | 4.3 | -0.3 | 3.4 | -0.5 |
ANTIGEN IDENTIFIED BY MONOCLONAL ANTIBODY KI-67 | AA004872 | -1.5 | -0.5 | -0.4 | -2.3 | -1.3 |
IMMUNOGLOBULIN LAMBDA LIGHT CHAIN | T63045 | 2.7 | 3.8 | 0.0 | 2.7 | -0.8 |
LEUKOCYTE COMMON ANTIGEN (CD45) | H71488 | 1.8 | 2.7 | 0.0 | 1.9 | -1.0 |
(2'-5') OLIGOADENYLATE SYNTHETASE E | AA100116 | 0.2 | -0.4 | 0.8 | -0.2 | -0.4 |
Les gènes qui correspondenx à ces ESTs: AA004872 - B48666 (Ki-67, cluster A), AA100116 - D00068 ((2'-5')oligoadenylate synthetase, cluster B), R11480 - X17115 (immunoglobin [Ig], cluster C), T63045 - X57809 (immnunoglobin, cluster C), H71488 (leukocyte antigen CD45, cluster C), N92646 (immunoglobin, cluster C), H46600 - J04599 (biglycan, cluster D). Les entrées de GenBank pour H71488 et N92646 ne correspondent pas aux noms donnés par les fichiers - sans doute ce sont des erreurs dans les fichiers de l'expérience originale.
Tableau initial de correlations.
Échantillon | BC4 | BC21 | -EGF | BC14 |
---|---|---|---|---|
BC21 | 0.79 | |||
-EGF | 0.00 | -0.26 | ||
BC14 | 0.829 | 0.830 | -0.02 | |
+Matrigel | 0.40 | 0.26 | 0.53 | 0.45 |
Cluster 1: (BC14, BC21). La vecteur qui y correspond est la moyenne des vecteurs (de logarithmes) pour BC14 et BC21. Taux de correlation entre cluster 1 et les autres: BC4 0.85, -EGF -0.15, +Matrigel 0.32. Cluster 2: (BC4, cluster 1). La vecteur qui y correspond est la moyenne des vecteurs pour BC4, BC14 et BC21. Taux de correlation: -EGF -0.11, +Matrigel 0.35. Cluster 3: (+Matrigel, -EGF). Arbre:
+-- BC14 +--+ | +-- BC21 +---+ +---+ +----- BC4 | | +--------- +Matrigel +---+ +--------- -EGF
Cet arbre ne correspond pas à celui de l'article:
là les échantillons BC4 et BC21 forment un cluster.
La raison pour laquelle on a un arbre différent est qu'on
a construit l'arbre sur un très petit ensemble de gènes
(l'article utilisait 1247 gènes dont l'expression
changeait beacoup
parmi les échantillons).
Alors que ce sont des gènes d'importance particulière
pour l'étude, ils déterminons des
relations différentes entre les échantillons.
3. Chevauchement exact (15 points)
Soit S et T deux séquences de longueurs n et m. Donnez un algortihme qui calcule le plus grand chevauchement exact entre les deux en temps O(n+m). On a un chevauchement exact de longueur t si S[1 ... t]=T[m-t+1 ... m] ou si T[1 ... t]=S[n-t+1 ... n]. L'algorithme doit trouver la valeur maximale de t, ainsi que spécifier si le chevauchement et entre le préfixe de S et le suffixe de T ou vice versa. (Les techniques de recherche exacte peuvent vous aider à trouver un algo.)
Solution.
Pour trouver le maximum t
avec
S[1 ... t]=T[m-t+1 ... m],
on utilise l'algorithme Z
sur le mot X=S$Tr
où Tr est T renversé.
Maintenant une boîte Z
de X après le $
correspond au cas où un préfixe de Tr
(donc un suffixe de T)
est identique à un préfixe de S.
Le meilleur t est le maximum des Zi(X)
avec i=n+1,...,n+m+1
calculés par l'algorithme.
On peut faire la même chose pour
trouver le plus long alignement entre un préfixe de T et
un suffixe de S en utilisant Z
sur le mot X'=T$Sr=Xr.