Accessibility Tools

Web Accessibility plugin by DJ-Extensions.com

- Le blog participatif de bioinformatique francophone depuis 2012 -

Les éléments répétés du génome humain : aperçu rapide avec R et le tidyverse

-

par

Guillaume Devailly
dans

Dans un pré­cé­dent article, nous avions regar­dé le fichier d'annotation des gènes du génome humain d’après Gen­code. J'avais uti­li­sé pour cela la puis­sante com­bi­nai­son dplyr + ggplot2 (packages cen­traux du tidy­verse), par­ti­cu­liè­re­ment adap­tée à tout ce qui est mani­pu­la­tion et visua­li­sa­tion de don­nées tabu­laires.

Mais notre génome n'est pas consti­tué que de gènes, loin s'en faut ! Les élé­ments répé­tés sont en fait bien plus majo­ri­taires. Je ne vais pas me ris­quer à don­ner ici une défi­ni­tion pré­cise de ce qu'est un élé­ment répé­té, je me conten­te­rai de rap­pe­ler que si les élé­ments trans­po­sables sont des élé­ments répé­tés, tout les élé­ments répé­tés ne sont pas trans­po­sables ! Comme sou­vent en bio-infor­ma­tique, je vais me conten­ter de la défi­ni­tion prag­ma­tique d'élément répé­té : un élé­ment répé­té est un élé­ment décrit dans ma table d'annotation des élé­ments répé­tés. :-p

Les sources d'annotation des élé­ments répé­tés du génome humain sont bien plus rares que pour ce qui concerne les gènes. Je vous pro­pose d'utiliser le temps de cet article une table dis­po­nible sur le UCSC table brow­ser. Alors oui, l'interface a mal vieilli, mais le UCSC table brow­ser reste une for­mi­dable col­lec­tion de fichiers d'annotation du génome. Pour obte­nir la table en ques­tion, il suf­fit nor­ma­le­ment de chan­ger le champ group sur Repeats et de lais­ser le reste par défaut.

Com­ment obte­nir une table d'annotation des élé­ments répé­tés du génome humain. Vous pou­vez cli­quer sur le bou­ton des­cribe table sche­ma pour une des­crip­tion des colonnes de la table.

J'ai per­son­nel­le­ment télé­char­gé cette table le 4 avril 2017. Peut-être la vôtre sera-t-elle plus récente, et donc légè­re­ment dif­fé­rente ? En effet, les anno­ta­tions du génome humain, gènes comme élé­ments répé­tés, ne sont pas encore par­faites et sont tou­jours acti­ve­ment amé­lio­rées. Cette table a été géné­rée à l'aide de l'outil Repeat­Mas­ker, outil qui per­met de mas­quer (en rem­pla­çant par des N) les nucléo­tides d'un fichier fas­ta qui sont inclus dans des élé­ments répé­tés. Je trouve assez iro­nique qu'une des meilleures sources d'annotation des élé­ments répé­tés soit issue d'un logi­ciel visant à s'en débar­ras­ser. ^^ Ce logi­ciel de plus de 20 ans sert notam­ment à faci­li­ter l'annotation des gènes des génomes en mas­quant les séquences répé­tées.

Si vous sou­hai­tez repro­duire les ana­lyses ci-des­sous, je vous laisse donc télé­char­ger la table, la mettre dans un réper­toire de tra­vail, et lan­cer R. Si vous n'en avez rien à faire de R, vous pou­vez tout à fait sau­ter les blocs de code et autres expli­ca­tions pour vous conten­ter de regar­der les jolies images. 🙂 Je détaille cepen­dant ma démarche, en espé­rant qu'au moins l'une ou l'un d'entre vous puisse en reti­rer une astuce utile, au prix d'un alour­dis­se­ment assez consé­quent de ce billet.

Import et toilettage des données

Après avoir lan­cé R et défi­ni un réper­toire de tra­vail appro­prié (via la com­mande

1
setwd()

 ), je com­mence par char­ger quelques packages que j'aime bien :

