Suivez l'guide :
Fréquences des dinucléotides dans le génome d'organismes modèles

L'analyse de séquences est au cœur de nombreux domaines de la bio-informatique. Le billet du jour s'intéressera aux séquences ADN, en se proposant de compter la fréquence en dinucléotides dans quelques génomes d'organismes modèles (avec une petite arrière-pensée derrière la tête).

Qu'est-ce qu'un dinucléotide ?

L'ADN double brins est classiquement structuré sous forme de double hélice, avec deux brins de direction opposée. À chaque position il y a deux acides nucléiques formant une paire, les appariements classiques étant A en face de T et C en face de G. Ce qu'on appelle dinucléotide est une séquence de deux nucléotides successifs, à ne pas confondre avec une paire de nucléotide. Un petit schéma sera peut-être plus clair :

Figure 1: Séquence double brin GATTACA
Figure 1: Séquence double brin GATTACA. En bleu : un dinucléotide. En rouge : une paire de base.

On peut aussi dire qu'un dinucléotide est un k-mer de longueur (k) égale à 2.

Qu'est-ce que la fréquence en dinucléotide d'une séquence ?

Prenons par exemple cette séquence :

Figure 2: Séquence TGTG
Figure 2: Séquence TGTG

Quels dinucléotides dénombrez vous ? Je vous laisse deux phrases pour réfléchir. C'est bon pour vous ?

On trouve deux fois le nucléotide TG (facile), et une fois le nucléotide GT (au milieu, peut-être l'avez-vous oublié ?). Donc pour une séquence de longueur n, il y aura n-1 dinucléotides. Ce qui donne une fréquence de 2/3 pour le dinucléotide TG et une fréquence de 1/3 pour le dinucléotide GT.

Mais est-ce bien tout ?

Et bien ça dépend. Si c'est un ADN simple brin, oui. Mais s'il s'agit d'un ADN double brin, surprise! il y a d'autres dinucléotides cachés en dessous !

Figure 3: Il y a un ACAC caché sous le TGTG

Sous la séquence du haut, on trouvera donc deux dinucléotides CA et un dinucléotide AC (les dinucléotides sont lus et écrits dans le sens 5' vers 3' par convention). Une séquence double brin de longueur n aura donc 2n-2 dinucléotides.

Cela nous donnerait donc une fréquence de 2/6 pour les dinucléotides TG et CA, et un fréquence de 1/6 pour les dinucléotides GT et CA.

Pour une petite séquence, compter à la main c'est faisable, mais dès que la séquence s’agrandit, autant le faire en utilisant la *bio-informatique* !

Calcul bio-informatique de la fréquence en dinucléotides avec R / Bioconductor

Je vous propose d'utiliser des fonctions existantes du package R / Bioconductor Biostrings. La séquence nucléique doit être présentée sous la forme d'un objet DNAString qui se crée via la fonction DNAString(). On peut ensuite appliquer la fonction dinucleotideFrequency() avec ou sans le paramètre as.prob = TRUE, en fonction de si l'on souhaite des comptages ou des fréquences.

La fonction dinucleotideFrequency() ne calcule que la fréquence du brin rentré, et non la fréquence du brin opposé. Pour la calculer, deux méthodes.

1) La méthode peu efficace mais simple : calculer le contenu en dinucléotide du complément inverse de la séquence et faire la moyenne des deux fréquences :

2) La méthode efficace, mais demandant de coder un peu plus : réassigner les fréquences du brin + aux dinucléotides complémentaires inverses, ce qui évite de tout recalculer.

Voici une petite fonction pour calculer la fréquence en dinucléotide d'une séquence ADN double brins. J'en profite pour faire un format de sortie sous forme de tableau tidy pour nous simplifier la vie un peu plus tard.

Fréquence en dinucléotides du génome de référence de la drosophile

