Vendredi 30 janvier 2004
| Nouvelles |
- Les cours 5 et 7 en vidéo sont disponibles
- la date de l'examen est fixée au 28 Janvier
- Le cours 8 est supprimé
- mes meilleurs voeux de bonheur de de réussite pour cette nouvelle année 2004
(cliquez-ici).
|
|
Le devoir 3 : Alignement de grandes séquences biologiques
Table des matières
Ce devoir est à rendre avant le Lundi 2 Février 2004 avant 10h, début du cours sur les Systèmes distribués.
(Je vous propose un delai jusqu'au Mercredi 4 Février (avant 23h59) -- cependant si l'enseignant concerné se plaint de l'absence des élèves en
Cours/TD/TP, alors je diviserai la note par 2 pour ceux qui ont rendu le devoir3 après le 2 Février. Sinon c'est bon. A vous de choisir.)
1 Présentation
Description de l'objectif du devoir
Il s'agit de programmer l'alignement de deux séquences et de l'appliquer à quelques exemples issus de la Biologie.
Cette année je me suis intéressé à l'alignement de grandes séquences biologiques. On dispose maintenant des génomes
complets de nombreuses espèces. Cela permet de deviner la fonction de certains gènes de nouvelles espèces en les comparant
à des espèces mieux connues. On est donc amené à comparer les séquences de chromosomes complets (ou de grandes parties).
Regarder par exemple l'article suivant :
Alignment of Whole Genomes. A.L. Delcher, S. Kasif, R.D. Fleischmann, J. Peterson, O. White, and S.L. Salzberg.
Nucleic Acids Research, 27:11 (1999), 2369-2376.
(l'article pour ceux qui n'arriveraient pas à le télécharger l'article
est ici)
L'objectif de ce devoir est de faire un programme qui permet de comparer deux séquences d'ADN (celles qui sont décrites à
la fin de l'article précédent : "human chromosome 12p13 and its syntenic region in mouse chromosome 6"). Plus précisement :
we chose a 222 930 bp subsequence of human
chromosome 12p13 (accession no. U47924) that is syntenic to a
227 538 bp contiguous subsequence of mouse chromosome 6
(accession no. AC002397). These sequences were the subject of
a recent study by AnsariLari et al. (PDF)
Cela ne se fait pas simplement : l'algorithme classique construit autour de la programmation dynamique ne peut pas
marcher dans de tels cas. En effet il utilise un capacité mémoire de l'ordre de O(n2) si n est la taille des chaînes à
comparer. Dans le cas qui nous intéresse, n vaut environ 230000. Il faut donc une mémoire de 230000*230000*4 =
211600 Moctets = 211 Go. (c'est assez rare de trouver un ordinateur avec une telle capacité). Il faut donc changer
d'algorithme. J'ai donc choisit d'utiliser un algorithme basé sur la méthode de Hirschberg qui utilise un espace mémoire en
O(n) pour un temps de calcul un peu plus lent que l'algorithme classique (ATTENTION, ce n'est pas la méthode proposée dans l'article précédent).
Remarque : Pour une utilisation traditionnelle des algorithmes d'alignement je vous engage à consulter le
devoir de l'année dernière disponibles dans
les archives 2003-2003.
2 Ce qu'il y a à faire
Le vrai sujet du devoir
-
Programmer l'algorithme classique d'alignement global de séquence avec coût de gap constant à l'aide de la programmation dynamique.
- La même chose avec l'algorithme de Hirschberg.
- Testez vos deux méthodes avec de petits exemples. Utilisez la valeur de gap et la matrice de similarité donnée dans le paragraphe suivant.
- Comparez sur des exemples tirés au hasard, les performances en mémoire et en temps de ces deux méthodes.
- Traitez à titre de vérification le premier exemple (avec la méthode de Hirschberg) -- vous pouvez construire le dotplot
- Passez le deuxième exemple -- attention cela est très long (1h20 chez moi) -- vous pouvez construire le dotplot avec "dotter", mais
c'est assez lent. Question subsidiaire : quel est la plus longue sous-chaîne qui s'aligne sans aucun gap ?
Votre document doit comporter les points suivant
-
Rédigez 2 ou 3 pages sur l'algorithme ainsi que sur la démonstration des complexité en mémoire et en temps.
- Environ 1 ou 2 pages de résultats de test
- 1 tableau de comparaison des 2 méthodes
- Les résultats de la première question (je ne veux pas des pages et des pages d'alignement)
- Les résultats de la deuxième question (Même remarque)
- Conclusions et commentaires (est-ce que cela vous a intéressé?). Pouvez-vous proposer une méthode pour améliorer l'efficacité de l'algorithme.
Remarque1 : L'année dernière, un certain nombre d'élèves ont eu des problèmes pour reconstruire l'alignement (il s'agit
d'un backtracking simple). Commencez ce devoir le plus tôt possible et n'hésitez pas à me poser des questions
Remarque2 : Je veux un document d'au plus une dizaine de pages avec des phrases en français. N'oubliez pas que c'est à
vous de me prouver par des exemples que vôtre programme est juste. Je ne souhaite le listing que pour vérifier que
vous n'avez pas copié!!
3 Le matériel
3.1 La description de l'algorithme de Hirschberg
Il y a une très bonne description de cet algorithme dans le livre de Meidanis and Setubal [Pages 58-61] (voir la référence à la fin de la page).
3.2 Le cout du gap et la matrice de similarité
J'ai choisit comme coût fixe du gap la valeur -8. La matrice est la suivante:
Matrice de similarité
|
A C G T N
A 5 -4 -4 -4 -2
C -4 5 -4 -4 -2
G -4 -4 5 -4 -2
T -4 -4 -4 5 -2
N -2 -2 -2 -2 -1
| |
Remarque : Attention : "N" n'est pas une nouvelle base de l'ADN, mais c'est ce qui est mis dans les séquences d'ADN quand la base n'est pas connue.
3.3 Le premier exemple
Il s'agit du même que celui de l'année dernière.
Pour vous familiariser avec les outils d'alignements, vous allez commencer par travailler avec les séquences
nucléiques des gènes pl6 chez l'homme et la souris. Allez chercher les séquences sur EMBL : [AC=U09584](Human PL6 protein (PL6) mRNA, complete cds) et
[AC=AF134238] (Mus musculus PL6 protein (Pl6) mRNA, complete cds) et sauvegardez les dans deux fichiers.
Pour ceux qui n'arrivent pas a trouver ces séquences, les voici :
pl6Human.adn,
pl6Mouse.adn.
3.4 Le deuxième exemple
Les deux séquences à comparer sont : AC002397 et
U47924. Comme d'habitude voici les séquences que j'ai récupéré :
AC002397.adn,
U47924.adn.
4 Mes résultats
4.1 Comparaison classique/hirschberg
Voici ce que j'ai trouvé en comparant des séquences aléatoires d'ADN de longueur 5000. Encore une fois je vous donne ces résultats à titre exemple, vous
pouvez trouver des choses différentes.
Comparaison classique/hirschberg
|
Normal : Nb Tests = 10, Lg chaînes = 3000
Temps total = 7385ms, Temps moyen = 738500µs, Temps moyen/n^2 = 82ps
Hirshberg : Nb Tests = 10, Lg chaînes = 5000
Temps total = 20832ms, Temps moyen = 2083200µs, Temps moyen/n^2 = 83ps
| |
Ces résultats sont bien meilleurs que ce que nous avons démontré en TD, mais c'est comme cela... (cela est peut-être du à la grande consommation mémoire de
la première méthode). [si vous utilisez Java et que vos pile débordent, utilisez l'option : java -Xmx250M ...]
Remarque : 1ps = 10-9s
4.2 Le premier exemple
Voici ce que j'ai trouvé (même remarque que précédemment, avec en plus le fait que si le score doit être identique, les alignements et le nombre de gaps,
peuvent être différents)
premier exemple
|
>embl|U09584|HS095841 Human PL6 protein (PL6) mRNA, complete cds.
longueur = 1860
Nombre de A = 340
Nombre de C = 605
Nombre de G = 544
Nombre de T = 371
Autres = 0
>embl|AF134238|AF134238 Mus musculus PL6 protein (Pl6) mRNA, complete cds.
longueur = 2177
Nombre de A = 423
Nombre de C = 667
Nombre de G = 639
Nombre de T = 447
Autres = 1
score = 3962
Longueur de l'alignement = 2208
Nombre d'identités = 1590 (72%)
Nombre de gaps = 379 (17%)
60
1: ggcga---gg------gg-c-c--ta-c---gc-tgc-g-gc--ccgg----caa-c---
: ::: :: :: : : :: : :: ::: : : :::: ::: :
2: gtcgactaggtcccaaggactccgtatcccagcatgccgagaagccggaaggcaagcgct
1: -a-agg-c---c--c-g---------------------a--ctcggc-ccc-t------c
: ::: : : : : : : :::: ::: : :
2: cagagggcgtactgccgcggtcgccggngggggcgcgcaggcgcggcgcccctgtttgtc
.................................
| |
4.3 Le deuxième exemple
ATTENTION : malgré tout le temps de calcul est très long : Temps total = 4756532ms = 1h20...Il faut être patient. C'était assez prévisible, en effet
le temps de calcul est en O(mn). Ce qui donne dans ce cas (en utilisant les mesures du début): 228*223/1000*83 s = 4220s (prévision pas trop mauvaise!)
Voici mes résultats.
deuxième
|
>AC002397 Sequence 227538 BP; 55803 A; 55982 C; 56282 G; 58635 T; 836 other;
longueur = 227538
Nombre de A = 55803
Nombre de C = 55982
Nombre de G = 56282
Nombre de T = 58635
Autres = 836
>U47924 Sequence 222930 BP; 53286 A; 56778 C; 57221 G; 55643 T; 2 other;
longueur = 222930
Nombre de A = 53286
Nombre de C = 56778
Nombre de G = 57221
Nombre de T = 55643
Autres = 2
score = 86124
Longueur de l'alignement = 245381
Nombre d'identités = 136350 (55%)
Nombre de gaps = 40294 (16%)
60
1: atattttccttcctcaacctttgttataaaaaaacactcttgaggcct-agcatggtggt
: : : :: :: : : : : : : :::::: :: :: :: :
2: ---tgtagc--cc-ca-gc-t-act-tg---g--gaggc-tgaggcaggagaat--tgct
1: gcacacctttgagtccagtcctttggaggccggggcaggctgatctccgtgagtttgagg
: ::: ::: ::: ::: : : : : :: :::: :: : : :
2: tga-acccgggag-gcag-aggttgca-g-tgagcca-a--gatca-cgccac-t-g--c
...........................
| |
5 Comment chercher une séquence d'ADN sur Internet?
Ceci est extrait de la page de l'année dernière. Prévenez-moi si il y a des liens faux (à revoir).
Si vous voulez quelque chose de plus juste et plus à jour, consultez
une page web de Pascal Hingamp.
Je vous conseille de vous connecter à Infobiogen. Si possible ouvrez cette page dans un autre
fenêtre. La suivez le mode d'emploi suivant :
-
Clicquez sur Start
- Cochez la base de donnée d'EMBL (European Molecular Biology Laboratory).
- Appuyez sur le bouton "Standard" de "Query Forms" (partie gauche de l'écran).
- Entrez alors dans la requête : "Xenopus laevis rhodopsin gene complete" (ne me demandez pas ce que c'est !!);
puis appuyez sur le bouton "Submit Query" et attendez.
- vous devez voir apparaitre deux références :
EMBL:XL23808 U23808 Xenopus laevis rhodopsin gene, complete cds. 4734
EMBL:XLRHODOP L07770 Xenopus laevis rhodopsin mRNA, complete cds. 1684
- Cliquez sur la première référence. Vous avez plein d'information : il est important de noter l'accession number [AC=U23808],
c'est une
référence unique, qui permet d'identifier de manière unique la séquence. Vers le bas vous allez trouver la séquence d'ADN correspondante.
- Vous pouvez sauver juste la séquence au format FASTA de la manière suivante : cliquez sur le bouton "Save", puis choisisssez dans
"use view" : FastaSeqs, puis recliquez sur le bouton "Save". Vous avez alors un fichier que vous pouvez sauvegarder au format texte.
Appelons le "xl23808.seq".
- IMPORTANT Pour retrouver une séquence par son AC (accession number). Dans la liste déroulante ou il y a écrit "AllText", choisissez "AccNumber". C'est
beaucoup plus sûr (selon les spécialistes)
 |
|
Juste pour votre Information : Xenopus laevis
est une grenouille africaine qui est utilisé
par la recherche car la femelle pond beaucoup
d'oeufs et peut-être pour d'autres raisons que j'ignore!!! |
6 Comment construire un "dotplot" ?
Comme on l'a vu dans le cours, c'est un moyen très simple de construire un dessin qui permet de visualiser l'alignement de deux
séquences. Hélas, cela ne donne que des informations qualitatives, cependant c'est très utile pour savoir ou l'on va. Il y a de nombreux
outils sur Internet qui permettent de visualiser des dotplots, en voici deux qui fonctionnent bien :
Voici un résultat obtenu pour le deuxième, c'est très lent!!!. On voit bien qu'il y a ressemblance : une diagonale semble se dessiner
7 Quelques liens fort utiles
Il y a vraiment beaucoup d'informations sur Internet, voici une sélection qui me parait particulièrement
intéressante :
Je pense que vous en avez assez!!! Il y a aussi le livre que j'ai utilisé pour la nouvelle version du cours :
Introduction to Computational Molecular Biology, Meidanis and Setubal, PWS Publishing Company
|