1
2
3
4
5
6
7
8
9
libra­ry(purrr)   # pro­gram­ma­tion fonc­tion­nelle
libra­ry(dplyr)   # mani­pu­la­tion de tableaux de données
libra­ry(rea­dr)   # import de fichiers txt
libra­ry(ggplot2) # des jolis plots
libra­ry(cow­plot) # thème ggplot2 + outils pour figures mul­ti-pan­neaux
libra­ry(for­cats) # mani­pu­la­tions de fac­tor
libra­ry(svglite) # export en svg
libra­ry(viri­dis) # une palette de cou­leur sym­pa
libra­ry(xtable)  # export des tableaux en HTML

J'importe la table dans R à l'aide d'une fonc­tion du package rea­dr, fonc­tion qui est plus rapide et qui a des valeurs par défaut de ses para­mètres plus adap­tées que la fonc­tion

1
read.table()

  de R base :

1
repeats_​dirty &lt ;- read_​tsv("repeats_rmsk_hg38.txt.gz")

Je vais ensuite ne gar­der que les colonnes qui m'intéressent, que je renomme. J'en pro­fite aus­si pour ne gar­der que les lignes concer­nant les chro­mo­somes stan­dards, en fil­trant les haplo­types alter­na­tifs qui ne feraient qu'alourdir cer­taines figures par la suite.

1
2
3
4
5
standard_​chromosomes &lt ;- c(paste0("chr", 1 :22), "chrX", "chrY")
 
repeats &lt ;- select(repeats_​dirty, geno­Name, genoS­tart, genoEnd, strand, rep­Name, rep­Fa­mi­ly, rep­Class) %&gt ;%
    dplyr::rename(chr = geno­Name, start = genoS­tart, end = genoEnd, name = rep­Name, fami­ly = rep­Fa­mi­ly, class = rep­Class) %&gt ;%
    fil­ter(chr %in% standard_​chromosomes)

Ce qui me donne cette table-ci :

1
2
3
head(repeats) %&gt ;%
    xtable() %&gt ;%
    print(type = "html", include.row­names = FALSE)
chrstartendstrandnamefami­lyclass
chr16710875367109046+L1P5L1LINE
chr183883158388618-AluYAluSINE
chr12516580325166380+L1MB5L1LINE
chr13355418533554483-AluScAluSINE
chr14194289441943205-AluYAluSINE
chr15033133650332274+HAL1L1LINE

Nous avons donc une clas­si­fi­ca­tions des élé­ments répé­tés en trois niveaux hié­rar­chiques, dans l'ordre : class > fami­ly > name.

Avant de regar­der plus en détail cette clas­si­fi­ca­tion, j'en pro­fite pour fil­trer les quelques lignes conte­nant un "?", qui cor­res­pondent à des clas­si­fi­ca­tions incer­taines. Je pour­rais les gar­der, mais il y en a rela­ti­ve­ment peu, et elles com­plexi­fie­raient l'analyse et alour­di­raient les figures.

1
2
3
4
5
6
7
repeats &lt ;- fil­ter(
    repeats,
    !(gre­pl("?", repeats$class, fixed = TRUE) |
        gre­pl("?", repeats$fami­ly, fixed = TRUE) |
        gre­pl("?", repeats$name, fixed = TRUE)
    )
)

Les classes d'éléments répétés

La pre­mière figure que nous allons géné­rer s’intéresse au 1er niveau hié­rar­chique de la clas­si­fi­ca­tion : les classes d'éléments répé­tés. Com­bien y en a‑t-il ? (divul­gâ­chis : 16) Quel est l’effectif de cha­cune des classes ? Quelle frac­tion du génome chaque classe couvre-t-elle ? Quelle est la dis­tri­bu­tion des lon­gueurs des élé­ments au sein de chaque classe ?

Je vais trier les classes par effec­tif décrois­sant pour rendre la figure plus jolie. Pour cela, comme j'utiliserai ggplot2, il me faut modi­fier l'ordre des levels de la colonne class après l'avoir trans­for­mée en fac­tor. J'utilise quelques fonc­tions du package for­cats.

1
repeats$class &lt ;- fac­tor(repeats$class) %&gt ;% fct_​infreq %&gt ;% fct_​rev

