UNIVERSITE DE LA MEDITERRANEE
MASTER INFORMATIQUE 1ère année
Faculté des Sciences de Luminy
Année 2004--2005
Département d'Informatique
M02: Complexité
 
Michel VAN CANEGHEM ( Michel.Van.Caneghem@lidil.univ-mrs.fr)
MAISON TURING COMPLEXITE BIO I BIO II ARO VIDEO


Nouvelles
- Le devoir 1 est à rendre le 29 Novembre 2004

- Le devoir 2 est à rendre le 10 Décembre 2004


Université de
la Méditerranée
Faculté
des Sciences
de Luminy
Département
d'Informatique
Laboratoire
d'Informatique
Fondamentale
Ce document a été traduit de LATEX par HEVEA.

Le devoir 2 : Alignement de grandes séquences biologiques

Table des matières

Attention, ce devoir est assez difficile : commencez assez tôt.

1  Présentation

1.1  Motivations 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 Ansari­Lari et al. (PDF)
ATTENTION : ces articles présentent une méthode qui n'est pas celle que nous allons utiliser

1.2  Description de l'objectif du devoir

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.

Une version plus élaborée de ce que je vous demande de programmer, se trouve dans un ensemble de programmes pour la biologie : EMBOSS. Il s'agit du programme stretcher. Voici l'article correspondant à ce programme : Eugene W. Myers, Webb Miller: Optimal alignments in linear space. Computer Applications in the Biosciences 4(1): 11-17 (1988) . La version PDF de cet article.

2  Ce qu'il y a à faire

Le vrai sujet du devoir
  1. Programmer l'algorithme classique d'alignement global de séquence avec coût de gap constant à l'aide de la programmation dynamique.
  2. La même chose avec l'algorithme de Hirschberg.
  3. 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.
  4. Comparez sur des exemples tirés au hasard, les performances en mémoire et en temps de ces deux méthodes. Le but est de vérifier que le temps de calcul n'est pas plus du double que dans la méthode classique et que la mémoire utilisée est effectivement linéaire.
  5. Traitez à titre de vérification le premier exemple (avec la méthode de Hirschberg) -- vous pouvez construire le dotplot
  6. 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 : quelle est la plus longue sous-chaîne qui s'aligne sans aucun gap ?
Votre document doit comporter les points suivant
  1. Rédigez 2 ou 3 pages sur l'algorithme de Hirschberg, ainsi que sur la démonstration des complexité en mémoire et en temps.
  2. Environ 1 ou 2 pages de résultats de test
  3. 1 tableau de comparaison des 2 méthodes
  4. Les résultats de la première question (je ne veux pas des pages et des pages d'alignement)
  5. Les résultats de la deuxième question (Même remarque)
  6. 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 classique

Je vous donne les Pages 52-53 du livre Introduction to Computational Molecular Biology (voir référence à la fin) qui décrivent le programme d'alignement. Voici mon programme en Java align.java qui fait la même chose (pour aider ceux qui n'y arrivent pas), mais ce programme est donné sans aucune garantie.

3.2  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. Vous pouvez aussi consulter les transparents : Linear Alignment

3.3  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.4  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.fasta, pl6Mouse.fasta.

3.5  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.fasta, U47924.fasta.

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 à embl. Si possible ouvrez cette page dans un autre fenêtre. La suivez le mode d'emploi suivant :
  1. Vous regardez en haut à droite, et vous devez trouver quelque chose du genre ;
    Get nucleotide sequence for xxxxx GO
  2. Vous pouvez entrer directement votre accession number. Si vous voulez des requètes plus complètes cliquez sur "Database query", juste en dessous de la ligne précédente. Vous vous retrouvez sur une page, repérez le tableau "Nucleotide sequence". Et la dans la colonne "Query Engine(s)" vous devez cliquer sur "SRS".
  3. Vous pouvez alors rentrer votre requête dans le champ "AllText" : "Xenopus laevis rhodopsin gene complete" (ne me demandez pas ce que c'est !!);
  4. Au bout d'un certain temps, vous devez voir apparaitre une liste de références. Cochez la dernière référence et modifier dans "Display option", pour obtenir "* complete entry * puis, cliquez sur . 
    EMBL:XL23808  U23808  U23808 Xenopus laevis rhodopsin gene, complete cds.  4734
  5. 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.
  6. Vous pouvez sauver juste la séquence au format FASTA de la manière suivante : Dans "Display option", vous choisissez : "FastaSeqs". Vous pouvez cocher également "Printer friendly view". Vous avez alors un fichier que vous pouvez sauvegarder au format texte. Appelons le "xl23808.seq".
  7. 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




Lundi 22 novembre 2004
©2004 Michel Van Caneghem