Maintenant que nous savons calculer à la fois pratiquement et théoriquement la fréquence en dinucléotides d'une séquence, attaquons nous à un organisme modèle, au hasard le génome de la dropsophile.

Le génome de la plupart des espèces modèles est déjà disponible sous forme d'objets BSgenome dans Bioconductor, mais pour être plus générique ici je fais un import directement depuis le fichier fasta du génome de référence Ensembl.

On tombe là sur un premier os. Il n'y a pas une séquence dans le génome de la drosophile, mais 1870 séquences ! Cinq autosomes, 2 chromosomes sexuels, un chromosome mitochondrial, et tout un tas de contigs non assemblés. Il nous faut donc modifier légèrement notre fonction d'origine pour qu'elle fonctionne bien avec ce genre d'objets (des DNAStringSet et non de simple DNAString) :

Figure 4 : Fréquences dinucléotidiques du génome de référence de la Drosophile
Figure 4 : Fréquences dinucléotidiques du génome de référence de la Drosophile

Comme on peut le voir, tout les dinucléotides n'ont pas la même fréquence. Mais est-ce un résultat attendu ? Ou pour reformuler la question, quel est la fréquence attendue pour chaque dinucléotide ?

Fréquence observée, fréquence théorique

Comme il y a 16 dinucléotides possibles, on pourrait s'attendre à ce que chaque dinucléotide soit observé avec une fréquence de 1/16 (soit 6,25%). Ce serait le cas si les nucléotides simples avaient une fréquence d'un quart. Hors ce n'est pas le cas en pratique ! S'il y a bien autant de A que de T et autant de C que de G, le contenu en G+C d'un génome peut varier énormément d'une espèce à l'autre (C'est d'ailleurs utilisé pour vérifier les puretés d'assemblages, voir par exemple les BlobPlots) !

On peut donc calculer de manière un peu plus fine la fréquence théorique d'un dinucléotide (TA par exemple, noté p(TA) ) comme la probabilité d'avoir un T en premier ( p(T) ) multiplié par la probabilité d'avoir un A en second ( p(A) ), où p(T) et p(A) correspondent à la fréquence en nucléotide A et T de la séquence étudiée :

p(AT) = p(A) x p(T)

Ou de manière plus générique :

p(N1N2) = p(N1) x p(N2)

On peut utiliser la fonction alphabetFrequency() de Biostrings pour obtenir les fréquences de chaque nucléotide, et ensuite il s'agit de tout faire fonctionner ensemble. Par exemple, et sans détailler :

Ce qui donne :

Figure 5 : Fréquences dinucléotidiques observées et attendues du génome de référence de la Drosophile

Il y a donc plus de AA et TT qu'attendu dans le génome de la drosophile ! On pourrait tester statistiquement cette différence, avec par exemple un test du chi2, mais il y a tellement de nucléotides que la p-value retournée est scandaleusement petite.

Fréquence en dinucléotides du génome de plusieurs espèces modèles

Il ne reste plus qu'à faire tourner ces fonctions sur les génomes de quelques espèces modèles populaires (en étant raisonnable avec l'API d'Ensembl lors de la récupération des génomes de références), et voilà :

Figure 6 : Fréquences en dinucléotides des génomes d’espèces modèles en biologie.

Les plus observateurs et observatrices d’entre vous remarqueront peut-être la disparition des dinucléotides CG dans les génomes des vertébrés. Ce mystère pourra faire l'objet d'un prochain billet si la demande se fait pressante.

Remerciements

Mille mercis à mes relecteurs et relectrices Sébastien Gradit, Sarah Guinchard, Virginie J, et à nos admins chéries !

  • À propos de
  • Après une thèse en cancérologie à Lyon et un postdoc en bioinformatique à Édimbourg, je suis chercheur à INRAE Toulouse depuis fin 2017. Régulation transcriptionnelle et épigénétique. Twitter: @G_Devailly

Catégorie: Suivez l'guide | Tags: , ,

Laisser un commentaire