Pour le pre­mier pan­neau, un dia­gramme en barres, j'utilise des astuces vues dans le pré­cé­dent billet. Si vous décou­vrez ggplot2, pour­quoi ne pas jeter un coup d’œil sur cet article star du blog ? Comme tou­jours, ce qui prend le plus de lignes, ce n'est pas la figure en elle-même, mais tous les petits ajus­tages néces­saires pour la rendre plus jolie.

1
2
3
4
5
6
7
plot_​Class &lt ;- ggplot(repeats, aes(x = class)) +
    geom_​bar(stat = "count", fill = "indianred1") +
    geom_​text(aes(label = ..count..), y = 10000, hjust = 0, stat = "count") +
    labs(x = "Classe", y = "Nombre d'éléments") +
    scale_​y_​continuous(sec.axis = dup_​axis()) +
    coord_​flip() +
    background_​grid(major = "xy", minor = "none")

Pour le deuxième pan­neau, j'ai envie de voir la lon­gueur totale cou­verte par chaque classe d'éléments répé­tés. Pour aider à la lec­ture, ce ne serait pas mal d'indiquer aus­si la frac­tion du génome cou­verte par chaque classe. Voi­là une excel­lente occa­sion d'utiliser une fea­ture récem­ment ajou­tée à ggplot2 : le second axe ! Avant toute chose, et comme la der­nière fois, je récu­père auprès de UCSC la lon­gueur totale de chaque chro­mo­some :

1
2
3
4
5
chr_​length &lt ;- read_​tsv("http://​hgdown​load​-test​.cse​.ucsc​.edu/​g​o​l​d​e​n​P​a​t​h​/​h​g​3​8​/​b​i​g​Z​i​p​s​/​h​g​3​8​.​c​h​r​o​m​.​s​i​zes", col_​names = FALSE) %&gt ;%
    dplyr::rename(seq­names = X1, length = X2) %&gt ;%
    fil­ter(seq­names %in% standard_​chromosomes) %&gt ;%
    mutate(seq­names = fac­tor(seq­names, levels = standard_​chromosomes))
genome_​length &lt ;- sum(as.nume­ric(chr_​length$length)) # 3 088 269 832

Pour faci­li­ter les pro­chains cal­culs, je rajoute une colonne conte­nant la lar­geur de chaque élé­ment :

1
repeats &lt ;- mutate(repeats, width = end - start)

Et c'est par­ti pour un peu de magie dplyr ! Je groupe mon tableau par classe d'éléments répé­tés (avec

1
group_​by()

), je cal­cule ensuite la lon­gueur totale cou­verte par chaque classe (avec

1
sum­ma­rise()

 ), et je lance le tout dans ggplot2 ! Je spé­ci­fie bien que je sou­haite un axe secon­daire, qui est une trans­for­ma­tion linéaire de l'axe prin­ci­pal (

1
sec.axis = sec_​axis(~100 * . /​ genome_​length, name = "% du génome")

 ).

1
2
3
4
5
6
7
8
9
10
11
12
plot_​genomeProp &lt ;- repeats %&gt ;%
    group_​by(class) %&gt ;%
    sum­ma­rise(total_​length = sum(width)) %&gt ;%
    arrange(total_​length) %&gt ;%
    ggplot(aes(x = class, y = total_​length)) +
    geom_​bar(stat = "iden­ti­ty", fill = "corn­flo­wer­blue") +
    geom_​text(aes(label = paste0(round(100 * total_​length /​ genome_​length, digits = 1), "%")), y = 10000, hjust = 0) +
    labs(y = "Lon­gueur cumu­lée (pb)") +
    scale_​y_​continuous(sec.axis = sec_​axis(~100 * . /​ genome_​length, name = "% du génome")) +
    coord_​flip() +
    theme(axis.text.y = element_​blank(), axis.title.y=element_​blank()) +
    background_​grid(major = "xy", minor = "none")

Le troi­sième et der­nier pan­neau sera un aper­çu de la dis­tri­bu­tion des lar­geurs pour chaque classe d'éléments répé­tés. Des boites à mous­taches géné­rées avec ggplot2 suf­fisent ici :

