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

—

par

Guillaume Devailly
dans

L'analyse de séquences est au cœur de nom­breux domaines de la bio-infor­ma­tique. Le billet du jour s'intéressera aux séquences ADN, en se pro­po­sant de comp­ter la fré­quence en dinu­cléo­tides dans quelques génomes d'organismes modèles (avec une petite arrière-pen­sée der­rière la tête).

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

L'ADN double brins est clas­si­que­ment struc­tu­ré sous forme de double hélice, avec deux brins de direc­tion oppo­sée. À chaque posi­tion il y a deux acides nucléiques for­mant une paire, les appa­rie­ments clas­siques étant A en face de T et C en face de G. Ce qu'on appelle dinu­cléo­tide est une séquence de deux nucléo­tides suc­ces­sifs, à ne pas confondre avec une paire de nucléo­tide. 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 dinu­cléo­tide. En rouge : une paire de base.

On peut aus­si dire qu'un dinu­cléo­tide est un k‑mer de lon­gueur (k) égale à 2.

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

Pre­nons par exemple cette séquence :

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

Quels dinu­cléo­tides dénom­brez vous ? Je vous laisse deux phrases pour réflé­chir. C'est bon pour vous ?

On trouve deux fois le nucléo­tide TG (facile), et une fois le nucléo­tide GT (au milieu, peut-être l'avez-vous oublié ?). Donc pour une séquence de lon­gueur n, il y aura n‑1 dinu­cléo­tides. Ce qui donne une fré­quence de 2/​3 pour le dinu­cléo­tide TG et une fré­quence de 1/​3 pour le dinu­cléo­tide 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, sur­prise ! il y a d'autres dinu­cléo­tides cachés en des­sous !

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

Sous la séquence du haut, on trou­ve­ra donc deux dinu­cléo­tides CA et un dinu­cléo­tide AC (les dinu­cléo­tides sont lus et écrits dans le sens 5' vers 3' par conven­tion). Une séquence double brin de lon­gueur n aura donc 2n‑2 dinu­cléo­tides.

Cela nous don­ne­rait donc une fré­quence de 2/​6 pour les dinu­cléo­tides TG et CA, et un fré­quence de 1/​6 pour les dinu­cléo­tides GT et CA.

Pour une petite séquence, comp­ter à la main c'est fai­sable, mais dès que la séquence s’agrandit, autant le faire en uti­li­sant la *bio-infor­ma­tique* !

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

Je vous pro­pose d'utiliser des fonc­tions exis­tantes du package R /​ Bio­con­duc­tor Bios­trings. La séquence nucléique doit être pré­sen­tée sous la forme d'un objet DNAS­tring qui se crée via la fonc­tion DNAString(). On peut ensuite appli­quer la fonc­tion dinucleotideFrequency() avec ou sans le para­mètre as.prob = TRUE, en fonc­tion de si l'on sou­haite des comp­tages ou des fré­quences.

1
2
3
4
5
6
7
dinu­cleo­ti­de­Fre­quen­cy(DNAS­tring("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(dinu­cleo­ti­de­Fre­quen­cy(DNAS­tring("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 fonc­tion dinucleotideFrequency() ne cal­cule que la fré­quence du brin ren­tré, et non la fré­quence du brin oppo­sé. Pour la cal­cu­ler, deux méthodes.

1) La méthode peu effi­cace mais simple : cal­cu­ler le conte­nu en dinu­cléo­tide du com­plé­ment inverse de la séquence et faire la moyenne des deux fré­quences :

1
2
3
4
5
my_​seq <- DNAS­tring("GTGT")
(dinu­cleo­ti­de­Fre­quen­cy(my_​seq, as.prob = TRUE) +
    dinu­cleo­ti­de­Fre­quen­cy(rever­se­Com­ple­ment(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 effi­cace, mais deman­dant de coder un peu plus : réas­si­gner les fré­quences du brin + aux dinu­cléo­tides com­plé­men­taires inverses, ce qui évite de tout recal­cu­ler.

Voi­ci une petite fonc­tion pour cal­cu­ler la fré­quence en dinu­cléo­tide d'une séquence ADN double brins. J'en pro­fite pour faire un for­mat de sor­tie sous forme de tableau tidy pour nous sim­pli­fier 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
dsdin­freq <- func­tion(dnas­tring, ...) {
    top_​strand_​freq <- dinu­cleo­ti­de­Fre­quen­cy(dnas­tring, ...)
    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(
        dinu­cleo­tide = names(freq),
        freq = freq
    )
}
 
dsdin­freq(DNAS­tring("GTGT"), as.prob = TRUE)
#>    dinu­cleo­tide      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

Main­te­nant que nous savons cal­cu­ler à la fois pra­ti­que­ment et théo­ri­que­ment la fré­quence en dinu­cléo­tides d'une séquence, atta­quons nous à un orga­nisme modèle, au hasard le génome de la drop­so­phile.

Le génome de la plu­part des espèces modèles est déjà dis­po­nible sous forme d'objets BSge­nome dans Bio­con­duc­tor, mais pour être plus géné­rique ici je fais un import direc­te­ment depuis le fichier fas­ta du génome de réfé­rence Ensem­bl.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
droso_​genome_​path <- "http://​ftp​.ensem​bl​ge​nomes​.org/​p​u​b​/​m​e​t​a​z​o​a​/​r​e​l​e​a​s​e​-​5​2​/​f​a​s​t​a​/​d​r​o​s​o​p​h​i​l​a​_​m​e​l​a​n​o​g​a​s​t​e​r​/​d​n​a​/​D​r​o​s​o​p​h​i​l​a​_​m​e​l​a​n​o​g​a​s​t​e​r​.​B​D​G​P​6​.​3​2​.​d​n​a​.​t​o​p​l​e​v​e​l​.​f​a​.gz"
 
droso_​genome <- readD­NAS­tring­Set(droso_​genome_​path)
 
droso_​genome
DNAS­tring­Set 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 pre­mier os. Il n'y a pas une séquence dans le génome de la dro­so­phile, mais 1870 séquences ! Cinq auto­somes, 2 chro­mo­somes sexuels, un chro­mo­some mito­chon­drial, et tout un tas de contigs non assem­blés. Il nous faut donc modi­fier légè­re­ment notre fonc­tion d'origine pour qu'elle fonc­tionne bien avec ce genre d'objets (des DNAS­tring­Set et non de simple DNAS­tring) :

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
dsdin­freq <- func­tion(dnas­tring) {
    top_​strand_​freq <- dinu­cleo­ti­de­Fre­quen­cy(dnas­tring) |> col­Sums()
    bottom_​strand_​freq <- top_​strand_​freq[c("TT", "GT", "CT", "AT", "TG", "GG", "CG", "AG", "TC", "GC", "CC", "AC", "TA", "GA", "CA", "AA")]
    freq <- pro­por­tions(top_​strand_​freq + bottom_​strand_​freq)
    data.frame(
        dinu­cleo­tide = names(freq),
        freq = freq
    )
}
 
dsdin­freq(droso_​genome)
#> dinu­cleo­tide       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
Figure 4 : Fréquences dinucléotidiques du génome de référence de la Drosophile
Figure 4 : Fré­quences dinu­cléo­ti­diques du génome de réfé­rence de la Dro­so­phile

Comme on peut le voir, tout les dinu­cléo­tides n'ont pas la même fré­quence. Mais est-ce un résul­tat atten­du ? Ou pour refor­mu­ler la ques­tion, quel est la fré­quence atten­due pour chaque dinu­cléo­tide ?

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

Comme il y a 16 dinu­cléo­tides pos­sibles, on pour­rait s'attendre à ce que chaque dinu­cléo­tide soit obser­vé avec une fré­quence de 1/​16 (soit 6,25%). Ce serait le cas si les nucléo­tides simples avaient une fré­quence d'un quart. Hors ce n'est pas le cas en pra­tique ! S'il y a bien autant de A que de T et autant de C que de G, le conte­nu en G+C d'un génome peut varier énor­mé­ment d'une espèce à l'autre (C'est d'ailleurs uti­li­sé pour véri­fier les pure­tés d'assemblages, voir par exemple les Blob­Plots) !

On peut donc cal­cu­ler de manière un peu plus fine la fré­quence théo­rique d'un dinu­cléo­tide (TA par exemple, noté p(TA) ) comme la pro­ba­bi­li­té d'avoir un T en pre­mier ( p(T) ) mul­ti­plié par la pro­ba­bi­li­té d'avoir un A en second ( p(A) ), où p(T) et p(A) cor­res­pondent à la fré­quence en nucléo­tide A et T de la séquence étu­dié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 uti­li­ser la fonc­tion alphabetFrequency() de Bios­trings pour obte­nir les fré­quences de chaque nucléo­tide, et ensuite il s'agit de tout faire fonc­tion­ner 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 <- func­tion(dnas­tring) {
    top_​strand_​n <- alpha­bet­Fre­quen­cy(dnas­tring, baseOn­ly = TRUE)[, 1 :4] |> col­Sums()
    bottom_​strand_​n <- top_​strand_​n[c("T", "G", "C", "A")]
    pro­por­tions(top_​strand_​n + bottom_​strand_​n)
}
 
double_​strand_​dinucleotide_​frequency <- func­tion(dnas­tring) {
 
    nucleotide_​frequency <- double_​strand_​nucleotide_​frequency(dnas­tring)
 
    top_​strand_​n <- dinu­cleo­ti­de­Fre­quen­cy(dnas­tring) |> col­Sums()
    bottom_​strand_​n <- top_​strand_​n[c("TT", "GT", "CT", "AT", "TG", "GG", "CG", "AG", "TC", "GC", "CC", "AC", "TA", "GA", "CA", "AA")]
 
    obser­ved <- tibble(
        dinu­cleo­tide = names(top_​strand_​n),
        nucleotide_​1 = sub­str(dinu­cleo­tide, 1, 1),
        nucleotide_​2 = sub­str(dinu­cleo­tide, 2, 2),
        type = "obser­ved",
        n = top_​strand_​n + bottom_​strand_​n,
        freq = pro­por­tions(n),
    )
 
    theo­ric <- tibble(
        dinu­cleo­tide = names(top_​strand_​n),
        nucleotide_​1 = sub­str(dinu­cleo­tide, 1, 1),
        nucleotide_​2 = sub­str(dinu­cleo­tide, 2, 2),
        type = "theo­ric",
        freq = nucleotide_​frequency[nucleotide_​1] * nucleotide_​frequency[nucleotide_​2],
        n = freq * sum(obser­ved$n)
    )
 
    bind_​rows(obser­ved, theo­ric)
}

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
#>   dinu­cleo­tide nucleotide_​1 nucleotide_​2 type            n   freq
#>   <chr>        <chr>        <chr>        <chr>       <dbl>  <dbl>
#> 1 AA           A            A            obser­ved 28995037 0.102
#> 2 AC           A            C            obser­ved 15004688 0.0526
#> 3 AG           A            G            obser­ved 15512583 0.0544
#> 4 AT           A            T            obser­ved 23169804 0.0813
#> 5 CA           C            A            obser­ved 19372988 0.0679
#> 6 CC           C            C            obser­ved 13215774 0.0463
#> 7 CG           C            G            obser­ved 11787346 0.0413
#> 8 CT           C            T            obser­ved 15512583 0.0544
#> 9 GA           G            A            obser­ved 15883894 0.0557
#>10 GC           G            C            obser­ved 15783960 0.0554
# … with 22 more rows
Figure 5 : Fré­quences dinu­cléo­ti­diques obser­vées et atten­dues du génome de réfé­rence de la Dro­so­phile

Il y a donc plus de AA et TT qu'attendu dans le génome de la dro­so­phile ! On pour­rait tes­ter sta­tis­ti­que­ment cette dif­fé­rence, avec par exemple un test du chi2, mais il y a tel­le­ment de nucléo­tides que la p‑value retour­née est scan­da­leu­se­ment petite.

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

Il ne reste plus qu'à faire tour­ner ces fonc­tions sur les génomes de quelques espèces modèles popu­laires (en étant rai­son­nable avec l'API d'Ensembl lors de la récu­pé­ra­tion des génomes de réfé­rences), et voi­là :

Figure 6 : Fré­quences en dinu­cléo­tides des génomes d’espèces modèles en bio­lo­gie.

Les plus obser­va­teurs et obser­va­trices d’entre vous remar­que­ront peut-être la dis­pa­ri­tion des dinu­cléo­tides CG dans les génomes des ver­té­brés. Ce mys­tère pour­ra faire l'objet d'un pro­chain billet si la demande se fait pres­sante.

Remerciements

Mille mer­cis à mes relec­teurs et relec­trices Sébas­tien Gra­dit, Sarah Guin­chard, Vir­gi­nie J, et à nos admins ché­ries !

Guillaume Devailly

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

Pour continuer la lecture :

Commentaires

Laisser un commentaire