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/x3x5 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 MUR11480 2.04.40.13.9-0.7
BIGLYCANH46600 2.22.60.22.5-0.8
BIGLYCANH46600 2.33.80.02.6-0.9
IMMUNOGLOBULIN CONSTANT HEAVY CHAIN REGIONN92646 2.44.3-0.33.4-0.5
ANTIGEN IDENTIFIED BY MONOCLONAL ANTIBODY KI-67AA004872 -1.5-0.5-0.4-2.3-1.3
IMMUNOGLOBULIN LAMBDA LIGHT CHAINT63045 2.73.80.02.7-0.8
LEUKOCYTE COMMON ANTIGEN (CD45)H71488 1.82.70.01.9-1.0
(2'-5') OLIGOADENYLATE SYNTHETASE EAA100116 0.2-0.40.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.530.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$TrTr 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.