1
2
3
4
5
6
7
plot_​sizeDistrib &lt ;- ggplot(repeats, aes(x = class, y = width)) +
    geom_​boxplot(fill = "mediu­mor­chid", out­lier.shape = NA) +
    labs(x = "Classe", y = "Taille des élé­ments (pb)") +
    scale_​y_​continuous(sec.axis = dup_​axis()) +
    coord_​flip(ylim = c(0, 2500)) +
    theme(axis.text.y = element_​blank(), axis.title.y=element_​blank()) +
    background_​grid(major = "xy", minor = "none")

Enfin, j'arrange labo­rieu­se­ment les pan­neaux à l'aide du package cow­plot et de quelques nombres magiques qui vont bien pour rendre la figure plai­sante à l’œil :

1
2
3
4
5
6
7
8
9
10
myoff­set &lt ;- 0.008
first­plot­size &lt ;- 0.44
svglite("plots/classeER.svg", width = 10, height = 5)
ggdraw() +
    draw_​plot(plot_​Class      , x = 0.0, y = myoff­set, w = first­plot­size, h = 0.96 - 2*myoff­set) +
    draw_​plot(plot_​genomeProp , x = first­plot­size, y = 0.0, w = (1-first­plot­size)/​2, h = 0.96) +
    draw_​plot(plot_​sizeDistrib, x = first­plot­size + (1-first­plot­size)/​2, y = 0.0, w = (1-first­plot­size)/​2, h = 0.96) +
    draw_​plot_​label(LETTERS[1 :3], x = c(0.14, first­plot­size - 0.01, first­plot­size + (1-first­plot­size)/​2 - 0.01), y = 0.92) +
    draw_​label("Classes d'éléments répé­tés", size = 15, x = 0.5, y = 0.97)
dev.off()

alt-text Figure 1 : Les classes d'éléments répé­tés du génome humain. A. Nombre d'éléments répé­tés pour chaque classe. B. Frac­tion du génome cou­vert par chaque classe. C. Dis­tri­bu­tion des tailles d'éléments répé­tés pour chaque classe.

Les plus obser­va­teurs d'entre vous auront peut être réa­li­sé, avec stu­peur, qu'en effet Repeat­Mas­ker caté­go­rise les gènes d'ARN ribo­so­maux (rRNA) et de trans­ferts (tRNA) comme étant des élé­ments répé­tés ! Ce qui est tech­ni­que­ment exact, mais m'a un peu sur­pris au début (ça va mieux main­te­nant, mer­ci). Je me suis amu­sé à com­pa­rer le nombre de copies de gènes d'ARN ribo­so­maux recen­sé par GENCODE, vu la der­nière fois (544) avec ceux repé­rés par Repeat­Mas­ker (1 751). Peut-être la dif­fé­rence est-elle due aux copies non fonc­tion­nelles, incluses dans la liste Repeat­Mas­ker mais pas dans celle de GENCODE ? Une telle dif­fé­rence se retrouve pour d'autres caté­go­ries de gènes ARN. Par exemple GENCODE recense 1 900 snR­NA et Repeat­Mas­ker 4 285.

Si les SINE (short inter­sper­sed nuclear ele­ments) sont plus nom­breux que les LINE (Long inter­sper­sed nuclear ele­ments), ils sont en géné­ral plus courts, et donc consti­tuent une frac­tion moindre de notre génome. La troi­sième classe la plus abon­dante, à la fois en effec­tif et en frac­tion du génome, est celle des élé­ments à LTR (long ter­mi­nal repeat). Il s'agit donc d'éléments issus de rétro­vi­rus endo­gènes.

Notez que la figure 1C ne montre pas les points ous­ti­ders. En effet, les plus longs élé­ments répé­tés le sont tel­le­ment que les mon­trer écra­se­raient le reste de la figure. Voyez plu­tôt :

1
2
3
4
arrange(repeats, desc(width)) %&gt ;%
    top_​n(5) %&gt ;%
    xtable() %&gt ;%
    print(type = "html", include.row­names = FALSE)
