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 :
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 :
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 !
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.
|
1 2 3 4 5 6 7 |
dinucleotideFrequency(DNAString("GTGT")) #> AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT #> 0 0 0 0 0 0 0 0 0 0 0 2 0 0 1 0 print(dinucleotideFrequency(DNAString("GTGT"), as.prob = TRUE), digits = 3) #> AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT #> 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.667 0.000 0.000 0.333 0.000 |
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 :
|
1 2 3 4 5 |
my_seq <- DNAString("GTGT") (dinucleotideFrequency(my_seq, as.prob = TRUE) + dinucleotideFrequency(reverseComplement(my_seq), as.prob = TRUE)) / 2 #> AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT #> 0.000 0.333 0.000 0.000 0.167 0.000 0.000 0.000 0.000 0.000 0.000 0.333 0.000 0.000 0.167 0.000 |
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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
dsdinfreq <- function(dnastring, ...) { top_strand_freq <- dinucleotideFrequency(dnastring, ...) bottom_strand_freq <- top_strand_freq[c("TT", "GT", "CT", "AT", "TG", "GG", "CG", "AG", "TC", "GC", "CC", "AC", "TA", "GA", "CA", "AA")] freq <- (top_strand_freq + bottom_strand_freq) / 2 data.frame( dinucleotide = names(freq), freq = freq ) } dsdinfreq(DNAString("GTGT"), as.prob = TRUE) #> dinucleotide freq #> AA AA 0.0000000 #> AC AC 0.3333333 #> AG AG 0.0000000 #> AT AT 0.0000000 #> CA CA 0.1666667 #> CC CC 0.0000000 #> CG CG 0.0000000 #> CT CT 0.0000000 #> GA GA 0.0000000 #> GC GC 0.0000000 #> GG GG 0.0000000 #> GT GT 0.3333333 #> TA TA 0.0000000 #> TC TC 0.0000000 #> TG TG 0.1666667 #> TT TT 0.0000000 |
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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
droso_genome_path <- "http://ftp.ensemblgenomes.org/pub/metazoa/release-52/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.32.dna.toplevel.fa.gz" droso_genome <- readDNAStringSet(droso_genome_path) droso_genome DNAStringSet object of length 1870: #> width seq names #> [1] 23513712 CGACAATGCACGACA...GAGAACAGAGAAGAG 2L dna:primary_as... #> [2] 25286936 CTCAAGATACCTTCT...AAACAGAATGAATTC 2R dna:primary_as... #> [3] 28110227 TAGGGAGAAATATGA...TCCATGTTATTCTAT 3L dna:primary_as... #> [4] 32079331 ACGGGACCGAGTATA...GCATTCTAGGAATTC 3R dna:primary_as... #> [5] 1348131 TTATTATATTATTAT...TTTGAGATATATGAA 4 dna:primary_ass... #> ... ... ... #> [1866] 1001 TATACATATATACAT...ATACGGATATACATA 211000022278576 d... #> [1867] 1001 AGCACGACGCGACTG...AACTAACAAAGACCG 211000022279016 d... #> [1868] 1001 TAATACGATTCACCT...CAGGACAAGGCACAC 211000022278495 d... #> [1869] 564 AAGAGAAGAGAAGAG...GTCGACCTGCAGGCA 211000022280427 d... #> [1870] 544 GTGTGCTGGAAAGGT...CGACCTACTCACAGC 211000022279089 d... |
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) :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
dsdinfreq <- function(dnastring) { top_strand_freq <- dinucleotideFrequency(dnastring) |> colSums() bottom_strand_freq <- top_strand_freq[c("TT", "GT", "CT", "AT", "TG", "GG", "CG", "AG", "TC", "GC", "CC", "AC", "TA", "GA", "CA", "AA")] freq <- proportions(top_strand_freq + bottom_strand_freq) data.frame( dinucleotide = names(freq), freq = freq ) } dsdinfreq(droso_genome) #> dinucleotide freq #> AA AA 0.10168661 #> AC AC 0.05262196 #> AG AG 0.05440317 #> AT AT 0.08125731 #> CA CA 0.06794174 #> CC CC 0.04634818 #> CG CG 0.04133863 #> CT CT 0.05440317 #> GA GA 0.05570537 #> GC GC 0.05535490 #> GG GG 0.04634818 #> GT GT 0.05262196 #> TA TA 0.06463509 #> TC TC 0.05570537 #> TG TG 0.06794174 #> TT TT 0.10168661 |
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 :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
double_strand_nucleotide_frequency <- function(dnastring) { top_strand_n <- alphabetFrequency(dnastring, baseOnly = TRUE)[, 1:4] |> colSums() bottom_strand_n <- top_strand_n[c("T", "G", "C", "A")] proportions(top_strand_n + bottom_strand_n) } double_strand_dinucleotide_frequency <- function(dnastring) { nucleotide_frequency <- double_strand_nucleotide_frequency(dnastring) top_strand_n <- dinucleotideFrequency(dnastring) |> colSums() bottom_strand_n <- top_strand_n[c("TT", "GT", "CT", "AT", "TG", "GG", "CG", "AG", "TC", "GC", "CC", "AC", "TA", "GA", "CA", "AA")] observed <- tibble( dinucleotide = names(top_strand_n), nucleotide_1 = substr(dinucleotide, 1, 1), nucleotide_2 = substr(dinucleotide, 2, 2), type = "observed", n = top_strand_n + bottom_strand_n, freq = proportions(n), ) theoric <- tibble( dinucleotide = names(top_strand_n), nucleotide_1 = substr(dinucleotide, 1, 1), nucleotide_2 = substr(dinucleotide, 2, 2), type = "theoric", freq = nucleotide_frequency[nucleotide_1] * nucleotide_frequency[nucleotide_2], n = freq * sum(observed$n) ) bind_rows(observed, theoric) } |
Ce qui donne :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
double_strand_dinucleotide_frequency(droso_genome) # A tibble: 32 × 6 #> dinucleotide nucleotide_1 nucleotide_2 type n freq #> <chr> <chr> <chr> <chr> <dbl> <dbl> #> 1 AA A A observed 28995037 0.102 #> 2 AC A C observed 15004688 0.0526 #> 3 AG A G observed 15512583 0.0544 #> 4 AT A T observed 23169804 0.0813 #> 5 CA C A observed 19372988 0.0679 #> 6 CC C C observed 13215774 0.0463 #> 7 CG C G observed 11787346 0.0413 #> 8 CT C T observed 15512583 0.0544 #> 9 GA G A observed 15883894 0.0557 #>10 GC G C observed 15783960 0.0554 # … with 22 more rows |
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à :
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 !
