Back to PBIL home page
Pôle Bioinformatique Lyonnais
Laboratoire de Biométrie et Biologie Evolutive
UMR CNRS 5558
perriere@biomserv.univ-lyon1.fr
Ce TP a pour objet l'étude de l'usage des codons chez certaines bactéries. Pour ce faire, nous allons utiliser les ressources du Pôle Bioinformatique Lyonnais (PBIL) accessibles via l'environnement R ou le site Web.
Le code génétique est dit dégénéré, c'est-à-dire que plusieurs codons peuvent correspondre à un même acide aminé, de tels codons étant dits synonymes. Depuis longtemps, on sait que l'usage des codons synonymes n'est pas le fait du hasard : il n'y a pas équirépartition de leur utilisation à l'intérieur d'une espèce, voire à l'intérieur d'un gène. Chez les bactéries, la composition en codons est dépendante de plusieurs facteurs (expression des gènes, hydrophobicité des protéines, localisation sur le chromosome, etc.) dont l'importance relative varie en fonction de l'espèce considérée (cf. cours).
Des banques de données publiques collectant les séquences génomiques existent depuis le début des années 80 (c'est-à-dire en gros depuis que des méthodes de séquençage rapide de l'ADN sont disponibles). Jusque vers 1994, l'accès à ces banques nécessitait de pouvoir disposer de comptes sur des centres serveurs. Avec le développement d'Internet, il est désormais possible d'interroger ces serveurs en utilisant des logiciels clients. Au PBIL, nous avons développé une interface d'accès aux banques utilisant l'environnement R. Cette interface est disponible dans la bibliothèque seqinr.
De très nombreuses banques de séquences sont disponibles, il existe des banques dédiées à des organismes particuliers (bactéries, levures, vertébrés, etc.), d'autres dédiées à des séquences particulières (immunoglobulines, ARN ribosomiques, etc.) Dans le cadre de ce TP, nous accèderons à une banque dédiée aux génomes procaryotes (bactéries et archées) complètement séquencés : EMGLib.
Tout d'abord, lancez R puis chargez les bibliothèques ade4 et
seqinr :
library(ade4) library(seqinr)
Ensuite, ouvrez la banque EMGLib et composez une première requête qui va nous permettre de déterminer quelles sont les différentes souches d'E. coli dont le génome est disponible :
em <- choosebank("emglib")
query(listname = "ecoli", query = "sp=escherichia coli")
ecoli$req
Quatre séquences correspondant à quatre souches sont disponibles. Dans le cadre de ce TP, nous ne nous intéresserons qu'à la souche K12. Afin de déterminer à quelle séquence correspond K12, il faut accéder aux annotations :
sapply(ecoli$req, getAnnot, nbl = 8)
On voit que c'est l'entrée de nom ECOLICG qui contient la séquence
d'E. coli K12. Afin de récuperer l'intégralité des séquences
codantes figurant dans cette entrée, il suffit alors de composer la
requête :
query(listname = "ecoli", query = "n=ECOLICG et t=cds")
La liste ecoli contient désormais les noms des 4290 CDS
d'E. coli K12. C'est le critère t=cds qui permet de
récuperer les parties codantes. En effet, en anglais, partie codante se dit
Coding DNA Sequences, d'ou l'abréviation CDS utilisée dans la
requête.
Nous allons maintenant sélectionner un sous-ensemble de gènes d'E. coli K12 dont on sait qu'ils sont hautement exprimés, ceci afin d'établir la table du facteur d'adaptation des codons (cf. cours). Une fois que ceci sera fait, il sera possible de calculer les valeurs du Codon Adaptation Index (CAI) pour l'ensemble des gènes. Ces valeurs constitueront un estimateur du niveau d'expression des gènes.
Pour composer la requête, il est possible de réutiliser la liste complète des CDS que nous avons construite auparavant. Par ailleurs, on sait que les gènes codant pour les protéines ribosomales sont toujours fortement exprimés chez les bactéries. Ce sont donc ces gènes que nous utiliserons. Pour les récuperer, la première étape consiste en la composition de la requête suivante :
query(listname = "rib", query = "ecoli et k=@ribosomal@protein@")
Le critère k=@ribosomal@protein@ permet de ne récupérer dans la
liste complète des CDS que ceux qui sont associés aux mots-clés contenant la
chaîne de caractères @ribosomal@protein@. Le caractère
@ est un joker, c'est-à-dire qu'il peut remplacer n'importe
quelle chaîne. La question qui se pose maintenant est de savoir si tous les
CDS retournés par cette requête codent bien pour des protéines ribosomales.
Pour le savoir, il faut une fois de plus accéder aux annotations :
sapply(rib$req, getAnnot, nbl = 10)
Tous les gènes codant pour des protéines ribosomales ont un nom qui commence par
les lettres rp. On voit donc qu'un certain nombre de CDS de la
liste ne codent pas pour de telles protéines mais leur sont associés, comme par
exemple les gènes prmA ou rimI :
getAnnot(rib$req[[20]], nbl = 10) getAnnot(rib$req[[60]], nbl = 10)
Il convient donc d'éliminer ces gènes de notre liste puisque rien ne garantit qu'ils soient effectivement hautement exprimés. Pour ce faire, on peut utiliser la requête suivante :
query(listname = "rib", query = "ecoli et k=@ribosomal@protein@ et no k=rim@ et no k=@methyltransferase@")
Lorsque l'on veut récupérer des gènes appartenant à une famille bien précise, il faut toujours – après une requête basée sur l'emploi de mots-clés – veiller à élimininer les entrées pour lesquelles les annotations montrent que l'on est pas en présence d'un représentant de la dite famille.
En utilisant les séquences des gènes de protéines ribosomales obtenues par la
requête précédente nous allons calculer les valeurs du facteur d'adaptation
des codons. Pour cela, il est nécessaire de définir la fonction
adaptf qui va récupérer les séquences proprement dites, calculer
les effectifs en codons puis le logarithme du rapport entre l'effectif d'un
synonyme et l'effectif du synonyme le plus abondant pour un acide aminé donné.
Dans le cas ou l'effectif pour un codon est égal à 0, on remplace cette
valeur par 0.5 afin de pouvoir calculer le logarithme :
adaptf <- function(q)
{
l <- 0
sum <- vector(mode = "numeric", length = 64)
for (i in seq(from = 1, to = q$nelem))
{
s <- getSequence(q$req[[i]])
l <- l + length(s)/3
sequence <- splitseq(seq = s, frame = 0, word = 3)
eff <- table(factor(sequence, levels = SEQINR.UTIL$CODON.AA$CODON))
sum <- sum + eff
}
for (i in seq(from = 1, to = 64))
{
if (sum[i] == 0)
sum[i] <- 0.5
}
T <- split(sum, SEQINR.UTIL$CODON.AA$AA)
w <- lapply(T, function(x)
{
return(log(x/max(x)))
})
names(w) <- NULL
w <- unlist(w)[as.character(SEQINR.UTIL$CODON.AA$CODON)]
return(w)
}
Une fois cette fonction définie, on l'applique à la liste contenant les séquences codant pour les protéines ribosomales :
w <- adaptf(rib)
Il va être nécessaire de calculer les effectifs en codons pour chaque gène d'E. coli K12. Pour ce faire, vous utiliserez la fonction suivante :
seqecoli <- lapply(ecoli$req, function(x)
{
tmp <- getSequence(x)
tab <- uco(tmp, index = "eff")
rm(tmp)
return(as.vector(tab))
})
De plus, il faut éliminer du jeu de données les séquences contenant des codons
stop en phase (au cas où cela se produise) ainsi que les séquences trop courtes
(moins de 50 codons). La nécessité d'une taille minimum des CDS est liée
au fait qu'il peut exister des variations stochastiques importantes sur des gènes
de petite taille (par exemple ceux codant pour les peptides leader de
certains opérons). Pour ce faire, vous pourrez utiliser la fonction
cleanup :
cleanup <- function(tab, n)
{
aa <- translate(sapply(rownames(tab), s2c), numcode = n)
tab <- tab[,colSums(tab[which(aa == "*"),]) == 1]
tab <- tab[,colSums(tab) > 50]
return(tab)
}
La construction du tableau d'usage des codons se fait à partir des commandes suivantes :
tecoli <- as.data.frame(seqecoli, row.names = SEQINR.UTIL$CODON.AA$CODON) names(tecoli) <- ecoli$req tecoli <- as.data.frame(t(cleanup(tecoli, 1))) length(ecoli$req)
Les séquences éliminées par la fonction cleanup doivent ensuite
être retirées de la liste ecoli$req :
ecoli$req <- ecoli$req[which(as.vector(ecoli$req) %in% row.names(tecoli))] length(ecoli$req)
Le nombre total de séquences figurant dans la liste doit alors être égal à 4248.
Il est maintenant possible de récupérer un certain nombre de valeurs d'indices d'usage du code comme le taux de G+C des gènes :
gc <- sapply(ecoli$req, function(x)
{
tmp <- getSequence(x)
tab <- GC(tmp)
rm(tmp)
return(tab)
})
Le taux de G+C en position 3 des codons :
gc3 <- sapply(ecoli$req, function(x)
{
tmp <- getSequence(x)
tab <- GC3(tmp)
rm(tmp)
return(tab)
})
L'indice de Kyte et Doolittle pour les protéines :
kd <- sapply(ecoli$req, function(x)
{
tmp <- getSequence(x)
tab <- uco(tmp, index = "freq") %*% EXP$KD
rm(tmp)
return(tab)
})
Le CAI, calculé à partir des valeurs du facteur d'adaptation que nous avons précedemment établies :
CAI <- function(q, w)
{
sapply(q$req, function(x)
{
tmp <- getSequence(x)
tab <- uco(tmp, index = "freq") %*% w
rm(tmp)
return(exp(tab))
})
}
cai <- CAI(ecoli, w)
Questions : Une fois que ces différents indices d'usage du code ont été calculés, examinez leur distribution sous la forme d'histogrammes. Eventuellement sauvegardez les graphiques car ces distributions seront utilisées plus tard. Ecrivez ensuite le code R vous permettant de calculer les valeurs cumulées de G+C3 et de CAI centrées par leurs moyennes respectives (calculées sur l'ensemble des gènes). Faites un graphe superposant les deux courbes montrant la variation de ces valeurs centrées-cumulées en fonction de la position des gènes sur le chromosome, ceci en sachant que :
query a
retourné la liste des CDS dans l'ordre de leur position sur le chromosome.
par(new = TRUE) autorise la superposition d'un graphique
à un autre précedemment tracé.
Nous allons calculer une Analyse Factorielle des Correspondances (AFC) sur le
tableau contenant les effectifs des codons. Ne conservez dans cette analyse que
les facteurs qui vous semblent pertinents en fonction de la décroissance des
valeurs propres. Une fois les calculs effectués, tracez les graphes croisant les
différents facteurs en utilisant la fonction scatter.coa (qui permet
la superposition du plan des individus à celui des variuables). Les arguments
clab permettent d'associer chaque point à un label qui correspond au
nom du gène (clab.row) ou au nom du codon
(clab.col) :
tecoli.coa <- dudi.coa(tecoli, scannf = TRUE) tecoli.coa scatter.coa(tecoli.coa, xax = 1, yax = 2, clab.col = 1, clab.row = 0, posieig = "none") scatter.coa(tecoli.coa, xax = 1, yax = 3, clab.col = 1, clab.row = 0, posieig = "none")
Questions : D'après la nature des codons correspondants, pensez-vous que les gènes ayant les valeurs les plus fortment négatives sur le premier facteur (F1) puissent être fortement exprimés ? Faiblement exprimés ? Obtenus par transfert horizontal ? Par ailleurs, que constatez-vous concernant le troisième facteur (F3) ?
Maintenant tracez le graphe croisant les valeurs de F1 et de CAI :
plot(tecoli.coa$li[,1], cai, xlab = "F1", ylab = "CAI")
Questions : Existe-t-il une relation entre les scores sur F1 et les
valeurs de CAI ? Cette relation est-elle linéaire ? Si la réponse à
cette deuxième question est non, quelle fonction faudrait-il à votre avis
appliquer aux valeurs de CAI pour obtenir une relation linéaire ? En
supposant que vous ayez déterminé cette fonction, appliquez-là aux valeurs
de CAI et faites afficher le graphe croisant les valeurs de F1 et les valeurs
de ƒ(CAI) rangées dans la variable cai2. Calculez ensuite
une régression linéaire entre ces deux valeurs puis superposez la droite de
régression au nuage de points :
plot(tecoli.coa$li[,1], cai2, xlab = "F1", ylab = "f(CAI)") r1 <- lsfit(tecoli.coa$li[,1], cai2) ls.print(r1) abline(r1)
Questions : Qu'en concluez-vous vis-à-vis du biais principal qui affecte la composition en codons chez E. coli K12 ? Vous pouvez vérifier cette hypothèse en regardant les annotations des séquences possèdant des scores élevés sur F1. Pour ce faire, les commandes suivantes vont vous permettre de récuperer la liste des 50 gènes ayant les scores les plus élevés sur ce facteur :
mnemo <- sapply(ecoli$req[1:length(ecoli$req)], getName) x <- data.frame(tecoli.coa$li[,1], mnemo) x <- x[order(x[,1], decreasing = TRUE),] write.table(x[1:50,2], file = "", row.names = FALSE, col.names = FALSE, quote = FALSE)
Pour accéder aux annotations de ces séquences, utilisez le formulaire de
soumission de
liste sur le site du PBIL. Une fois la liste des mnémoniques collée, sélectionnez
la banque EMGLib et cliquez sur le bouton SUBMIT. Une fois la page de
résultats affichée, vous pouvez accéder aux annotations des CDS en cliquant sur les
liens ad hoc. A quelles catégories fonctionnelles appartiennent ces gènes ?
Réiterez cette opération pour les 50 gènes possédant les scores les plus faibles sur
F1 et répondez à cette même question.
Recommencez l'opération de calcul de régression linéaire avec affichage du graphe pour les couples croisant les valeurs F2×G+C3 et F3×KD.
plot(tecoli.coa$li[,2], gc3, xlab = "F2", ylab = "G+C3") r2 <- lsfit(tecoli.coa$li[,2], gc3) ls.print(r2) abline(r2) plot(tecoli.coa$li[,3], kd, xlab = "F3", ylab = "KD") r3 <- lsfit(tecoli.coa$li[,3], kd) ls.print(r3) abline(r3)
Questions : Que pouvez-vous conclure de ces différentes observations quant aux différents biais affectant la composition en codons des gènes d'E. coli K12 ?
Le problème de l'AFC classique est qu'elle mélange les effets portants directement
sur la composition en codons et les effets relatifs à la composition en acides
aminés des protéines. Un des moyens de s'affranchir de cet effet est d'utiliser les
fréquences relatives ou les valeurs du Relative Synonymous Codon Usage
(RSCU). Cependant, comme nous l'avons vu en cours, ces valeurs ne sont pas adaptées
à la méthode et leur usage peut provoquer l'apparition de biais importants. Il est
cependant possible de s'affranchir de l'effet de la composition en acides aminés en
utilisant une variante de l'AFC : l'AFC intra-classe. Cette méthode est
également implémentée dans la bibliothèque ade4. Cette fois, nous
extrairons les trois premiers facteurs de l'analyse sans vérifier la
décroissance des valeurs propres :
facaa <- as.factor(aaa(translate(sapply(names(tecoli), s2c)))) tecoli.coa2 <- dudi.coa(as.data.frame(t(tecoli)), scannf = FALSE, nf = 3) with.dudi <- within(tecoli.coa2, facaa, scannf = FALSE, nf = 3) with.dudi
Questions : Une fois le calcul effectué, regardez la corrélation entre les scores sur le premier facteur de l'AFC intra-classe et les valeurs de CAI ainsi que la corrélation entre les scores sur le troisième facteur avec les valeurs de l'indice de Kyte et Doolittle. En particulier, comparez les résultats avec ceux obtenus avec l'AFC classique. Maintenant, si vous regardez le graphe de la décroissance des valeurs propres, avait-on raison de sélectionner trois facteurs ?
Comme exercice d'application vous allez réitérer toutes les analyses précédemment
effectuées (y compris celles relatives aux variations de G+C3 et de CAI le long du
chromosome) sur les gènes de M. genitalium G37. Par ailleurs, vous
pouvez en plus calculer une AFC basée non plus sur les effectifs en codons, mais
sur les valeurs de RSCU. Le code de la fonction permettant de calculer les valeurs
de RSCU pour une séquence sur figure ci-dessous. Cette fonction s'applique aux
séquences provenant de requêtes seqinr de la même
façon que uco :
RSCU <- function(seq, frame = 0, NA.rscu = 1)
{
sequence <- splitseq(seq = seq, frame = frame, word = 3)
eff <- table(factor(sequence, levels = SEQINR.UTIL$CODON.AA$CODON))
freq <- eff/(floor(length(seq)/3))
T <- split(freq, SEQINR.UTIL$CODON.AA$AA)
rscu <- lapply(T, function(x)
{
return(x/((1/length(x)) * sum(x)))
})
names(rscu) <- NULL
rscu <- unlist(rscu)[as.character(SEQINR.UTIL$CODON.AA$CODON)]
is.na(rscu[!is.finite(rscu)]) <- TRUE
rscu[is.na(rscu)] <- NA.rscu
return(rscu)
}
Remarque : Les bactéries du genre Mycoplasma n'utilisent pas
le code génétique standard, en effet chez ces organismes le tryptophane est codé
par deux codons au lieu d'un, avec UGA en plus de UGG.
Cette particularité devra être prise en compte lorsque vous utiliserez la
fonction cleanup (le numéro pour le code génétique de
Mycoplasma est 4). De même, il ne vous sera pas
possible d'utiliser EXP$KD pour calculer les valeurs de l'indice de
Kyte et Doolittle car ce vecteur est basé sur le code génétique standard. Vous
devrez donc le modifier avant de l'utiliser.
Questions : Une fois que vous aurez effectués tous les calculs sur les gènes de M. genitalium G37, comparez les résultats avec ceux obtenus chez E. coli K12. Quelles sont les différences les plus notables entre les deux espèces ? Y a-t-il des différences au niveau de la distribution des valeurs de G+C3 et du CAI ? Si oui, comment les expliquez-vous ? Les biais d'usage des codons sont-ils identiques chez ces deux bactéries ? Quelles interprétations biologiques peut-on donner de ces différences ?
Comparez également les résultats de l'AFC sur les effectifs, l'AFC sur les RSCU et l'AFC intra-classe, ceci relativement à la corrélation avec les valeurs de CAI. Comment expliquez-vous les différences constatées ? Qu'en concluez-vous relativement à l'utilisation des valeurs de RSCU ?
S'il vous reste du temps, utilisez les outils de ce TP pour identifier les gènes de N. meningitidis Z2491 présentant à la fois des faibles valeurs de CAI et de GC. En accédant aux annotations de ces séquences dans la banque, déterminez quelles sont les catégories fonctionnelles les plus fréquemment rencontrées pour les gènes ayant de faibles valeur de CAI et de G+C3.
Back to PBIL home page