chrstartendstrandnameclassfami­lywidth
chr1123500000124000000+ALR/​AlphaSatel­litecen­tr500000
chr1123000000123500000+ALR/​AlphaSatel­litecen­tr500000
chr54800000048500000+ALR/​AlphaSatel­litecen­tr500000
chr75900000059500000+ALR/​AlphaSatel­litecen­tr500000
chr84450000045000000+ALR/​AlphaSatel­litecen­tr500000
chr123500000035500000+ALR/​AlphaSatel­litecen­tr500000

D’après la table d'annotation, les élé­ments répé­tés les plus longs sont donc les cen­tro­mères, fai­sant tous exac­te­ment 500 000 paires de base. Quelle coïn­ci­dence ! En fait, à l'heure d’écriture de cet article, les cen­tro­mères du génome humain ne sont tou­jours pas assem­blés… Parce que figu­rez-vous qu'assembler 23 ou 24 régions d'environ 500 kb très hau­te­ment répé­tées, ce n'est pas de la tarte ! En atten­dant, les cen­tro­mères sont donc anno­tés avec une lon­gueur esti­mée arbi­traire. Mais avec le rapide déve­lop­pe­ment des tech­no­lo­gies de séquen­çage de frag­ments longs, il est pos­sible que les cen­tro­mères humains soient assem­blés pro­chai­ne­ment. Les plus longs reads séquen­cés par la tech­no­lo­gie Nano­pore se rap­prochent de la méga-base !

Notre génome est en tout cas consti­tué par envi­ron :

1
sum(repeats$width) /​ genome_​length

49,4% d'éléments répé­tés ! Sont-ils homo­gè­ne­ment répar­tis entre les chro­mo­somes ? C'est ce que je vous pro­pose de décou­vrir ensuite.

Distributions des éléments répétés entre chromosomes

Tout d'abord, sou­hai­tant mettre en évi­dence les trois plus grandes caté­go­ries d'éléments répé­tés (LINE, SINE et LTR), je crée une nou­velle colonne  via un

1
mutate()

  et un

1
if_​else()

 . Je regroupe ensuite le tableau par chro­mo­some (

1
group_​by()

 ) et par classe et somme les lar­geurs d'éléments répé­tés (

1
mutate(sum(width))

 ). Je joins le tableau à celui conte­nant la lon­gueur des chro­mo­somes (

1
left_​join()

 ) pour pou­voir cal­cu­ler la frac­tion de chaque chro­mo­some conte­nant des élé­ments répé­tés (le second

1
mutate()

 ). J'en pro­fite pour réorien­ter les levels de fac­tors pour ordon­ner les dif­fé­rentes colonnes dans la figure. Et enfin j'envoie les don­nées dans ggplot2, en ajus­tant tout un tas de micro-détails pour avoir une figure exac­te­ment comme j'aime :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
svglite("plots/ERparChrom.svg", width = 5, height = 3)
repeats %&gt ;%
    mutate(simple_​class = if_​else(
        !(class %in% c("SINE", "LINE", "LTR")),
        "other",
        as.cha­rac­ter(class)
    )) %&gt ;%
    group_​by(chr, simple_​class) %&gt ;%
    sum­ma­rize(cum_​size = sum(width)) %&gt ;%
    left_​join(chr_​length, by = c("chr" = "seq­names")) %&gt ;%
    ungroup %&gt ;%
    mutate(
        frac_​repeat = cum_​size/​length,
        simple_​class = fac­tor(simple_​class, levels = c("other", "LTR","SINE", "LINE")),
        chr = fac­tor(chr, levels = standard_​chromosomes)
    ) %&gt ;%
    ggplot(aes(x = chr, y = frac_​repeat, fill = simple_​class)) +
    geom_​bar(stat = "iden­ti­ty") +
    theme(axis.text.x = element_​text(angle = 90, hjust = 0, vjust = 0.5)) +
    scale_​y_​continuous(labels = scales::percent) +
    scale_​fill_​manual(values = c("gray", viri­dis(3, begin = 0.25, end = 0.9))) +
    labs(x = NULL, y = "% du chro­mo­some", title = "Conte­nu répé­té de chaque chro­mo­some", fill = "Classe") +
    background_​grid(major = "xy", minor = "none")
dev.off()

