These (Classification ARN Codants & ARN Non Codants)
These (Classification ARN Codants & ARN Non Codants)
d’ARN non-codants
THÈSE
pour l’obtention du
par
Arnaud Fontaine
Composition du jury
Examinateurs : Fabrice Leclerc, C.R. CNRS MAEM, Université Henri Poincaré – Nancy 1
Nouredine Melab, Professeur LIFL, Université des Sciences et Technologies de Lille
Fariza Tahi, Maı̂tre de conférences IBISC, Université d’Evry-Val d’Essonne
Directeur : Hélène Touzet, C.R. CNRS LIFL, Université des Sciences et Technologies de Lille
Introduction 1
i
Table des matières
ii
4 Deux exemples d’intégration de Protea et caRNAc 113
4.1 L’alignement multiple de séquences nucléiques . . . . . . . . . . . . . . . . . . 113
4.1.1 L’alignement multiple de séquences codantes homologues . . . . . . . 114
4.1.2 L’alignement multiple de séquences partageant une structure commune 114
4.1.3 Magnolia, alignement de séquences fonctionnelles homologues . . . . 115
4.1.4 Les résultats expérimentaux de Magnolia . . . . . . . . . . . . . . . 120
4.2 L’annotation par génomique comparative . . . . . . . . . . . . . . . . . . . . . 121
4.2.1 Le pipeline d’annotation . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4.2.2 Résultats expérimentaux du pipeline . . . . . . . . . . . . . . . . . . . 133
Conclusion 137
Bibliographie 141
iii
Table des matières
iv
Introduction
Tous les organismes vivants, des plus simples aux plus complexes, sont composés de cellules
qui présentent des caractéristiques communes en terme de structure mais également de fonc-
tionnement. Trois types de macromolécules fondamentales sont impliquées dans cette unité
cellulaire et moléculaire du vivant : les ADN, les ARN et les protéines. Schématiquement,
la séquence des éléments qui composent ces molécules constitue la représentation minimale
permettant de décrire l’information qu’elles contiennent.
La séquence des ADN est responsable du stockage de l’information génétique qui détermine
le patrimoine génétique d’un organisme. L’information génétique est segmentée en gènes dont
l’expression est à l’origine de la synthèse d’ARN puis de protéines. La séquence d’une protéine
ne constitue qu’un premier niveau dans sa description. Les protéines se replient en effet dans
l’espace pour former une structure tridimensionnelle qui détermine leur fonction. Les ARN
sont quant à eux des acteurs plus polyvalents. Les ARN codants contiennent l’information
nécessaire à la synthèse d’une protéine, tandis que les ARN non-codants se comportent som-
mairement comme des protéines en se repliant sur eux-mêmes pour adopter une conformation
spatiale qui détermine leur fonction.
Alors que l’on dispose de plus en plus de séquences d’ADN et d’ARN provenant de la
génomique et de la transcriptomique, la signification de la plupart de ces séquences reste
encore à élucider. Les premiers travaux d’annotation automatique de séquences remontent au
début des années 80, suite au premier séquençage complet d’un génome. A l’heure actuelle,
près de 1 000 génomes d’organismes différents sont complètement séquencés et disponibles
publiquement, et plus de 4 000 projets de séquençage sont en cours 1 . De prime abord,
l’analyse systématique de ces séquences par des techniques expérimentales en vue de leur
annotation n’est plus envisagée à cause de leurs coûts humain et financier trop importants.
L’analyse automatique de séquences par des moyens informatiques est donc plus que jamais
un challenge majeur.
L’aboutissement de nombreux projets de séquençage ces dernières années a toutefois
quelque peu changé la donne en contribuant à l’émergence de la génomique comparative.
En effet, bien que portée par les mêmes supports et exprimée selon des mécanismes com-
muns, l’information génétique est également la source de la diversité des organismes vivants.
L’unicité cellulaire du vivant suggère ainsi l’évolution à partir d’ancêtres communs plus ou
moins éloignés durant laquelle les séquences génomiques se transforment, tout en préservant
certaines fonctions.
La génomique comparative consiste à étudier et analyser les ressemblances et les différences
apparues durant l’évolution entre des séquences génomiques. Ces séquences ne sont ainsi plus
considérées de manière individuelle mais reliées entre elles par l’évolution. Dans ce contexte,
1
http://www.genomesonline.org
1
Introduction
Plan de lecture
Ce document est organisé en quatre chapitres.
Le premier chapitre est consacré à l’introduction des notions biologiques nécessaires à la
bonne compréhension des méthodes présentées dans ce document. Ce chapitre débute par la
présentation des mécanismes en lien direct avec le stockage et l’expression de l’information
génétique communs à tous les organismes vivants. Ensuite, nous détaillons la transcription,
première étape de l’expression de l’information génétique, qui permet de synthétiser les ARN,
et la maturation des transcrits. Puis, nous nous intéressons plus en détails aux deux types
d’ARN existants : les ARN codants et les ARN non-codants. Enfin, la dernière partie de ce
chapitre est consacrée à la mise en place de notre fil conducteur : l’analyse comparative de
séquences avec notamment la définition de méta-séquence qui revient de manière récurrente
dans les chapitres suivants.
Le second chapitre porte sur la détection de séquences codantes. Tout d’abord, nous
présentons trois types de méthodes existantes : l’approche ab initio, l’approche par similarité
de séquence et enfin l’approche comparative. La littérature foisonne de méthodes dédiées à
ce problème. Nous n’en présentons qu’une sélection éclairée représentative de la diversité des
méthodes existantes. Ensuite, nous présentons Protea, notre contribution au problème. Par-
tant d’observations concrètes sur deux séquences, nous formalisons dans un premier temps
la détection d’une séquence d’acides aminés conservée sur deux séquences par la comparai-
son de leurs traductions potentielles. Puis, ce principe est étendu à des ensembles de taille
quelconque de séquences. Enfin, la dernière partie de ce chapitre est consacrée aux résultats
expérimentaux de Protea. Les performances de Protea sont évaluées sur un large éventail
de séquences, puis nous passons à un cas pratique en appliquant Protea à l’annotation de
2
nouvelles séquences codantes sur le génome humain.
Le troisième chapitre porte sur la prédiction de structures d’ARN et la prédiction de gènes
à ARN. Nous commençons par poser le problème de la prédiction de structures en introduisant
les deux grandes familles d’approches : l’approche thermodynamique et l’approche compara-
tive. Nous refermons ce rapide état de l’art sur la prédiction de structures par l’évaluation
de référence des méthodes existantes, BRAliBase I proposée par Gardner en 2004 [GG04].
Ensuite, nous nous intéressons à un problème étroitement lié à la prédiction de structures
d’ARN, la prédiction de gènes à ARN. Pour présenter ce problème, nous dressons un parallèle
avec les approches mises en œuvre pour la prédiction de gènes codants en partant des biais
de composition en séquence pour terminer par l’approche comparative. La troisième partie
qui compose ce chapitre décrit notre contribution au problème de la prédiction de structures.
Nous y présentons ainsi les évolutions apportées au logiciel caRNAc, dédié à la prédiction de
structures secondaires conservées, notamment l’introduction des méta-séquences. Enfin, les
résultats de caRNAc par rapport à ceux des méthodes existantes sur le jeu de données de
référence, BRAliBase I, sont exposés dans la dernière partie de ce chapitre. Avant de refer-
mer ce chapitre, nous décrivons les travaux que nous avons menés dans l’optique de définir
une méthode de prédiction de gènes à ARN basée sur l’existence de structure conservée signi-
ficative prédite par caRNAc.
Le quatrième et dernier chapitre de ce document est dédié à la présentation de deux tra-
vaux collaboratifs menés au sein de l’équipe qui intègrent Protea et caRNAc. La première
partie de ce chapitre est consacrée à l’alignement multiple de séquences codantes homologues
et de séquences qui partagent une structure commune, avec le logiciel Magnolia. Magnolia
est un logiciel d’alignement multiple qui résulte de la combinaison de Protea, caRNAc et
Gardenia, un logiciel d’alignement de structures d’ARN développé dans l’équipe. La seconde
partie de ce chapitre est dédiée à l’application de Protea et caRNAc pour l’annotation de
séquences génomiques. Nous présentons le pipeline logiciel développé dans l’équipe avant de
fournir quelques résultats expérimentaux.
3
Introduction
4
Chapitre 1
Les travaux décrits dans ce manuscrit portent sur l’analyse des acides ribonucléiques
(ARN) codants et non-codants par des approches de bio-informatique. Dans ce premier cha-
pitre, nous présentons le contexte biologique général et les mécanismes moléculaires sur les-
quels s’appuient nos travaux.
Nous commençons par situer les ARN, leur synthèse et leurs fonctions au sein de la cellule,
en section 1.1. Nous nous intéressons ensuite plus précisément aux caractéristiques des ARN
codants dans la section 1.2, des ARN non-codants dans la section 1.3 et à leur évolution
dans la section 1.4. Etant informaticien, toute cette présentation n’est pas rédigée par un
spécialiste, et s’adresse en priorité à des non-spécialistes. Le but est de fournir les notions de
base en biologie moléculaire nécessaires à la compréhension du document et de légitimer les
choix de modèles que nous avons faits dans la suite du travail.
Enfin, dans la dernière section du chapitre, nous abordons les premiers formalismes et
méthodes de bio-informatique avec la génomique comparative, l’alignement de séquences et
l’introduction du concept de méta-séquence. Les méta-séquences et les ensembles de méta-
séquences nous serviront tout au long de ce document pour représenter des ensembles de
séquences d’ARN de distance évolutive hétérogène.
5
Chapitre 1. Les Acides RiboNucléiques
1. Nucléole
2. Noyau
3. Ribosome
4. Vésicule
5. Réticulum endoplasmique rugueux
6. Appareil de Golgi
7. Microtubule
8. Réticulum endoplasmique lisse
9. Mitochondrie
10. Vacuole
11. Cytoplasme
12. Lysosome
13. Centrosome
Source http://en.wikipedia.org/wiki/File:Biological_cell.svg
Parmi les eucaryotes, on retrouve aussi bien des organismes unicellulaires, tels que les
levures, que des organismes pluricellulaires tels que les plantes et les animaux. La majorité des
procaryotes sont quant à eux des organismes unicellulaires microscopiques. Les procaryotes
sont divisés en deux règnes : les bactéries et les archées. Bien que la taille et la forme des
archées sont similaires à celles des bactéries, les archées s’en distinguent par des caractères
plus similaires à ceux des eucaryotes tels que la structure des gènes et les mécanismes relatifs
à leur expression.
Toute cellule est régie essentiellement par cinq types de macromolécules : les lipides, les glu-
cides, les acides désoxyribonucléiques (ADN), les acides ribonucléiques (ARN) et les protéines.
L’ADN assure le stockage de l’information génétique et sa transmission au fil des générations.
L’ensemble de l’ADN qui définit un organisme est appelé son génome. Chez les eucaryotes,
l’ADN est stocké dans le noyau. Les lipides sont les principaux constituants des membranes
cellulaires. Les glucides sont le soutien de la vie : ils servent de source et de stockage d’énergie,
ils participent aux paroies cellulaires, . . . Les protéines sont des molécules indispensables à la
structuration et au fonctionnement des cellules, et résultent de l’expression de l’information
génétique. Les acides ribonucléiques (ARN) sont quant à eux des protagonistes plus ambigus
qui peuvent endosser plusieurs rôles auxquels nous nous intéressons par la suite.
6
1.1. L’ARN au sein de la cellule
Fig. 1.2 – Le dogme central présentant la transcription de l’ADN en ARN messagers traduits
par la suite en protéine, et la transcription de l’ADN en ARN non-codants.
Le dogme central mentionne également deux autres types d’ARN : les ARN ribosomiques
et les ARN de transfert. Contrairement aux ARN messagers, ces ARN sont des molécules
fonctionnelles non traduites en protéine et que l’on regroupe sous le terme d’ARN non-codants.
Au moment de leurs découvertes, ces deux types d’ARN non-codants apparaissent comme
des exceptions au dogme central. Depuis, de nombreux autres ARN non-codants ont été
découverts, portant à plus de 600 le nombre de familles d’ARN non-codants connues à ce
jour. Ces ARN sont impliqués dans de nombreux processus essentiels des cellules tels que la
synthèse des protéines, la maturation des ARN messagers, les processus de régulation pré- et
post-transcriptionnelle, . . .
Afin de clarifier le discours, les gènes à l’origine d’ARN messagers seront par la suite
appelés des gènes codants, et par opposition, les gènes à l’origine d’ARN non-codants, des
gènes à ARN.
7
Chapitre 1. Les Acides RiboNucléiques
Fig. 1.3 – Structure des acides nucléiques. Chaque groupe phosphate (P) est lié à un sucre
(dR ou R) lui-même lié à une base azotée (A, C, G, T ou U).
8
1.1. L’ARN au sein de la cellule
Source http://openclipart.org/media/files/hs/1771
L’organisation de l’information génétique stockée dans l’ADN varie selon les organismes.
Le génome des eucaryotes est généralement organisé en plusieurs molécules d’ADN empa-
quetées, les chromosomes. Le génome des procaryotes et des archées n’est en général constitué
que d’un seul chromosome qui se présente sous forme circulaire, c’est-à-dire que les extrémités
5′ et 3′ de la molécule d’ADN sont liées ce qui a pour effet de fermer la molécule.
9
Chapitre 1. Les Acides RiboNucléiques
(a) Initiation.
(b) Elongation.
(c) Terminaison.
Source http://en.wikipedia.org/wiki/File:Simple_transcription_initiation1.svg
10
1.2. Les ARN codants
11
Chapitre 1. Les Acides RiboNucléiques
12
1.2. Les ARN codants
Source http://fr.wikipedia.org/wiki/Fichier:Acides_amin%E9s_propri%E9t%E9s_diagramme_Venn.svg
Fig. 1.7 – Diagramme de Venn des acides aminés selon leur propriétés.
13
Chapitre 1. Les Acides RiboNucléiques
MQIFVKTLTGKTITLEVEPSDTIENVKAKIQDKEGIPPDQQRLIFAGKQLEDGRTLSDYNIQKESTLHLVLRLRGG
(a) Structure primaire.
Fig. 1.8 – Structure de l’ubiquitine d’Homo sapiens. Les hélices alpha sont en rouge, les
feuillets bêta en jaune.
Fig. 1.10 – Traduction du début d’un gène codant pour la poly-ubiquitine C. Source : Ref-
Seq, Homo sapiens, NM 021009.
14
1.2. Les ARN codants
traduits correpond à la séquence codante d’un ARN messager et est appelé le cadre ouvert
de lecture.
Même si l’immense majorité des organismes vivants utilisent le code génétique standard,
on note toutefois quelques exceptions à cette règle chez certains organismes pour lesquels
les acides aminés codés ne sont pas les mêmes. Par exemple, chez les champignons Candida,
le codon CUG habituellement traduit par la leucine correspond à la sérine, ou encore chez
certains procaryotes où le codon STOP UAG code parfois pour un acide aminé supplémentaire,
la pyrrolysine.
La région 3′ non traduite d’un ARN messager contient le signal de poly-adénylation utilisé
lors de la maturation, mais également des sites de fixation pour des protéines dont le rôle est
d’orienter vers sa destination finale dans la cellule l’ARN messager puis la protéine produite.
Les régions 5′ et 3′ non traduites peuvent en sus comporter d’autres signaux pour la régulation
traductionnelle, c’est-à-dire des sites destinées à la fixation d’autres molécules venant activer
ou au contraire éteindre la traduction.
La traduction est un processus séquentiel et linéaire. Certains éléments dans l’ARN mes-
sager peuvent parfois perturber la lecture des triplets par le ribosome allant jusqu’à entraı̂ner
un glissement du ribosome d’une ou plusieurs bases “en avant” ou “en arrière”. Le glisse-
ment du ribosome induit un changement de cadre de lecture, également appelé frameshift.
Ce processus est induit par certaines répétitions ou motifs dans la séquence souvent accom-
pagnés d’une petite structure locale formée par l’ARN messager. Bien qu’il s’agisse le plus
souvent d’une erreur, le glissement du ribosome peut être une action programmée. Plusieurs
15
Chapitre 1. Les Acides RiboNucléiques
études mettent en évidence des changements de cadre de lecture programmés chez les virus,
les éléments transposables (section 1.4.2) [Jac88, GA96].
tige
A
A C U C G
G
U
U G A G C A G boucle
Une structure est décrite par une classification en quatre niveaux hiérarchiques :
– la structure primaire est simplement la séquence, orientée de 5′ en 3′ , des bases qui
composent la molécule ;
16
1.3. Les ARN non-codants
– la structure secondaire est l’ensemble des appariements sans croisement, formant des
tiges emboı̂tées ou juxtaposées comme illustré sur les schémas de la figure 1.13 ;
– la structure tertiaire est l’ensemble de tous les appariements. En plus des appariements
de la structure secondaire, les appariements suivants sont donc autorisés : les pseudo-
nœuds (appariements chevauchants), les triplets (appariements à trois), les quadruplets
(à quatre) et les appariements isolés ;
– la structure spatiale désigne la configuration de la molécule dans l’espace, généralement
en interaction avec d’autres molécules.
La figure 1.14 présente les différentes manières usuelles de représenter les structures d’ARN
selon cette hiérarchie. La figure 1.15 montre la structure tertiaire d’un ARN ribosomique 18S
de levure.
La stabilité d’une molécule d’ARN est mesurée par son énergie libre qui est issue des
principes de la thermodynamique. Plus l’énergie libre d’une structure est faible, plus celle-
ci est stable. Les tiges stabilisent une structure, tandis que les boucles la déstabilisent. La
stabilité apportée par une tige est fonction de sa longueur et de la nature de ses appariements :
les appariements canoniques (G≡C, A=U et G=U) sont plus stables que les appariements non
canoniques (G−A, C−U, . . .). Toutes ces caractéristiques sont reprises dans le modèle d’énergie
libre de Turner [TSF88, MSZT99].
La description d’une structure secondaire ou tertiaire d’un ARN par la simple énumération
de ses appariements peut s’avérer insuffisante pour décrire les interactions avec d’autres
molécules. Ces interactions font intervenir des motifs particuliers dont la description requiert le
plus souvent l’utilisation d’une représentation plus fine. A cet effet, la classification de Leontis-
Westhof [LW01] a largement été adoptée par la communauté pour la description d’intéractions
tridimensionnelles caractéristiques. L’observation de structures tridimensionnelles réelles a no-
tamment permis l’élaboration de la base de données Scor [KTHB02] contenant plus de 8 000
motifs récurrents.
17
Chapitre 1. Les Acides RiboNucléiques
GUAUCUCGCACGGUUACACCCCAUCCUUCGGGAGGGCUUUAGGGGUGUGCUG
GAAGCACCACCGGACAGCCGGAACAUUGCCGAAAGGCAGCCC
(a) Structure primaire.
60
C A C C A
40 G C
A C
U U A G
U G G
C
A G A
G
G 50 U C
G G G
C G U C A
G C C G U
A A C A G G
U C A U C 70
G C 20
G U C
C
30 G G G
U
C U G G
C A
A A
10 C C
G A
C U
U G U 80
C C A
U C G
A U G C C
90 C
G C
5’ G G
A A
A
GUAUCUCGCACGGUUACACCCCAUCCUUCGGGAGGGCUUUAGGGGUGUGCUGGAAGCACCACCGGACAGCCGGAACAUUGCCGAAAGGCAGCCC
(c) Structure secondaire sous forme arc-annotée.
Fig. 1.14 – Structures de l’ARN. Exemple d’un élément non traduit structuré de l’ARN
génomique du Tombusvirus. Source : Rfam, RF00176.
18
1.3. Les ARN non-codants
19
Chapitre 1. Les Acides RiboNucléiques
pliquée dans la régulation de la traduction, les micro ARN, a été découverte grâce à une
approche bio-informatique, l’analyse comparative de génomes, confirmée par des méthodes
plus classiques de biochimie [Ruv01].
1.4.1 Généralités
L’évolution est causée par la présence de variations parmi les traits héréditaires d’une
population, et par divers mécanismes qui favorisent la propagation de certains traits plutôt que
d’autres. Par la sélection naturelle, les traits héréditaires favorisant la survie et la reproduction
des individus d’une population voient leurs fréquences croı̂tre d’une génération à l’autre. La
pression de sélection correspond à un ensemble de contraintes environnementales auxquelles
est assujettie une population d’individus, comme par exemple la composition chimique d’un
milieu.
D’un point de vue génétique, l’évolution des espèces, guidée par la sélection naturelle et
la pression de sélection, est engendrée par une évolution des séquences génomiques. Plusieurs
mécanismes sont impliqués dans l’apparition de mutations dans les séquences génomiques,
c’est-à-dire des altérations induites ou spontanées des acides nucléiques.
A l’échelle génomique
Essentiellement trois familles de mécanismes interviennent dans l’évolution des acides
nucléiques à l’échelle génomique : les transferts horizontaux, les éléments transposables et
les recombinaisons chromosomiques. Ces mécanismes sont étroitement liés à la plasticité du
génome et produisent des insertions et/ou des délétions de séquences plus ou moins longues
de nucléotides dans l’ADN, soit délibérément, soit par erreur.
Un transfert horizontal est une intégration d’un fragment d’ADN au sein du matériel
génétique d’un organisme. Il peut prendre trois formes : la transformation, la conjugaison et
la transduction. La transformation consiste pour un organisme à insérer dans son génome du
matériel génétique présent dans son environnement proche. Ce matériel peut, par exemple,
20
1.4. L’évolution des acides nucléiques
provenir d’un autre organisme mort dont l’ADN s’est retrouvé “libéré”. Les transferts ho-
rizontaux concernent principalement les organismes procaryotes mais ont déjà été observés
chez des eucaryotes unicellulaires comme les levures et les champignons. Comme son nom l’in-
dique, la conjugaison bactérienne ne s’applique qu’aux bactéries. Ce processus est un échange
unidirectionnel de matériel génétique d’une bactérie vers une autre bactérie. Après avoir re-
copié une portion de son ADN, appelée plasmide, une bactérie transfère cette copie à une
autre bactérie par le biais d’un canal chimique. Une fois transmise, cette copie peut ou non
être intégrée dans le génome de la bactérie receveuse. Le dernier type de transfert horizontal
est la transduction. L’intégration de matériel génétique étranger est ici réalisée par un virus.
Tous les virus ont un cycle lytique, ou infectieux, pendant lequel ils injectent leur matériel
génétique. La machinerie cellulaire de l’hôte est alors détournée pour produire des copies du
virus jusqu’à ce que la cellule hôte éclate. Chaque copie produite peut alors infecter une autre
cellule et soit suivre un nouveau cycle lytique, soit suivre un cycle lysogénique. Lors d’un
cycle lysogénique, le matériel génétique du virus s’intègre au matériel génétique de l’hôte qui
le transmet à ses descendants. Durant le cycle lytique du virus, une partie de l’ADN de l’hôte
peut être “emportée” avec une copie de l’ADN du virus. S’il s’en suit un cycle lysogénique,
alors le fragment d’ADN emporté est intégré au matériel génétique de l’hôte en même temps
que l’ADN du virus, ce qui constitue la transduction.
Un élément transposable est une séquence capable de se déplacer et/ou de se multiplier de
manière autonome dans un génome. On distingue couramment deux types d’éléments trans-
posables selon leur mode de transposition : les éléments à ARN, ou de classe I, et les éléments
à ADN, ou transposons de classe II. Les éléments de classe I, présents uniquement chez les
eucaryotes, fonctionnent selon le principe du “copier-coller”, c’est-à-dire qu’ils ne se déplacent
pas dans un génome mais qu’une copie de l’élément est insérée ailleurs dans le génome. Ce
processus, appelé transposition réplicative, fonctionne en trois temps : la séquence génomique
de l’élément est transcrite en ARN, puis l’ARN synthétisé est rétro-transcrit en ADN, et ce
fragment d’ADN est finalement inséré dans le génome. Les éléments de classe I sont ainsi
appelés rétro-transposons ou rétro-posons. Les éléments de classe II, également appelés trans-
posons, peuvent être sujets à une transposition réplicative ou conservative. La transposition
conservative suit alors le principe du “couper-coller”, c’est-à-dire que l’élément transposable
est excisé du génome puis réinséré ailleurs. Quelque soit la nature de la transposition, celle-
ci ne se fait pas de manière aléatoire, mais au niveau d’un court motif spécifique dans la
séquence génomique. Toutefois, l’identification du site peut être approximative provoquant
l’excision ou l’insertion d’un fragment de séquence voisin de l’élément transposable lors d’une
transposition.
Les recombinaisons chromosomiques interviennent exclusivement chez les organismes eu-
caryotes durant la reproduction sexuée des espèces, et plus particulièrement au cours de la
méı̈ose, c’est-à-dire la division d’une cellule permettant d’obtenir quatre cellules sexuelles.
Les recombinaisons assurent le brassage génétique par la formation de nouvelles combinai-
sons génétiques, essentielles à la diversité d’une population et à l’évolution des espèces. Il
existe deux principes complémentaires de recombinaisons chromosomiques : les recombinai-
sons inter-chromosomiques et les recombinaisons intra-chromosomiques. La recombinaison
inter-chromosomique désigne la séparation aléatoire des chromosomes homologues d’une cel-
lule au cours de sa méı̈ose. Etant donné une cellule comportant n paires de chromosomes
homologues, la méı̈ose de cette cellule produit quatre cellules comportant chacune un exem-
plaire de chaque paire. Le nombre de combinaisons de chromosomes possibles est donc de
21
Chapitre 1. Les Acides RiboNucléiques
2n ce qui, pour une cellule humaine composée de 23 paires de chromosomes, représente plus
de huit millions de combinaisons. La recombinaison intra-chromosomique, également appelée
enjambement, est un échange de segments entre deux chromosomes homologues au niveau de
sites précis des chromosomes, appelés chiasmas. On dénombre en moyenne entre un et cinq
sites possibles entre deux chromosomes homologues. La recombinaison intra-chromosomique
peut-être déséquilibrée, c’est-à-dire qu’un fragment d’ADN peut être inséré ou au contraire
délété dans les chromosomes. Dans ce cas, les conséquences varient selon la longueur du
fragment inséré ou délété, mais surtout de la région affectée.
A l’échelle nucléotidique
Les mécanismes décrits précédemment ont des conséquences à l’échelle des génomes en
provoquant des insertions ou des délétions de fragments de séquences entiers relativement
longs. Les mécanismes auxquels nous allons maintenant nous intéressés sont plus discrets et
portent sur l’insertion, la délétion et la substitution, c’est-à-dire le remplacement, de quelques
nucléotides.
Les insertions de nucléotides peuvent être le résultat de duplications de certains fragments,
à l’image des transpositions réplicatives à l’échelle du génome. Certains facteurs extérieurs
peuvent également à l’origine d’insertion ou la délétion de nucléotides, telle que l’exposition
à des rayons ultraviolets par exemple.
Les substitutions de nucléotides peuvent être classées en deux groupes : les transitions et
les transversions. Une transition est une substitution d’une purine par une autre purine, ou
d’une pyrimidine par une autre pyrimidine, tandis qu’une transversion est une substitution
d’une pyrimidine par une purine ou inversement. La transition d’un C en U est un phénomène
fréquent résultant de la dégradation spontanée par désamination de la cytosine. Cette modi-
fication est réversible grâce à un processus qui détecte l’uracile dans l’ADN. Cependant si la
réparation n’est pas effectuée avant la prochaine réplication de l’ADN, la guanine appariée à
la cytosine d’origine sur le brin opposé est substituée par une adénine lors de la réplication,
et l’uracile remplacée par une thymine. Le second type de mutation spontanée est lié au di-
nucléotide 5′ –CG, appelé un “point chaud”, car il est l’objet de fréquentes mutations lorsque
la cytosine en 5′ est sous sa forme méthylée. La désamination de la méthylcytosine produit
en effet une thymine, qui ne peut alors être reconnue par aucun mécanisme de réparation. Le
troisième type de mutation spontanée est l’oxydation des nucléotides par les radicaux libres
de l’oxygène, sous-produits du métabolisme oxydatif normal des cellules. Une guanine oxydée,
par exemple, s’apparie à tort avec une adénosine induisant une transversion de G-C en T-A.
Les probabilités d’apparition, de conservation et de transmission d’une mutation sont les
objets de nombreuses études [Rus93, Cro97, DCCC98]. Dans le cadre de nos travaux, nous
nous sommes plus particulièrement intéressés aux mutations dans les séquences des gènes
codants et des gènes à ARN car sous la pression de sélection, ces mutations suivent certains
schémas à l’origine de biais locaux.
22
1.4. L’évolution des acides nucléiques
Les substitutions dans une séquence codante sont classées en trois catégories en fonction
de leur impact sur le codon modifié :
– une mutation faux-sens : le codon affecté ne code plus pour le même acide aminé. L’im-
pact de ce type de mutation sur la protéine produite dépend du rôle de l’acide aminé
original dans la protéine et de l’acide aminé qui lui a été substitué. On parle parfois
de mutation synonyme lorsque le nouvel acide aminé codé a des propriétés physico-
chimiques proches de l’acide aminé codé avant la mutation (section 1.2.1) ;
– une mutation non-sens : le codon affecté ne code plus pour un acide aminé mais pour
un codon STOP. La protéine produite est alors tronquée comme cela peut se produire
avec une insertion ou une délétion décalante ;
– une mutation silencieuse : l’acide aminé codé reste le même, donc cette mutation n’a
aucune conséquence sur la protéine codée. Ce type de substitution est rendu possible
grâce aux nombreuses redondances dans le code génétique. Cette propriété est également
appelée dégénérescence du code génétique.
Une insertion ou une délétion dans une séquence codante peut allonger ou réduire la
longueur de la protéine codée. En particulier, on parle de mutation décalante si la séquence
insérée ou supprimée provoque un changement de cadre de lecture, c’est-à-dire lorsque la
longueur de la séquence en question n’est pas multiple de trois. Le plus souvent, une mutation
décalante entraı̂ne l’apparition d’un codon STOP prématuré dans le cadre de lecture, ce qui a
pour effet de produire une protéine tronquée lors de la traduction. Dans la plupart des cas, la
protéine ainsi tronquée n’est plus capable d’assurer sa fonction. Il existe des exemples où ce
type de mutation n’est pas létal et peut même conférer un avantage significatif aux individus
qui en sont porteurs, comme par exemple le variant CCR5-∆32 du gène CCR5 [GS03]. Le
gène CCR5 code pour une protéine qui se trouve à la surface de certaines cellules, notamment
des cellules immunitaires, servant de récepteur à chemokine. Le variant ∆32 du gène CCR5
est issu de la délétion décalante de 32 nucléotides dans le gène et dont la traduction génère
une protéine plus courte incapable d’assurer son rôle. Cependant, cette mutation prodigue
aux individus porteurs une immunité naturelle à certains virus comme la petite variole et le
HIV car ces virus sont alors incapables d’infecter les cellules immunitaires dont les récepteurs
à chemokine sont absents ou non fonctionnels. Le variant CCR5-∆32 est largement répandu
de nos jours en Europe où il touche entre 5 et 14% de la population, mais est beaucoup plus
rare en Afrique ou en Asie. L’hypothèse la plus vraissemblable est que ce variant a fait l’objet
d’une sélection naturelle au cours de la pandémie de peste noire qui a décimé plus d’un tiers
de la population européenne au XIVème siècle, les porteurs du variant CCR5-∆32 étant alors
plus résistants que les autres à la maladie.
23
Chapitre 1. Les Acides RiboNucléiques
Toutefois, contrairement aux mutations silencieuses dans les régions codantes qui n’ont
aucun effet sur la protéine codée, les substitutions qui préservent les appariements modifient
la stabilité locale et globale de la molécule. Par exemple, pour deux bases appariées G et C, la
transition de C en T peut préserver l’appariement puisque G et T peuvent s’apparier mais ce
nouvel appariement est moins stable que l’appariement original. Au delà de la stabilité des
structures d’ARN c’est la conformation spatiale locale et globale de la structure de l’ARN qui
est soumise à la pression de sélection. Au voisinage de certains motifs structuraux, certains
nucléotides, appariés ou non, sont en effet exposés d’une manière particulière. La conformation
de ces nucléotides est le plus souvent primordiale à la fonction de la molécule. Les séquences
de ces motifs structuraux, au même titre que les régions non appariées, sont donc contraintes
d’être conservées.
24
1.5. L’analyse comparative de séquences nucléiques
AAN33049 1 CGAATGCCAGGCCCAGCCCTCA---CCTCTCGCTCCGCAGGGGGGAGTCG 47
||| ||||..|| ..||||||||| ||.||
AAA31576 1 ATG--------AGCCGGCAGAGTATCTCGCTCC--------GATTC- 30
AAN33049 48 CCTGCACCGGTGGCCGCTGCTCCTGCTGCTGCTGCTGCTGC-TCCC---- 92
||||||||.||.|||||||||..|| | ||||
AAA31576 31 -------------CCGCTGCTTCTCCTGCTGCTGTCGC--CATCCCCCGT 65
Fig. 1.16 – Exemple d’alignement semi-global entre deux fragments de séquences homologues
de gènes codants pour la prostaglandine dont le pourcentage d’identité est de 44,9%.
25
Chapitre 1. Les Acides RiboNucléiques
(a) Alignement correct de deux fragments de séquences codant pour une prostaglandine dont le pourcentage
d’identité est 77,8%.
AE008837.1 CGCGGGGUGGAGCAGCCUGGUAGCUCGUCGGGCUCAUAACCCGAAGGUCGUCGGUUCAAAUCCGGCCCCCGCAA
||||||||.|||||| .||||||||||.|||||||||||||||.||||||..|||||.|.|||.||||||||||
X16759.1 CGCGGGGUAGAGCAG-UUGGUAGCUCGCCGGGCUCAUAACCCGGAGGUCGCAGGUUCGAGUCCUGCCCCCGCAA
(b) Alignement correct de deux séquences d’ARN de transfert dont le pourcentage d’identité est de 86,5%.
Fig. 1.17 – Deux alignements semi-globaux optimaux corrects de séquences codantes homo-
logues (a) et de séquences partageant une structure commune (b).
26
1.5. L’analyse comparative de séquences nucléiques
27
Chapitre 1. Les Acides RiboNucléiques
(a) Alignement incorrect de deux fragments de séquences codant pour une intégrine dont le pourcentage d’iden-
tité est 37,8%.
M22657.1 ACUUUUAAAGGAUAGAAGU-----------AAUCCAUU--GGCCUUAGG-----AGCCAAAAAAUUGGUG--------CAACUCCAAAUAAAAGUA--
|| ||||..|| .||.||.|| |.||.|||.||||||| |.||.| |.|.|
K01921.1 -----------------GUCUCUGUGGCGCAAUCGGUUAGCGCGUUCGGCUGUUAACCGAAAGAUUGGUGGUUCGAGCCCACCC-------AGGGACG
(b) Alignement erroné de deux séquences d’ARN de transfert dont le pourcentage d’identité est de 34,7%.
Fig. 1.18 – Deux alignements semi-globaux optimaux incorrects de séquences codantes ho-
mologues (a) et de séquences partageant une structure commune (b).
28
1.5. L’analyse comparative de séquences nucléiques
Les méta-séquences
Afin d’analyser convenablement les ensembles de séquences hétérogènes en terme de
conservation, nous proposons de créer des sous-ensembles de séquences suffisamment simi-
laires pour pouvoir les aligner de manière fiable. Pour la construction des sous-ensembles de
séquences similaires, nous définissons une représentation des données sous forme de graphe.
Soit S = {s1 , . . . , sn } un ensemble de n séquences nucléiques et id(si , sj ) le pourcentage
d’identité entre deux séquences si et sj . Nous créons ensuite une partition P de S compor-
tant m parties {P1 , . . . , Pm }. Chaque partie Pk est composée d’un sous-ensemble de séquences
de S connectées par la relation de similarité id(si , sj ) ≥ α, où α est un seuil sur le pourcen-
tage d’identité au delà duquel on considère que les séquences présentent un degré de similarité
suffisant. Chaque partie de P est appelée une méta-séquence. Une méta-séquence fait ainsi
référence soit à une seule séquence, soit à un ensemble de plusieurs séquences qui seront
représentées par leur alignement. La figure 1.19 présente un exemple de quatre séquences
regroupées en deux méta-séquences, symbolisées par les parties grisées. Sur cet exemple, le
couple de séquences s1 et s2 d’une part, et le couple s1 et s3 d’autre part, présentent un
pourcentage d’identité supérieur à α. Les séquences s1 , s2 et s3 sont donc regroupées pour
former la partie P1 , et la séquence s4 forme à elle seule la partie P2 .
P1
s1 s2
s4 s3
P2
Fig. 1.19 – Exemple de création de deux méta-séquences à partir de quatre séquences. Chaque
sommet noir correspond à une séquence, les arêtes relient les séquences dont le pourcentage
d’identité est supérieur au seuil α. Les régions grisées correspondent aux méta-séquences
correspondantes.
29
Chapitre 1. Les Acides RiboNucléiques
30
Chapitre 2
Les approches ab initio ont pour objectif de prédire l’ensemble des gènes présents dans
une séquence nucléique sans autre connaissance extérieure. Pour cela, elles tirent parti des
signaux présents dans la séquence et des biais de composition des séquences codantes.
31
Chapitre 2. Recherche de gènes et régions codantes
32
2.1. Les méthodes ab initio
Source http://www.pnas.org/cgi/doi/10.1073/pnas.0703773104
Fig. 2.1 – Représentation graphique de deux PWM décrivant les sites d’épissage. La fréquence
de chaque nucléotide à une position est proportionnelle à la taille dans la lettre qui le
représente. Sur la représentation du site donneur, à gauche, on voit clairement apparaı̂tre une
séquence consensus qui caractérise le début de l’intron GT. On remarque le même phénomène
pour le site accepteur qui marque la fin des introns et qui est caractérisé par la séquence
consensus AG.
obtenir ces probabilités, on procède par comptage sur un ensemble d’apprentissage constitué
de séquences contenant le motif à caractériser. Cet ensemble doit contenir suffisamment de
séquences pour que les observations réalisées fassent du sens. De plus, comme la forme des
signaux peut être propre à chaque espèce, il est nécessaire de disposer des séquences de même
provenance, ou dont on est sûr qu’elles sont biaisées de la même manière. Un modèle construit
sur une espèce ne pourra donc pas être directement appliqué à une autre espèce. Toutefois,
certains signaux comme les signaux liés à la transcription sont relativement bien conservés
entre les espèces. Un modèle descriptif d’un signal peut donc sous certaines conditions être
transféré d’une espèce à une autre. Plus l’ordre du modèle de Markov est important, plus
la quantité de données d’apprentissage nécessaire est importante. Etant donné que quatre
nucléotides différents sont possibles à une position donnée, il est nécessaire de déterminer
4k+1 probabilités pour un modèle d’ordre k. Une fois le modèle de Markov construit, son
exploitation est assez simple : étant donnés une séquence et un modèle de Markov, on peut
calculer la probabilité que la séquence ait été générée selon ce modèle avec, par exemple,
l’algorithme de Viterbi [Vit67]. A l’heure actuelle, les WAM sont très largement utilisés pour
décrire et détecter les sites d’épissage [TMM07, Sto90, Gel95].
33
Chapitre 2. Recherche de gènes et régions codantes
34
2.2. Les approches par homologie de séquence
de HMM la détection des régions codantes et non codantes est confiée à des modèles Marko-
viens dont les caractéristiques varient selon la nature de la région modélisée. Historiquement,
EcoParse [KMH94] est le premier programme à avoir intégré un modèle de Markov caché.
D’autres ont suivi comme Genie [KHRE96], Fgenesh [SS00], Augustus [SW03, Sta03],
Genie [KHRE96], GeneMark.hmm [LB98, LTHCB05] et Genscan [BK97]. Genscan et
GeneMark.hmm sont les programmes les plus utilisés actuellement pour la richesse de leur
modélisation de la structure des gènes et parce qu’ils ont été paramétrés sur plus d’une cen-
taine d’organismes.
35
Chapitre 2. Recherche de gènes et régions codantes
que les décalages négatifs (-1 ou -2). Contrairement à ses homologues, BlastX ne gère pas
les décalages potentiels du cadre de lecture. Enfin, il existe une version enrichie de BlastX,
nommée BlastC [SG94], qui prend en compte dans son système de score les éventuels biais
d’usage des codons connus et avérés. On estime qu’environ 50% des gènes peuvent être iden-
tifiés à partir des séquences présentes dans SwissProt et Pir [MSSR02] à l’aide de ces
outils.
Même lorsque l’on dispose d’une protéine similaire, il est difficile de déterminer la structure
complète d’un gène, particulièrement les bornes des extrémités 5′ et 3′ non traduites. Qui plus
est, l’identification précise de la séquence codante peut également s’avérer incomplète : tous
les domaines protéiques ne sont pas nécessairement partagés, en cas d’épissage alternatif
notamment, et certains exons peuvent ne pas être attrapés à cause de leur petite taille.
La recherche de similarité peut ensuite donner lieu à des traitement plus sophistiqués.
Historiquement, Procrustes [GMP96] est la première méthode de prédiction de gènes à
exploiter la similarité de séquences au niveau peptidique. Etant donnée une protéine simi-
laire à la séquence nucléique d’intérêt, trouvée manuellement par FastX ou BlastX par
exemple, la méthode déployée par Procrustes consiste à aligner la séquence nucléique et la
séquence protéique. Dans un premier temps, Procrustes identifie tous les exons potentiels
en recherchant les bornes exons/introns selon leurs séquences consensus strictes. Les exons
potentiels sont ensuite traduits puis alignés sur la séquence peptidique fournie. Enfin, Pro-
crustes assemble par programmation dynamique un enchaı̂nement d’exons qui maximise le
score de similarité avec la protéine tout en respectant la structure minimale du gène, c’est-
à-dire une séquence codante qui commence par un codon START, ne contient pas de codon
STOP dans le cadre de lecture et se termine par un codon STOP. GeneWise [BCD04] ou
encore GenomeScan [YLB01] effectuent la même tâche que Procrustes. Ces derniers sont
globalement plus tolérants que Procrustes : ils n’interdisent pas strictement les décalages
de cadres de lecture ainsi que les codons STOP dans le cadre de lecture mais ces deux éléments
sont fortement pénalisés. De plus, les sites d’épissage dans GeneWise et GenomeScan sont
détectés à l’aide d’un HMM, fortement inspiré de celui de Genscan. Dans GenomeThrea-
der, les sites d’épissage sont décris par un représentation ad hoc équivalente aux HMM
précédents. La différence entre ces trois programmes réside essentiellement dans l’évaluation
de la similarité entre les séquences : GenomeScan utilise les résultats de BlastX, Gene-
Wise évalue leur similarité à l’aide d’un modèle pair-Markov caché (pair-HMM), et Ge-
nomeThreader utilise un algorithme d’alignement local par programmation dynamique.
ORFGene2 [RMK96] et PredictGenes [GHKB00] sont deux méthodes équivalentes à Ge-
neWise dans la modélisation adoptée. Toutefois elles intègrent la recherche de protéines
similaires en interrogeant la banque SwissProt. Toutes ces méthodes sont destinées à la
prédiction de gènes eucaryotes, en majorité entraı̂nées et testées chez l’Homme, la souris ou des
plantes. Elle ont été conçues pour travailler sur des séquences provenant d’organismes proches.
Elles fournissent d’ailleurs d’excellents résultats sur des séquences très conservées, mais leurs
performances se révèlent moyennes sur des séquences relativement peu conservées [TPP99].
36
2.2. Les approches par homologie de séquence
complémentaires. Ces ADN complémentaires, notés ADNc, sont obtenus par transcription in-
verse d’ARN matures. Il existent plusieurs protocoles pour obtenir les séquences d’ADNc rétro-
transcrits. Le séquençage “classique” d’un ADNc permet d’en obtenir la séquence complète
et ce de manière fiable. Avant que soit mis au point ce protocole, le séquençage des ARN
messagers se faisait par un protocole à haut débit moins fiable. Ce protocole consiste à
séquencer quelques centaines de nucléotides en une seule fois à chaque extrémité d’un ADNc.
Ces fragments nommés des EST, acronyme pour Expressed Sequence Tags, ne représentent
donc qu’une information partielle par rapport à la taille de certains ADNc qui peuvent at-
teindre plusieurs milliers de nucléotides. La figure 2.2 illustre de manière schématique des
séquences d’EST obtenues pour un ARN messager mature.
5’ EST 3’ EST
AAAAAAAAAA
5’ UTR Séquence codante 3’ UTR
Fig. 2.2 – Les EST sont issus du séquençage partiel des extrémités d’un ARN mature, ici un
ARN messager.
Les ADNc et les EST représentent les informations les plus pertinentes dont on peut
disposer pour établir la structure précise des gènes surtout s’ils sont issus du même orga-
nisme que la séquence nucléique à annoter [FSY+ 99]. En effet, comme les ADNc proviennent
d’ARN transcrits, ils contiennent en plus de la séquence codante les extrémités 5′ et 3′ non
traduites. Les données d’ADNc et d’EST disponibles dans les banques dépendent des condi-
tions expérimentales dans lesquelles les ARN ont été extraits : dans quel type de tissus ? A
quel stade de développement ? . . . Les données disponibles pour un organisme ne couvrent
donc qu’une partie de son transcriptome, elles ne sont pas nécessairement représentatives de
tous ses gènes. De plus, bien que les ADNc et les EST correspondent à des séquences trans-
crites, ces séquences ne sont pas obligatoirement traduites en protéines. Ces caractéristiques
des ADNc et des EST en font des indices précieux pour confirmer des prédictions effectuées
par ailleurs, notamment pour déterminer la structure des gènes. A ce titre, les ADNc sont,
par nature, une source d’information plus complète que les EST qui ne représentent qu’une
source d’information partielle.
Les banques de données d’EST sont fortement redondantes du fait de la technique em-
ployée pour les obtenir. Pour traiter ce problème de redondance, EbEST [JJ98] et Pa-
gan [KS01] procèdent à un regroupement des EST trouvés dans les banques. Les EST
chevauchants sont regroupés puis un EST représentatif de chaque groupe est sélectionné.
Les EST représentatifs sont enfin réalignés avec la séquence nucléique. Les séquences d’EST
comportent un certain nombre d’erreurs introduites durant leur séquençage. La procédure
d’alignement est donc paramétrée pour tolérer les substitutions, les insertions et les délétions
d’un nucléotide.
Quelques méthodes se servent des séquences d’ADN complémentaires présentes dans les
banques de données. Ces programmes, tels que Aat [HAZK97] ou Sim4 [FHZ+ 98] interrogent
les banques de données d’ADNc à la recherche d’ADNc similaires puis réalignent les ADNc
trouvés avec la séquence d’intérêt afin de localiser plus précisément les exons. Lorsque l’on
dispose d’ADNc complémentaires provenant du même organisme que la séquence à annoter,
l’identification des exons peut se révéler particulièrement fiable et précise [FSY+ 99].
37
Chapitre 2. Recherche de gènes et régions codantes
38
2.4. Protea
séquences codantes. Cependant, Qrna emploie trois modèles de Markov cachés différents pour
évaluer la probabilité de trois hypothèses : le modèle COD permet d’évaluer la probabilité
que les séquences alignées soient des séquences codantes homologues, le modèle RNA évalue
la probabilité que les séquences alignées soient des séquences non-codantes structurées homo-
logues, et enfin le modèle OTH qui évalue la probabilité que les mutations entre les séquences
soient fortuites. L’avantage certain de Qrna par rapport à Syncod est de distinguer dans son
modèle COD les mutations faux-sens qui produisent des acides aminés aux propriétés physico-
chimiques proches. Cette caractéristique lui permet de considérer des séquences séparées par
une distance évolutive plus importante.
Alors que Syncod et Qrna travaillent sur des séquences sorties de leur contexte
génomique, Rosetta [BPM+ 00a, BPM+ 00b] propose un point de vue différent en tra-
vaillant sur un alignement complet de deux génomes. Rosetta est cependant exclusive-
ment paramétré et destiné à la comparaison des génomes de l’Homme et de la souris.
Rosetta utilise deux critères essentiels pour effectuer ses prédictions : la synténie d’une
part, et la similarité au niveau nucléique et peptidique d’autre part. Les régions codantes
des deux organismes sont en effet exceptionnellement bien conservées, approximativement
85% d’identité au niveau nucléique, en comparaison de leurs séquences introniques, envi-
ron 35% d’identité [MZB96, MB98, LMS+ 95, KH94]. D’autres méthodes telles que TwinS-
can [KFDB01], AGenDA [TRG+ 03], Utopia [BRS03], Pro-Gen [NGM01], CEM [BH00]
et Sgp-1 [WGJMOG01] fonctionnent de manière analogue à Rosetta. Plus récemment, Pro-
jector [MD04] et GeneMapper [CP06] constituent un compromis entre Rosetta et Qrna.
Ces méthodes proposent une adaptation de la technique déployée dans Qrna en intégrant la
synténie comme le fait Rosetta pour la comparaison de deux génomes complets alignés.
Enfin, un dernier ensemble de méthodes travaille non plus sur deux séquences, mais sur un
ensemble de séquences alignées. ExoniPhy [SH04] et EvoGene [PH03] notamment étendent
et affinent l’approche de Qrna. Pour cela elles requièrent en plus de l’alignement multiple
un arbre phylogénétique correspondant aux séquences alignées. La connaissance des distances
évolutives permet d’estimer avec plus précision les taux de mutations attendus entre séquences
codantes, notamment le taux de mutations silencieuses. Ces deux méthodes utilisent des
modèles phylogénétiques de Markov cachés. Le principe est identique à celui d’un modèle
de Markov caché classique à la différence que les probabilités conditionnelles du modèle sont
obtenues par des fonctions paramétrées par l’arbre phylogénétique. Dans le cadre de la compa-
raison de génomes complets alignés, N-Scan [GB06] propose une adaptation de TwinScan
pour traiter des alignements multiples. L’approche de N-Scan, plus simple que celle de Exo-
niPhy et EvoGene, ne requiert pas d’arbre phylogénétique mais est limitée à la prédiction
de gènes pour l’Homme et la drosophile.
2.4 Protea
Dans cette section, nous présentons la méthode Protea, que nous avons développée pour
l’identification de séquences codantes. Par rapport aux approches existantes que nous venons
de présenter, Protea présente au moins deux caractéristiques qui font son originalité. En
premier lieu, Protea peut traiter un nombre quelconque de séquences sans nécessiter d’ali-
gnement multiple préalable entre ces séquences. Disposer d’un alignement multiple correct
est en effet une tâche délicate quand la distance évolutive entre les séquences est impor-
tante [Mar08], et peut se révéler une source d’erreurs comme nous l’avons mis en évidence
39
Chapitre 2. Recherche de gènes et régions codantes
dans la section 1.5. En second lieu, Protea repose sur les principes de l’analyse comparative
avec une idée générale qui à notre connaissance n’a jamais été utilisée : il s’agit d’exploiter le
schéma évolutif observé entre séquences codantes à travers la conservation du cadre de lecture.
Les séquences non codantes ne sont pas censées présenter ce schéma évolutif particulier.
Nous commençons par illustrer et valider cet argument pour un couple de séquences. Nous
expliquons ensuite comment étendre ce principe à une famille comprenant un nombre quel-
conque de séquences, sans passer par un alignement multiple, mais en utilisant une structure
de graphe.
40
2.4. Protea
u1 = TACCATAGTCGACATGA v1 = TACCACTGTCAACATGA
U V
u2 = TAGCACTGTCAACAGGA v2 = TACCATTGCCAGCACGA
Cadre de u2 Cadre de v2
-3 -2 -1 1 2 3 -3 -2 -1 1 2 3
-3 10 2 -3 10
Cadre de u1
Cadre de v1
-2 1 12 -2 2 6 1
-1 2 1 18 -1 25
1 15 1 1 37
2 8 1 2 6
3 2 8 3 5 6
(a) Comparaison de toutes les paires de traductions (b) Comparaison de toutes les paires de traductions
de deux séquences similaires. de deux séquences codantes homologues.
41
Chapitre 2. Recherche de gènes et régions codantes
100 100
Prédictions correctes (%)
80 80
60 60
40 40
20 20
0 0
30 40 50 60 70 80 90 100 100 200 300 400 500
Identité moyenne (%) Longueur moyenne
Fig. 2.4 – Proportion de couples de cadres de lecture correctement prédits parmi tous les
couples de séquences de chaque famille de Pandit. Les résultats sont exprimés en fonction
du pourcentage d’identité des séquences (gauche) et de leur longueur moyenne (droite).
42
2.4. Protea
il faudrait alors produire un alignement multiple des n traductions potentielles afin de les
comparer. Cette solution n’est évidemment pas praticable car elle requiert de générer une
quantité d’alignements multiples exponentielle par rapport au nombre de séquences à traiter.
Pour surmonter cette difficulté, il faut trouver une autre manière de généraliser l’approche
pour deux séquences à un ensemble de séquences de taille quelconque. Si l’on suppose qu’il
existe un cadre de lecture globalement conservé entre toutes les séquences homologues, alors on
doit pouvoir détecter cette conservation à partir des comparaisons des traductions potentielles
uniquement entre paires de séquences.
Soit S = {s1 , . . . , sn } un ensemble de n séquences. Nous utilisons la notion de méta-
séquence introduite dans la section 1.5.3. Selon ce processus, on crée une partition P des
séquences de S, où chaque partie correspond à une méta-séquence. On construit ensuite le
graphe des cadres de lecture à partir des m parties de P. Le but de cette structure est de
combiner les comparaisons deux à deux des traductions potentielles.
Définition 1 (Graphe des Cadres de Lecture (GCL)). Le graphe des cadres de lecture G =
(P, E, φ) est un 36-graphe valué non orienté tel que
– P est l’ensemble des m sommets de G correspondants aux m méta-séquences ;
– E ⊂ P × P × C × C est l’ensemble des arêtes du graphe où C = {−3, −2, −1, 1, 2, 3}
est l’ensemble des cadres de lecture possibles ;
– φ(Pi , Pj , c, c′ ) associe à chaque arête de E une valuation issue de la comparaison entre
la méta-séquence Pi traduite selon le cadre c et la méta-séquence Pj traduite selon c′ .
43
Chapitre 2. Recherche de gènes et régions codantes
Score de consistance
Etant donnés deux sommets, les arêtes qui les relient sont triées par valuation décroissante.
On attribue alors à chacune un score qui dépend de son rang : l’arête dont la valuation est
la plus élevée se voit attribuée un score de 6, la seconde 5, jusqu’à 1. Seules les six meilleures
arêtes sont retenues, car comme il est mentionné en section 2.4.1, seuls six couples de cadres
de lecture sont censés émerger de la conservation au niveau nucléique.
On définit une affectation globale de cadres de lecture comme une fonction de P dans C =
{−3, −2, −1, 1, 2, 3} qui attribue un cadre de lecture à chaque sommet d’un GCL. Pour chaque
affectation globale A = (c1 , . . . , cm ), on définit le score de consistance noté consistance(A),
comme la somme des rangs des arêtes induites par A.
X
consistance(A) = rang(Pi , Pj , ci , cj )
1≤i<j≤m
La valeur de ce score est comprise entre M et 6M , où M est le nombre de sommets
connectés dans le GCL. Dans la plupart des cas, M = m(m−1) 2 car chaque couple de sommets
est au moins relié par une arête. La valeur optimale 6M est atteinte lorsqu’une affectation
globale couvre toutes les meilleures comparaisons deux à deux.
La qualité d’un GCL est donnée par le score de consistance le plus élevé obtenu par une
affectation globale. En pratique, il est rarement nécessaire d’énumérer toutes les affectations
pour trouver celle dont le score de consistance est maximal. Les scores des arêtes étant bornés,
on applique les techniques standards de rebroussement et d’élagage pour éviter des calculs
inutiles.
44
2.4. Protea
6M ⌊(i−M )/6⌋
X 1 X M i − 6j − 1
Pr[consistency score(A) ≥ s] = (−1)j
6M j M −1
i=s j=0
f (x) = (x + x2 + · · · + x6 )N ,
où chaque arrangement possible correspond à un terme. f (x) peut également s’écrire comme
une série multinômiale
5
!N
X
N i
f (x) = x x
i=0
N
1 − x6
= xN
1−x
= xN (1 − x6 )N (1 − x)−N
D’où l’expression,
N ∞
N
X
k N 6k
X N +l−1
x (−1) x xl ,
k l
k=0 l=0
donc le coefficient c de xp inclus tous les termes où
p = N + 6k + l.
45
Chapitre 2. Recherche de gènes et régions codantes
p − 6k − 1 p − 6k − 1
= ,
p − 6k − N N −1
donc
⌊(p−N )/6⌋
X
k N p − 6k − 1
c= (−1) .
k N −1
k=0
46
2.4. Protea
ainsi uniquement d’un motif évolutif particulier et non d’une similarité au niveau nucléique.
On exclut toutefois une autre affectation des affectations compatibles, l’affectation sur le
brin opposé où la troisième position des codons coı̈ncide avec celle de l’affectation à évaluer.
Le score de cette affectation est en effet clairement biaisé à cause d’une propriété du code
génétique. En effet, les mutations silencieuses apparaissent le plus souvent sur la troisième
base des codons (section 1.4.2). Lorsque l’on observe un grand nombre de mutations silen-
cieuses entre deux cadre de lecture, on peut donc naturellement obtenir deux cadres de lecture
sur les brins opposés qui présentent ce même biais, comme illustré sur la figure 2.5 où les trois
mutations silencieuses entre les cadres +1 produisent sur les brins opposés deux mutations
silencieuses fortuites. Sur l’exemple de la figure 2.3, le couple de cadres de lecture correct
pour les séquences de V est (1, 1). On constate un certain nombre de mutations silencieuses
entre ces deux cadres, c’est pourquoi le cadre (−1, −1) obtient également un bon score. Pour
interpréter le score d’alignement d’une affectation, on ne garde donc finalement que quatre
des cinq affectations compatibles. La significativité de la déviation du score d’alignement
de la meilleure affectation globale d’un GCL est mesuré par un z-score, même si le nombre
d’observations est faible.
Dans le cas des petits GCL, la détection d’ensembles de séquences codantes homologues se
fait au moyen d’un seuil sur le z-score du score de la meilleure affectation globale. A l’instar de
la classification réalisée sur les GCL de plus grande taille, la valeur de ce seuil a été déterminée
de manière empirique au cours de la validation de la méthode.
Cadre +1 V G N Cadre +1 V G N
GTAGGCAACCGA GTCGGTAATCG
P L R Cadre -2 P F R Cadre -1
Fig. 2.5 – Du fait de la redondance du code génétique, les mutations dites silencieuses (posi-
tions grisées) n’ont aucun effet sur les acides aminés codés. Ces mutations apparaissent plus
fréquemment à la troisième position des codons.
47
Chapitre 2. Recherche de gènes et régions codantes
Trois jeux de données ont été construits pour mesurer les performances de Protea :
des familles de séquences codantes, des ARN non-codants et des séquences aléatoires. Pour
chaque jeu de données, des sous-ensembles ont été construits aléatoirement de façon à ce que
chaque famille ne soit représentée qu’une seule fois et que les sous-ensembles soient les plus
équivalents possible en terme de pourcentage d’identité et de longueur.
48
2.5. Résultats expérimentaux de Protea
VGJ_BPG4 AAAAAATCAATTCGCCGCTCTGGT---------------------------------------GGCAAATCTAAGGGTGCCCGTCTCTGGTATGTAGGCGGAACACAATAC
Q9G087_BPS13 TCTAAAGGTAAAAAACGTTTTGGCGCTCGCTCCGGTCGTCCACAGCCGTTGCGAGGCACTAAAGGCAAGCGTAAAGGCGCTCGTCTTTGGTATGTAGGCGGTCAACAATTT
VGJ_BPAL3 ATGAAGAAAGCACGTCGTTCTCCT---------------------------------------AGTCGTCGTAAAGGTGCTCGCCTCTGGTATGTAGGCGGTTCTCAGTTT
** ** * * * *** ** ** ** ** ************** ** *
seq0 GACATGTTCAACAAATATGCGAAT---------------------------------------ATCAGATCTGACTACCGGTTAGGTAGTCAGGCGGTCGCACCGAGTCAC
seq1 ACGATGAGTACTTTTTACGAGAACCCGGGGTTCTGAATCTGGAGTGCCCCGAACCGCACCACTATCAGGAGTGATTATCGATTAGGCAGTCAGGCGGTCGCTAAGAGTCTT
seq2 GTAATACACGTTAAGTACGCGTGT---------------------------------------CTTGCTAGTGATTACCGATTGGGTAGTCAGGCGGTCGCTTCCAGCCTT
** ** * * * *** ** ** ** ** ************** ** *
Fig. 2.6 – En haut, l’alignement multiple de la famille PF04726 du jeu de données CO-
DANT ; en bas, un alignement multiple obtenu par la procédure de mélange conservative. Les
caractères ’*’ marquent les positions parfaitement conservées.
49
Chapitre 2. Recherche de gènes et régions codantes
Mesure de performances. Pour évaluer la qualité des prédictions effectuées par un clas-
sifieur binaire, deux mesures sont couramment utilisées : la sensibilité et la spécificité. La
sensibilité de la méthode correspond ici à la proportion de séquences codantes homologues
classifiées “codant”, tandis que la spécificité correspond à la proportion de séquences d’un
autre type classifiées “autre”. La sensibilité Sn sur un jeu de données est donnée par
TP
Sn =
TP + FN
où T P désigne la quantité de vrais positifs, c’est-à-dire les prédictions “codant” correctes, et
F N la quantité de faux négatifs, c’est-à-dire les prédictions “autre” incorrectes. De manière
analogue, la spécificité Sp sur un jeu de données est donnée par
TN
Sp =
TN + FP
où T N désigne la quantité de vrais négatifs, c’est-à-dire les prédictions “autre” correctes, et
F P , c’est-à-dire les prédictions “codant” incorrectes.
Résultats généraux. Les résultats de Protea sur les jeux de données CODANT, NON-
CODANT et ALEATOIRE pour des ensembles de 3, 5 et 11 séquences sont répertoriés dans
la table 2.1 selon le pourcentage d’identité moyen et la longueur moyenne des séquences. Ces
résultats sont également présentés graphiquement sur la figure 2.7.
(a) Répartition des prédictions “codant” sur le jeu (b) Répartition des prédictions “codant” sur le jeu
de données CODANT. de données NONCODANT.
Fig. 2.7 – Répartition des prédictions “codant” de Protea sur les ensembles de 11 séquences
des jeux de données CODANT (a) et NONCODANT (b). Les histogrammes tracés en poin-
tillés représentent les ensembles de données initiaux, tandis que les histogrammes en trait
plein représentent les prédictions “codant” de Protea.
Dans la plupart des cas, la sensibilité et la spécificité sont supérieures à 80%. Comme
on pouvait l’espérer, les performances des Protea augmentent avec le nombre de séquences
ainsi que leur longueur. En dessous de 50% d’identité moyenne, Protea est particulièrement
performant et affiche une spécificité supérieure à 90%. En revanche, les performances de
Protea se dégradent au delà de 90% d’identité. Sur de telles séquences très conservées, la
50
2.5. Résultats expérimentaux de Protea
quantité de mutations est insuffisante pour détecter un schéma de substitution en lien avec
la conservation d’une séquence d’acides aminés particulière. Le comportement de Protea
sur des courtes séquences est soumis à un biais causé par la présence d’acides aminés rares,
tel que le tryptophane par exemple, conservés entre des traductions potentielles. Les scores
de substitution de ces acides aminés dans les matrices BLOSUM sont relativement élevés
et la conservation fortuite d’un seul acide aminé de ce type sur des séquences de moins de
cinquante nucléotides entraı̂ne une élévation mécanique des scores lors de leur alignement.
Comparaison avec Qrna. Protea a été conçu pour traiter des ensembles de plusieurs
séquences non alignées. Néanmoins, il peut quand même être utilisé sur des paires de
séquences, bien qu’il ne soit pas spécialement conçu pour ce cas de figure.
Comme Qrna nécessite en entrée des séquences alignées, nous avons testés plusieurs
méthodes d’alignement : ClustalW, T-Coffee, Dialign2-2 et Blast. Les performances
de Qrna étant meilleures sur les alignements produits par ClustalW, seuls les résultats de
Qrna obtenus avec cette méthode sont présentés.
Les résultats de Qrna et de Protea sur les couples de séquences des jeux de données
CODANT et NONCODANT sont reportés dans la table 2.2. Globalement, Qrna est plus
spécifique que Protea, tandis que Protea est plus sensible sur les mêmes données. Si l’on
considère le compromis entre sensibilité et spécificité, Protea se montre plus performant
que Qrna qui dispose de plus, rappelons le, d’un modèle pour détecter les séquences non-
codantes homologues qui l’avantage sur le jeu de données NONCODANT. De plus, Protea
est clairement plus performant que Qrna en dessous de 50% d’identité : 38% de sensibilité et
100% de spécificité pour Qrna, contre 81,3% et 95,5% respectivement pour Protea. Cette
observation permet de mettre en évidence les limites intrinsèques des méthodes basées sur
des alignements en présence de séquences divergentes. En terme de temps de calculs, Qrna
s’avère beaucoup plus gourmand que Protea notamment à cause de l’évaluation du modèle
non-codant très coûteuse. Pour réaliser cette expérience, il a fallu moins d’une heure à Protea
contre plus de 40 heures à Qrna.
51
Chapitre 2. Recherche de gènes et régions codantes
Tab. 2.1 – Les résultats de Protea sur les jeux de données CODANT, NONCODANT
et ALEATOIRE de 3, 5 et 11 séquences. Les résultats avec 11 séquences pour le jeu de
données NONCODANT ne sont pas mentionnés car très peu de familles de ce jeu de données
contiennent autant de séquences. Le tableau (a) contient la proportion de données correcte-
ment classifiées “codant”. Les tableaux (b) et (c) les proportions de données correctement
classifiées “autre” par Protea.
52
2.5. Résultats expérimentaux de Protea
Tab. 2.2 – Les résultats de Qrna et Protea sur les couples de séquences.
ensembles de séquences sont construits à partir d’alignements multiples générés par Mul-
tiz [BKR+ 04] puis filtrés par phastCons [SH05]. Dix huit génomes d’eucaryotes supérieurs
dont l’Homme, la souris, le chien et le poulet ont ainsi été comparés.
L’objectif de notre expérience est de découvrir de nouvelles séquences codantes chez
l’Homme, c’est pourquoi nous avons retiré de ce jeu de données toute ensemble contenant une
séquence déjà identifiée comme telle de manière expérimentale. A cet effet, nous avons utilisé
les ressources fournies par l’UCSC TableBrowser pour filtrer les séquences chevauchantes
ou incluses dans les pistes KnownGene [HKC+ 06] ou MGC (Mammalian Gene Collection).
Etant donné les performances de Protea sur les séquences courtes ou trop conservées, les
éléments de moins de cinquante nucléotides ou dont le pourcentage d’identité est supérieur à
90% ont été écartés. Au total, 97 956 ensembles d’au moins douze séquences ont été soumis à
Protea. Ce jeu de données peut être séparé en deux groupes : les séquences annotées dont la
fonction putative a déjà été prédite par d’autres méthodes bio-informatiques, et les séquences
non annotées. Parmi les 97 956 ensembles analysés, on compte 37 318 ensembles contenant
des séquences annotées et 60 638 sans annotation. Les résultats obtenus sont schématisés sur
la figure 2.8.
Une partie des séquences traitées comportent des annotations réalisées par des méthodes
de prédiction ab initio ou par analyse comparative. L’UCSC Table Browser nous a permis
de récupérer les annotations réalisées par Augustus, ExoniPhy, ExonWalk, GeneID,
Genscan et N-Scan. Nous avons également réalisé des annotations en utilisant une approche
classique par homologie de séquences au niveau peptidique grâce à BlastX sur la base
SwissProt en filtrant les alignements dont la E-valeur est inférieure à 10−4 .
Globalement, 23 220 ensembles de séquences annotées sont confirmés par Protea, soit
62% des séquences annotées par d’autres méthodes. La table 2.3 explicite la répartition des
prédictions de Protea en fonction des autres méthodes. Cette table contient également le
cœfficient de corrélation entre les prédictions de chaque méthode et celles de Protea. Les
valeurs de ce cœfficient comprises entre 0,1 et 0,26 sont relativement faibles. Cette expérience
montre que l’approche comparative de Protea est complémentaire des approches existantes,
notamment des approches ab initio et par analyse comparative.
53
Chapitre 2. Recherche de gènes et régions codantes
Familles
97 956
23 200 6 023
232 2 950
Fig. 2.8 – Découpage des résultats de Protea sur les groupes de séquences similaires de
l’UCSC.
Tab. 2.3 – Les résultats de Protea sur les éléments conservés comportant des annotations
putatives.
54
2.5. Résultats expérimentaux de Protea
2.5.3 Conclusions
Les expériences présentées dans la section 2.5.1 permettent d’apprécier en pratique les
forces et les faiblesses de Protea. Protea est une méthode efficace et performante, capable
de traiter des séquences non alignées très faiblement conservées. Les performances de Protea
sont, qui plus est, cohérentes avec le comportement attendu pour une méthode à base d’ana-
lyse comparative : ses performances croissent avec le nombre et la longueur des séquences
comparées. Au cours de ces expériences, nous avons noté que Protea n’était pas à l’aise sur
des séquences courtes ou très bien conservées. Ces faiblesses proviennent du principe même de
l’analyse comparative de séquences. Comparer des séquences quasiment identiques n’apporte
pas plus d’information qu’une seule séquence. Les séquences trop courtes ne contiennent pas
suffisamment d’information pour réaliser des observations significatives et pertinentes. L’appli-
cation de Protea à l’annotation de séquences codantes sur le génome humain (section 2.5.2)
a permis d’une part de mettre en évidence la complémentarité de Protea par rapport aux
méthodes de prédictions existantes, et d’autre part de détecter de nouveaux fragments de
séquences codantes putatives avec un degré de confiance élevé.
55
Chapitre 2. Recherche de gènes et régions codantes
56
Chapitre 3
Dans le chapitre précédent, nous avons abordé le problème de l’identification des régions
codantes. Nous nous intéressons maintenant aux ARN non-codants. Comme pour les ARN
codants, plusieurs types d’informations peuvent être pris en compte tels que l’homologie et les
biais de composition. Toutefois, la majorité des ARN non-codants présentent une particularité
supplémentaire qui est la formation d’une structure spatiale stable (section 1.3.1), que l’on
peut capturer partiellement à travers la structure secondaire. C’est un signal important qui
s’avère fort utile pour la prédiction de gènes à ARN. De ce fait, nous commençons par présenter
dans la première section les approches principales pour la prédiction de structures secondaires.
En section 3.2, nous expliquons ensuite comment ces méthodes s’appliquent à la prédiction
de gènes à ARN. Enfin, dans la section 3.3, nous présentons notre contribution au problème,
avec une évolution du logiciel caRNAc et des premiers résultats pour la prédiction de gènes.
57
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
séquences selon le schéma évolutif des ARN non-codants. Nous détaillons dans la suite de
cette section les méthodes thermodynamiques, comparatives et hybrides.
58
3.1. La prédiction de structures secondaires, état de l’art
les deux cas, le problème se résout par programmation dynamique. Le calcul de la struc-
ture d’énergie libre minimale se décompose ainsi en deux étapes : le remplissage de la table
de programmation dynamique afin de calculer l’énergie libre minimale atteignable, puis la
reconstruction d’une structure optimale par remontée dans la matrice.
Soit une séquence d’ARN s = s[1..n], c’est-à-dire un mot de longueur n sur l’alphabet
{A, C, G, U}. On définit une matrice carrée E de n × n cellules. Chaque cellule E(i, j) de la
matrice E correspond à l’énergie libre de la structure d’énergie libre minimale de la sous-
séquence s[i..j], avec i ≤ j, du mot s. Le remplissage de la matrice E s’effectue en suivant la
relation suivante, illustrée de manière schématique sur la figure 3.1
(
E(i + 1, j)
E(i, j) = min min {E(i + 1, k − 1) + E(k + 1, j) + α(i, k)}
i<k<j
(a) La base i est libre, l’énergie libre de la (b) La base i est appariée avec une certaine
séquence s[i..j] est alors celle de la séquence base k, l’énergie libre de la séquence s[i..j]
s[i + 1..j]. est alors la somme des énergies associées à
l’appariement de i avec k et aux séquences
s[i + 1..k − 1] et s[k + 1..j].
Par construction, l’énergie libre de la structure d’énergie libre minimale de s se trouve dans
la cellule E(1, n). Pour reconstruire une structure optimale associée à cette valeur d’énergie
libre, on remonte la matrice en partant de la cellule E(1, n) afin de retracer le chemin suivi
dans la matrice pour obtenir cette valeur. La complexité spatiale de l’algorithme est en O(n2 )
à cause du stockage de la matrice carrée E. Chaque cellule de la matrice nécessite un calcul
en temps linéairement proportionnel à longueur de la sous-séquence correspondante. Cet
algorithme a donc une complexité temporelle en O(n3 ).
Bien que cette modélisation soit fortement limitée, l’algorithme défini par Nussinov et
Jacobson est la base de la plupart des algorithmes de prédiction de structures d’ARN qui
visent à déterminer une structure d’énergie libre minimale.
L’algorithme de Zuker
La modélisation adoptée par l’algorithme de Nussinov et Jacobson ne prend pas en compte
de nombreux éléments qui contribuent à stabiliser une structure, comme les empilements
d’appariements, ou à la déstabiliser, comme les régions non appariées.
L’algorithme proposé par Zuker [ZS81, JTZ89, JTZ90] est une extension de l’algorithme
de Nussinov et Jacobson pour le modèle d’énergie plus réaliste de Freier-Turner [TSF88,
MSZT99, MDC+ 04] : empilements d’appariements pour former des tiges, boucles terminant
59
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
une tige, branchements multiples,. . . La figure 3.2 montre une partie des éléments pris en
compte dans ce modèle d’énergie.
Quatre matrices F, C, M etM 1 sont nécessaires au découpage d’une structure dans ce
modèle, comme illustré en figure 3.3. Pour une séquence s = s[1..n],
– F (i, j) correspond à l’énergie libre de la structure d’énergie libre minimale de la sous-
séquence s[i..j] ;
– C(i, j) correspond à l’énergie libre de la structure d’énergie libre minimale de la sous-
séquence s[i..j] où l’appariement entre les nucléotides aux positions i et j est forcé ;
– M (i, j) correspond à l’énergie libre de la structure d’énergie libre minimale de la sous-
séquence s[i..j] sachant que cette sous-séquence fait partie d’un embranchement multiple
comportant au moins une composante, c’est-à-dire une structure quelconque fermée un
appariement ;
– M 1 (i, j) correspond à l’énergie libre de la structure d’énergie libre minimale de la sous-
séquence s[i..j] sachant que cette sous-séquence fait partie d’un embranchement multiple
comportant exactement une composante.
Les récurrences établies par Zuker sont les suivantes.
(
F (i + 1, j)
F (i, j) = min min {C(i, k) + F (k + 1, j)}
i<k<j
H(i, j)
min {C(k, l) + I(i, j, k, l)}
C(i, j) = min i<k<l<j
min M (i + 1, u) + M 1 (u + 1, j − 1) + a
i<u<j
min {(u − i + 1).c + C(u + 1, j) + b}
i<u<j
M (i, j) = min min {M (i, u) + C(u + 1, j) + b}
i<u<j
M (i, j − 1) + c
M 1 (i, j − 1) + c
M 1 (i, j) = min
C(i, j) + b
où H(i, j) est l’énergie d’une boucle terminale fermée par l’appariement entre les
nucléotides en position i et j, I(i, j, k, l) est l’énergie d’une boucle interne formée des deux
sous-séquences s[i..j] et s[k..l] et où les variables a, b et c sont des constantes qui proviennent
du modèle d’énergie linéaire des embranchements multiples, à savoir que l’énergie d’un em-
branchement multiple de degré degree et de longueur size est E = a+b.degree+c.size. Le degré
d’un embranchement multiple est le nombre de sous-séquences non appariées qui séparent ses
composantes. La longueur d’un embranchement multiple est la somme des longueurs de ses
régions non appariées, comme illustré sur la figure 3.2.
L’algorithme de Zuker a une complexité temporelle en O(n4 ) et spatiale en O(n2 ). Histori-
quement, on compte deux implémentations strictes de l’algorithme de Zuker : Mfold [Zuk89]
et RNAfold [HFS+ 94]. Ces deux logiciels sont les plus utilisés pour la prédiction de struc-
tures secondaires. Toutefois, une complexité temporelle en O(n3 ) de l’algorithme a pu être
atteinte grâce à un traitement différent des boucles internes basé sur une fonction de coût
concave ou convexe [LZP99]. Des travaux récents de Roytberg et al ont permis d’améliorer
60
3.1. La prédiction de structures secondaires, état de l’art
5′ 3′ 5′ 3′
i j i j
i+1 j−1
5′ 3′ 5′ 3′
5′ 3′
i j i j
i j
5′ 3′
i j
i+k j − k′
5′ 3′
i i+k j − k′ j
5′ 3′
i j
5′ 3′
i j
61
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Source http://www.zbit.uni-tuebingen.de/pas/EMBO-RNACourse/handouts/HandoutBook.pdf
Fig. 3.3 – Schématisation des règles de décomposition appliquées dans l’algorithme de Zuker.
Les arcs pleins correspondent à des appariements entre des bases reliées. Les traits discontinus
indiquent des régions qui ne contiennent aucun appariement.
62
3.1. La prédiction de structures secondaires, état de l’art
Shape GGGCCCAUAGCUCAGUGGUAGAGUGCCUCCUUUGCAAGGAGGAUGCCCUGGGUUCGAAUCCCAGUGGGUCCA
[] (((((((((((((((.((((.....(((((((...))))))).))))))))))).........)))))))). -35.9 kcal/mol
[[][]] ((((((((.....((.((((.....(((((((...))))))).))))))(((.......))).)))))))). -32.2 kcal/mol
[[][][]] ((((((...((((.......)))).(((((((...))))))).....(((((.......))))).)))))). -31.7 kcal/mol
GGGCCCAUAGCUCAGUGGUAGAGUGCCUCCUUUGCAAGGAGGAUGCCCUGGGUUCGAAUCCCAGUGGGUCCA
(a) Structure primaire d’un ARN de transfert associé à l’alanine provenant du génome de Natronobac-
terium pharaonis (AB003409.1)
A A
GC
GC GC
GU
CG
GC
CG GU
CG CG
A UA C
UG C
C
CG
AU
CG
G UC
CG G AA
U
UGA UA U UAA
CUCG A C C C G
UG
CG
G GA
AU G GAGU C U G G GU C
UG C U A C U
GC G C
GC
UG C GAUG
G
A
AG U CG
U
GC A UA
G
CG CG
UA
CG
CG
CG UA
U U
U A A
A U C
U C
G G
(b) Structure secondaire d’énergie libre minimale (c) Structure secondaire sous-optimale prédite par
(−35.9 kcal/mol) prédite par RNAfold RNAfold (−31.7 kcal/mol) qui correspond à la
structure secondaire réelle de l’ARN de transfert
Fig. 3.4 – La structure secondaire de gauche correspond à la structure d’énergie libre minimale
prédite pour un ARN de transfert par l’algorithme de Zuker et le modèle d’énergie de Freier-
Turner. La structure de droite est la structure secondaire réelle des ARN de transfert.
L’algorithme de Zuker a été étendu par Eddy et Rivas [RE99] afin de permettre la
prédiction de structures tertiaires d’ARN incluant des pseudonœuds. Ces derniers ont ainsi
pu montrer que la prédiction de structures tertiaires est un problème dont la complexité
temporelle en O(n6 ) est quasiment impraticable sur des séquences de plus d’une centaine
de bases. Toutefois, en restreignant l’investigation à certaines classes de pseudonœuds, des
solutions algorithmiques exactes et efficaces ont été proposées telles que Pknots [RE99] et
Pknots-RG [RSG07]. Le modèle énergétique sous-jacent à la formation des pseudonœuds
reste néanmoins trop flou pour permettre d’établir des prédictions aussi pertinentes que pour
les structures secondaires.
Bon nombre d’approches se sont inspirées de l’algorithme de Zuker en utilisant d’autres
modèles énergétiques pour la prédiction de structures secondaires. CONTRAfold [DWB06]
est une méthode qui adopte un modèle d’énergie obtenu par apprentissage sur un ensemble
de séquences annotées par leur structure connue et vérifiée. Une méthode plus récente, MC-
fold [PM08], adopte un schéma différent du modèle d’énergie de Turner en utilisant une
base de motifs identifiés in silico sur des structures tertiaires d’ARN, les NCM, acronyme
pour Nucleotide Cyclic Motifs. Plusieurs structures jusqu’ici incorrectement prédites grâce
63
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
au modèle de Turner ont ainsi pu être prédites par ce logiciel. Cependant, la recherche de
ces motifs particuliers et leur assemblage pour former une structure secondaire a un coût
algorithmique non négligeable qui limite son utilisation à des séquences relativement courtes.
La fonction de partition
L’approche thermodynamique peut être abordée sous un autre angle avec la fonction de
partition. Le but n’est alors plus de minimiser l’énergie libre mais de maximiser la probabilité
d’une structure donnée d’ARN connaissant l’ensemble des structures que la séquence peut
adopter et la probabilité de formation d’un appariement dans ce contexte. Selon les prin-
cipes de la thermodynamique, la probabilité d’une structure Ψ dans un système équilibré est
proportionnelle à son facteur de Boltzmann
−E(Ψ)
exp
RT
où E(Ψ) est l’énergie libre de la structure Ψ, R est la constante du gaz parfait (en
Joules/(Kelvin mol)), T est la température absolue (en Kelvin). L’ensemble des structures
est déterminée par la fonction de partition notée Z. Cette fonction est une grandeur fonda-
mentale qui englobe les propriétés statistiques d’un système à l’équilibre thermodynamique.
Le système considéré ici étant l’ensemble des structures secondaires possibles, la fonction de
partition est définie par
X −E (Ψ)
Z= exp
RT
Ψ
Grâce à cette fonction, on peut déterminer la probabilité d’une structure Ψ dans l’ensemble
des structures possibles considérées :
exp −E(Ψ)
RT
p(Ψ) =
Z
Cette approche pour calculer la probabilité d’une structure nécessite de calculer la fonc-
tion de partition complète. Directement, ce calcul est impraticable car il demande de cal-
culer l’énergie libre de toutes les structures secondaires possibles dont le nombre croı̂t de
manière exponentielle en fonction de la longueur de la séquence [WFHS99]. Grâce aux tra-
vaux de McCaskill [McC90], la fonction de partition peut être calculée de manière partielle et
récursivement par programmation dynamique avec une complexité spatiale en O(n2 ) et tem-
porelle en O(n3 ). Pour expliquer l’idée mise en œuvre par McCaskill, on se place dans le cas
du modèle énergétique simple utilisé dans l’algorithme de Nussinov et Jacobson où l’énergie
d’un structure est obtenue en sommant les contribution des appariements individuels. Soit
Z(i, j) la fonction de partition pour toutes les structures de la séquence s[i..j]
X −α(i, j)
Z(i, j) = Z(i + 1, j) + Z(i + 1, k − 1)Z(k + 1, j) exp
RT
k
Cette formule peut être obtenue en transformant l’équation de récurrence utilisée dans
l’algorithme de Nussinov et Jacobson présentée à la page 58 en remplaçant les opérations
de minimisation par des sommes, les sommes par des multiplications, et les énergies par les
64
3.1. La prédiction de structures secondaires, état de l’art
65
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Fig. 3.5 – Alignement de sept séquences d’ARN de transfert, avec représentation de la struc-
ture commune. Les deux colonnes isolées c1 et c2 participent à un même appariement et font
apparaı̂tre des mutations compensatoires.
S(i, j, k, l) = min
S(i + 1, j, k, l)
S(i, j, k + 1, l)
min {S(i + 1, p − 1, k, l) + S(p + 1, j, 0, 0) + α(s1 [i], s1 [p])}
i<p<j
min {S(i + 1, p − 1, 0, 0) + S(p + 1, j, k, l) + α(s1 [i], s1 [p])}
i<p<j
min {S(i, j, k + 1, q − 1) + S(0, 0, q + 1, l) + α(s2 [k], s2 [q])}
k<q<l
min {S(0, 0, k + 1, q − 1) + S(i, j, q + 1, l) + α(s2 [k], s2 [q])}
k<q<l
min {S(i + 1, p − 1, k + 1, q − 1) + S(p + 1, j, q + 1, l) + α′ (s1 [i], s1 [p], s2 [q], s2 [k])}
i<p<j
k<q<l
où α(s1 [i], s1 [j]) correspond à la contribution énergétique apportée par l’appariement
entre s1 [i] et s1 [j] telle qu’elle est utilisée dans l’algorithme de Nussinov et Jacobson, et
66
3.1. La prédiction de structures secondaires, état de l’art
α′ (s1 [i], s1 [j], s2 [k], s2 [l]) correspond à la contribution énergétique apportée par l’appariement
conjoint entre s1 [i] et s1 [j] d’une part, et s2 [k] et s2 [l] d’autre part. La définition la plus na-
turelle pour α′ consiste à prendre la somme des contributions énergétiques des appariements
individuels
α′ (s1 [i], s1 [j], s2 [k], s2 [l]) = α(s1 [i], s1 [j]) + α(s2 [k], s2 [l])
On peut introduire un facteur bonifiant les appariements conjoints afin de favoriser les
corepliements en présence de mutations, particulièrement en présence de mutations compen-
satoires. Un malus peut également être considéré en cas d’introduction d’insertion ou de
délétion, c’est-à-dire dans les deux premières règles.
L’algorithme de Sankoff a une complexité temporelle en O(n3 m3 ) et une complexité
spatiale en O(n2 m2 ). Cet algorithme peut être étendu à plus de deux séquences. Pour N
séquences, sa complexité temporelle est alors en O(l3N ) et sa complexité spatiale en O(l2N ),
où l est la longueur de la plus longue des séquences traitées.
La complexité élevée de l’algorithme de Sankoff, même sur deux séquences, le rend im-
praticable sur des séquences qui dépassent la centaine de nucléotides. Il existe cependant de
nombreuses déclinaisons de cet algorithme qui tentent de traiter ce problème. Seule une partie
d’entre elles sont présentées, les plus originales. FoldAlign [HLG05, HTG07] implante une
version où les embranchements multiples sont interdits et seules les tiges terminées par une
boucle terminale sont considérées. La restriction appliquée dans Dynalign [MT02, HSM07]
est une borne maximale sur la distance qui sépare des nucléotides alignés, ce qui permet
de restreindre l’exploration de la matrice S à son hyperdiagonale. Dans Consan [DE06] et
Stemloc [Hol05], des régions fortement conservées sont identifiées entre les séquences pour
produire un alignement grossier des séquences, ce qui permet de contraindre la formation
d’appariements respectueux de cet alignement.
Toujours inspirées de l’algorithme de Sankoff, d’autres heuristiques s’attachent à la
prédiction de structures communes à plus de deux séquences. FoldAlignM [THG07] réalise
un alignement global de manière progressive à partir des alignements deux à deux générés
par FoldAlign pour produire une structure globalement conservée. Cette manière de passer
de deux à plus de séquences est fortement inspirée de PMcomp/PMmulti [HBS04], nouvel-
lement LocaRNA [WRH+ 07], également repris dans StrAl [DWMS06]. Dans ces logiciels,
les repliements deux à deux sont réalisés par comparaison des matrices de probabilité d’appa-
riement produite à l’aide de la fonction de partition. Murlet [KTKA07] réalise de manière
itérative une alignement multiple des séquences à partir des repliements deux à deux calculés
par l’algorithme de Sankoff restreint à la manière de Consan, c’est-à-dire en établissant un
alignement préliminaire entre les séquences. Cette idée d’alignement préliminaire est reprise
dans MxscaRNA [TTKA06, TKKA08] où cette fois seules les parties ouvrantes et fermantes
des tiges sont alignées. L’alignement par morceaux ainsi obtenu est ensuite utilisé pour inférer
une structure globale.
67
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
La valeur de Mij varie entre 0 et 2 bits et mesure le degré de corrélation des deux co-
lonnes. Elle est maximale lorsque les deux colonnes sont parfaitement corrélées, c’est-à-dire
qu’un appariement est totalement absent ou conservé, et que leur contenu individuel est
pourtant totalement aléatoire, c’est-à-dire que toutes les nucléotides apparaissent de manière
équiprobable. Mij est nulle en l’absence de mutation dans les deux colonnes ou lorsque les
colonnes varient de façon indépendante, c’est-à-dire fij (x, y) = fi (x)fj (y). Des corrections
peuvent être apportées à cette mesure pour prendre en compte la composition globale des
séquences ou encore un arbre phylogénétique pour prendre en compte le taux de mutations
attendues par colonne [KH99, GHH+ 94]. Sur l’exemple de la figure 3.5, l’information mu-
tuelle des deux colonnes appariées c1 et c2 qui présentent des mutations compensatoires est
de 1,37 bit, alors qu’entre les colonnes non appariées m1 et m2 cette valeur est de 0,52 bit.
L’information mutuelle constitue la noyau de bon nombre de méthodes de prédiction de
structure à partir d’un alignement. La plupart de ces méthodes utilisent cette information
comme score soit à la place de l’information énergétique pour les méthodes à base de l’algo-
rithme de Zuker, soit comme bonus ou malus complémentaire à une approche énergétique.
RNAalifold [HFS02] est actuellement la méthode plus utilisée pour la prédiction de
structure secondaire à partir d’un alignement, car cette méthode est celle qui exploite le
modèle énergétique de Turner. Son fonctionnement repose sur l’algorithme de Zuker généralisé
à un alignement. La contribution énergétique d’un appariement entre deux colonnes est sim-
plement calculée en moyennant les contributions individuelles. Un bonus est appliqué en fonc-
tion de la corrélation entre les colonnes appariées. Ilm [RSZ04a, RSZ04b] est une méthode
strictement analogue à RNAalifold où l’algorithme de repliement adapté à l’alignement
multiple est celui de Nussinov et Jacobson. Plus récemment, RNAlishapes [Vos06] est en
fait une variante de RNAalifold où l’algorithme de repliement est une version modifiée
de RNAshapes : le repliement est effectué entre des représentations abstraites des struc-
tures optimales et sous-optimales prédites individuellement. L’algorithme travaille ainsi au
niveau des tiges en tentant de faire correspondre les parties ouvrantes d’une même tige
conservée d’une part, et les parties fermantes correspondantes d’autre part. Cove [ED94],
Pfold [KH99, KH03] et CMFinder [YWR06] sont trois méthodes à base de grammaire sto-
chastiques hors contexte entraı̂nées sur des séquences exemples et où l’information mutuelle
mesurée entre les colonnes est utilisée comme pondération des informations apprises. Pfold
présente toutefois une originalité supplémentaire : la mesure de l’information mutuelle peut
68
3.1. La prédiction de structures secondaires, état de l’art
être affinée pour tenir compte de la distance évolutive qui sépare les espèces dont sont issues
les séquences.
P-DCfold [TGR02, TER03, TER05, Eng06, ET07] propose une alternative originale
aux méthodes citées précédemment. Le logiciel commence par chercher des séquences palin-
dromiques conservées et alignées qui exhibent une ou plusieurs mutations compensatoires.
Puis en appliquant une heuristique gloutonne, il construit successivement des ensembles de
palindromes tous compatibles entre eux, c’est-à-dire sans croisement ni chevauchement. Cette
approche de type “diviser pour régner” permet d’une part de réduire de manière drastique la
complexité du repliement et d’autre part d’autoriser la formation de pseudonœuds.
TP
Sp =
T P + (F P − ξ)
69
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
T P × F N − (F P − ξ) × F N
M CC = p
(T P + (F P − ξ))(T P + F N )(T N + (F P − ξ))(T N + F N )
70
3.2. La prédiction de gènes à ARN
(a) Résultats des approches qui suivent le para- (b) Résultats des approches qui suivent le para-
digme “Aligner puis replier” comparés aux ap- digme “Aligner et replier simultanément” com-
proches thermodynamiques. parés aux approches thermodynamiques.
banques de données se sont développées. Par la suite, les approches comparatives ont fait leur
apparition, intégrant à la fois des informations de similarité et des informations intrinsèques
aux séquences.
La prédiction de gènes à ARN est un problème plus complexe que la prédiction de gènes
codants pour plusieurs raisons. On ne dispose pas des signaux forts présents dans les gènes
codants : l’existence d’un cadre ouvert de lecture et les biais dans l’usage des codons. De plus,
contrairement aux gènes codants, la production de certains ARN non-codants ne suit pas
le schéma classique “un gène pour une molécule” : certains ARN non-codants sont localisés
dans les introns d’autres gènes ou encore dans les régions non traduites des ARN messagers.
A cause de ce type d’ARN, il devient plus compliqué d’exploiter les signaux qui balisent
traditionnellement les gènes. Enfin, les propriétés à l’origine de la fonction des ARN varient
d’une famille d’ARN à une autre : certaines familles sont caractérisées par la seule conserva-
tion d’un motif de séquence, d’autres par une structure commune. Pour toutes ces raisons,
des méthodes de prédiction d’ARN non-codants ad hoc ont été développées, c’est-à-dire des
méthodes qui ciblent une seule famille d’ARN non-codants. L’idée est alors de rechercher des
séquences ou de vérifier si des séquences respectent un ensemble de contraintes qui décrivent
les propriétés conservées au sein d’une famille particulière.
L’organisation de cet état de l’art est calquée sur celui des méthodes de prédiction de
séquences codantes du chapitre 2. Nous envisageons dans un premier temps les approches ab
initio, c’est-à-dire l’exploitation de signaux exclusivement présents dans les séquences d’ARN
non-codants. Dans la section 3.2.1, nous nous intéressons donc à l’analyse de différents biais
de composition des séquences d’ARN non-codants liés à la formation d’une structure. Cette
analyse nous conduit naturellement vers l’analyse de la stabilité des structures d’ARN non-
codants présentée dans la section 3.2.2. Ensuite, nous nous intéressons aux méthodes de
prédiction d’ARN non-codants par homologie, c’est-à-dire la recherche de séquences homo-
logues dans les banques de données. Dans la section 3.2.3, deux types de similarité sont ainsi
envisagées : au niveau nucléique et au niveau structural. Enfin, nous clôturons cet état de
l’art par les approches comparatives de prédiction d’ARN non-codants (section 3.2.4).
71
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Tab. 3.2 – Résultats des mesures effectuées dans [Sch02]. (G+C)% correspond à la moyenne du
pourcentage en G et en C observé. ρ(CG) correspond à la moyenne de la fréquence normalisée
du di-nucléotide CG. Les valeurs entre parenthèses sont les écarts-types associés.
72
3.2. La prédiction de gènes à ARN
73
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
d’une région à former une structure complète fonctionnelle. Les programmes de prédiction
de structures qui adoptent une approche thermodynamique (section 3.1.1) sont conçus pour
fournir la structure d’énergie libre minimale qui peut se former à partir d’une séquence.
Quelque soit la séquence choisie, ces programmes prédisent toujours une structure. Seule,
l’existence d’une structure prédite n’est donc pas informative. C’est pourquoi il est nécessaire
de s’intéresser plus précisément à la qualité, en terme de stabilité thermodynamique, des
structures prédites.
Pour qu’une structure se forme, de nombreux appariements se font puis se défont jusqu’à
ce qu’un état stable soit atteint. On peut donc s’attendre à ce que les structures des ARN non-
codants soient remarquablement stables et caractérisées par une énergie libre particulièrement
faible. Evaluer la significativité de la stabilité d’une structure nécessite de disposer d’une
distribution de l’énergie libre avec laquelle effectuer la comparaison. Il n’existe cependant
aucune théorie pour la construire, elle est donc établie de manière empirique grâce à de
nombreuses séquences aléatoires équivalentes. Le protocole suivi pour évaluer la significativité
de la stabilité d’une structure est le suivant. A partir d’une séquence s,
1. calculer E, l’énergie libre minimale de la structure optimale prédite pour s ;
2. construire la distribution de l’énergie libre, c’est-à-dire inférer les structures d’un grand
nombre de séquences aléatoires obtenues par mélange de s ou en utilisant un processus
Markovien, de telle sorte que la composition en mono-nucléotides et/ou en di-nucléotides
de s soit conservée ;
3. évaluer la significativité de E à partir de la distribution obtenue grâce au z-score ou à
la P-valeur de E ; ces mesures sont équivalentes en terme d’information apportée.
Le z-score de E mesure l’écart de E par rapport à la distribution. On l’obtient en calculant
le rapport
E−µ
σ
où µ et σ sont respectivement la moyenne et l’écart-type de la distribution. Comme l’énergie
libre est à valeur négative, plus le z-score de E est faible, plus la structure optimale de s est
stable. La seconde mesure que l’on utilise est la P-valeur de E, qui correspond à la probabilité
d’obtenir une valeur d’énergie libre inférieure ou égale à E dans la distribution. Plus la P-
valeur de E est proche de 0, plus le nombre de structures prédites ayant une énergie libre
inférieure est faible. Par conséquent, plus la P-valeur de E est faible, plus la stabilité de la
structure optimale de s est significative.
Toute la difficulté dans ce processus d’évaluation réside dans la constitution de la distribu-
tion, et donc dans le choix de la composition des séquences. La conservation de la composition
mono-nucléotidique de s permet de tenir compte d’un éventuel biais de composition en G et
en C de s dû à l’existence d’une structure. La conservation de la composition di-nucléotidique
de s a une propriété supplémentaire : prendre en considération la formation éventuelle d’em-
pilements d’appariements.
74
3.2. La prédiction de gènes à ARN
Fig. 3.7 – Distribution de la négation des z-scores de l’énergie libre des structures de 243
ARN non-codants par rapport à des structures optimales de séquences aléatoires de même
composition en mono-nucléotides.
L’observation réalisée sur les ARN de transfert n’est pas valable pour toutes les familles
d’ARN non-codants : en moyenne, les structures d’ARN non-codants sont plus stables que
les structures optimales de séquences aléatoires équivalentes. Cependant, la stabilité moyenne
des ARN non-codants n’est pas assez significative pour constituer un signal suffisant pour
les détecter lorsque la distribution de référence est construite avec des séquences de même
longueur et de même composition en mono-nucléotides. Cette étude a tout de même donné lieu
au développement d’un logiciel, ncRNAScan, qui réalise des prédictions d’ARN non-codants
structurés selon ce procédé.
Les expériences précédentes ont été reprises dans [BWRVdP04] en considérant des
séquences de même composition en di-nucléotides pour 500 ARN de transfert, 581 ARN
ribosomiques et 506 micro ARN. Les résultats révèlent que les micro ARN possèdent
systématiquement des structures plus stables que des structures de séquences aléatoires de
même composition en di-nucléotides. Les structures des ARN ribosomiques et des ARN de
transfert ne sont pas systématiquement plus stables, mais en moyenne plus stables.
75
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Ces investigations ont été étendues à 300 familles d’ARN non-codants [CFKK05]. Cette
étude ouvre des perspectives intéressantes quant à la conservation de la composition di-
nucléotidique pour mesurer la significativité de la stabilité des structures d’ARN non-codants.
Le tableau 3.3 reprend une partie de leurs résultats.
Tab. 3.3 – Extraits des résultats de Clote et al [CFKK05]. Les z-scores et P-valeurs sont ceux
de l’énergie libre des structures.
Les structures des ARN non-codants sont en moyenne plus stables que ce qui est at-
tendu par hasard. Pour certaines familles, comme les ARN de transfert et les petits ARN
nucléolaires U1, les résultats sont plus modérés : la stabilité moyenne des structures n’est pas
aussi significative que pour les autres familles d’ARN. A l’instar de l’étude de Rivas et al,
cette étude a donné lieu au développement d’un logiciel, RANDfold, qui évalue la stabilité
thermodynamique d’une structure par rapport à une distribution d’énergie libre construite à
partir de séquences aléatoires de même longueur et composition en di-nucléotides.
76
3.2. La prédiction de gènes à ARN
pour les différentes sous-familles d’ARN ribosomiques. De manière générale, une similarité
significative entre des régions non codantes peut également constituer une information perti-
nente sur la présence de séquences soumises à une contrainte fonctionnelle.
77
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
matrices RIBOSUM. HomoStRscan [LMZ04] est une méthode en tout point équivalente
à Rsearch adaptée au traitement du génome humain : les matrices de substitutions sont
construites à partir de séquences d’ARN non-codants provenant du génome humain.
Rsearch et FastR ont la propriété d’être génériques au regard des structures recherchées
grâce à l’utilisation des matrices RIBOSUM. Toutefois, cette approche générale ne permet pas
de prendre en compte certaines caractéristiques propres à un ensemble de séquences struc-
turées homologues, telle que la conservation locale de motifs de séquence ou de structure
(section 1.3.1). Si l’on souhaite cibler les recherches à un ensemble de séquences homologues
particulier, les modèles de covariance sont alors plus appropriés [ED94, DEKM99]. Sans entrer
dans les détails, un modèle de covariance est une grammaire stochastique hors contexte pro-
filée par plusieurs matrices de substitutions de nucléotides et d’appariements. A partir d’un
alignement multiple de séquences homologues annoté par la structure commune partagée par
les séquences, on construit une matrice de substitutions d’appariements pour chaque couple
de colonnes appariées et une matrice de substitutions de nucléotides pour chaque colonne
non appariée. Ainsi, un modèle de covariance capture à la fois des informations locales sur la
structure, mais également sur la séquence dans les régions appariées et non appariées.
Plusieurs méthodes s’appuient sur les modèles de covariance, notamment Infer-
nal [GJBM+ 03] utilisée pour maintenir la banque Rfam. Un modèle de covariance est en effet
disponible pour chaque famille de Rfam, construit à partir d’un alignement multiple vérifié
manuellement et annoté par la structure conservée. Infernal est une implémentation stricte
et rigoureuse des modèles de covariance, c’est-à-dire qu’aucun traitement ne précède l’aligne-
ment. L’alignement à l’aide d’un modèle de covariance est une tâche plus gourmande en temps
de calculs et en espace mémoire qu’un alignement basé sur une grammaire “générique”, telle
que celle utilisée dans Rsearch. L’exploitation des modèles de covariance est particulièrement
gourmande en ressources, d’autant plus que la taille d’un modèle de covariance, c’est-à-dire
le nombre de règles dans la grammaire associée, est proportionnelle à la longueur de l’aligne-
ment multiple ayant servi à sa construction et au nombre d’appariements impliqués dans la
structure commune.
Plusieurs manières de contourner ce problème ont été envisagées. Dans RNAcad [Bro99],
un modèle de Markov caché sert à détecter des régions non appariées précises de la séquence
à aligner. Ce marquage permet d’imposer des contraintes lors de l’évaluation du modèle de co-
variance, sous forme de points de passage forcés dans la dérivation des règles de la grammaire
associée au modèle. Cette solution s’avère extrêmement efficace, bien qu’elle nécessite beau-
coup de données supplémentaires pour entraı̂ner le modèle de Markov. Dans le même esprit,
RaveNnA [WR04, WR06] construit un modèle de Markov particulier, un profile HMM, à par-
tir d’un modèle de covariance. Un profile HMM est un modèle de Markov caché adapté pour
modéliser un alignement multiple et non une simple séquence comme les modèles de Markov
cachés classiques. Comme les modèles de Markov classiques, ce type de modèle ne permet pas
de décrire la formation d’appariements entre nucléotides. Dans RaveNnA, le profile HMM est
construit à l’aide d’une heuristique basée sur le principe du maximum de vraisemblance pour
approximer au mieux le modèle de covariance original. Le filtrage à l’aide du profile HMM
permet de cibler les régions dont l’évaluation à l’aide du modèle de covariance est susceptible
de donner un score significatif. En moyenne, 90% des alignements trouvés avec un modèle de
covariance peuvent ainsi être retrouvés avec le profile HMM correspondant dans RaveNnA,
et ce environ 600 fois plus rapidement. Jusqu’ici nous n’avons qu’une seule facette des modèles
de covariance : construire un alignement structure/séquence. CMFinder [YWR06] propose
78
3.2. La prédiction de gènes à ARN
d’utiliser ce type de modèles pour tenter d’améliorer un alignement classique réalisé unique-
ment sur la séquence primaire. A partir d’un ensemble d’alignements classiques, CMFinder
sélectionne des bons candidats potentiels sur la base de critères énergétiques, puis affinent les
alignements retenus grâce au modèle de covariance selon le principe de maximum de vraisem-
blance. Dans sa première version, Infernal était relativement gourmand en ressources car il
ne procédait à aucune optimisation spécifique ou filtrage particulier en amont. Récemment, In-
fernal a évolué [NE07] et intègre maintenant un filtre qui permet de déterminer de manière
exacte les dérivations de la grammaire qui n’aboutiront pas à un alignement significatif étant
donné un modèle de covariance. Cette optimisation permet de diminuer de manière dras-
tique les ressources physiques et temporelles nécessaires à Infernal, et par extension, aux
approches à base de modèles de covariance.
Comme nous venons de le voir, le modèle de covariance est un outil puissant pour
modéliser des séquences homologues partageant une structure commune. L’exploitation de
cette modélisation s’avère cependant généralement coûteuse sans un filtrage approprié de
l’espace des dérivations à explorer. Erpin [GL01, LFL+ 04, LLFG05] propose une version
simplifiée des modèles de covariance : une matrice de score pour la substitution d’apparie-
ments est calculée pour chaque tige, et non pour chaque appariement, et une matrice pour
la substitution de nucléotides non appariées pour chaque sous-séquence non appariée. A la
manière de FastR, un ensemble de contraintes est associé à chaque tige, notamment l’inter-
valle des distances autorisées entre la partie ouvrante et la partie fermante. Les tiges sont
détectées de manière indépendante, puis assemblées pour former un alignement par program-
mation dynamique : il n’y a donc pas de notion de structure globale dans l’algorithme, ce qui
réduit drastiquement sa complexité. Le facteur déterminant de la complexité de Erpin est la
distance maximale autorisée qui sépare les parties ouvrante et fermante des tiges.
79
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
et des boucles internes symétriques ou non, et peuvent former des pseudo-œuds. Les mots sont
systématiquement décrits par une séquence, contrairement aux espaceurs et aux tiges dont la
description du contenu en nucléotides est facultative. Quelque soit le type d’élément décrit, le
contenu en nucléotides peut être approximatif, c’est-à-dire que les erreurs sont tolérées dans
une certaine mesure variable d’un élément à un autre. De plus, pour les tiges il est possible
d’autoriser certaines liaisons bancales et d’en préciser la quantité maximale autorisée. Pour
trouver les séquences qui satisfont un descripteur, RNAMot ordonne les éléments du des-
cripteur en fonction de leur probabilité marginale d’apparition estimée de manière empirique,
puis tente de placer les éléments récursivement, du moins probable au plus probable. Lorsque
deux occurrences chevauchantes satisfont un descripteur, RNAMot calcule un score qui per-
met de déterminer la meilleure occurrence, c’est-à-dire celle qui remplie au mieux les éléments
décrits : les tiges plus longues sont favorisées, les mésappariements et les erreurs défavorisées,
...
Inspiré de RNAMot, Palingol [BKV96] reprend les mêmes types d’éléments descriptifs.
Toutefois, Palingol est beaucoup plus expressif car il offre la possibilité de définir des expres-
sions logiques et de construire des branchements. Par exemple, il devient possible de décrire :
« si cette tige n’est pas présente ou qu’elle contient plus d’un mésappariement, alors appliquer
une pénalité ». Le programme s’avère également plus souple dans la gestion des mots ; un mot
peut être décrit par une matrice poids-position, comme celle utilisée pour modéliser les sites
d’épissage (section 2.1.2). Pour trouver les séquences qui satisfont un descripteur, l’algorithme
de Palingol comporte une phase préliminaire d’indexation de toutes les tiges présentent dans
la séquence afin de ne pas rejouer plusieurs fois inutilement des calculs coûteux.
PatScan [DLO97] offrent les mêmes possibilités que Palingol, bien que la gestion des
erreurs diffère : les substitutions, les insertions et les délétions sont considérées comme des
erreurs différentes, sans aucun moyen simple de les banaliser. Cette caractéristique peut rendre
l’élaboration du descripteur particulièrement difficile puisque qu’il faut alors explicitement
écrire les expressions conditionnelles qui simulent ce comportement. PatSearch [PLD00,
GLL+ 03] est le descendant de PatScan. Outre certains aspects syntaxiques différents dans la
définition des descripteurs, PatSearch estime le nombre d’occurrences attendues par hasard
avec un descripteur. Cette valeur est obtenue par simulation, en exécutant PatSearch avec
le même descripteur sur des séquences aléatoires de même composition en mono-nucléotides.
L’espérance du nombre d’occurrences attendues donne une indication sur la “qualité” d’un
descripteur afin de relativiser le nombre d’occurrences réellement observées.
RNAMotif [MEG+ 01], le descendant de RNAMot, est actuellement le logiciel le plus
utilisé. La principale nouveauté apportée dans RNAMotif est la totale liberté laissée à
l’utilisateur pour définir son propre système de score. La définition de ce système peut très
sophistiquée et faire appel à des structures algorithmiques relativement évoluées, telles que
des expressions conditionnelles et des boucles.
MilPat [TdGSG06] est une méthode récente qui se distingue des autres à deux points de
vue : la possibilité de décrire des interactions inter-séquences, et l’algorithme d’énumération
des occurrences. La fonction de bon nombre d’ARN non-codants implique une interaction
avec une autre séquence nucléique (section 1.3) qui se traduit par la formation d’apparie-
ments. MilPat permet de décrire ce type d’interactions comme la formation d’une “tige”
dont une partie se trouve sur la séquence cible et l’autre partie sur une séquence fournie
en plus du descripteur. L’autre originalité de MilPat est l’algorithme d’énumération des
occurrences. Le problème de trouver les occurrences qui satisfont un descripteur est ici for-
80
3.2. La prédiction de gènes à ARN
81
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Tab. 3.4 – Résultats de BRAliBase III sur les ensembles de cinq séquences. La sensibilité
et la spécificité sont exprimées en pourcentage.
Tab. 3.5 – Résultats de BRAliBase III sur les ensembles de vingt séquences. La sensibilité
et la spécificité sont exprimées en pourcentage.
82
3.2. La prédiction de gènes à ARN
les alignements de cinq séquences. A l’opposé des ARN de transfert, les petits ARN U5 sont
plutôt bien prédits par les deux types d’approches, à l’exception de Erpin dont la sensibilité ne
dépasse pas 41%. Dans les faits, les ARN U5 comportent des sites de fixation très conservés
qui, à eux seuls, constituent un signal suffisant pour les prédire. La figure 3.8 montre la
structure d’un ARN U5 où les bases sont colorées en fonction de leur degré de conservation
au sein de la famille. Trois sites sont particulièrement conservés : la boucle terminale et une
des boucles internes de la tige en 3′ ainsi que la multiboucle.
Fig. 3.8 – Structure d’un petit ARN U5 composée de deux tiges juxtaposées. Les couleurs
indiquent le degré conservation de chaque base au sein des séquences connues de la famille.
Le nombre de séquences utilisées pour réaliser les prédictions influe particulièrement sur le
comportement des méthodes à base d’homologie de séquences : leur sensibilité double lorsque
l’on passe de cinq à vingt séquences. Bien qu’on ne dispose pas des résultats en fonction du
pourcentage d’identité, les auteurs de BRAliBase III mentionnent que la faible sensibilité
des méthodes par pure homologie de séquences provient de la divergence entre les séquences
requêtes et la séquence ciblée. En pratique, ils notent une nette dégradation de la sensibilité
de ce type d’approche lorsque l’ensemble des séquences utilisé ne comporte aucune séquence
dont le pourcentage d’identité avec la séquence à prédire est supérieur à 65%.
Les résultats en terme d’efficacité sont donnés dans la table 3.6. Sans surprise, les méthodes
les plus rapides sont celles qui n’intègrent pas ou peu d’information sur la structure à savoir
Blast, Fasta et Erpin. Si l’on omet les temps d’initialisation, les méthodes les plus per-
formantes Infernal et Rsearch sont également les plus lentes. A titre de comparaison, In-
fernal et Rsearch traitent en moyenne entre mille et deux milles fois moins de nucléotides
par seconde que Erpin, Blast et Fasta. Bien que négligeable pour évaluer des banques de
données conséquentes, le temps d’initialisation qui précède la phase de recherche des méthodes
83
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
les plus sophistiquées comme RaveNnA, Infernal et Rsearch est relativement élevé.
84
3.2. La prédiction de gènes à ARN
les probabilités du modèle, EvoFold s’appuie sur un arbre phylogénétique contenant les
distances évolutives qui sépare les organismes dont sont extraites les séquences alignées.
Qrna et EvoFold procèdent à une analyse comparative plutôt fine des mutations entre
les séquences. MSARi [CKB04] adopte la même démarche à base de modèles probabilistes
mais de manière heuristique. MSARi procède entre trois temps. Premièrement, les probabi-
lités d’appariement de tous les couples de nucléotides sont calculées individuellement pour
chaque séquence grâce à la fonction de partition (section 3.1.1). A partir des résultats obtenus,
MSARi recherche des tiges conservées grossièrement cohérentes avec l’alignement multiple,
c’est-à-dire que les tiges mises en correspondance ne respectent pas nécessairement stricte-
ment l’alignement multiple. MSARi suppose en effet que l’alignement peut contenir quelques
erreurs et que les bases appariées ne sont pas nécessairement correctement alignées. Enfin,
MSARi sélectionne les tiges conservées de manière gloutonne, par nombre de mutations com-
pensatoires décroissant, pour former une structure secondaire commune. La classification est
enfin réalisée en fonction de la significativité de la structure obtenue, évaluée en fonction du
nombre de mutations compensatoires globalement observées.
Contrairement aux approches précédentes, une autre classe de méthodes s’attache à la sta-
bilité thermodynamique d’une éventuelle structure secondaire commune. Cette approche suit
le même schéma que celle présentée dans la section 3.2.2 où l’on évalue la stabilité d’une struc-
ture d’énergie minimale prédite sur une seule séquence. La difficulté supplémentaire ici est de
construire une distribution de l’énergie libre de structures obtenues non plus sur des séquences
individuelles de composition équivalente, mais sur des ensembles de séquences alignées ou non.
AlifoldZ [WH04] et ddbRNA [DBDH03] procèdent ainsi à l’évaluation de la stabilité d’une
structure commune par rapport à une distribution d’énergie libre construite à partir d’aligne-
ments générés en mélangeant les positions de l’alignement multiple original. Par construction,
un alignement obtenu par mélange de ses positions respecte deux propriétés : la composition
en mono-nucléotides de chaque séquence est préservée et la conservation globale des séquences.
Dans ddbRNA, la structure commune est calculée en assemblant de manière gloutonne des
tiges conservées, et la procédure de mélange s’efforce en plus de détruire au moins partiel-
lement cette structure commune. Dans AlifoldZ en revanche, la prédiction de la structure
commune est déléguée à RNAalifold, et la procédure de mélange est plus complexe car elle
respecte le degré de conservation locale, c’est-à-dire que le degré de conservation de chaque
position est préservé entre tous les couples de séquences. Cette propriété est particulièrement
importante pour préserver les longues insertions/délétions qui pourraient alors être éclatées
en plusieurs petites régions et empêcher la formation de tiges lors du repliement commun.
La figure 3.9 montre la distribution du z-score de l’énergie libre de la structure commune
prédite par RNAalifold sur des alignements d’ARN de transfert comportant de une à quatre
séquences. La distribution de l’énergie libre tracée en trait plein, est comparée à la distribu-
tion de l’énergie libre de la structure commune prédite par RNAalifold à partir du même
alignement mélangé par la procédure de AlifoldZ. D’après ces résultats, la structure secon-
daire commune à plusieurs séquences semble significativement plus stable qu’une structure
d’ARN non-codant prédite à partir d’une seule séquence. En fait, plus l’alignement comporte
de séquences, plus l’énergie libre de la structure commune prédite est faible. RNAalifold
intègre en effet dans son calcul de l’énergie libre de la structure commune un bonus négatif
pour chaque covariation. Le nombre de covariations observables augmente avec le nombre
de séquences homologues différentes, par conséquent l’énergie libre de la structure commune
prédite diminue. En revanche, on ne s’attend pas à trouver de covariations sur les alignements
85
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
mélangés, quelque soit le nombre de séquences de l’alignement. En fixant un seuil sur la valeur
du z-score de l’énergie de la structure prédite à −4, il devient alors possible de distinguer clai-
rement les alignements d’ARN de transfert des alignements de séquences aléatoires. Plus le
nombre de séquences alignées est élevé, plus cette classification s’avère efficace : pour quatre
séquences, 98,36% des alignements d’ARN de transfert peuvent ainsi être discriminés sans
prédire à tort un alignement de séquences aléatoires.
Fig. 3.9 – Distribution du z-score de l’énergie libre de la structure commune prédite par
RNAalifold sur des alignements de plusieurs familles d’ARN non-codants évaluée par rap-
port à une distribution empirique de l’énergie libre des structures prédites par RNAali-
fold sur des alignements générés par mélange des séquences originales selon la procédure
d’AlifoldZ. N est le nombre de séquences présentent dans les alignements. Pour N = 1, les
structures sont prédites par RNAfold.
RNAz [WHS05] est une amélioration de AlifoldZ qui intègre une mesure supplémentaire
de la stabilité de la structure commune : le SCI (Structure Conservation Index). Le SCI évalue
la stabilité de la structure commune par rapport aux structures prédites individuellement. Il
s’obtient en calculant le rapport EA /Ē, où EA est l’énergie libre de la structure commune
prédite par RNAalifold, et Ē est la moyenne de l’énergie libre des structures individuelles
prédites par RNAfold. Lorsque le SCI est proche de 0, la structure trouvée par RNAalifold
a une énergie libre plus faible que la moyenne de l’énergie libre des structures individuelles : la
structure trouvée pour l’alignement n’est pas significative ; les structures sont mal conservées.
Un SCI proche de 1 indique au contraire que les structures sont parfaitement conservées. Un
SCI plus grand que 1 indique non seulement que les structures sont parfaitement conservées,
86
3.3. Evolution et enrichissement du logiciel caRNAc
mais qu’il existe en plus des mutations compensatoires. Afin d’éviter la construction em-
pirique coûteuse d’une distribution d’énergie libre, les paramètres de cette distribution sont
approximés au moyen d’un processus d’apprentissage supervisé, les SVM (Support Vector Ma-
chine). Ce même type de processus est utilisé pour effectuer la classification de l’alignement
en fonction du SCI et du z-score de l’énergie libre de la structure commune. Actuellement,
RNAz est la méthode la plus utilisée pour la prédiction d’ARN non-codants.
Comme nous l’avons déjà remarqué dans la section 3.2.2, évaluer la stabilité d’une struc-
ture par rapport à une distribution d’énergie libre établie à partir de séquences de même
composition en di-nucléotides apporte de meilleurs résultats qu’en ne préservant que la compo-
sition en mono-nucléotides. Récemment, Sissiz [GW08] reprend le protocole employé jusqu’ici
mais avec une procédure de génération d’alignements multiples qui préserve une composition
en di-nucléotides donnée en plus de toutes les propriétés déjà évoquées, notamment la conser-
vation locale. La distribution d’énergie libre obtenue à partir des alignements générés par
cette procédure est selon les auteurs plus proche de la distribution réelle.
3.3.1 L’existant
caRNAc admet en entrée n séquences d’ARN non alignées et produit pour chaque
séquence un structure secondaire sous forme d’une liste de tiges conservées. Cela se fait en
deux temps. La première phase de l’algorithme consiste à procéder à tous les repliements deux
à deux suivant une adaptation de l’algorithme de Sankoff. Ensuite, ces prédictions sont com-
binées à l’aide d’une structure de graphe pour obtenir une structure secondaire pour chaque
séquence.
87
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
2
tiges
potentielles
paires 3 paires de structure
3
de tiges secondaire
motifs tiges compatibles 4 commune
conservés
tiges
potentielles
2
covariations
1 cohérence
séquence B des motifs
conservés
Fig. 3.10 – Déroulement de la prédiction d’une structure secondaire conservée entre deux
séquences dans caRNAc
La recherche de tiges Les tiges potentielles énumérées lors de cette étape contiennent
aux moins trois appariements canoniques consécutifs, peuvent contenir des mésappariements
et sont systématiquement fermées par un appariement canonique A-U, C-G ou G-U. Ces tiges
sont dites maximales car elles ne peuvent être étendues pour obtenir des tiges d’énergie libre
inférieure selon ces règles. L’énergie associée à une tige est calculée en utilisant le modèle de
Turner restreint aux empilements d’appariements, aux motifs des boucles terminales et aux
mésappariements symétriques. Comme les tiges sont prédites indépendamment les unes des
autres, les règles relatives aux boucles internes et aux embranchements ne peuvent pas être
appliquées à cette étape. De même, seules les boucles internes d’une longueur inférieure à huit
88
3.3. Evolution et enrichissement du logiciel caRNAc
nucléotides sont évaluées car on peut déjà affirmer à cette étape qu’aucune tige ne pourra être
prédite dans cette région non appariée. Toutes les tiges sont énumérées par programmation
dynamique avec une complexité spatiale et temporelle quadratique par rapport à la longueur
de la séquence, puis filtrées selon leur valeur d’énergie libre grâce à une fonction de seuil.
Cette fonction, établie de manière empirique, admet deux paramètres : la longueur de la tige
et le taux en G et en C de la séquence pour tenir compte d’un éventuel biais favorisant la
formation de tiges particulièrement stables.
Les points d’ancrage. caRNAc s’appuie sur des points d’ancrage entre les séquences
pour guider et accélérer le repliement et l’alignement des séquences. Ces points d’ancrage sont
des régions significativement conservées entre les deux séquences, sans insertion ni délétion.
L’algorithme se déroule de la manière suivante :
1. la recherche de tous les blocs maximaux conservés entre les deux séquences ;
2. le tri des blocs trouvés par score décroissant et filtrage des blocs chevauchants ;
3. la sélection gloutonne des blocs compatibles pour former des points d’ancrage.
Le système de score utilisé pour l’alignement est +1 en cas d’identité, −2 en cas de
substitution. Un bloc maximal est un bloc qui ne peut être étendu pour atteindre un score
plus élevé sans que ce score prenne une valeur négative ou nulle durant l’extension. Tout bloc
maximal dont la taille et le score sont supérieurs à 8 est ainsi conservé. Les blocs conservés
sont ensuite triés et filtrés pour éliminer les blocs qui impliquent au moins une même base
d’une des deux séquences. Bien que drastique, ce critère permet d’éviter de trancher entre
deux blocs qui pourraient a priori être corrects mais qui pourraient introduire une contrainte
erronée dans la suite du déroulement de l’algorithme. Une fois triés, les blocs sont sélectionnés
de manière gloutonne par score décroissant pour former des points d’ancrage. Un bloc est ainsi
sélectionné s’il est compatible avec l’alignement local déjà construit. Sur l’exemple suivant,
les blocs conservés rouges et bleus ne peuvent sélectionnés simultanément comme points
d’ancrage car ils introduiraient une incohérence dans l’alignement des deux séquences.
Séquence 1
Séquence 2
Le filtrage des tiges En fonction des points d’ancrage déterminés à l’étape précédente,
les couples de tiges copliables sont énumérés. Deux tiges sont copliables si elles présentent au
moins une covariation et si elles respectent les contraintes introduites par les points d’ancrage.
Le terme de covariation est ici à prendre au sens fort du terme, c’est-à-dire en présence d’une
mutation compensée : une covariation est comptée lorsque les deux bases d’un appariement
sont mutées d’une tige à l’autre. La recherche de covariations s’effectue sur les deux tiges
alignées sur leur structure primaire. Lorsqu’un couple de tiges présente au moins une cova-
riation, les deux tiges sont dites copliables si elles sont compatibles avec les points d’ancrage,
c’est-à-dire si les replier simultanément ne contredit pas l’alignement local des séquences selon
les points d’ancrage. Deux tiges ne sont donc pas copliables si elles correspondent à l’un de
ces trois cas :
89
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
1. violation d’un point d’ancrage : si (t1 , t1 ) et (t2 , t2 ) étaient repliées simultanément, alors
l’alignement de t1 et t2 contredirait l’alignement local au niveau du point d’ancrage.
t1
t1
t2
t2
2. décalage trop large à l’extérieur d’un point d’ancrage : lorsque l’ouverture ou la ferme-
ture d’une tige tombe entre deux points d’ancrage, un décalage borné est autorisé. Ce
décalage est variable selon les zones. Il dépend de la différence de longueur des fragments
de séquences entre les points d’ancrage. Sur cet exemple, le décalage entre t1 et t2 est
trop large.
t1
t1
t2 t2
t1
t2
t2
A ce niveau, toute tige copliable avec une autre tige est conservée pour l’étape suivante.
Une tige qui ne peut est copliée avec aucune tige est supprimée sauf si elle satisfait certains
critères :
– il s’agit d’une tige-boucle, c’est-à-dire si la taille de la boucle est d’une longueur
inférieure ou égale à huit nucléotides ;
– son énergie libre est relativement faible, en pratique le seuil est fixé à -1500 cal/mol ;
– elle se situe dans une région d’insertion potentielle, c’est-à-dire une région située
entre deux points d’ancrage consécutifs significativement plus longue que dans l’autre
séquence.
90
3.3. Evolution et enrichissement du logiciel caRNAc
Le corepliement des tiges copliables Le corepliement des tiges est le cœur algorith-
mique de la prédiction de structure secondaire commune de caRNAc, c’est aussi son origina-
lité. L’algorithme de repliement est une adaptation des récurrences de Sankoff, normalement
appliquées au niveau nucléotidique, au repliement de tiges complètes. La complexité de l’algo-
rithme de Sankoff n’est alors plus fonction de la taille des séquences mais du nombre de tiges
potentielles. De plus, comme seuls les couples de tiges copliables sont considérés, la taille du
problème se retrouve alors encore considérablement réduite.
Une tige t est caractérisée par les positions des extrémités de sa partie ouvrante, t.lef topen
et t.lef tclose, et de sa partie fermante, t.rightopen et t.rightclose, comme illustré sur la
figure 3.11. Les récurrences de caRNAc reposent sur trois applications next, last et prev
permettant à l’algorithme de naviguer entre les tiges, comme illustré sur la figure 3.12. A
partir de l’ensemble A des n tiges potentielles d’une séquence sa , on définit deux listes A→
et A← ordonnées des tiges de A.
AU ... CG CG ... AU
↑ ↑ ↑ ↑
t.lef topen t.lef tclose t.rightopen t.rightclose
Fig. 3.11 – Chaque tige t est décrite par les positions des extrémités de sa partie ou-
vrante (gauche), t.lef topen et t.lef tclose, et de sa partie fermante (droite), t.rightopen et
t.rightclose.
a1 a2 a3 a2 a4 a3 a4
Fig. 3.12 –
A→ = (a1 , a2 , . . . , an ) désigne la liste des tiges potentielles ordonnées par ordre crois-
sant d’ouverture : ai ≤ aj si et seulement si ai .lef topen ≤ aj .lef topen. L’application
next : [1..n] −→ [1..n] permet d’obtenir pour une tige t la prochaine tige dont la partie
ouvrante ne chevauche pas celle de t :
91
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
L’application prev : [1..n] −→ [1..n] permet d’obtenir pour une tige t la tige précédente
dont la partie fermante ne chevauche pas la partie ouvrante de t :
S(i, j, k, l) = min
S(i, j, k, l − 1)
S(i, j − 1, k, l)
min {S(i, prev(x), k, l) + S(next(x), last(j), 0, 0) + bind(ax = aj , −)}
1≤x≤n
min {S(i, prev(x), 0, 0) + S(next(x), last(j), k, l) + bind(ax = aj , −)}
1≤x≤n
min S(i, j, k, prev(y)) + S(0, 0, next(y), last(l)) + bind(−, by = bl )
1≤y≤m
min S(0, 0, k, prev(y)) + S(i, j, next(y), last(l)) + bind(−, by = bl )
1≤y≤m
min S(i, prev(x), k, prev(y)) + S(next(x), last(j), next(y), last(l)) + bind(ax = aj , by = bl )
1≤x≤n
1≤y≤m
Construction du graphe des tiges Le graphe des tiges est un graphe non dirigé où
chaque nœud correspond à une tige apparaissant dans au moins un corepliement, et une
arête entre deux nœuds indique que les tiges associées aux nœuds reliés ont été copliées. La
figure 3.13 montre un exemple de graphe des tiges obtenu sur cinq ARN de transfert.
92
3.3. Evolution et enrichissement du logiciel caRNAc
graphe associé
2 3 z }| {
4
a
2 3 4
1 e
a a a
a e e
b e
c → b b
b
d d
d
e c c c
Fig. 3.13 – Graphe des tiges construit après les corepliements de cinq ARN de transfert.
b
b1
a b2
a
a
1 2 2 1
b
c2
1 2 2 1 c
c1
c
Fig. 3.14 – Regroupement de tiges lorsqu’elles sont emboı̂tées. La tige de la séquence a s’est
repliée avec la tige b1 et la tige c2 tandis que les tiges b2 et c1 ont été copliées. Toutes ces tiges
sont correctes, mais celles des séquences b et c sont considérées comme deux tiges distinctes
emboı̂tées. Dans le graphe, elles sont fusionnées et les arêtes correspondantes sont regroupées.
93
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Pour qu’un couple de tiges soit copliable, il est nécessaire que les tiges présentent au moins
une covariation. Lorsque qu’un jeu de données comporte deux séquences proches partageant
une structure commune, il est fort probable qu’une partie des tiges communes ne présentent
aucune covariation et n’aient donc pas été copliées. Toutefois, ces tiges peuvent avoir été
copliées par ailleurs et se retrouvent donc dans le graphe des tiges. Un deuxième type d’arête
est donc introduit pour identifier les couples de tiges qui ne présentent pas de covariation :
les arêtes étiquetées identité. Les arêtes qui correspondent à un corepliement de deux tiges
seront étiquetées coplié. Cette modification permet d’améliorer la connexité du graphe en
présence de tiges qui n’ont pas pu être copliées car elles ne présentaient pas de mutations
compensatoires.
La sélection des tiges séquence par séquence Comme chaque tige appartient à une
et une seule composante connexe, on attribue à une tige l’indice de la composante qui la
contient. Pour chaque séquence, les tiges sont ensuite incorporées de manière gloutonne par
indice décroissant jusqu’à un certain seuil. L’incorporation se fait également sous la contrainte
de ne pas entrer de conflit avec la structure secondaire en cours de construction : les croise-
ments d’appariements et les chevauchements de tiges sont ainsi interdits. Toutefois, un léger
chevauchement est autorisé entre les tiges, et résolu en tronquant la tige la plus longue. Cette
liberté par rapport aux contraintes initiales permet de récupérer a posteriori des tiges maxi-
males légèrement chevauchantes qui n’auraient pas pu être repliées simultanément.
94
3.3. Evolution et enrichissement du logiciel caRNAc
qui ont guidé notre démarche sont que les approches “aligner puis replier” sont très perfor-
mantes quand les séquences sont proches, alors que les approches “aligner et replier simul-
tanément” sont plus robustes quand les séquences sont plus éloignées. Cela a été clairement
établi dans [GG04]. L’idéal serait donc d’avoir une approche tout terrain, qui permette de
traiter correctement des jeux de données hétérogènes, contenant des séquences à des distances
évolutives quelconques. Pour cela, nous proposons une solution basée sur les méta-séquences,
à l’image de ce que nous avons fait dans Protea.
Dans l’algorithme original de caRNAc, il existe une contrainte forte au repliement si-
multané de deux tiges : pour pouvoir être copliées, deux tiges doivent présenter une cova-
riation. Ce critère permet de prédire des tiges pour lesquelles il existe une réelle évidence
d’une évolution sous contrainte fonctionnelle. Cependant, ce critère de sélection peut poser
problème face à un jeu de données qui comporte un sous ensemble de séquences fortement
conservées. La redondance d’information apportée par des séquences proches perturbe ainsi le
fonctionnement de l’algorithme. Pour traiter ce problème, on propose de regrouper en amont
les séquences ressemblantes sous forme d’un alignement multiple, et d’adapter caRNAc pour
ne plus travailler uniquement sur des séquences individuelles, mais sur des méta-séquences,
c’est-à-dire des séquences individuelles et/ou des ensembles de séquences représentées par des
alignements multiples.
L’énergie associée à une méta-tige est définie comme la moyenne des énergies des tiges
individuelles qu’elle contient. Toutefois, on bonifie cette énergie pour chaque mutation qui
préserve un appariement.
95
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
AAUGGCCGUGUCAUCGGCCGGG
−−AGGCCGAGUCAUCGGCC−GG
ACUGGCCGAGUCAUCGGCCGGG
Les points d’ancrage entre méta-séquences La recherche de points d’ancrage est elle
aussi adaptée pour traiter les méta-séquences et s’effectue directement entre les alignements
multiples. La comparaison intra-méta-séquence n’est pas nécessaire car les séquences sont déjà
alignées. Etant données deux méta-séquences représentées par deux alignements multiples
U = {u1 , . . . , um } composé de m séquences alignées u1 , . . . , um et V = {v1 , . . . , vn } composé
de n séquences alignées v1 , . . . , vn , le score attribué à la comparaison de deux colonnes i et
j respectives de U et V est obtenu en sommant les scores de toutes les comparaisons deux à
deux entre la position i d’une séquence alignée de U et la position j d’une séquence alignée
de V . Plus formellement, ce score est calculé par la relation suivante
X X
s(i, j) = score(uk [i], vl [j])
1≤k≤m 1≤l≤n
où score(uk [i], vl [j]) est le score attribué dans la version originale de caRNAc lors de la
recherche de blocs conservés, c’est-à-dire +1 en cas d’identité entre uk [i] et vl [j] et −2 dans
le cas contraire. On assimile la comparaison entre un gap et un nucléotide à une substitution.
La procédure de sélection gloutonne des points d’ancrage reste identique à ceci près que le
seuil sur le score est corrigé. Dans la version originale, le seuil minimal pour retenir un bloc
est de 8. Comme le nombre de comparaisons réalisées est ici de m.n, ce seuil est maintenant
de 8m.n.
96
3.3. Evolution et enrichissement du logiciel caRNAc
Le filtrage des méta-tiges Pour que deux tiges soient décrétées copliables dans caRNAc,
il est nécessaire qu’elles présentent au moins une covariation et que leur repliement simul-
tané soit compatible avec les points d’ancrage. Pour un couple de méta-tiges, ces définitions
s’adaptent naturellement. Aucune covariation n’est attendue entre les tiges contenues dans
une même méta-tige. Entre deux méta-tiges T1 et T2 , on considère donc qu’une covariation
existe si au moins une tige de T1 présente une covariation avec au moins une tige de T2 . On
pourrait exiger que chaque tige de T1 présente au moins une covariation avec une tige de
T2 , cependant, dans les faits le premier critère est suffisamment sélectif pour retenir les bons
couples de méta-tiges. Pour la compatibilité avec les points d’ancrage en revanche, la modifi-
cation obéit aux mêmes contraintes que si on avait à faire à des séquences individuelles : on
impose que tous les couples de tiges (t1 , t2 ) issus d’un couple de méta-tiges (T1 , T2 ) soient
compatibles avec les points d’ancrage.
Le chevauchement de tiges maximales L’algorithme tel qu’il est conçu ne permet pas
de prédire simultanément deux tiges qui se chevauchent ne serait-ce que d’une base. Cette res-
triction peut poser un problème lorsque les vraies tiges d’une structure ne sont pas maximales
et que leur extension entraı̂ne un chevauchement comme illustré sur la figure 3.16.
G G G C C C A U A G ... ... U C A U G G G A U C G A A U C C C A U G G G C C C
G G G C C C A U A G ... ... U C A U G G G A U C G A A U C C C A U G G G C C C
Fig. 3.16 – Les tiges réelles ne sont pas nécessairement maximales. Les tiges maximales (b)
qui correspondent aux tiges de l’exemple (a) se chevauchent de deux bases et sont donc
incompatibles entre elles.
Les tiges potentielles considérées dans caRNAc à l’issue de la première étape sont
systématiquement des tiges maximales. Par conséquent, sur l’exemple de la figure 3.16 une
seule des deux tiges pourrait donc être prédite.
La gestion des chevauchements est introduite dans l’algorithme en modifiant les définitions
des applications next, last et prev (page 91). Soit δ le nombre de bases autorisées à se che-
vaucher, on redéfinit ces applications de la manière suivante
En pratique, on fixe δ = 2, ce qui est suffisant pour rattraper les vraies tiges à partir de
tiges maximales correspondantes, sans pour autant introduire de fausses tiges.
97
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
98
3.4. Résultats expérimentaux
optimale du repliement simultané des tiges a′j et a′l , c’est-à-dire l’énergie du corepliement de b′j
et b′l ajoutée à l’énergie optimale du repliement des deux séquences entre les parties ouvrantes
et fermantes de ces tiges. La table T est une table de travail pour les calculs intermédiaires.
Etant donné que seules les parties fermantes des tiges sont indexées dans les tables S et ST ,
on défini l’application open(i) : [1..n] −→ [1..n] qui permet de localiser la partie ouvrante
d’une tige ai dans la liste A→ .
open(i) = {k ∈ [1..n]|ak = ai }
Cette application est également définie pour l’ensemble des tiges B de la seconde séquence.
La table S est remplie par position d’ouverture de tige décroissante selon la règle suivante
où la table ST est partiellement recalculée pour chaque couple (j, l). Les règles de remplis-
sage de ST dans le couple d’intervalles ([prev(open(j)); last(l)], [prev(open(j)); last(l)]) sont
les suivantes
ST (i − 1, k)
ST (i, k − 1)
ST (i, k) = min
si open(i) ≥ next(open(j))
ST (prev(open(i)), prev(open(k))) + S(i, k)
et open(k) ≥ next(open(l))
99
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Les résultats obtenus sur BRAliBase I par les deux versions de caRNAc sont synthétisés
dans les tables 3.7 et 3.8. Globalement, on constate qu’on obtient de meilleurs résultats avec
la nouvelle version de caRNAc, quelque soit le jeu de données et le degré de conservation.
Le M CC est en effet toujours supérieur ou égal à celui atteint par la version 2004. En terme
d’efficacité, la version 2008 de caRNAc s’avère beaucoup plus rapide, en particulier sur les
séquences longues où le temps de calcul est au minimum divisé par 7.
Pour les ARN de transfert et les RNAse P, la sensibilité et la spécificité des struc-
tures prédites sont systématiquement supérieures ou égales aux anciennes valeurs. Pour les
ARN ribosomiques en revanche, bien que le compromis sensibilité/spécificité soit globalement
amélioré, on constate une légère perte de spécificité. Si l’on observe plus finement les struc-
tures prédites pour ces séquences, on remarque que les faux positifs supplémentaires sont des
appariements qui appartiennent à la structure tertiaire. caRNAc prédit en réalité une tige
de la structure tertiaire au détriment d’une autre tige moins stable de la structure secondaire.
La tige prédite par caRNAc n’est pas considérée comme correcte dans le benchmark, bien
qu’elle existe dans la structure réelle.
L’évolution des performances varie en fonction du degré moyen de conservation du jeu de
données. Sur les jeux de données très conservés en moyenne, les résultats n’évoluent quasiment
pas. C’est ici la faible quantité de mutations qui est en cause : une tige ne peut en effet être
prédite que si elle présente au moins une covariation. En revanche sur les jeux de données
moyennement conservés, les résultats sont nettement meilleurs grâce à l’utilisation des méta-
séquences. Ces jeux de données comportent en effet des sous-ensembles de séquences très
conservées qui, lorsqu’ils ne font pas l’objet d’un traitement particulier, perturbent la méthode.
D’une part ces séquences présentent peu, voire pas du tout, de covariations et d’autre part
introduisent une redondance d’information qui fausse les statistiques du graphe des tiges.
Tab. 3.7 – Résultats de caRNAc version 2004 et version 2008 sur BRAliBase I.
100
3.4. Résultats expérimentaux
caRNAc caRNAc
Famille Conservation Longueur Accélération
2004 2008
medium 73 0,125 s 0,052 s 2,40
tRNA
high 73 0,515 s 0,043 s 11,98
medium 377 48 s 0,941 s 51,00
RNaseP
high 377 0,831 s 0,500 s 1,67
medium 1542 153 s 20 s 7,65
SSU
high 1542 1149 s 22 s 52,23
medium 2904 2916 s 116 s 25,14
LSU
high 2904 1394 s 97 s 14,37
Les résultats obtenus avant et après modification de caRNAc sont donnés dans la
table 3.9. Le repliement complémentaire par RNAfold permet d’améliorer les résultats glo-
baux. RNAfold a en effet tendance à compléter les structures prédites par caRNAc par
plus d’appariements corrects que de mauvais appariements. Cela se traduit par un meilleur
M CC, c’est-à-dire un meilleur compromis sensibilité/spécificité, provenant d’une sensibilité
qui augmente en moyenne de 22% alors que la spécificité ne diminue que de 15% en moyenne.
Par rapport à la version 2004 de caRNAc, tous les résultats vont dans le sens de ces obser-
vations. Toutefois, pour les petites sous-unités ribosomiques, la tige de la structure tertiaire
prédite par caRNAc induit en erreur le repliement thermodynamique en structure secondaire
de RNAfold, ce qui diminue de facto la sensibilité et la spécificité.
Tab. 3.9 – Résultats de caRNAc dont les structures prédites sont complétées par RNAfold.
101
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
(a) Résultats sur les ARN de transfert (b) Résultats sur les ARN de RNase P
Méthode Conserv. Sens. Spé. MCC Méthode Conserv. Sens. Spé. MCC
medium 77,8 100,0 0,880 medium 57,4 57,4 0,571
RNAalifold RNAalifold
high 90,5 100,0 0,950 high 78,9 77,8 0,782
medium 100,0 75,0 0,863 medium 70,4 55,1 0,620
Ilm Ilm
high 76,2 69,6 0,722 high 43,7 36,5 0,395
medium 100,0 100,0 1,000 medium 87,0 92,2 0,895
Pfold Pfold
high 95,2 100,0 0,975 high 66,2 88,7 0,765
medium 100,0 100,0 1,000 medium 61,5 96,7 0,770
caRNAc caRNAc
high 76,2 100,0 0,871 high 51,4 100,0 0,716
medium 94,3 95,0 0,945 medium 32,0 32,8 0,321
Dynalign Dynalign
high 54,8 54,5 0,535 high 40,3 39,6 0,397
medium 23,8 33,3 0,268 medium 5,2 22,7 0,107
FoldAlign FoldAlign
high 23,8 31,2 0,259 high 19,7 35,9 0,265
Structures prédites complétées par RNAfold Structures prédites complétées par RNAfold
medium 100,0 100,0 1,000 medium 61,1 67,3 0,639
RNAalifold RNAalifold
high 100,0 100,0 1,000 high 77,5 77,5 0,773
medium 100,0 100,0 1,000 medium 89,6 89,6 0,895
caRNAc caRNAc
high 100,0 100,0 1,000 high 70,8 66,2 0,684
Tab. 3.10 – Résultats de BRAliBase I sur les ARN de transfert et sur les ARN de RNase
P.
(a) Résultats sur les petites sous-unités ribosomiques (b) Résultats sur les grosses sous-unités ribosomiques
Méthode Conserv. Sens. Spé. MCC Méthode Conserv. Sens. Spé. MCC
medium 84,4 92,1 0,881 medium 75,0 92,1 0,831
RNAalifold RNAalifold
high 59,8 60,6 0,601 high 79,0 76,3 0,776
medium 59,9 51,5 0,554 medium 68,4 58,0 0,630
Ilm Ilm
high 51,3 43,0 0,469 high 49,0 39,3 0,438
medium - - - medium - - -
Pfold Pfold
high 70,9 92,6 0,810 high - - -
medium 53,4 91,5 0,699 medium 50,0 95,8 0,692
caRNAc caRNAc
high 41,9 94,1 0,628 high 50,9 94,1 0,692
medium - - - medium - - -
Dynalign Dynalign
high - - - high - - -
medium - - - medium - - -
FoldAlign FoldAlign
high - - - high - - -
Structures prédites complétées par RNAfold Structures prédites complétées par RNAfold
medium 88,0 89,8 0,889 medium 84,4 89,9 0,871
RNAalifold RNAalifold
high 59,3 58,3 0,588 high 79,2 77,3 0,782
medium 71,3 71,8 0,715 medium 86,3 85,6 0,859
caRNAc caRNAc
high 72,7 74,8 0,737 high 83,9 82,5 0,832
Tab. 3.11 – Résultats de BRAliBase I sur les ARN des petites et grosses sous-unités ribo-
somiques.
102
3.4. Résultats expérimentaux
résultats sur les séquences courtes, c’est-à-dire les ARN de transfert, et de longueur moyenne,
les RNase P, il s’avère incapable de traiter les séquences plus longues pour des raisons de
complexité et de applications numériques. En effet, pour effectuer une prédiction Pfold
calcule des probabilités qui peuvent être très faibles jusqu’à descendre sous la capacité du
type primitif utilisé dans l’implémentation de Pfold.
Sur les ARN de transfert, RNAalifold et Pfold sont globalement meilleurs que caR-
NAc. Pour Pfold, ses excellents résultats sur ce jeu de données sont biaisés car ces séquences
ont été utilisés pour entraı̂ner la méthode. Sur les ARN de RNase P, caRNAc et Pfold ob-
tiennent les meilleurs résultats sur le jeu de données medium avec un M CC respectif de 0,77 et
0,895. caRNAc est légèrement plus spécifique que Pfold avec une spécificité de 96,7% contre
92,2% pour Pfold, mais Pfold se montre beaucoup plus sensible. Si l’on compare les struc-
tures de caRNAc repliées par RNAfold à Pfold, les compromis sensibilité/spécificité des
deux méthodes sont strictement équivalents. Sur les ARN de RNase P très conservés (high),
les meilleures prédictions sont produites par RNAalifold avec un M CC de 0,782. On note
cependant sur ces données que caRNAc est le seul à ne prédire aucun appariement incorrect
puisque sa spécificité est de 100%, tout en prédisant plus d’un appariement sur deux de la
structure réelle. De plus, sur cette structure prédite par caRNAc les appariements ajoutés
par RNAfold n’améliore pas le compromis sensibilité/spécificité puisque le M CC passe de
0,716 à 0,684.
Les ARN ribosomiques sont des séquences particulièrement longues qui posent des
problèmes de complexité à Dynalign et FoldAlign. La mémoire requise par ces méthodes
pour replier ces séquences dépasse largement les capacités offertes par la machine utilisée
pour le benchmark, c’est à dire 1Go. A l’exception du jeu de données ssu high où caRNAc
s’avère la méthode la plus performante avec un M CC à 0,628, RNAalifold est la méthode
dont le compromis sensibilité/spécificité est meilleur sur les ARN ribosomiques. caRNAc
tient ses objectifs puisque la spécificité des structures qu’il prédit ne descend jamais en des-
sous de 90%. Cette spécificité accrue est un atout important : les appariements prédits par
caRNAc constituent des contraintes sûres pour guider un repliement purement thermodyna-
mique. Sur les données très conservées, cette stratégie permet d’atteindre un M CC supérieur
à RNAalifold. Notamment dans le cas des petites sous-unités très conservées, le M CC de
caRNAc+RNAfold vaut 0,737 contre 0,601 pour RNAalifold et 0,588 pour RNAali-
fold+RNAfold .
103
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
Cette expérience nous amène à confronter nos résultats à ceux des méthodes existantes qui
procèdent à la prédiction d’ARN non-codants par analyse comparative.
Nous avons dû adapter le protocole d’AlifoldZ pour deux raisons : caRNAc travaille
sur des séquences non alignées, et caRNAc ne calcule pas une structure commune consensus
mais une structure par séquence globalement partagée par toutes les séquences. Le protocole
que nous utilisons est donc le suivant pour un ensemble S de n séquences :
– prédiction par caRNAc de la structure commune des séquences de S qui se traduit par
l’obtention de n structures ;
– production d’un alignement A des séquences de S avec ClustalW ;
– mélange des positions de A à l’aide du script shuffle aln.pl d’AlifoldZ afin d’obtenir
100 alignements mélangés ;
– prédiction par caRNAc de la structure commune des séquences de chacun des aligne-
ments mélangés ;
– calcul d’un z-score individuel des structures prédites pour les séquences de S.
A l’aide des n z-scores individuels obtenus pour chaque séquence de S, on réalise enfin
une prédiction suivant un vote majoritaire : si plus de la moitié des n z-scores associés aux
n structures prédites par caRNAc sont inférieurs à un certain seuil α, alors on prédit un
ensemble d’ARN non-codants homologues. Dans le cas contraire, on considère que l’on ne se
trouve pas en présence d’ARN non-codants homologues.
Afin d’évaluer le potentiel de cette approche, nous avons sélectionné de manière aléatoire
quatre séquences de chaque famille présente dans Rfam. Pour chaque famille, nous avons
calculé le z-score moyen des structures prédites par caRNAc par rapport à cinquante aligne-
ments mélangés de ces séquences. Pour apporter un contrôle négatif, nous avons constitué un
jeu de données composé des “familles” de séquences extraites de dix alignements aléatoires
générés par famille en mélangeant l’alignement structural des séquences originales. Les
résultats obtenus sur ces deux jeux de données sont présentés sur la figure 3.17. Sur les
familles d’ARN non-codants homologues, le z-score moyen des structures prédites par caR-
NAc est en moyenne inférieur à celui calculé sur les ensembles de séquences aléatoires. Les
deux distributions du z-score moyen observées ne se détachent toutefois pas complètement,
et leur intersection est relativement importante. En fixant à −1 le seuil sur le z-score calculé,
environ 20% des séquences aléatoires sont prédits ARN non-codants à tort, et un peu moins
de 70% ARN non-codants sont correctement détectés.
Nous avons également envisagé l’utilisation de Sissiz pour obtenir des alignements mul-
tiples de même composition en di-nucléotides (section 3.2.4). Toutefois, l’implémentation de
Sissiz pose plusieurs problèmes. Sissiz approxime la distribution en di-nucléotides de l’ali-
gnement original de manière asymptotique, ce qui ne garantit pas d’obtenir strictement la
même distribution, en particulier sur des séquences courtes comme c’est souvent le cas pour
les ARN non-codants. D’autre part, Sissiz rencontre quelques problèmes numériques gênants
pour une utilisation systématique. L’absence complète d’un di-nuléotide provoque sous cer-
taines conditions une erreur fatale, de même qu’une répartition totalement équiprobable de
tous les di-nucléotides qui est considérée comme une “absence de signature” significative. Pour
toutes ces raisons, nous n’avons pas poussé nos investigations avec ce logiciel.
104
3.4. Résultats expérimentaux
0.18
ARN
0.16 Autre
0.14
0.12
Fréquence
0.1
0.08
0.06
0.04
0.02
0
-15 -10 -5 0 5 10
z-score
Fig. 3.17 – Répartition du z-score moyen observé de l’énergie libre des structures prédites
par caRNAc sur les familles d’ARN non-codants de Rfam (trait plein), et sur des familles
de séquences aléatoires (trait discontinu).
Nous avons décidé de faire varier plusieurs propriétés susceptibles d’influencer les per-
formances des méthodes : la méthode d’alignement employée, la degré de conservation des
séquences, le nombre de séquences et la qualité des structures des ARN non-codants. Les
choix des séquences a donc été réalisé avec l’objectif de pouvoir construire des ensembles de
séquences dont le pourcentage d’identité moyen varie de 40 à plus de 95%. Afin de mesurer
le gain d’information apporté par l’utilisation de plusieurs séquences similaires, nous avons
réalisé des alignements comportant de deux à cinq séquences. Nous n’avons pas construit d’ali-
gnement de plus de cinq séquences pour une raison simple : lorsque l’on dispose d’autant de
séquences similaires, les outils d’inférence de structures communes peuvent suffire à détecter
des ARN non-codants homologues. En effet, l’existence d’une structure commune à plus de
dix séquences constitue un signal fort pour identifier des ARN non-codants homologues. Dans
la section 3.2, nous avons vu que la détection des ARN non-codants à partir d’une séquence
dépend de la qualité des structures des ARN à détecter. Pour évaluer l’influence de la qualité
des structures sur la détection à partir de plusieurs séquences, nous avons sélectionné des
familles d’ARN non-codants dont les structures communes ont une stabilité variable. Nous
avons également retenu quelques familles dont les structures communes comportent des pseu-
donœuds pouvant gêner la prédiction.
A partir de ces propriétés, nous avons recueilli trois types de données provenant de
plusieurs organismes : des ARN non-codants homologues pour évaluer la sensibilité, des
fragments codants homologues d’ARN messagers et des séquence aléatoires pour évaluer
la spécificité. La répartition des données est la suivante : 21 familles d’ARN non-codants,
15 familles d’ARN messagers et 21 “familles” de séquences aléatoires. La majorité des fa-
105
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
milles d’ARN non-codants retenues proviennent de Rfam. Afin d’assurer la variabilité des
paramètres définies précédemment nous avons ajouté deux familles de micro ARN provenant
de miRBase. Dans la section 1.3, nous avons vu que des fragments d’ARN messagers, les
introns et les extrémités 3′ et 5′ non traduites, sont susceptibles de contenir des structures.
Pour produire des séquences qui ne contiennent a priori pas de structure, nous avons re-
tiré ces fragments des ARN messagers que nous avons utilisés. Les ARN messagers sont en
général composés de plusieurs milliers de bases, contrairement aux ARN non-codants dont
la longueur dépasse rarement les 300 bases. Pour constituer un jeu de données comparable
au jeu de données positif d’ARN non-codants, certains alignements ont par conséquent été
tronqués. Le second jeu de données négatives, les shuffles, est composé d’alignements positifs
mélangés par la procédure employée dans AlifoldZ [WH04] (section 3.2.4).
Le protocole expérimental
Trois des méthodes les plus récentes ont été testées sur les jeux de données ainsi constitués :
Qrna [RE01], ddbRNA [DBDH03] et RNAz [WHS05]. Contrairement à ddbRNA et RNAz
qui effectuent une classification binaire “ARN non-codants homologues”/“autre”, Qrna ef-
fectue une classification en trois classes : “codant”, “non-codant” et “autre”. L’objectif de
nos tests est d’évaluer les performances de détection des ARN non-codants. De notre point
de vue, les classes “codant” et “autre” sont équivalentes car elles ne correspondent pas à
la prédiction d’ARN non-codants homologues. Ces deux classes sont donc fusionnées. Pour
évaluer les performances des méthodes nous utilisons les trois notions classiques, à savoir la
sensibilité, la spécificité et le cœfficient de corrélation de Matthews. La sensibilité désigne ici
la proportion d’ARN non-codants homologues détectés, la spécificité la proportion d’ARN
messagers et d’alignements mélangés non prédits comme des ARN non-codants homologues.
Pour chaque famille, des ensembles de deux, trois et cinq séquences sont créés
aléatoirement. Chaque ensemble de séquences est ensuite aligné et soumis aux différentes
méthodes, caRNAc étant la seule méthode qui ne nécessite pas d’alignement préalable. Au
total, plus de 80 000 alignements ont ainsi été constitués.
L’influence de l’alignement
Nous avons fait appel à cinq méthodes d’alignement largement utilisées :
Blast [WBB+ 08] et Needleman&Wunsch [NW70] pour les alignements deux à deux,
ClustalW [THG94], T-Coffee [NHH00] et Dialign2-2 [Mor99] pour les alignements
de deux séquences et plus. Selon la méthode utilisée, les résultats varient sensiblement. En
moyenne, les meilleurs résultats sont obtenus sur les alignements produits par ClustalW
(figure 3.18). Cette observation globale se vérifie sur les résultats moyens de chaque méthode
de détection, quelque soit le nombre de séquences utilisées.
RNAz a été entraı̂né à reconnaı̂tre les ARN non-codants sur des alignements produits par
ClustalW. Il est donc normal qu’il obtienne de meilleurs résultats sur ce type d’alignements.
Les performances de ddbRNA sur les alignements de ClustalW s’expliquent par un nombre
moyen de gaps moins important dans ces alignements que dans les alignements de Dialign2-
2 et T-Coffee. En effet, les gaps sont une entrave à la recherche de tiges pratiquée dans
ddbRNA : les positions qui contiennent au moins un gap sont ignorées et ne contribuent
donc pas à la formation des tiges. Les résultats obtenus sur les alignements produits par
Blast et Needleman&Wunsch sont en moyenne inférieurs aux alignements de deux séquences
106
3.4. Résultats expérimentaux
0.8
0.6
MCC
0.4
0.2
ClustalW
Dialign2-2
T-Coffee
0
50 55 60 65 70 75 80 85 90 95
Pourcentage d’identité moyen
Fig. 3.18 – Résultats moyens des méthodes de détection selon la méthode d’alignement mul-
tiple utilisée. Les résultats sont exprimés en fonction du pourcentage d’identité des aligne-
ments. Ces résultats sont calculés à partir de tous les alignements de deux, trois et cinq
séquences.
produits par ClustalW (table 3.12). Cette observation est notamment valable pour Qrna
qui a pourtant été entraı̂né à reconnaı̂tre les ARN non-codants sur des alignements produits
par Blast. Ces résultats nous amènent à nous focaliser sur les alignements générés par
ClustalW.
Tab. 3.12 – Résultats moyens obtenus sur des alignements produits par Blast, Needle-
man&Wunsch et ClustalW avec deux séquences.
Le table 3.13 et la figure 3.19 donnent les résultats obtenus selon le nombre de séquences
utilisées. Les performances de ddbRNA et de RNAz sont en moyenne meilleures sur des
alignements de trois séquences. Toutefois, pour RNAz, la sensibilité est bien plus élevée en
utilisant cinq séquences. Quant à caRNAc, ses résultats croissent strictement avec le nombre
de séquences utilisées.
107
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
0.8
0.6
caRNAc (shuffle-aln.pl)
MCC
ddbRNA
RNAz
0.4 Qrna
0.2
0
50 55 60 65 70 75 80 85 90 95
Pourcentage d’identité moyen
0.8
0.6
caRNAc (shuffle-aln.pl)
MCC
ddbRNA
RNAz
0.4
0.2
0
50 55 60 65 70 75 80 85 90 95
Pourcentage d’identité moyen
0.8
0.6
caRNAc (shuffle-aln.pl)
MCC
ddbRNA
RNAz
0.4
0.2
0
50 55 60 65 70 75 80 85 90 95
Pourcentage d’identité moyen
108
3.4. Résultats expérimentaux
Tab. 3.13 – Résultats selon le nombre de séquences utilisées. Ces résultats sont établis sur
les alignements réalisés avec ClustalW, et exprimés en pourcentage pour la sensibilité et la
spécificité. La spécificité moyenne est calculée sur les ARN messagers et les shuffles.
L’influence de la conservation
La conservation entre les séquences est un paramètre dont l’influence varie suivant la
méthode employée. En moyenne, les meilleurs résultats proviennent des alignements dont le
pourcentage d’identité moyen est compris entre 60% et 85%. caRNAc est la seule méthode
capable de traiter convenablement des jeux de données dont la conservation moyenne est
inférieure à 60% d’identité. Entre 60% et 85% d’identité, toutes les méthodes ont une
spécificité moyenne supérieure à 80%. Néanmoins, hors de cet intervalle la spécificité moyenne
reste élevée et ne descend pas en dessous de 75%. Par contre, la sensibilité de Qrna et de
ddbRNA se dégrade rapidement lorsque l’on dépasse 90% d’identité moyenne.
La spécificité moyenne sur les shuffles est à peu près équivalente à la spécificité moyenne
sur les ARN messagers. Lorsque le pourcentage d’identité est inférieur à 80%, la spécificité
moyenne sur les shuffles est très légèrement inférieure à celle des ARN messagers ; la situation
s’inverse au delà de 80% d’identité.
L’influence de la conservation est en réalité étroitement liée à la qualité d’alignement.
En effet, des séquences homologues mal conservées partagent une structure commune que
toutes les méthodes ont dû mal à prédire à partir d’un alignement qui n’est pas correct. C’est
également la raison pour laquelle caRNAc est la seule méthode qui produise des résultats
probants en dessous de 60% d’identité moyenne. Sur la figure 3.20 est représentée à gauche
la structure secondaire d’un ARN non-codant présent dans la partie 3′ de l’ARN des virus
de la famille des pomovirus. Sur cette figure est également présenté un alignement de trois
séquences homologues de ce type d’ARN non-codant produit ClustalW et extrait de notre
jeu de données. Bien que plus de 90% des positions de cet alignement soient correctes, seul
caRNAc prédit des séquences d’ARN non-codants homologues. La structure prédite par
caRNAc sur ces séquences est présentée en partie droite de la figure 3.20.
109
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
40
50 G U A
A A A
A
C A
C U
G C
U U G
A 50
G C G
A
C A C U A
G C A G G C A C
30 U U
U C G C
G C 60 C C
C C
G C C
U G U C
40 A U
30
A C U C
A C
C G G 60
A U C
A A A C 70
C U
A G U
C U U G U C
C C U C C C A C
C C G A C C
G G G
G G C
20
G U
20 U U
G A G A C
A 80
A C C
C C A C C G 70
G G C G C
C U G C G U A A G
C U U 80 G
U G A A C C A A
G U U C
G G U C
G C U 5’
C G
A C G A
5’ U
A C C
10 10
T91413.1 UUAGCUCGC-CAGUGCGAGGCCUCUUCCUACACAAGAGGUAU---UGG-GGUGCGACUCCCCCGUCUAUCCUGAACGUCAUCAGGACCA
X54354.1 UUAGCUCGC-CAGUGCGAGGCUCGUUCCCACACAACAAGUAA---UGGUGGUGCAACUCCCCCGUCC-UCCCGAACGUCAUCGGGACCA
Y16104.1 UAAUUGAGGACAGUUCCUCUCCCUCUAGCACACAGA-GGUCAAACUGGGUG--CAACUCCCCCC-CCUUCCGUGG-GUAACGGAAACC-
(c) Alignement extrait de Rfam.
T91413.1 UUAGCUC--GCCAGUGCGAGGCCUCUUCCUACACAAGAGGUAUUGG-GGUGCGACUCCCCCGUCUAUCCUGAACGUCAUCAGGACCA
X54354.1 UUAGCUC--GCCAGUGCGAGGCUCGUUCCCACACAACAAGUAAUGGUGGUGCAACUCCCCCGUCC-UCCCGAACGUCAUCGGGACCA
Y16104.1 -UAAUUGAGGACAGUUCCUCUCCCUCUAGCACACAGAGGUCAAACUGGGUGCAACUCCCCC--CCCUUCCGUGGGUAACGGAAACC-
************************** ********** *************
(d) Alignement produit par ClustalW. Le symbole * marque les positions correctes.
Fig. 3.20 – Structure secondaire réelle (à gauche) et structure prédite par caRNAc (à droite)
d’un ARN non-codant présent dans l’ARN des pomovirus, ici celle du Cacao yellow mo-
saic virus. En partie inférieure, l’alignement produit par ClustalW de ladite séquence et
de deux séquences homologues ainsi que l’alignement structural correspondant extrait de
Rfam (RF00233). Ces séquences ont un pourcentage d’identité moyen est de 66%.
110
3.4. Résultats expérimentaux
Les conclusions
Face à ses concurrents, caRNAc tire son épingle du jeu sur les séquences faiblement
conservées où il est le seul à fournir des résultats pertinents. Sur ce type de données, les autres
méthodes sont induites en erreur par un alignement incorrect. Sur les séquences relativement
bien conservées en revanche, RNAz se dégage nettement de toutes les méthodes existantes
en terme de sensibilité. Au niveau spécificité, Qrna, ddbRNA et RNAz sont à peu près
équivalentes, bien que Qrna soit nettement plus spécifique sur les ARN messagers grâce à
son modèle pour détecter les séquences codantes homologues.
Contrairement à caRNAc, la conception de RNAz repose sur un système d’apprentissage
très sophistiqué entraı̂né sur un très grand nombre de séquences issues d’un large éventail de
familles d’ARN non-codants. Cette caractéristique est un point fort pour RNAz lorsque le
jeu de données qu’on lui soumet répond aux critères pour lesquels il a été entraı̂né. Cet atout
peut toutefois s’avérer limitant, notamment en ce qui concerne le nombre de séquences. En
effet, RNAz ne peut pas traiter de jeux de données comportant plus de dix séquences car son
processus d’apprentissage a été limité à des jeux de données contenant au plus dix séquences.
caRNAc en revanche n’est pas limité en nombre de séquences. Qui plus est, les différentes
expériences que nous avons présentées montre que les performances de caRNAc augmentent
avec le nombre de séquences. Cette propriété lui confère un net avantage par rapport aux
méthodes existantes dont les performances sont limitées par les difficultés à produire un
alignement multiple d’un grand nombre de séquences.
111
Chapitre 3. Prédiction de structures communes d’ARN non-codants homologues
112
Chapitre 4
Au cours des chapitres 2 et 3 nous avons présenté deux méthodes que nous avons mis au
point pour la prédiction de séquences codantes homologues et de séquences non codantes qui
partagent une structure. L’originalité de ces méthodes réside dans le traitement d’ensembles de
séquences non alignées par analyse comparative qui permet d’obtenir des résultats significatifs
sur des séquences faiblement conservées. De plus, ces méthodes tirent parti du concept que
nous avons introduit, les méta-séquences (section 1.5.3), qui permet d’éliminer les redondances
de séquences au sein d’un jeu de données et donc de traiter des ensembles de séquences
hétérogènes en terme de conservation.
Dans ce chapitre, nous nous intéressons à l’intégration de ces méthodes dans deux projets
collaboratifs réalisés au sein de l’équipe. Dans la section 4.1, nous présentons Magnolia, une
méthode d’alignement multiple de séquences nucléiques fonctionnelles basée sur les prédictions
de Protea et de caRNAc. Dans la section 4.2, nous présentons un pipeline d’annotation
par génomique comparative.
113
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
commune que nous avons développée. Enfin, nous terminons cette section par une évaluation
des performances de Magnolia.
114
4.1. L’alignement multiple de séquences nucléiques
ClustalW,
Dialign2-2,
Protea T-Coffee
d’acides aminés
séquences couples de
nucléiques séquences
structures
2 séquences toutes les séquences primaire et secondaire
caRNAc Gardenia
115
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
N
X N
X
Si = pijk
j=1,j6=k k=1
1 si Si = N (N − 1)
Ci =
0
PM
i=1 Si
SP S = PM r
i=1 Sri
PM
i=1 Ci
CS =
M
où Mr est le nombre de colonnes de l’alignement de référence et Sri est le score de la ième
colonne de l’alignement de référence. Le SPS et le CS prennent leurs valeurs dans l’intervalle
[0; 1], 1 étant la valeur maximale où l’alignement A correspond exactement à l’alignement de
référence. Sur l’exemple de la figure 4.2, l’alignement par Protea +ClustalW est signifi-
cativement plus proche de l’alignement de référence de Pandit que les alignements générés
par ClustalW, Dialign2-2 et MultAlin.
116
4.1. L’alignement multiple de séquences nucléiques
O97702_CANFA TGCAGCCCCCGGGAGGGCCAGCCCGCCTGCAGCCAGCGGGGCGAGTGCCTG------TGTGGCCAATGTGTCTGCCATAGCAGTGACTTTGGCAAGATCACGGGCAAGTACTGC
Q86G85_PSEIC TGCCGGTCACCTGAAAACAACGAAATCTGCAGTGGAAACGGACAATGTGTA------TGTGGACAATGTATGTGTAACTCTGACGATGACCGCCACTATAGTGGCAAATACTGC
Q19267_CAEEL TGTTTTGGAAAAGGATCC---------TGTCATGGAGATGGAAGCCGCGAAGGCAGT---GGAAAGTGTAAATGTGAGACTGGA------------TATACTGGAAATCTATGC
** * * ** ** * ** * *** ** * * ** ** ***
O97702_CANFA TGCAGCCCCCGGGAGGGCCAGCCCGCCTGCAGCCAGCGGGGC---GAGTGCCTGTGTGGCCAATGTGTCTGCCATAGCAGTGACTTTGGCAAGATCACGGGCAAGTACTGC
Q86G85_PSEIC TGCCGGTCACCTGAAAACAACGAAATCTGCAGTGGAAACGGA---CAATGTGTATGTGGACAATGTATGTGTAACTCTGACGATGACCGCCACTATAGTGGCAAATACTGC
Q19267_CAEEL TGTTTT---------GGAAAAGGATCCTGTCATGGAGATGGAAGCCGCGAAGGCAGTGGAAAGTGTAAATGTGAGACTGGA------------TATACTGGAAATCTATGC
** * *** ** **** * *** ** * * ** ** ***
(b) Alignement de Protea. L’alignement au niveau peptidique a été confié à ClustalW. Le SPS de cet
alignement vaut 0,83 et son CS 0,75.
O97702_CANFA TGCAGCCCCCGGGAGGGCCAGCCCGCCTGCAG-CCAGCGGGGCGAGTGCCTGTGTGGCCAATGTGTCTGCCATAGCAGTGACTTTGGCAAGATCACGGGCAAGTACTGC
Q86G85_PSEIC TGCCGGTCACCTGAAAACAACGAAATCTGCAG-TGGAAACGGACAATGTGTATGTGGACAATGTATGTGTAACTCTGACGATGACCGCCACTATAGTGGCAAATACTGC
Q19267_CAEEL ----TGTTTTGGAAAAGGATCCTGTCATGGAGATGGAAGCCGCGAAGGCA---GTGGAAAGTGTAAATGTGAGACTGGATATACTGGAAATCTATGC------------
* ** ** * * * **** * *** ** * * * *
(c) Alignement de ClustalW des séquences nucléiques. Le SPS de cet alignement vaut 0,524 et son CS 0,286.
O97702_CANFA TGCAGCCCCCGGGAGGGCCAGCCCGCCTGCAGCCAGCGGGGCGAGTGCCTGTGTGGCCAATGTGTCTGCCATAGCAGTGACTTTGGCAAGATCACGGGCAAGTACTGC
Q86G85_PSEIC TGCCGGTCACCTGAAAACAACGAAATCTGCAGTGGAAACGGACAATGTGTATGTGGACAATGTATGTGTAACTCTGACGATGACCGCCACTATAGTGGCAAATACTGC
Q19267_CAEEL ----TGTTTTGGAAAAGGATCCTGTCATGGAGATGGAAGCCGCGAAGGCAGTGGAAAGTGTAAATGTGAGACTGGA--------------TATACTGGAAATCTATGC
* ***** * ** * **** * * ** ** ***
(d) Alignement de MultAlin des séquences nucléiques en utilisant des informations au niveau peptidique. Le
SPS de cet alignement vaut 0,476 et son CS 0,210.
O97702_CANFA TGCAGCCCCCGGGAGGGCCAGCCCGCCTGCAGCCAGCGGGGCGAGTGCCTGTGTGGCCAATGTGTCTGCCATAGCAGTGACTTTGGCAAG--------------------
Q86G85_PSEIC TGCCGGTCACCTGAAAACAACGAAATCTGCAGTGGAAACGGACAATGTGTATGTGGACAATGTATGTGTAACTCTGACGATGACCGCCAC--------------------
Q19267_CAEEL ---------------------------TGTTTTGGAAAAGGATCCTGTca------------------------TGGAGATGGAAGCCGcgaaggcagtggaaagtgtaa
** ** ** ** **
O97702_CANFA -------------ATCACGGGCAAGTACTGC
Q86G85_PSEIC -------------TATAGTGGCAAATACTGC
Q19267_CAEEL atgtgagactggaTATACTGGAAATCTATGC
* ** ** ***
(e) Alignement de Dialign2-2 des séquences nucléiques en utilisant des informations au niveau peptidique. Le
SPS de cet alignement vaut 0,476 et son CS 0,210.
Fig. 4.2 – Un exemple de trois séquences de la famille PF07974 de Pandit alignées par
Protea +ClustalW, ClustalW, Dialign2-2 et MultAlin. Le pourcentage d’identité
moyen de ces séquences est de 44,3%. Le pourcentage d’identité moyen au niveau peptidique
est de 30,3%. Le SPS et le CS de chaque alignement est calculé en utilisant l’alignement fourni
dans Pandit comme référence.
117
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
Implantation de Magnolia
Magnolia est développé sous forme d’un site Web qui fait appel aux différentes com-
posantes et synthétise les résultats. Lorsqu’une prédiction “codant” est réalisée par Pro-
tea, deux alignements sont produits : l’alignement multiple des séquences d’acides aminées
prédites, et la rétro-transcription de cet alignement. Plusieurs méthodes sont proposées à
l’utilisateur pour l’alignement au niveau peptidique : ClustalW, Dialign2-2 et T-Coffee.
La mise en couleur des acides aminés du premier alignement et des codons correspondants
dans le second alignement est inspirée des couleurs de RasMol2 . La figure 4.3 montre un
exemple des alignements produits par Magnolia pour une famille de séquences codantes
homologues.
Fig. 4.3 – L’alignement par Magnolia (Protea) du domaine Zn-finger des protéines Ran
(Pfam PF00641). La longueur moyenne des séquences est de 92 nucléotides et leur pourcen-
tage d’identité moyen de 45,1%. Les triplets de bases sont coloriés en fonction de l’acide aminé
codé. L’alignement de référence est quasiment identique à celui fourni dans Pandit.
118
4.1. L’alignement multiple de séquences nucléiques
multiple des séquences annoté par la structure est généré. Chaque tige prédite par caRNAc
est coloriée. La figure 4.4 montre un exemple d’alignement pour une famille de séquences
non codantes homologues partageant une structure commune. Les structures secondaires de
chaque séquences sont également représentées sous forme graphique à l’aide du programme
NaView [BH88].
119
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
Pour chaque famille de Pandit, un sous-ensemble de quatre séquences ont été choisies
aléatoirement. Sur les 6 491 ensembles ainsi construits, 6 122 (94,3%) sont correctement
prédites “codant” par Magnolia, et pour plus de 99% d’entre elles les cadres de lecture
prédits sont corrects. Moins de 3% des familles sont prédites “structurés” par Magnolia.
Pour estimer la qualité des alignements multiples produits, nous utilisons le SPS décrit à la
page 116. Comme Magnolia s’appuie sur une méthode d’alignement multiple externe pour
aligner les séquences d’acides aminés prédites, nous avons testé trois méthodes différentes :
ClustalW, Dialign2-2 et T-Coffee. Les alignements de Magnolia sont comparés aux
120
4.2. L’annotation par génomique comparative
alignements produits par ces mêmes méthodes utilisées sur les séquences nucléiques initiales.
Les résultats sont présentés en figure 4.5 en fonction du pourcentage d’identité moyen des
séquences nucléiques. Quelque soit le degré de conservation des séquences, les alignements de
Magnolia sont plus proches des alignements de référence de Pandit que les autres méthodes
d’alignement multiple testées. Plus les séquences sont divergentes, plus l’écart se creuse entre
les alignements de Magnolia et les autres, quelque soit la méthode d’alignement sous-jacente
utilisée.
121
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
SCI
0.4 0.6
Magnolia 0.4
0.2 ClustalW
Poa 0.2
ProAlign
Pcma
0 0
32 34 36 38 40 42 44 46 48 32 34 36 38 40 42 44 46 48
Pourcentage d’identité moyen Pourcentage d’identité moyen
SCI
0.4 0.6
Magnolia 0.4
0.2 ClustalW
Poa 0.2
ProAlign
Pcma
0 0
34 36 38 40 42 44 46 48 34 36 38 40 42 44 46 48
Pourcentage d’identité moyen Pourcentage d’identité moyen
SCI
0.4 0.6
Magnolia 0.4
0.2 ClustalW
Poa 0.2
ProAlign
Pcma
0 0
38 40 42 44 46 48 38 40 42 44 46 48
Pourcentage d’identité moyen Pourcentage d’identité moyen
Fig. 4.6 – Résultats de Magnolia sur BRAliBase 2.1 comparés aux résultats des méthodes
d’alignement multiple traditionnelles. Le SPS et le SCI des alignements sont présentés en
fonction du pourcentage d’identité moyen des séquences alignées et du nombre de séquences
utilisées.
122
4.2. L’annotation par génomique comparative
1.2
0.8
1
0.6 0.8
SPS
SCI
0.4 Magnolia 0.6
FoldAlignM
Lara
Mlocarna 0.4
0.2 StrAl
Marna
MxscaRNA 0.2
PMmulti
R-Coffee
0 0
32 34 36 38 40 42 44 46 48 32 34 36 38 40 42 44 46 48
Pourcentage d’identité moyen Pourcentage d’identité moyen
1.2
0.8
1
0.6 0.8
SPS
SCI
1.2
0.8
1
0.6 0.8
SPS
SCI
Fig. 4.7 – Résultats de Magnolia sur BRAliBase 2.1 comparés aux résultats des méthodes
qui construisent un alignement multiple structural. Le SPS et le SCI des alignements sont
présentés en fonction du pourcentage d’identité moyen des séquences alignées et du nombre
de séquences utilisées.
123
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
Séquence à Banque de
annoter séquences
Yass Comparaisons 2 à 2
Régions conservées
caRNAc Protea
R S Y L
Structure conservée Séquence d’acides
GGGGGTAACCCC aminés conservés CGATCCTATTTA
CGCGGCAACGCG CGCAGTTACTTG
TGGGGTAACTCG AGAAGCTACCTA
Fig. 4.8 – Pipeline d’annotation automatique de séquences codantes par Protea et d’ARN
non-codants par caRNAc dans une séquence génomique.
Le choix des génomes à comparer est cruciale car la qualité des prédictions qui peuvent être
réalisées dépend pleinement de ces séquences. De manière générale, le facteur déterminant est
les distances évolutives qui séparent les organismes dont elles sont issues de l’organisme dont
provient la séquence à annoter. Prenons l’exemple du génome d’Escherichia coli dans lequel
on cherche à identifier de nouvelles séquences codantes ou d’ARN non-codants. La séquence à
124
4.2. L’annotation par génomique comparative
annoter sera alors la séquence génomique d’une souche d’Escherichia coli, et la banque pourra
alors être constituée de séquences génomiques d’autres bactéries plus ou moins éloignées en
terme d’évolution. Avec des séquences génomiques très proches de celle de l’organisme ciblé,
le risque majeur est que les séquences homologues exhibent trop peu de mutations pour
qu’elles puissent à elles seules faire l’objet d’une analyse comparative pertinente. A l’inverse,
si l’on choisit des séquences d’organismes trop éloignés, on risque de ne pas être en mesure
d’identifier les séquences homologues trop peu conservées, ou pire, que les organismes choisis
ne possèdent pas d’homologue pour une séquence fonctionnelle putative du génome à annoter.
Bien qu’il n’existe pas de critère absolu pour déterminer les séquences génomiques à choisir,
certaines situations font naturellement émerger des contraintes. Par exemple, si on cherche
à identifier une séquence liée à un phénotype particulier tel que la production d’un agent
pathogène, il apparaı̂t alors nécessaire de considérer d’autres souches du même organisme qui
partagent ce même phénotype.
Dans le cas où l’on souhaite découvrir de nouvelles séquences codantes ou d’ARN non-
codants, il apparaı̂t naturel de vouloir masquer les régions qui comportent déjà des annotations
de la séquence à annoter, ou les régions susceptibles de parasiter les comparaisons telles que
des régions hautement conservées (plus de 95%) ou des éléments répétés. Le pipeline offre ainsi
la possibilité de masquer automatiquement des régions à ignorer durant la comparaison avec
la banque. Le masquage des régions déjà annotées permet de diminuer de manière radicale
le nombre de régions à comparer et par conséquent le nombre de prédictions à traiter en
sortie du pipeline. Toutefois, cette mesure drastique peut également conduire à manquer des
séquences codantes ou d’ARN non-codants inconnues qui chevaucheraient ou seraient incluses
dans des régions déjà annotées pour une autre fonction.
125
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
(a) Identification des k-mers puis des (b) Evaluation puis filtrage des 10
enchaı̂nements de k-mers compatibles meilleures chaı̂nes de k-mers selon leur
entre eux, c’est-à-dire les diagonales score
qui apparaissent sur le dotplot
Fig. 4.9 – Schéma des étapes principales de Fasta pour deux séquences A et B.
Chronologiquement, Fasta [LP85] est la première heuristique d’alignement local qui uti-
lise le principe des graines contiguës. Une graine contiguë est une séquence d’une certaine
longueur k, un k-mer, parfaitement conservée entre deux séquences. La valeur de k consti-
tue un paramètre crucial de la méthode qui affecte à la fois sa sensibilité et son efficacité.
126
4.2. L’annotation par génomique comparative
Plus k est grand, plus la méthode est rapide au détriment de sa sensibilité. L’heuristique de
Fasta se décompose en quatre grandes étapes illustrées par les schémas de la figure 4.9. La
première étape consiste à rechercher à l’aide d’une graine contiguë les k-mers conservés. Fasta
cherche ensuite à créer des régions locales similaires, c’est-à-dire à regrouper des k-mers afin
de détecter des séquences conservées plus longues, pouvant contenir des substitutions mais
n’introduisant aucun gap. Pour chaque région, Fasta calcule ensuite un score puis filtre les
meilleures régions afin de ne garder que les dix meilleures. Le calcul de ce score fait intervenir
une matrice de substitution triviale, la matrice identité. Les régions trouvées sont ensuite
incorporées de manière gloutonne par score décroissant afin de ne garder que les régions com-
patibles entre elles, c’est-à-dire des régions qui peuvent faire partie d’un même alignement.
L’alignement local est enfin produit par l’algorithme de Smith&Waterman où la matrice de
programmation dynamique est réduite à la zone qui englobe les régions retenues à l’étape
précédente.
L’heuristique de Blast [AGM+ 90] est globalement identique à celle de Fasta. La
différence majeure entre Blast et Fasta se situe au niveau du passage des k-mers à un
alignement. Pour chaque k-mer trouvé, Blast procède à son extension de part et d’autre
de façon à trouver une région conservée la plus longue possible. Blast continue d’étendre la
région tant que le score cumulé calculé au fur et à mesure de l’extension ne descend pas en
dessous d’un certain seuil. Contrairement à Fasta, Blast applique des coûts variables aux
substitutions : 5 pour un match et −4 pour un mismatch. Les régions ainsi obtenues, appelées
des HSP (High-scoring Segment Pairs), sont ensuite filtrées en fonction de l’espérance statis-
tique de leurs scores. Dans la version “avec gap” de Blast les k-mers trouvés peuvent être
groupés pour former une HSP avant leur extension si la distance qui les sépare ne dépasse
pas un certain seuil.
PatternHunter [MTL02] et Yass [NK05] sont des heuristiques qui fonctionnent selon le
même schéma global que Blast. Leur originalité tient à l’utilisation des graines espacées, c’est-
à-dire des sous-séquences conservées. Les graines espacées apportent une meilleure sensibilité
que les graines contiguës, permettant ainsi de trouver des régions moins conservées, sans pour
autant dégrader ni la spécificité ni l’efficacité apportées par les graines contiguës. Le graphique
de la figure 4.10 montre le gain de sensibilité théorique apporté par l’utilisation d’une graine
espacée par rapport à une graine contiguë de même contenu informationnel, c’est-à-dire où
le nombre de nucléotides comparés est identique et appelé poids d’une graine. La figure 4.11
montre un exemple d’alignement qui ne pourrait être obtenu avec une graine contiguë de
même poids que la graine espacée utilisée. En effet, bien que ces séquences soient similaires,
celles-ci ne contiennent pas de mot de longueur 6 qui soit parfaitement conservé.
Quelque soit la méthode d’alignement utilisée, il est nécessaire d’évaluer la significativité
des alignements trouvés. Le score d’un alignement n’est en soi pas un critère de décision pour
plusieurs raisons : il dépend du système de score utilisé mais également de la longueur des
séquences comparées. Pour évaluer la significativité d’un alignement, il est donc nécessaire
d’évaluer sa probabilité d’occurrence afin de répondre à la question suivante : quelle était la
probabilité de trouver cet alignement “par hasard” en comparant des séquences ne possédant
aucune homologie a priori ? Cette question suppose que les séquences ont été générées selon un
modèle aléatoire. Les travaux de Karlin et Altschul [KO87, KO88, KA90, KB92] ont permis
d’établir le modèle actuellement utilisé par toutes les méthodes et selon lequel, pour un
système de score fixé, la distribution des scores d’alignements locaux suit une loi de Gumbel
[Gum58]. Dès lors, on est en mesure d’évaluer la significativité d’un alignement de score s en
127
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
1
Graine contiguë (#########)
Graine espacée (###---#-#-##-##)
0.8
0.6
Sensibilité
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Identité
GACTGAACTCAT TAGACTCGACGA
|.||.||||.|| |||.||.|||.|
GGCTAAACTAAT TAGGCTAGACTA
Graine contiguë ####
Graine espacée ##-## ##-##
Fig. 4.11 – Les deux alignements présentent une identité de 9/12 = 75%. On considère deux
graines de poids 4 : la graine contiguë #### et la graine espacée ##-##. La graine contiguë ne
détecte que le premier alignement, alors que la graine espacée détecte les deux. Le symbole #
correspond à une position d’identité et le symbole - à une position quelconque.
128
4.2. L’annotation par génomique comparative
Afin d’évaluer le gain pratique que peuvent apporter les heuristiques à base de graines
espacées, nous avons cherché à comparer systématiquement les résultats de Blast et de
Yass. Nous nous sommes restreints à ces deux logiciels pour plusieurs raisons. Blast reporte
toutes les similitudes locales détectées sous forme de plusieurs alignements, contrairement à
Fasta qui n’en sélectionne qu’une partie pour former un seul et même alignement “global”.
Fasta peu donc être amené à écarter certaines séquences localement similaires qu’il n’arrive
pas à regrouper pour former son alignement, telles que des séquences distantes, répétées,
permutées ou inversées. Des différences minimes séparent Yass et PatternHunter. Pour
comparer les performances pratiques en terme de sensibilité des approches à base de graines
contiguës et espacées, nous avons choisi de travailler sur la comparaison de séquences d’ARN
non-codants car leurs séquences tendent à être moins bien conservées que celles des régions
codantes. A cet effet, nous avons comparé les performances de Blast et de Yass sur deux
jeux de données : les 574 familles d’ARN non-codants de Rfam, et les trois familles d’ARN
non-codants de BRAliBase III (page 81). Le protocole expérimental est identique à celui de
BRAliBase III : chaque séquence de chaque famille est utilisée comme séquence “requête”
pour retrouver ses homologues, soit dans Rfam, soit dans BRAliBase III selon le jeu de
données dont elle est issue. Le compromis sensibilité/spécificité de Yass et de Blast sur Rfam
est présenté en figure 4.12, toutes familles confondues. Ce compromis est calculé en faisant
varier le seuil sur la E-valeur des alignements produits par les deux logiciels. Cette expérience
fait clairement apparaı̂tre que Yass est plus sensible que Blast à spécificité équivalente. La
table 4.1 présente la sensibilité de Blast et de Yass sur chaque famille de BRAliBase III
obtenue en fixant à seuil à 10−4 sur la E-valeur des alignements produits. Cette seconde
expérience confirme les résultats précédents. A spécificité équivalente, Yass détecte quatre
fois plus d’ARN de transfert que Blast, presque deux fois plus d’ARN ribosomiques 5S, et 8%
d’ARN U5 supplémentaires. En terme de temps d’exécution, Blast et Yass sont équivalents.
Cependant, Blast est entre 5 et 10% plus rapide que Yass dans ces expériences simplement
car Yass produit plus d’alignements que Blast. Les résultats obtenus par Yass dans ces
expériences nous ont conduit à le choisir pour la première étape du pipeline.
129
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
Fig. 4.12 – Compromis sensibilité/spécificité de Blast (croix) et de Yass (ronds) sur Rfam,
toutes familles confondues, avec des graines de poids 9. Chaque rond correspond à une
exécution de Yass avec une graine espacée différente de même poids.
Tab. 4.1 – Sensibilité de Blast et de Yass sur les familles d’ARN non-codants de BRAli-
Base III avec des graines de poids 9 . Le seuil sur la E-valeur est ici fixé à 10−4 .
130
4.2. L’annotation par génomique comparative
effectuer les prédictions par analyse comparative sur la séquence à annoter. Chaque groupe de
séquences doit donc comporter un fragment de la séquence à annoter et un nombre suffisant
de séquences similaires. Les alignements deux à deux produits à l’étape précédente impliquent
chacun un fragment de la séquence à annoter et une fragment de séquence provenant de la
banque de données. La constitution des groupes de séquences à partir de ces alignements
revient à identifier sur la séquence à annoter des régions particulièrement conservées, c’est-à-
dire des régions pour lesquelles un nombre significatif de séquences similaires ont été trouvées.
L’identification de ces régions et la génération des groupes de séquences similaires est réalisée
par un logiciel développé dans l’équipe selon un protocole en trois temps. Premièrement,
l’ensemble des alignements deux à deux produits à l’étape précédente sont repositionnés sur
la séquence à annoter, comme illustré sur la figure 4.13. Ensuite, on calcule pour chaque base
de cette séquence son score cumulé de densité grâce auquel on délimite les régions les plus
conservées. Enfin, on extrait les séquences impliquées dans les régions conservées pour former
des groupes de séquences pertinents pour la suite du pipeline.
GACAACCGAAACTCG
ACGACAACCGAA
GATACGACAAC
Alignements 2 à 2 CCGCCCACATCTGCGAGGGTA TCGGA−−−GACAA
TCCGTCCGCACCTTCGCGGATATTC CGATCGGATACGAC
ATCCGGCCGCAT−−−−−−GGCTATTCTCAGC CCCGATCGGATACGA
Séquence S ...CGGATCCGACCACATCTGCGCGGGTATTCTCAGCGA...ACCCGATCGGATACGACAACCGAAACTCGACTCAGCTCACCC...
Fig. 4.13 – Repositionnement des alignements deux à deux sur la séquence à annoter.
Chaque position i de la séquence S à annoter est caractérisée par une densité en aligne-
ments, notée di , donné par la formule suivante
ni
di =
A
où ni est le nombre d’alignements présents à la position i de la séquence S et A est le
nombre moyen d’alignements par base de S. La valeur de di représente la quantité d’aligne-
ments qui impliquent la base en position i de la séquence S par rapport au nombre d’ali-
gnements moyens par base pour la séquence S. Par nature, la densité en alignements d’une
position ne dépend pas de celle des positions voisines. L’interprétation de sa valeur dépend
donc du contexte dans lequel on l’observe et ne représente donc pas un indicateur suffisant
pour repérer des régions denses en alignements. Afin de simplifier l’identification des régions
conservées, on calcule pour chaque position i de S son score cumulé de densité Di selon la
relation suivante
0
Di = max
Di−1 + log( dλi )
où λ est un paramètre strictement supérieur à 1 qui permet d’ajuster le seuil de détection
des régions conservées. La valeur de Di inclus non seulement la densité en alignements de
la position i mais également le score Di−1 de la position précédente. Ce dernier élément est
très important car il introduit une certaine inertie dans les variations du score qui va nous
permettre de détecter simplement les régions denses en alignements. La figure 4.14 illustre la
131
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
variation du score cumulé sur un exemple fictif. Cet exemple fait apparaı̂tre trois pics du score
cumulé qui correspondent chacun à une région particulièrement dense en alignements. On
définit une région conservée par un intervalle de positions [i; j] sur la séquence S. Cet intervalle
est déterminé de la manière suivante : Dj est un maximum global du score cumulé et i est la
plus grande valeur inférieure à j telle que Di−1 = 0. Sur l’exemple de la figure 4.14, les trois
zones grisées correspondent aux trois régions conservées identifiées définies par les intervalles
de positions [11; 26], [50; 83] et [122; 130] sur la séquence S. Pour cet exemple, le paramètre λ
est fixé à 1, 3. Dans les faits, ce paramètre permet de jouer sur la sensibilité de détection des
régions conservées en modifiant l’amplitude de la valeur du score cumulé de densité. Lorsque
λ augmente les valeurs du score diminuent, et inversement. Cependant, comme la valeur du
score est majorée par zéro, plus λ augmente moins on observe de maximums, et inversement.
Augmenter la valeur de λ permet donc de diminuer le nombre de régions conservées détectées ;
diminuer sa valeur permet au contraire d’augmenter le nombre de régions détectées.
Alignements 2 à 2
Séquence S
10
7
Score cumulé de densité
0
0 20 40 60 80 100 120 140 160
Position sur S
Lorsque les intervalles des régions conservées sont déterminés, les séquences correspon-
dantes des alignements qui intersectent ces intervalles sont extraites. Il arrive que les aligne-
ments soient incomplets, notamment lorsque les génomes contiennent des séquences similaires
132
4.2. L’annotation par génomique comparative
faiblement conservées. Sur l’exemple de la figure 4.14, la seconde région identifiée correspond
à une accumulation de plusieurs alignements dont la majorité ne couvrent que partiellement
le fragment de S en cause. Pour régler ce problème a posteriori, on étend de part et d’autre
chaque séquence de la région afin d’obtenir une séquence de longueur identique au fragment
de S. C’est étape n’est possible que si l’on dispose dans la banque de données du contexte
des séquences à étendre.
A ce niveau du pipeline, plusieurs traitements peuvent être appliqués aux groupes de
séquences détectés. Ces traitements dépendent essentiellement des méthodes de prédictions
auxquelles on souhaite les soumettre. Si l’on souhaite les soumettre à Protea et/ou à caR-
NAc pour prédire des séquences codantes ou d’ARN non-codants homologues, aucun trai-
tement particulier supplémentaire n’est requis. Ces deux logiciels travaillent en effet sur des
séquences non alignées et intègrent un pré-traitement pour éliminer les séquences redondantes.
Si toutefois l’on souhaite utiliser d’autres méthodes basées sur une analyse comparative, il
convient de vérifier certaines propriétés sur les groupes de séquences. Notamment pour les
méthodes qui s’appuient sur un alignement, il est nécessaire de s’assurer qu’il est possible de
construire un alignement fiable. Dans cette optique, on propose d’épurer chaque groupe de
séquences en éliminant les séquences trop divergentes selon un procédé strictement analogue
à celui mis en œuvre dans Protea et caRNAc pour construire les méta-séquences (sec-
tion 1.5.3). On propose également de filtrer les groupes de séquences en fonction du nombre
de séquences qu’ils contiennent.
133
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
Les régions inter-géniques du génome d’Escherichia coli sont déterminées à partir des 115
gènes à ARN et des 4 290 gènes codants annotés. Seules les régions dont la longueur dépasse
cinquante nucléotides sont retenues, ce qui représente au total 2 367 séquences couvrant 500
kilobases. La longueur moyenne de ces séquences est de 211 nucléotides, et la séquence la plus
longue fait 1 729 nucléotides. Quatre gènes à ARN n’ont volontairement pas été exclus pour
fournir un contrôle positif.
Les régions inter-géniques ainsi déterminées sont comparées aux séquences génomiques
complètes de quatre organismes proches d’Escherichia coli en terme d’évolution :
– Klebsiella pneumoniae, souche 342 ;
– Salmonella enterica enteridis, souche PT4 ;
– Salmonella enterica serovar Paratyphi A, souche AKU 12601 ;
– Salmonella enterica serovar Typhi, souche CT18.
Les comparaisons sont réalisées par Blast et les alignements obtenus filtrés pour
répondre à trois critères : une E-valeur inférieure à 0,01, une longueur supérieure à cinquante
nucléotides, et un pourcentage d’identité supérieur à 65%. Ces critères sont fixés pour fournir
des alignements pertinents à Qrna. Au total, 23 674 alignements sont produits par Blast,
dont plus de la moitié proviennent de Salmonella enterica serovar Typhi.
Pour optimiser le traitement par Qrna, les alignements dont la longueur dépasse deux
cents nucléotides sont découpés en fragments de deux cents nucléotides qui se chevauchent de
cinquante nucléotides. Parmi ces alignements, Qrna prédit 556 couples de séquences d’ARN
non-codants homologues. Ces 556 candidats contiennent les quatre gènes à ARN laissés vo-
lontairement. Parmi ces candidats, 281 correspondent à des éléments reconnus a posteriori
comme n’étant pas des ARN non-codants mais des éléments possédant une structure se-
condaire caractéristique tels que des terminateurs de gènes, des séquences répétées et des
séquences régulatrices. Parmi les 275 candidats restants, 49 ont été choisis manuellement en
fonction de l’aspect général de la structure secondaire prédite et de la proximité avec des
gènes connus par ailleurs. Après vérification expérimentale par Nothern Blot, il apparaı̂t que
11 des 49 candidats retenus sont effectivement transcrits en ARN de longueurs inférieures à
quatre cents nucléotides, et 6 semblent faire partie d’ARN plus longs supposés être des ARN
messagers. Les 32 candidats restants n’apparaissent pas être transcrits, au moins dans les
conditions expérimentales observées.
134
4.2. L’annotation par génomique comparative
taille des régions à 1 000 nucléotides. A l’issue de cette étape, 309 groupes de séquences sont
constitués, dont 113 contiennent au moins un fragment des 171 gènes présents à l’origine :
22 ARN ribosomiques, 80 des 89 ARN de transfert et 11 des 60 autres types de gène à ARN
restants.
Fig. 4.15 – Arbre phylogénétique reliant les deux groupes de quatre organismes utilisés dans
les comparaisons avec Escherichia coli.
Les résultats que nous avons obtenus sont les suivants. 1 844 alignements sont générés par
Yass à la première étape du pipeline, et 71 groupes de séquences conservées sont détectés. Au
sein de ces groupes de séquences on retrouve 87 des 171 gènes à ARN présents dans le génome
d’Escherichia coli : les 22 ARN ribosomiques, 64 des 89 ARN de transfert et seulement 1 des
60 autres types de gènes à ARN. Parmi les 71 groupes de séquences, 65 présentent au moins
une tige conservée par toutes les séquences détectée par caRNAc. Ces 65 groupes intersectent
exactement la même quantité de gènes à ARN qu’à l’étape précédente, soit 87 gènes à ARN.
Le volume de données circulant dans le pipeline est drastiquement réduit par rapport à
l’expérience précédente, sans pour autant que la quantité de gènes à ARN détectés en sortie
135
Chapitre 4. Deux exemples d’intégration de Protea et caRNAc
soit réduite dans les mêmes proportions. En effet, près de 27 fois moins d’alignements deux à
deux sont produits au cours de la première étape et quatre fois moins de régions conservées
sont identifiées, mais on retrouve néanmoins 92% des gènes à ARN déjà prédits à partir des
séquences proches.
136
Conclusion
137
Conclusion
Markov cachés pour réaliser l’alignement des séquences nucléiques en travaillant au niveau des
acides aminés codés et en supportant les changements éventuels du cadre de lecture. Il serait
également intéressant de travailler plus finement sur les mutations silencieuses et synonymes
attendues et observées entre les séquences. Par exemple, il serait intéressant de proposer
système de bonification pour les mutations silencieuses afin de prendre en compte les acides
aminés identiques produits par des codons différents.
138
Toujours selon le principe de fenêtre glissante, il serait particulièrement intéressant de
proposer un mode de fonctionnement incrémental pour Protea et caRNAc. Partant d’un
ensemble de séquences homologues, codantes ou structurées, les comparaisons de tous les
couples de séquences sont réalisées une seule fois, et le graphe obtenu stocké. A l’aide du
fenêtre glissante, on balaye ensuite un génome à la recherche d’une séquence homologue à
l’ensemble de séquences pré-traitées, sans refaire de calculs inutiles. Dans l’idée, ce mode de
fonctionnement est équivalent aux modèles de covariances utilisés dans des méthodes comme
Infernal ou Erpin pour la détection de structures conservées.
Concernant le concept de méta-séquence que nous avons introduit et mis en œuvre dans
Protea et caRNAc, celui-ci pourrait être affiné. Actuellement, les séquences fortement si-
milaires sont regroupées et représentées par un alignement multiple, tandis que les séquences
“uniques” restent inchangées. Ce regroupement réalisé de manière binaire en fonction du pour-
centage d’identité pourrait être complété en intégrant des informations phylogénétiques. Ces
informations permettraient de pondérer les comparaisons entre méta-séquences dans Pro-
tea et caRNAc. Ce procédé déjà appliqué dans des méthodes comme ExoniPhy pour
la prédiction de séquences codantes ou d’EvoFold pour la prédiction d’ARN non-codants
semble contribuer à améliorer les résultats notamment entre les séquences qui présentent peu
de mutations.
139
Conclusion
140
Bibliographie
[AGM+ 90] Stephen F. Altschul, Warren Gish, Webb Miller, Eugene W. Myers, and Da-
vid J. Lipman. Basic local alignment tool. Journal of Molecular Biology,
215(3) :403–410, October 1990.
[AS05] Julien Allali and Marie-France Sagot. A Multiple Graph Layers Model with
Application to RNA Secondary Structures Comparison. String Processing and
Information Retrieval, pages 348–359, 2005.
[BAB+ 04] Ewan Birney, T. Daniel Andrews, Paul Bevan, Mario Caccamo, Yuan Chen,
Laura Clarke, Guy Coates, James Cuff, Val Curwen, Tim Cutts, Thomas
Down, Eduardo Eyras, Xose M. Fernandez-Suarez, Paul Gane, Brian Gib-
bins, James Gilbert, Martin Hammond, Hans-Rudolf Hotz, Vivek Iyer, Kers-
tin Jekosch, Andreas Kahari, Arek Kasprzyk, Damian Keefe, Stephen Kee-
nan, Heikki Lehvaslaiho, Graham McVicker, Craig Melsopp, Patrick Meidl,
Emmanuel Mongin, Roger Pettett, Simon Potter, Glenn Proctor, Mark Rae,
Steve Searle, Guy Slater, Damian Smedley, James Smith, Will Spooner, Arne
Stabenau, James Stalker, Roy Storey, Abel Ureta-Vidal, K. Cara Woodwark,
Graham Cameron, Richard Durbin, Anthony Cox, Tim Hubbard, and Mi-
chele Clamp. An overview of Ensembl. Genome Research, 14(5) :925–928,
2004. doi:10.1101/gr.1860604.
[BAW+ 05] Amos Bairoch, Rolf Apweiler, Cathy H. Wu, Winona C. Barker, Brigitte
Boeckmann, Serenella Ferro, Elisabeth Gasteiger, Hongzhan Huang, Rodrigo
Lopez, Michele Magrane, Maria J. Martin, Darren A. Natale, Claire O’Do-
novan, Nicole Redaschi, and Lai-Su L. Yeh. The Universal Protein Re-
source (UniProt). Nucleic Acids Research, 33(Suppl 1) :D154–159, 2005.
doi:10.1093/nar/gki070.
[BCD04] Ewan Birney, Michele Clamp, and Richard Durbin. GeneWise and Genome-
Wise. Genome Research, 14(5) :988–995, 2004. doi:10.1101/gr.1865504.
[BG96] Moisès Burset and Roderic Guigo. Evaluation of gene struc-
ture prediction programs. Genomics, 34(3) :353–367, June 1996.
doi:10.1006/geno.1996.0298.
[BGH+ 98] Winona C. Barker, John S. Garavelli, Daniel H. Haft, Lois T. Hunt, Chris-
topher R. Marzec, Bruce C. Orcutt, Geetha Y. Srinivasarao, Lai-Su L. Yeh,
Robert S. Ledley, Hans-Werner Mewes, Friedhelm Pfeiffer, and Akira Tsugita.
The PIR-International Protein Sequence Database. Nucleic Acids Research,
26(1) :27–32, 1998. doi:10.1093/nar/26.1.27.
141
Bibliographie
[BH88] Robert E. Bruccoleri and Gerhard Heinrich. An improved algorithm for nu-
cleic acid secondary structure display. Computational Applications in Bios-
ciences, 4(1) :167–173, 1988. doi:10.1093/bioinformatics/4.1.167.
[BH00] Vineet Bafna and Daniel H. Huson. The conserved exon method for gene fin-
ding. Proceedings of the 8th International Conference on Intellignet Systems
for Molecular Biology ISMB, 8 :3–12, 2000.
[BHGP94] James W. Brown, Elizabeth S. Haas, Donald G. Gilbert, and Norman R. Pace.
The Ribonuclease P Database. Nucleic Acids Research, 22(17) :3660–3662,
1994.
[BK97] Christopher B. Burge and Samuel Karlin. Prediction of complete gene struc-
tures in human genomic DNA. Journal of Molecular Biology, 268(1) :78–94,
1997. doi:10.1006/jmbi.1997.0951.
[BK06] Rajnish Bharadwaj and Alex L. Kolodkin. Descrambling Dscam diversity.
Cell, 125(3) :421–424, May 2006. doi:10.1016/j.cell.2006.04.012.
[BKR+ 04] Mathieu Blanchette, W. James Kent, Cathy Riemer, Laura Elnitski,
Arian F.A. Smit, Krishna M. Roskin, Robert Baertsch, Kate Rosenbloom,
Hiram Clawson, Eric D. Green, David Haussler, and Webb Miller. Aligning
multiple genomic sequences with the threaded blockset aligner. Genome Re-
search, 14(4) :708–723, April 2004.
[BKR07] Markus Bauer, Gunnar W. Klau, and Knut Reinert. Accurate multiple
sequence-structure alignment of RNA sequences using combinatorial optimi-
zation. BMC Bioinformatics, 8 :271, 2007. doi:10.1186/1471-2105-8-271.
[BKV96] Bernard Billoud, Milutin Kontic, and Alain Viari. Palingol : a declarative
programming language to describe nucleic acids’ secondary structures and
to scan sequence database. Nucleic Acids Research, 24(8) :1395–1403, April
1996.
[BLT93] Mark S. Boguski, Todd M. Lowe, and Carolyn M. Tolstoshev. dbEST –
database for ”expressed sequence tags”. Nature Genetics, 4 :332–333, 1993.
doi:10.1038/ng0893-332.
[BPM+ 00a] Serafim Batzoglou, Lior Pachter, Jill P. Mesirov, Bonnie Berger, and Eric S.
Lander. Human and Mouse Gene Structure : Comparative Analysis and
Application to Exon Prediction. Genome Research, 10(7) :950–958, 2000.
doi:10.1101/gr.10.7.950.
[BPM+ 00b] Serafim Batzoglou, Lior Pachter, Jill P. Mesirov, Bonnie Berger, and Eric S.
Lander. Human and mouse gene structure : comparative analysis and ap-
plication to exon prediction. In Proceedings of the 4th Annual Internatio-
nal Conference on Computational Molecular Biology RECOMB, pages 46–53,
2000.
[Bro99] Michael P.S. Brown. RNA modeling using stochastic context-free grammars.
PhD thesis, University of California, Santa Cruz, 1999.
[BRS03] Philippe Blayo, Pierre Rouzé, and Marie-France Sagot. Orphan gene finding :
an exon assembly approach. Theoritical Computer Science, 290(3) :1407–1431,
January 2003. doi:10.1016/S0304-3975(02)00043-9.
142
[BT06] Guillaume Blin and Hélène Touzet. How to Compare Arc-Annotated Se-
quences : The Alignment Hierarchy. In String Processing and Information
Retrieval (SPIRE), volume 4209 of Lecture Notes in Computer Science, pages
291–303. Springer Berlin / Heidelberg, 2006. doi:10.1007/11880561_24.
[BWRVdP04] Eric Bonnet, Jan Wuyts, Pierre Rouzé, and Yves Van de Peer. Evidence that
microRNA precursors, unlike other non-coding RNAs, have lower folding free
energies than random sequences. Bioinformatics, 20(17) :2911–2917, 2004.
[BZ04] Vineet Bafna and Shaojie Zhang. FastR : fast database search tool for non-
coding RNA. Proceedings of the IEEE Computer Society Bioinformatics
Conference (CSB’04), pages 52–61, 2004.
[CBG+ 05] Liu Changning, Bai Baoyan, Skogerbø Geir, Cai Lun, Deng Wei, Zhang Yong,
Bu Dongbo, Zhao Yi, and Chen Runsheng. NONCODE : an integrated know-
ledge database of non-coding RNAs. Nucleic Acids Research, 33(Database
issue) :D112–115, 2005.
[CDH01] Richard J. Carter, Inna Dubchak, and Stephen R. Holbrook. A computational
approach to identify genes for functional RNAs in genomic sequences. Nucleic
Acids Research, 29(19) :3928–3938, 2001.
[CFKK05] Peter Clote, Fabrizio Ferré, Evangelos Kranakis, and Danny Krizanc. Structu-
ral RNA has lower folding energy than random RNA of the same dinucleotide
frequency. RNA, 11 :578–591, 2005.
[CK91] David K. Y. Chiu and Ted Kolodziejczak. Inferring consensus structure from
nucleic acid sequences. Computational Applications in Biosciences, 7(3) :347–
352, July 1991.
[CKB04] Alex Coventry, Daniel J. Kleitman, and Bonnie Berger. MSARi : Multiple se-
quence alignments for statistical detection of RNA secondary structure. Pro-
ceedings of the National Academy of Sciences of the United States of America,
101(33) :12102–12107, 2004.
[Con08] The UniProt Consortium. The universal protein resource (UniProt).
Nucleic Acids Research, 36(Database issue) :D190–195, January 2008.
doi:10.1093/nar/gkm895.
[Cor88] Florence Corpet. Multiple sequence alignment with hierarchical clustering.
Nucleic Acids Research, 16(22) :10881–10890, 1988.
[CP06] Sourav Chatterji and Lior Pachter. Reference based annotation with Gene-
Mapper. Genome Biology, 7(R29), 2006. doi:10.1186/gb-2006-7-4-r29.
[CPL+ 07] Ségolène Caboche, Maude Pupin, Valérie Leclère, Arnaud Fontaine, Philippe
Jacques, and Gregory Kucherov. NORINE : a database of nonribosomal pep-
tides. Nucleic Acids Research, 2007. doi:10.1093/nar/gkm792.
[Cri70] Francis Crick. Central Dogma of Molecular Biology. Nature, 227 :561–563,
August 1970. doi:10.1038/227561a0.
[Cro97] James F. Crow. The high spontaneous mutation rate : is it a health risk ? Pro-
ceedings of the National Academy of Sciences of the United States of America,
94(16) :8380–8386, August 1997.
143
Bibliographie
144
[DWMS06] Deniz Dalli, Andreas Wilm, Indra Mainz, and Gerhard Steger. STRAL :
progressive alignment of non-coding RNA using base pairing probability
vectors in quadratic time. Bioinformatics, 22(13) :1593–1599, July 2006.
doi:10.1093/bioinformatics/btl142.
[ED94] Sean R. Eddy and Richard Durbin. RNA sequence analysis using covariance
models. Nucleic Acids Research, 22(11) :2079–2088, June 1994.
[Edd01] Sean R. Eddy. Non-coding RNA genes and the modern RNA world. Nature
Reviews Genetics, 2(12) :919–929, 2001. doi:10.1038/35103511.
[Eng06] Stefan Engelen. Algorithmes pour la prédiction de structures secondaires
d’ARN. PhD thesis, Université d’Evry Val d’Essonne, 2006.
[ET07] Stefan Engelen and Fariza Tahi. Predicting RNA secondary structure by
the comparative approach : how to select the homologous sequences. BMC
Bioinformatics, 8 :464, 2007.
[FBG07] Eva K. Freyhult, Jonathan P. Bollback, and Paul P. Gardner. Exploring ge-
nomic dark matter : a critical assessment of the performance of homology
search methods on noncoding RNA. Genome Research, 17(1) :117–125, Ja-
nuary 2007. doi:10.1101/gr.5890907.
[FdMT08] Arnaud Fontaine, Antoine de Monte, and Hélène Touzet. MAGNO-
LIA : multiple alignment of protein-coding and structural RNA se-
quences. Nucleic Acids Research, 36(Web Server issue) :W14–W18, 2008.
doi:10.1093/nar/gkn321.
[FHL+ 07] Laurent Fousse, Guillaume Hanrot, Vincent Lefèvre, Patrick Pélissier, and
Paul Zimmermann. MPFR : A multiple-precision binary floating-point li-
brary with correct rounding. ACM Transactions on Mathematical Software,
33(2) :13, 2007. doi:http://doi.acm.org/10.1145/1236463.1236468.
[FHS00] Martin Fekete, Ivo L. Hofacker, and Peter F. Stadler. Prediction of RNA base
pairing probabilities on massively parallel computers. Journal of Computa-
tional Biology, 7(1-2) :171–182, 2000. doi:10.1089/10665270050081441.
[FHZ+ 98] Liliane Florea, George Hartzell, Zheng Zhang, Gerald M. Rubin, and Webb
Miller. A computer program for aligning a cDNA sequence with a genomic
DNA sequence. Genome Research, 8(9) :967–974, September 1998.
[Fic95] James W. Fickett. ORFs and genes : how strong a connection ? Journal of
Computational Biology, 2(1) :117–123, 1995.
[FMSB+ 06] Robert D. Finn, Jaina Mistry, Benjamin Schuster-Bockler, Sam Griffiths-
Jones, Volker Hollich, Timo Lassmann, Simon Moxon, Mhairi Marshall, Ajay
Khanna, Richard Durbin, Sean R Eddy, Erik L L Sonnhammer, and Alex Ba-
teman. Pfam : clans, web tools and services. Nucleic Acids Research, 34(Da-
tabase issue) :D247–D251, 2006. doi:10.1093/nar/gkj149.
[FSY+ 99] Yoshifumi Fukunishi, Harukazu Suzuki, Masayasu Yoshino, Hideaki Konno,
and Yoshihide Hayashizaki. Prediction of human cDNA from its homolo-
gous mouse full-length cDNA and human shotgun database. FEBS Letters,
464(3) :129–132, December 1999.
145
Bibliographie
146
[GJGvD+ 06] Sam Griffiths-Jones, Russell J. Grocock, Stijn van Dongen, Alex Bateman,
and Anton J. Enright. miRBase : microRNA sequences, targets and gene
nomenclature. Nucleic Acids Research, 34(Database issue) :D140–4, January
2006. doi:10.1093/nar/gkj112.
[GJMM+ 05] Sam Griffiths-Jones, Simon Moxon, Mhairi Marshall, Ajay Khanna, Sean R.
Eddy, and Alex Bateman. Rfam : annotating non-coding RNAs in com-
plete genomes. Nucleic Acids Research, 33(Database issue) :D121–D124, 2005.
doi:10.1093/nar/gki081.
[GL01] Daniel Gautheret and André Lambert. Direct RNA motif definition and
identification from multiple sequence alignments using secondary structure
profiles. Journal of Molecular Biology, 313(5) :1003–1011, November 2001.
doi:10.1006/jmbi.2001.5102.
[GLL+ 03] Giorgio Grillo, Flavio Licciulli, Sabino Liuni, Elisabetta Sbisa, and Graziano
Pesole. PatSearch : A program for the detection of patterns and structural
motifs in nucleotide sequences. Nucleic Acids Research, 31(13) :3608–3612,
July 2003.
[GMC90] Daniel Gautheret, Francois Major, and Robert Cedergren. Pattern sear-
ching/alignment with RNA primary and secondary structures : an effective
descriptor for tRNA. Computational Applications in Biosciences, 6(4) :325–
331, 1990.
[GMP96] Mikhail S. Gelfand, Andrey A. Mironov, and Pavel A. Pevzner. Gene recog-
nition via spliced sequence alignment. Proceedings of the National Academy
of Sciences of the United States of America, 93(17) :9061–9066, August 1996.
[GS03] Alison P. Galvani and Montgomery Slatkin. Evaluating plague and smallpox
as historical selective pressures for the CCR5-∆32 HIV-resistance allele. Pro-
ceedings of the National Academy of Sciences of the United States of America,
100(25) :15276–15279, 2003. doi:10.1073/pnas.2435085100.
[Gui98] Roderic Guigo. Assembling genes from predicted exons in linear time with dy-
namic programming. Journal of Computational Biology, 5(4) :681–702, 1998.
[Gum58] Emil J. Gumbel. Statistics of extremes. Columbia University Press, 1958.
[GW08] Tanja Gesell and Stefan Washietl. Dinucleotide controlled null models
for comparative RNA gene prediction. BMC Bioinformatics, 9 :248, 2008.
doi:10.1186/1471-2105-9-248.
[GWW05] Paul P. Gardner, Andreas Wilm, and Stefan Washietl. A benchmark of mul-
tiple sequence alignment programs upon structural RNAs. Nucleic Acids
Research, 33(8) :2433–2439, 2005. doi:10.1093/nar/gki541.
[HAZK97] Xiaoqiu Huang, Mark D. Adams, Hao Zhou, and Anthony R. Kerlavage. A
tool for analyzing and annotating genomic sequences. Genomics, 46(1) :37–45,
November 1997. doi:10.1006/geno.1997.4984.
[HBB+ 02] T. Hubbard, D. Barker, E. Birney, G. Cameron, Y. Chen, L. Clark, T. Cox,
J. Cuff, V. Curwen, T. Down, R. Durbin, E. Eyras, J. Gilbert, M. Hammond,
L. Huminiecki, A. Kasprzyk, H. Lehvaslaiho, P. Lijnzaad, C. Melsopp, E. Mon-
gin, R. Pettett, M. Pocock, S. Potter, A. Rust, E. Schmidt, S. Searle, G. Sla-
ter, J. Smith, W. Spooner, A. Stabenau, J. Stalker, E. Stupka, A. Ureta-Vidal,
147
Bibliographie
148
[Ike81a] Toshimichi Ikemura. Correlation between the abundance of Escherichia coli
transfer RNAs and the occurrence of the respective codons in its protein genes.
Journal of Molecular Biology, 146(1) :1–21, February 1981.
[Ike81b] Toshimichi Ikemura. Correlation between the abundance of Escherichia coli
transfer RNAs and the occurrence of the respective codons in its protein
genes : a proposal for a synonymous codon choice that is optimal for the
E. coli translational system. Journal of Molecular Biology, 151(3) :389–409,
September 1981.
[Ike82] Toshimichi Ikemura. Correlation between the abundance of yeast transfer
RNAs and the occurrence of the respective codons in protein genes. Diffe-
rences in synonymous codon choice patterns of yeast and Escherichia coli
with reference to the abundance of isoaccepting transfer RNAs. Journal of
Molecular Biology, 158(4) :573–597, July 1982.
[Jac88] Tyler E. Jacks. Ribosomal frameshifting in retroviral gene expression. PhD
thesis, University of California, 1988.
[JJ98] Jian Jiang and Howard J. Jacob. EbEST : an automated tool using expressed
sequence tags to delineate gene structure. Genome Research, 8(3) :268–275,
March 1998.
[JLMZ02] Tao Jiang, Guohui Lin, Bin Ma, and Kaizhong Zhang. A General Edit Dis-
tance between RNA Structures. Journal of Computational Biology, 9(2) :371–
388, 2002. doi:10.1089/10665270252935511.
[JTZ89] John A. Jaeger, Douglas H. Turner, and Michael Zuker. Improved predictions
of secondary structures for RNA. Proceedings of the National Academy of
Sciences of the United States of America, 86(20) :7706–7710, October 1989.
[JTZ90] John A. Jaeger, Douglas H. Turner, and Michael Zuker. Predicting opti-
mal and suboptimal secondary structure for RNA. Methods in Enzymology,
183 :281–306, 1990.
[KA90] Samuel Karlin and Stephen F. Altschul. Methods for assessing the statistical
significance of molecular sequence feature by using general scoring schemes.
Proceedings of the National Academy of Sciences of the United States of Ame-
rica, 87 :2264–2268, 1990.
[KB92] Samuel Karlin and Volker Brendel. Chance and significance in protein and
DNA sequence analysis. Science, 257 :39–49, 1992.
[KBD+ 03] D. Karolchik, R. Baertsch, M. Diekhans, T.S. Furey, A. Hinrichs, Y.T. Lu,
K.M. Roskin, M. Schwartz, C.W. Sugnet, D.J. Thomas, R.J. Weber, D. Hauss-
ler, and W.J. Kent. The UCSC Genome Browser Database. Nucleic Acids
Research, 31(1) :51–54, 2003.
[KC07] Keith Knapp and Yi-Ping Phoebe Chen. An evaluation of contemporary
hidden Markov model genefinders with a predicted exon taxonomy. Nucleic
Acids Research, 35(1) :317–324, 2007. doi:10.1093/nar/gkl1026.
[KE03] Robert J. Klein and Sean R. Eddy. RSEARCH : Finding homologs of single
structured RNA sequences. BMC Bioinformatics, 4(1) :44, 2003.
149
Bibliographie
[KFDB01] Ian Korf, Paul Flicek, Daniel Duan, and Michael R. Brent. Integrating ge-
nomic homology into gene structure prediction. Bioinformatics, 17(suppl
1) :S140–S148, 2001.
[KH94] Ben F. Koop and Leroy Hood. Striking sequence similarity over almost 100 ki-
lobases of human and mouse T-cell receptor DNA. Nature Genetics, 7(1) :48–
53, May 1994. doi:10.1038/ng0594-48.
[KH99] Bjarne Knudsen and Jotun Hein. RNA secondary structure prediction using
stochastic context-free grammars and evolutionary history. Bioinformatics,
15(6) :446–454, June 1999.
[KH03] Bjarne Knudsen and Jotun Hein. Pfold : RNA secondary structure prediction
using stochastic context-free grammars. Nucleic Acids Research, 31(13) :3423–
3428, July 2003.
[KHF+ 04] D. Karolchik, A.S. Hinrichs, T.S. Furey, K.M. Roskin, C.W. Sugnet, D. Hauss-
ler, and W.J. Kent. The UCSC Table Browser data retrieval tool. Nucleic
Acids Research, 32(Suppl 1) :D493–496, 2004.
[KHRE96] David Kulp, David Haussler, Martin G. Reese, and Frank H. Eeckman. A
generalized hidden Markov model for the recognition of human genes in DNA.
Proceedings of the 4th International Conference on Intellignet Systems for
Molecular Biology ISMB, 4 :134–142, 1996.
[KME02] Robert J. Klein, Ziva Misulovin, and Sean R. Eddy. Noncoding RNA genes
identified in AT-rich hyperthermophiles. Proceedings of the National Academy
of Sciences of the United States of America, 99(11) :7542–7547, May 2002.
doi:10.1073/pnas.112063799.
[KMH94] Anders Krogh, I. Saira Mian, and David Haussler. A hidden Markov model
that finds genes in E. coli DNA. Nucleic Acids Research, 22(22) :4768–4778,
November 1994.
[KO87] Samuel Karlin and Friedemann Ost. Counts of long aligned word matches
among random letter sequences. Advances in applied probability, 19 :293–351,
1987.
[KO88] Samuel Karlin and Friedemann Ost. Maximal length of common words among
random letter sequences. Annals of Probability, 16 :535–563, 1988.
[KS01] Sasivimol Kittivoravitkul and Marek Sergot. PAGAN : Predict and Annotate
Genes in genomic sequence based on ANalysis of EST Clusters. In Interna-
tional Conference on Intellignet Systems for Molecular Biology ISMB, 2001.
[KTHB02] Peter S. Klosterman, Makio Tamura, Stephen R. Holbrook, and Steven E.
Brenner. SCOR : a structural classification of RNA database. Nucleic Acids
Research, 30(1) :392–394, 2002.
[KTKA07] Hisanori Kiryu, Yasuo Tabei, Taishin Kin, and Kiyoshi Asai. Murlet : a prac-
tical multiple alignment tool for structural RNA sequences. Bioinformatics,
23(13) :1588–1598, July 2007. doi:10.1093/bioinformatics/btm146.
[KYT+ 07] Taishin Kin, Kouichirou Yamada, Goro Terai, Hiroaki Okida, Yasuhiko Yoshi-
nari, Yukiteru Ono, Aya Kojima, Yuki Kimura, Takashi Komori, and Kiyoshi
150
Asai. fRNAdb : a platform for mining/annotating functional RNA candi-
dates from non-coding RNA sequences. Nucleic Acids Research, 35(Database
issue) :D145–D148, January 2007. doi:10.1093/nar/gkl837.
[LB98] Alexander V. Lukashin and Mark Borodovsky. GeneMark.hmm : new solution
for gene finding. Nucleic Acids Research, 26(4) :1107–1115, 1998.
[LFL+ 04] André Lambert, Jean-Fred Fontaine, Matthieu Legendre, Fabrice Leclerc, Em-
manuelle Permal, François Major, Harald Putzer, Olivier Delfour, Bernard
Michot, and Daniel Gautheret. The ERPIN server : an interface to profile-
based RNA motif identification. Nucleic Acids Research, 32(Web Server is-
sue) :W160–5, July 2004. doi:10.1093/nar/gkh418.
[LGC94] Alain Laferriere, Daniel Gautheret, and Robert Cedergren. An RNA pattern
matching program with enhanced performance and portability. Computatio-
nal Applications in Biosciences, 10(2) :211–212, April 1994.
[LLFG05] Andre Lambert, Matthieu Legendre, Jean-Fred Fontaine, and Daniel Gauthe-
ret. Computing expectation values for RNA motifs using discrete convolutions.
BMC Bioinformatics, 6 :118, 2005. doi:10.1186/1471-2105-6-118.
[LMS+ 95] Jane E. Lamerdin, Mishcelle A. Montgomery, Stephanie A. Stilwagen, Lisa K.
Scheidecker, Robert S. Tebbs, Kerry W. Brookman, Larry H. Thompson, and
Anthony V. Carrano. Genomic sequence comparison of the human and mouse
XRCC1 DNA repair gene regions. Genomics, 25(2) :547–554, January 1995.
[LMZ04] Shu-Yun Le, Jacob V. Jr Maizel, and Kaizhong Zhang. An algorithm for
detecting homologues of known structured RNAs in genomes. Proceedings
of the IEEE Computer Society Bioinformatics Conference (CSB’04), pages
300–310, 2004.
[LP85] David J. Lipman and William R. Pearson. Rapid and sensi-
tive protein similarity searches. Science, 227(4693) :1435–1441, 1985.
doi:10.1126/science.2983426.
[LS79] Michael R. Lerner and Joan A. Steitz. Antibodies to small nuclear RNAs
complexed with proteins are produced by patients with systemic lupus ery-
thematosus. Proceedings of the National Academy of Sciences of the United
States of America, 76(11) :5495–5499, November 1979.
[LTHCB05] Alexandre Lomsadze, Vardges Ter-Hovhannisyan, Yury O. Chernoff, and
Mark Borodovsky. Gene identification in novel eukaryotic genomes by
self-training algorithm. Nucleic Acids Research, 33(20) :6494–6506, 2005.
doi:10.1093/nar/gki937.
[LW01] Neocles B. Leontis and Eric Westhof. Geometric nomenclature and classifica-
tion of RNA base pairs. RNA, 7(4) :499–512, April 2001.
[LWHT05] Jianghui Liu, Jason T. L. Wang, Jun Hu, and Bin Tian. A method for aligning
RNA secondary structures and its application to RNA motif detection. BMC
Bioinformatics, 6 :89, 2005. doi:10.1186/1471-2105-6-89.
[LZP99] Rune B. Lyngsø, Michael Zuker, and Christian N. S. Pedersen. Fast evalua-
tion of internal loops in RNA secondary structure prediction. Bioinformatics,
15(6) :440–445, June 1999.
151
Bibliographie
152
RNA sequences. Nucleic Acids Research, 36(Web Server issue) :W10–W13,
July 2008. doi:10.1093/nar/gkn278.
[MYH+ 08] Toutai Mituyama, Kouichirou Yamada, Emi Hattori, Hiroaki Okida, Yu-
kiteru Ono, Goro Terai, Aya Yoshizawa, Takashi Komori, and Kiyoshi
Asai. The Functional RNA Database 3.0 : databases to support mining
and annotation of functional RNAs. Nucleic Acids Research, October 2008.
doi:10.1093/nar/gkn805.
[MZB96] Wojciech Makalowski, Jinghui Zhang, and Mark S. Boguski. Comparative
analysis of 1196 orthologous mouse and human full-length mRNA and protein
sequences. Genome Research, 6(9) :846–857, September 1996.
[NE07] Eric P. Nawrocki and Sean R. Eddy. Query-dependent banding (QDB) for fas-
ter RNA similarity searches. PLoS Computational Biology, 3(3) :e56, March
2007. doi:10.1371/journal.pcbi.0030056.
[NGM01] Pavel S. Novichkov, Mikhail S. Gelfand, and Andrey A. Mironov. Gene recog-
nition in eukaryotic DNA by comparison of genomic sequences. Bioinforma-
tics, 17(11) :1011–1018, 2001.
[NHH00] Cédric Notredame, Desmond G. Higgins, and Jaap Heringa. T-Coffee : A
novel method for fast and accurate multiple sequence alignment. Journal of
Molecular Biology, 302 :205–217, 2000.
[NJ80] Ruth Nussinov and Ann B. Jacobson. Fast algorithm for predicting the secon-
dary structure of single-stranded RNA. Proceedings of the National Academy
of Sciences of the United States of America, 77(11) :6309–6313, 1980.
[NK05] Laurent Noé and Gregory Kucherov. YASS : enhancing the sensitivity of DNA
similarity search. Nucleic Acids Research, 33(suppl2) :W540–543, 2005.
[NPGK78] Ruth Nussinov, George Piecznik, Jerrold R. Grigg, and Daniel J. Kleitman.
Algorithms for loop matchings. SIAM Journal on Applied Mathematics,
35(1) :68–82, July 1978.
[NW70] Saul B. Needleman and Christian D. Wunsch. A general method applicable to
the search for similarities in the amino acid sequence of two proteins. Journal
of Molecular Biology, 48(3) :443–453, 1970.
[OSKR06] Aleksey Y. Ogurtsov, Svetlana A. Shabalina, Alexey S. Kondrashov, and Mi-
khail A. Roytberg. Analysis of internal loops within the RNA secondary
structure in almost quadratic time. Bioinformatics, 22(11) :1317–1324, June
2006. doi:10.1093/bioinformatics/btl083.
[PAA+ 03] Genis Parra, Pankaj Agarwal, Josep F. Abril, Thomas Wiehe, James W. Fi-
ckett, and Roderic Guigò. Comparative Gene Prediction in Human and Mouse.
Genome Research, 13(1) :108–117, 2003. doi:10.1101/gr.871403.
[PBG00] Genis Parra, Enrique Blanco, and Roderic Guigò. GeneID in Drosophilia.
Genome Research, 10(4) :511–515, April 2000. doi:10.1101/gr.10.4.511.
[PBS+ 06] Jakob Skou Pedersen, Gill Bejerano, Adam Siepel, Kate Rosenbloom, Kers-
tin Lindblad-Toh, Eric S Lander, Jim Kent, Webb Miller, and David Hauss-
ler. Identification and classification of conserved RNA secondary structures
in the human genome. PLoS Computational Biology, 2(4) :e33, April 2006.
doi:10.1371/journal.pcbi.0020033.
153
Bibliographie
154
[RKJE01] Elena Rivas, Robert J. Klein, Thomas A. Jones, and Sean R. Eddy. Compu-
tational identification of noncoding RNAs in E. coli by comparative genomics.
Current Biology, 11(17) :1369–1373, September 2001.
[RMK96] Igor B. Rogozin, Luciano Milanesi, and Nikolay A. Kolchanov. Gene structure
prediction using information on homologous protein sequence. Computational
Applications in Biosciences, 12(3) :161–170, June 1996.
[RMO01] Sanja Rogic, Alan K. Mackworth, and Francis B. F. Ouellette. Evalua-
tion of gene-finding programs on mammalian sequences. Genome Research,
11(5) :817–832, 2001. doi:10.1101/gr.147901.
[RSG07] Jens Reeder, Peter Steffen, and Robert Giegerich. pknotsRG : RNA
pseudoknot folding including near-optimal structures and sliding win-
dows. Nucleic Acids Research, 35(Web Server issue) :W320–4, July 2007.
doi:10.1093/nar/gkm258.
[RSZ04a] Jianhua Ruan, Gary D. Stormo, and Weixiong Zhang. An iterated loop mat-
ching approach to the prediction of RNA secondary structures with pseudok-
nots. Bioinformatics, 20(1) :58–66, January 2004.
[RSZ04b] Jianhua Ruan, Gary D. Stormo, and Weixiong Zhang. ILM : a web server for
predicting RNA secondary structures with pseudoknots. Nucleic Acids Re-
search, 32(Web Server issue) :W146–9, July 2004. doi:10.1093/nar/gkh444.
[Rus93] Peter J. Russell. Fundamentals of Genetics and the Biology Place. Pearson
Education, Limited, 1993.
[Ruv01] G Ruvkun. Molecular biology. Glimpses of a tiny RNA world. Science,
294(5543) :797–799, October 2001. doi:10.1126/science.1066315.
[San85] David Sankoff. Simultaneous solution of the RNA folding, alignment and
protosequence problems. SIAM Journal on Applied Mathematics, 45 :810–
825, 1985. doi:10.1137/0145048.
[SB05] Sven Siebert and Rolf Backofen. MARNA : multiple alignment and
consensus structure prediction of RNAs based on sequence struc-
ture comparisons. Bioinformatics, 21(16) :3352–3359, August 2005.
doi:10.1093/bioinformatics/bti550.
[SC09] Dietmar Schmucker and Brian Chen. Dscam and DSCAM : complex genes
in simple animals, complex animals yet simple genes. Genes & Development,
23(2) :147–156, January 2009. doi:10.1101/gad.1752909.
[Sch02] Peter Schattner. Searching for RNA genes using base-composition statistics.
Nucleic Acids Research, 30(9) :2076–2082, 2002.
[SDKW98] Steven L. Salzberg, Arthur L. Delcher, Simon Kasif, and Owen White. Mi-
crobial gene identification using interpolated Markov models. Nucleic Acids
Research, 26(2) :544–548, 1998.
[Sea92] David B. Searls. The Linguistics of DNA. American Scientist, 80 :579–591,
1992.
[SG94] David J. States and Warren Gish. Combined use of sequence similarity and
codon bias for coding region identification. Journal of Computational Biology,
1(1) :39–50, 1994.
155
Bibliographie
156
[TER03] Fariza Tahi, Stefan Engelen, and Mireille Régnier. A Fast Algorithm for
RNA Secondary Structure Prediction Including Pseudoknots. In Bioinforma-
tic and Bioengineering (BIBE), IEEE International Symposium on, pages 11–
17. IEEE Computer Society, March 2003. doi:10.1109/BIBE.2003.1188924.
[TER05] Fariza Tahi, Stefan Engelen, and Mireille Régnier. P-DCfold or How
to Predict all Kinds of Pseudoknots in RNA Secondary Structures. In-
ternational Journal on Artificial Intelligence Tools, 14(5) :703–716, 2005.
doi:10.1142/S021821300500234X.
[TGR02] Fariza Tahi, Manolo Gouy, and Mireille Régnier. Automatic RNA secondary
structure prediction with a comparative approach. Computers and Chemistry,
26(5) :521–530, July 2002.
[THG94] Julie D. Thompson, Desmond G. Higgins, and Toby J. Gibson. CLUSTAL W :
improving the sensitivity of progressive multiple sequence alignment through
sequence weighting, position-specific gap penalties and weight matrix choice.
Nucleic Acids Research, 22(22) :4673–4680, 1994.
[THG07] Elfar Torarinsson, Jakob H. Havgaard, and Jan Gorodkin. Multiple structural
alignment and clustering of RNA sequences. Bioinformatics, 23(8) :926–932,
April 2007. doi:10.1093/bioinformatics/btm049.
[TKKA08] Yasuo Tabei, Hisanori Kiryu, Taishin Kin, and Kiyoshi Asai. A fast structural
multiple alignment method for long RNA sequences. BMC Bioinformatics,
9 :33, 2008. doi:10.1186/1471-2105-9-33.
[TMM07] Leila Taher, Peter Meinicke, and Burkhard Morgenstern. On splice site pre-
diction using weight array models : a comparison of smoothing techniques.
Journal of Physics : Conference Series, 90 :012004 (8pp), 2007.
[Tou07] Hélène Touzet. Comparative analysis of RNA genes : the caRNAc software.
Methods in Molecular Biology, 395 :465–474, 2007.
[TP04] Hélène Touzet and Olivier Perriquet. CARNAC : folding families of related
RNAs. Nucleic Acids Research, 32(Web Server issue) :W142–5, July 2004.
doi:10.1093/nar/gkh415.
[TPP99] JD Thompson, F Plewniak, and O Poch. BAliBASE : a benchmark alignment
database for the evaluation of multiple alignment programs. Bioinformatics,
15(1) :87–88, 1999. doi:10.1093/bioinformatics/15.1.87.
[TRG+ 03] Leila Taher, Oliver Rinner, Saurabh Garg, Alexander Sczyrba, Michael
Brudno, Serafim Batzoglou, and Burkhard Morgenstern. AGenDA :
homology-based gene prediction. Bioinformatics, 19(12) :1575–1577, 2003.
doi:10.1093/bioinformatics/btg181.
[TSF88] Douglas H. Turner, Naoki Sugimoto, and Susan M. Freier. RNA structure
prediction. Annual Review of Biophysics and Biophysical Chemistry, 17 :167–
192, 1988.
[TTKA06] Yasuo Tabei, Koji Tsuda, Taishin Kin, and Kiyoshi Asai. SCARNA :
fast and accurate structural alignment of RNA sequences by matching
fixed-length stem fragments. Bioinformatics, 22(14) :1723–1729, July 2006.
doi:10.1093/bioinformatics/btl177.
157
Bibliographie
158
[WHN08] Andreas Wilm, Desmond G. Higgins, and Cédric Notredame. R-Coffee : a
method for multiple alignment of non-coding RNA. Nucleic Acids Research,
36(9) :e52, May 2008. doi:10.1093/nar/gkn174.
[WHS05] Stefan Washietl, Ivo L. Hofacker, and Peter F. Stadler. Fast and reliable pre-
diction of noncoding RNAs. Proceedings of the National Academy of Sciences
of the United States of America, 102(7) :2454–2459, 2005.
[WMS06] Andreas Wilm, Indra Mainz, and Gerhard Steger. An enhanced RNA align-
ment benchmark for sequence alignment programs. Algorithms for Molecular
Biology, 1 :19, 2006. doi:10.1186/1748-7188-1-19.
[WPK+ 07] Stefan Washietl, Jakob S. Pedersen, Jan O. Korbel, Claudia Stocsits, An-
dreas R. Gruber, Jorg Hackermuller, Jana Hertel, Manja Lindemeyer, Kris-
tin Reiche, Andrea Tanzer, Catherine Ucla, Carine Wyss, Stylianos E. An-
tonarakis, France Denoeud, Julien Lagarde, Jorg Drenkow, Philipp Kapra-
nov, Thomas R. Gingeras, Roderic Guigo, Michael Snyder, Mark B. Gerstein,
Alexandre Reymond, Ivo L. Hofacker, and Peter F. Stadler. Structured RNAs
in the ENCODE selected regions of the human genome. Genome Research,
17(6) :852–864, June 2007. doi:10.1101/gr.5650707.
[WPVdP04] Jan Wuyts, Guy Perriere, and Yves Van de Peer. The European riboso-
mal RNA database. Nucleic Acids Research, 32(Suppl 1) :D101–103, 2004.
doi:10.1093/nar/gkh065.
[WR04] Zasha Weinberg and Walter L. Ruzzo. Exploiting conserved structure for fas-
ter annotation of non-coding RNAs without loss of accuracy. Bioinformatics,
20 Suppl 1 :i334–41, August 2004. doi:10.1093/bioinformatics/bth925.
[WR06] Zasha Weinberg and Walter L. Ruzzo. Sequence-based heuristics for faster an-
notation of non-coding RNA families. Bioinformatics, 22(1) :35–39, January
2006. doi:10.1093/bioinformatics/bti743.
[WRH+ 07] Sebastian Will, Kristin Reiche, Ivo L Hofacker, Peter F Stadler, and Rolf
Backofen. Inferring noncoding RNA families and classes by means of genome-
scale structure-based clustering. PLoS Computational Biology, 3(4) :e65,
April 2007. doi:10.1371/journal.pcbi.0030065.
[XMU94] Ying Xu, Richard J. Mural, and Edward C. Uberbacher. Constructing gene
models from accurately predicted exons : an application of dynamic program-
ming. Computational Applications in Biosciences, 10(6) :613–623, December
1994.
[XU97] Ying Xu and Edward C. Uberbacher. Automated gene identification in large-
scale genomic sequences. Journal of Computational Biology, 4(3) :325–338,
1997.
[YLB01] Ru-Fang Yeh, Lee P. Lim, and Christopher B. Burge. Computational Inference
of Homologous Gene Structures in the Human Genome. Genome Research,
11(5) :803–816, 2001. doi:10.1101/gr.175701.
[YLLL04] Xiaomin Ying, Hong Luo, Jingchu Luo, and Wuju Li. RDfolder : a web server
for prediction of RNA secondary structure. Nucleic Acids Research, 32(Web
Server issue) :W150–3, July 2004. doi:10.1093/nar/gkh445.
159
Bibliographie
[YWR06] Zizhen Yao, Zasha Weinberg, and Walter L Ruzzo. CMfinder–a covariance
model based RNA motif finding algorithm. Bioinformatics, 22(4) :445–452,
February 2006. doi:10.1093/bioinformatics/btk008.
[ZGS08] Matthias Zytnicki, Christine Gaspin, and Thomas Schiex. DARN ! A Weigh-
ted Constraint Solver for RNA Motif Localization. Constraints, 13(1–2) :91–
109, February 2008. doi:10.1007/s10601-007-9033-9.
[ZS81] Michael Zuker and Patrick Stiegler. Optimal computer folding of large RNA
sequences using thermodynamic and auxiliary information RNA sequences
using thermodynamic and auxiliary information. Nucleic Acids Research,
9(1) :133–148, 1981.
[ZS84] Michael Zuker and David Sankoff. RNA secondary structures and their pre-
diction. Bulletin of Mathematical Biology, 46 :591–621, 1984.
[Zuk89] Michael Zuker. On finding all suboptimal foldings of an RNA molecule.
Science, 244 :48–52, 1989.
[Zyt07] Matthias Zytnicki. Localisation d’ARN non-codants par réseaux de contraintes
pondérées. PhD thesis, Université de Toulouse III - Paul Sabatier, 2007.
160