alt-text Figure 2 : Conte­nu en élé­ments répé­tés de chaque chro­mo­some.

En pre­mière approxi­ma­tion, il semble que chaque auto­some ait un conte­nu en élé­ments répé­tés à peu près équi­valent, oscil­lant entre envi­ron 42% pour le chro­mo­some 22, et 59% pour le chro­mo­some 19. Il est amu­sant de com­pa­rer cette figure avec celle du conte­nu en gènes de chaque chro­mo­some géné­rée la der­nière fois. Ain­si le chro­mo­some 19 est à la fois l'autosome le plus riche en gènes pro­téiques et le plus riche en élé­ments répé­tés ! L'énigmatique chro­mo­some 13 est rela­ti­ve­ment pauvre en élé­ments répé­tés, et en même temps pauvre en gènes. Les chro­mo­somes sexuels font ici leur malin, avec le chro­mo­some X ayant le plus fort taux en élé­ments répé­tés (62%), et le chro­mo­some Y le plus faible (28%). Éton­nam­ment (en tout cas pour moi), notre chro­mo­some Y est donc ni riche en pseu­do­gènes, ni riche en élé­ments répé­tés, il est juste… petit.

Les familles d'éléments répétés

Après avoir détaillé les classes d'éléments répé­tés, jetons un œil aux niveaux de clas­si­fi­ca­tions sui­vants, familles et sous-familles :

1
2
3
4
5
group_​by(repeats, class) %&gt ;%
    sum­ma­rise(n_​family = length(unique(fami­ly)), n_​subfamily = length(unique(name)), n_​element = n()) %&gt ;%
    arrange(desc(n_​family), desc(n_​subfamily)) %&gt ;%
    xtable() %&gt ;%
    print(type = "html", include.row­names = FALSE)
classn_​familyn_​subfamilyn_​element
DNA16226479941
LINE71711516226
LTR6567709475
SINE6601779233
Satel­lite4227018
Simple_​repeat114162678663
Unk­nown1715531
tRNA1621777
snR­NA1124285
Retro­po­son165397
scR­NA151334
Low_​complexity1498618
rRNA131751
RC131754
RNA11666
srpR­NA111595

Nous allons essayer de repré­sen­ter gra­phi­que­ment cette diver­si­té, en affi­chant des dia­grammes en barres d'effectif de chaque famille de répé­tés. Je vais colo­rier les barres par le nombre de sous-familles pour chaque famille. Les effec­tifs variant énor­mé­ment, je suis contraint d'utiliser une échelle loga­rith­mique. J'ai alors été sur­pris de décou­vrir que pour l'instant,

1
coord_​flip()

  et

1
annotation_​logticks()

  sont mutuel­le­ment exclu­sif !

Je com­mence par pré­pa­rer les don­nées :

1
2
3
effectif_​table &lt ;- group_​by(repeats, class, fami­ly) %&gt ;%
    sum­ma­rise(diff_​name = length(unique(name)), size = n()) %&gt ;%
    ungroup()

Géné­rons ensuite un pan­neau de figure par classe d'éléments répé­tés pos­sé­dant de mul­tiples sous-familles (classes Satel­lite, LTR, LINE, SINE et DNA). J'utilise pour cela la fonc­tion

1
map()

  du package purrr, une variante de

1
lap­ply()

 , en défi­nis­sant une fonc­tion ano­nyme via les nota­tions quelque peu éso­té­riques

1
~

  et

1
.x

 .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
mul­ti­Fa­mi­ly &lt ;- c("Satel­lite", "LTR", "LINE", "SINE", "DNA")
 
myFa­mi­ly­Plots &lt ;- map(
    mul­ti­Fa­mi­ly,
    ~fil­ter(effectif_​table, class == .x) %&gt ;%
        arrange(size) %&gt ;%
        mutate(fami­ly = fac­tor(fami­ly, levels = fami­ly)) %&gt ;%
        mutate(text_​dark = if_​else(diff_​name &gt ; 50, TRUE, FALSE)) %&gt ;%
        ggplot(aes(x = fami­ly, y = size, fill = diff_​name)) +
        geom_​bar(stat = "iden­ti­ty") +
        geom_​text(aes(label = size, color = text_​dark), y = 1, hjust = 0) +
        scale_​fill_​viridis(limits=c(0, 100),  oob = scales::squish) +
        scale_​color_​manual(values = c("grey80", "black"), guide = FALSE) +
        scale_​y_​log10() +
        coord_​flip(ylim = c(10, 1e7)) +
        background_​grid(major = "xy", minor = "none") +
        labs(x = NULL, y = NULL, title = .x) +
        theme(legend.posi­tion = "none")
)
names(myFa­mi­ly­Plots) &lt ;- mul­ti­Fa­mi­ly
myFa­mi­ly­Plots$DNA &lt ;- myFa­mi­ly­Plots$DNA + labs(y = "Nombre d'éléments")

Géné­rons ensuite le même type de pan­neau pour toutes les classes d'éléments ayant une seule famille :

1
2
3
4
5
6
7
8
9
10
11
12
plot_​other &lt ;- fil­ter(effectif_​table, !(class %in% c("Satel­lite", "LTR", "LINE", "SINE", "DNA"))) %&gt ;%
    arrange(size) %&gt ;%
    mutate(fami­ly = fac­tor(fami­ly, levels = fami­ly)) %&gt ;%
    ggplot(aes(x = fami­ly, y = size, fill = diff_​name)) +
    geom_​bar(stat = "iden­ti­ty") +
    geom_​text(aes(label = size), y = 1, hjust = 0, stat = "count") +
    scale_​fill_​viridis(limits=c(0, 100),  oob = scales::squish) +
    scale_​y_​log10() +
    coord_​flip(ylim = c(10, 1e7)) +
    background_​grid(major = "xy", minor = "none") +
    labs(x = NULL, y = NULL, title = "Autres", fill = "Nombre de\nsous-familles") +
    theme(legend.posi­tion = "bot­tom")

Je récu­père la légende pour l'afficher à part, à l'aide d'une fonc­tion de cow­plot :

1
2
myLe­gend &lt ;- get_​legend(plot_​other)
plot_​other &lt ;- plot_​other + theme(legend.posi­tion = "none")

J'utilise un peu de magie noire pour homo­gé­néi­ser les marges de mes dif­fé­rents pan­neaux et gérer l'alignement ver­ti­cal. Ne me deman­dez pas d'expliquer, j'ai juste copié-col­lé un bout de code depuis inter­net.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
myFa­mi­ly­Plots &lt ;- map(myFa­mi­ly­Plots, ggplot­Grob)
plot_​other &lt ;- ggplot­Grob(plot_​other)
 
myFa­mi­ly­Plots.widths &lt ;- map(myFa­mi­ly­Plots, ~.x$widths[1 :3])
plot_​other.widths &lt ;- plot_​other$widths[1 :3]
 
max.widths &lt ;- grid::unit.pmax(
    plot_​other.widths,
    do.call(grid::unit.pmax, myFa­mi­ly­Plots.widths)
)
 
plot_​other$widths[1 :3] &lt ;- max.widths
myFa­mi­ly­Plots &lt ;- map(myFa­mi­ly­Plots, func­tion(x) {
    x$widths[1 :3] &lt ;- max.widths
    return(x)
})

Enfin, j'arrange les dif­fé­rents pan­neaux, et j'exporte la figure dans un .SVG :

1
2
3
4
5
6
7
8
9
10
11
12
svglite("plots/familleER.svg", width = 10, height = 7)
ggdraw() +
    draw_​text("Les familles\nd'éléments répé­tés", x = 1/​6, y = 0.95, size = 18) +
    draw_​plot(myLe­gend               , x = 0.0, y = 0.75, w = 1/​3, h = 0.2 ) +
    draw_​plot(myFa­mi­ly­Plots$DNA      , x = 0.0, y = 0.0 , w = 1/​3, h = 0.8 ) +
    draw_​plot(myFa­mi­ly­Plots$LINE     , x = 1/​3, y = 0.62, w = 1/​3, h = 0.38) +
    draw_​plot(myFa­mi­ly­Plots$SINE     , x = 1/​3, y = 0.31, w = 1/​3, h = 0.31) +
    draw_​plot(myFa­mi­ly­Plots$LTR      , x = 1/​3, y = 0.0 , w = 1/​3, h = 0.31) +
    draw_​plot(myFa­mi­ly­Plots$LTR      , x = 1/​3, y = 0.0 , w = 1/​3, h = 0.31) +
    draw_​plot(myFa­mi­ly­Plots$Satel­lite, x = 2/​3, y = 0.7 , w = 1/​3, h = 0.3 ) +
    draw_​plot(plot_​other             , x = 2/​3, y = 0.0 , w = 1/​3, h = 0.7 )
dev.off()

alt-text Figure 3 : Les familles d'éléments répé­tés.

Recon­nais­sez-vous des noms fami­liers ? Par exemple, nous avons 2 118 inser­tions d'éléments Pig­gy­Bac dans notre génome. Ce trans­po­son est à l'origine d'une méthode de clo­nage de gènes dans des plas­mides assez popu­laire.

Ce que je remarque sur­tout, c'est que des ARN de trans­ferts (tRNA) se baladent dans la caté­go­rie des SINE. MAIS POURQUOI ! POURQUOI ON NE PEUT PAS AVOIR DES CLASSIFICATIONS COHÉRENTES EN BIOINFORMATIQUE !

Hum hum, par­don.

En fait tout va bien : la clas­si­fi­ca­tion est stric­te­ment non che­vau­chante au niveau des sous-familles : les tRNA de classe SINE ne contiennent pas les même sous-familles de tRNA que les tRNA de classe tRNA. Oui, je sais, ce n'est pas très clair.  Mais il se trouve qu'un cer­tain nombre de SINE dérivent de séquences d'ARN de trans­ferts. Je pense donc que cette clas­si­fi­ca­tion est tout à fait jus­ti­fiée.

Je pour­rais me perdre ensuite dans les détails des dif­fé­rentes sous-familles d'élément répé­tés, mais je pré­fère lais­ser les plus curieux d'entre vous se perdre dans ce fas­ci­nant tableau, et nous racon­ter leurs trou­vailles en com­men­taires. Et c'est donc sur cette abrupte conclu­sion que je conclus.

Un grand mer­ci aux super relec­teurs et relec­trice : Clé­mence, eorn, Mathu­rin et Max, sans qui cet article serait beau­coup moins bien.

Vous avez aimé ? Dites-le nous !

Moyenne : 5 /​ 5. Nb de votes : 1

Pas encore de vote pour cet article.

Partagez cet article

Guillaume Devailly

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

2 réponses à “Les éléments répétés du génome humain : aperçu rapide avec R et le tidyverse”

  1. Kumquat
    Kumquatum

    "Je trouve assez iro­nique qu'une des meilleures sources d'annotation des élé­ments répé­tés soit issue d'un logi­ciel visant à s'en débar­ras­ser" => Mieux connaître son enne­mi pour l'éliminer :p

    "divul­gâ­chis" <3

    Sinon très bon article, clair et péda­gogue.
    Ce serait amu­sant d'appliquer cette métho­do­lo­gie sur d'autres orga­nisme pour com­pa­rer 🙂

    Répondre
  2. Guillaume Devailly
    Guillaume Devailly

    Une des vignettes de cow­plot a dis­pa­ru de CRAN, et a été mise à jour, résul­tant un un lien mort.

    "Ne me deman­dez pas d'expliquer, j'ai juste copié-col­lé un bout de code depuis inter­net."

    La vignette en ques­tion est tou­jours visible entre autre via la way­back machine de l'internet archive :

    https://web.archive.org/web/20170711035139/https://cran.r‑project.org/web/packages/cowplot/vignettes/plot_grid.html

    Mais il y a main­te­nant des méthodes plus recom­man­dables d'alignement des plots :

    https://​wil​ke​lab​.org/​c​o​w​p​l​o​t​/​a​r​t​i​c​l​e​s​/​a​l​i​g​n​i​n​g​_​p​l​o​t​s​.​h​tml

    Répondre

Laisser un commentaire

Pour insérer du code dans vos commentaires, utilisez les balises <code> et <\code>.