Introduction
Non, ne fuyez pas tout de suite, chers lecteurs, tout va s'éclaircir : dplyr, c’est plyr pour les data.frame (les tableaux de données). Attendez, j’y viens, plyr, c’est un package R pour appliquer (apply) des fonctions. Donc, dplyr (prononcez “diplir”), c’est un package R, pour appliquer des fonctions à un tableau de données.
Et ça, mes amis, c'est super cool.
Rendez-vous service : imprimez immédiatement cette feuille d’astuces le concernant (pdf, 0.5Mo). En couleur. En fait imprimez en 10, et donnez en 8 à vos collègues, et mettez en une chez vous, et une sur votre lieu de travail. dplyr est centré autour du concept de données propres (ou tidy) : chaque ligne est une observation, chaque colonne un paramètre observé. Finissons notre petite digression en évoquant broom, un package qui permet de rendre propres (tidy) certains retour de fonctions sales (messy) de R, via la fonction tidy() .
Pour vous initier à dplyr, je vous propose une petite série d'exercices appliqués à notre génome à tous. Ici on ne parlera pas tant du génome dans sa définition moderne (ensemble des séquences d’ADN d’un organisme), mais dans une version plus étymologique, et plus désuète : ensemble des gènes d’un organisme (ici, l’humain).
Si vous n’en avez rien à faire de dplyr, ou a fortiori de R, ne fuyez pas pour autant : vous perdriez une occasion de mieux connaître votre génome !
GENCODE et l’annotation du génome humain
Nous allons utiliser les données d’annotations de GENCODE, un consortium pour l’annotation du génome humain et du génome murin, notamment utilisé par Ensembl.
Pour télécharger le dernier fichier d’annotation disponible, il suffit de se rendre sur leur site internet, et de sélectionner Data > Human > Current release. Lors de l’écriture de cet article, la dernière version est la 25, correspondant aux coordonnées génomiques GRCh38 (aussi appelé hg38). Plusieurs fichiers sont disponibles mais un seul suffira aujourd’hui, le premier : “Comprehensive gene annotation | CHR”, en version GFF3. Le lecteur curieux se demandera quelle est la différence entre GTF et GFF3 ? C’est compliqué. Le GTF est une variante mineure du GFF2. Le GFF3 est plus récent et mieux spécifié, nous allons donc utiliser cette version. Vous devez télécharger ce fichier (45Mo) si vous souhaitez reproduire chez vous l’analyse ci dessous (ce que je vous conseille !). Le fichier GFF3 est un fichier texte, délimité par des tabulations, dont les colonnes suivent un formatage défini.
Lecteur du futur, peut-être aurez-vous la chance de disposer d’un fichier d’annotation plus récent ? Deux choix s’offrent alors à vous : refaire l’analyse avec les jolies annotations du futur, vous pourrez alors comparer avec l’état des connaissances de 2017 et bien vous marrer ; ou bien télécharger la version que nous utiliserons ici (le lien devrait rester valide - 45Mo).
Pré-requis
Une connaissance basique de R, et notamment quelques notions de ggplot2, vous aidera à comprendre le code produit ici (pour s’initier à R, pourquoi ne pas commencer par ici, et pour ggplot2 par là ?). Je vais cependant commenter le code dans l'espoir de le rendre plus intelligible. Si vous souhaitez reproduire l’analyse en parallèle de votre lecture, vous aurez besoin d’une version de R un peu récente (j'utilise la 3.3.1 “Bug in your hair”), et de quelques packages. L'installation de ces derniers ne devrait pas être trop longue :
|
1 2 3 4 5 6 7 8 9 10 11 |
# packages à installer ---------------- # - depuis CRAN: install.packages("dplyr") # manipulation de données tabulaires install.packages("readr") # parseur de fichiers texte install.packages("broom") # nettoie après que R ait tout salit install.packages("cowplot") # figures multi-panaux, inclus ggplot2 install.packages("svglite") # pour exporter les figures en svg # - depuis Bioconductor: source("https://bioconductor.org/biocLite.R") biocLite("rtracklayer") # parseur, notamment de fichiers GFF |
Il vous faut aussi, malheureusement, pas mal de mémoire vive, je dirai au moins 3 ou 4Go pour R (j’ai écrit cet article sur un ordinateur ayant 8Go). Un des défaut majeur de R étant qu’il travaille essentiellement sur la mémoire vive… Oui, je sais, ça fait pas mal de défauts.
Analyse
Parsons le fichier GENCODE
Pour charger le fichier dans R, plutôt que d’écrire notre propre parseur, nous allons utiliser celui du package rtracklayer en utilisant la fonction import() . Cette fonction retourne un objet du type GRanges, que nous allons convertir en data_frame via la fonction as_data_frame() du package dplyr. Les data_frame, (aussi appelés tibble) (ou encore tbl_df) (oui, je sais) de dplyr sont légèrement différents des data.frame de R tout court : la méthode d’affichage (print) par défaut est plus sympa, et les row.names sont interdits !
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
setwd("F:/divers/moi/Bioinfo-fr") # changez votre chemin ! compiler::enableJIT(3) # Invocation de magie noire pour que R tourne plus vite compiler::setCompilerOptions(suppressAll = TRUE) # Cachons la magie noir dans le... noir. library(dplyr) library(readr) library(broom) library(rtracklayer) library(cowplot) library(svglite) # l'import peu prendre 5 minutes en fonction de votre configuration # notez qu'on importe directement la version compressé .gz gencode <- import("gencode.v25.annotation.gff3.gz", format = "GFF") %>% # import() retoune un objet GRanges, nous le convertissons en data_frame as_data_frame |
%>%: le pipe R
Petite digression à propos du pipe : %>% , inclus dans le package dplyr (et issus de magrtittr). Il s’agit de l’équivalent R du pipe | de shell/bash. Il passe le résultat de la fonction de gauche comme premier argument (par défaut) de la fonction de droite :
|
1 2 3 4 5 6 7 8 9 10 11 |
# ainsi, au lieu d’écrire fonction3( fonction2( fonction1(x) ) ) # on pourra écrire fonction1(x) %>% fonction2 %>% fonction3 |
Le raccourcis . permet de définir la position de l’objet dans le cas où il n’est pas le premier argument :
|
1 2 3 4 |
fonction2(x1 = myData1, x2 = fonction1(myData2)) # pourra s’écrire: fonction1(myData2) %>% fonction2(x1 = myData1, x2 = .) |
Mes exemples bidons ne sont pas trop flatteurs, mais ce petit su-sucre syntaxique devient vite indispensable dès qu’on l'essaye. Il permet d’accomplir des tâches complexes sans multiplier les variables intermédiaires ou les fonctions nichées, et comme nous allons le voir, il se marie admirablement bien avec les fonctions du package dplyr, qui prennent toujours comme premier paramètre un tableau de données, et retournent toujours un tableau de données.
Un bien beau jeu de données
Examinons un peu le jeu de données tout juste chargé.
|
1 2 3 4 5 |
> object.size(gencode) 1609890024 bytes > dim(gencode) [1] 2577236 29 > gencode |
| seqnames | start | end | width | strand | source | type | score | phase | ID | gene_id | gene_type | gene_status | gene_name | level | havana_gene | Parent | transcript_id | transcript_type | transcript_status | transcript_name | transcript_support_level | tag | havana_transcript | exon_number | exon_id | ont | protein_id | ccdsid |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| chr1 | 11869 | 14409 | 2541 | + | HAVANA | gene | NA | NA | ENSG00000223972.5 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | character(0) | NA | NA | NA | NA | NA | character(0) | NA | NA | NA | character(0) | NA | NA |
| chr1 | 11869 | 14409 | 2541 | + | HAVANA | transcript | NA | NA | ENST00000456328.2 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENSG00000223972.5 | ENST00000456328.2 | processed_transcript | KNOWN | DDX11L1-002 | 1 | basic | OTTHUMT00000362751.1 | NA | NA | character(0) | NA | NA |
| chr1 | 11869 | 12227 | 359 | + | HAVANA | exon | NA | NA | exon:ENST00000456328.2:1 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000456328.2 | ENST00000456328.2 | processed_transcript | KNOWN | DDX11L1-002 | 1 | basic | OTTHUMT00000362751.1 | 1 | ENSE00002234944.1 | character(0) | NA | NA |
| chr1 | 12613 | 12721 | 109 | + | HAVANA | exon | NA | NA | exon:ENST00000456328.2:2 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000456328.2 | ENST00000456328.2 | processed_transcript | KNOWN | DDX11L1-002 | 1 | basic | OTTHUMT00000362751.1 | 2 | ENSE00003582793.1 | character(0) | NA | NA |
| chr1 | 13221 | 14409 | 1189 | + | HAVANA | exon | NA | NA | exon:ENST00000456328.2:3 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000456328.2 | ENST00000456328.2 | processed_transcript | KNOWN | DDX11L1-002 | 1 | basic | OTTHUMT00000362751.1 | 3 | ENSE00002312635.1 | character(0) | NA | NA |
| chr1 | 12010 | 13670 | 1661 | + | HAVANA | transcript | NA | NA | ENST00000450305.2 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENSG00000223972.5 | ENST00000450305.2 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1-001 | NA | basic | OTTHUMT00000002844.2 | NA | NA | c("PGO:0000005", "PGO:0000019") | NA | NA |
| chr1 | 12010 | 12057 | 48 | + | HAVANA | exon | NA | NA | exon:ENST00000450305.2:1 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000450305.2 | ENST00000450305.2 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1-001 | NA | basic | OTTHUMT00000002844.2 | 1 | ENSE00001948541.1 | c("PGO:0000005", "PGO:0000019") | NA | NA |
| chr1 | 12179 | 12227 | 49 | + | HAVANA | exon | NA | NA | exon:ENST00000450305.2:2 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000450305.2 | ENST00000450305.2 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1-001 | NA | basic | OTTHUMT00000002844.2 | 2 | ENSE00001671638.2 | c("PGO:0000005", "PGO:0000019") | NA | NA |
| chr1 | 12613 | 12697 | 85 | + | HAVANA | exon | NA | NA | exon:ENST00000450305.2:3 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000450305.2 | ENST00000450305.2 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1-001 | NA | basic | OTTHUMT00000002844.2 | 3 | ENSE00001758273.2 | c("PGO:0000005", "PGO:0000019") | NA | NA |
| chr1 | 12975 | 13052 | 78 | + | HAVANA | exon | NA | NA | exon:ENST00000450305.2:4 | ENSG00000223972.5 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1 | 2 | OTTHUMG00000000961.2 | ENST00000450305.2 | ENST00000450305.2 | transcribed_unprocessed_pseudogene | KNOWN | DDX11L1-001 | NA | basic | OTTHUMT00000002844.2 | 4 | ENSE00001799933.2 | c("PGO:0000005", "PGO:0000019") | NA | NA |
1.6 Go ! O_O
2,5 millions de lignes !
29 colonnes !
Voilà un bien joli jeu de données. Presque un gros jeu de données, même !
Les types de gènes
Combien de gènes dans le génome humain : un par ligne ? 2.5 millions de gènes ? Sans doute pas : regardons d’un peu plus près la colonne type. Pour cela, nous allons produire un petit histogramme comptant l'occurrence de chaque valeur possible de type dans le tableau. Et là, surprise, pas besoin de dplyr : ggplot2 s’en sort très bien tout seul.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
# le code minimal pour produire la figure qui nous intéresse est assez simple : # on passe notre tableau de donnée en 1er paramètre de ggplot2 # puis on définit quelle colonne va donner quoi (axe x, y, z, contour, couleur, taille) dans aes() ggplot(gencode, aes(x = type)) + # enfin, on choisit un geom # stat = "count" est inutile: c'est la valeur par défaut de geom_bar # donc, ggplot2 va compter l’effectif de chaque type de gènes geom_bar(stat = "count") # mais il faut se rendre a l'évidence : la figure généré est moche: :( # pour la rendre plus jolie, il nous faut alourdir notre code: # trions les catégories dans un ordre a priori logique fig1 <- ggplot(gencode, aes(x = factor(type, levels = rev(unique(type))))) + geom_bar(fill = "darkorange") + # affichons l’effectif sur chaque bar, avec le mot magique de ggplot2 ..count.. geom_text(aes(label = ..count..), y = 10000, hjust = 0, stat = "count") + labs(x = "", y = "", title = "Gencode v25") + # renversons les axes coord_flip() # exportons le plot ggsave("figure1.svg", fig1 + theme(plot.margin = unit(c(5,8,1,1), "mm")), width = 18, height = 7, units = "cm", device = svglite) # et libérons de la mémoire vive rm(fig1) |
Figure 1: décomposition du tableau d'annotation GENCODE. Beaucoup de lignes, beaucoup d’éléments.
GENCODE compte donc 58 037 gènes. C’est beaucoup, on entend plutôt parler de 20 ou 25 000 gènes d’habitude. Une petite division nous donnera un nombre moyen de 3,4 transcrits différents par gènes, et de 5,9 exons par transcrits. Il y a un peu plus de codons start que de codons stop. Plus ou moins lié, il y a un peu plus de régions non traduites en 5’ qu’en 3’. Nous trouvons aussi 117 codons stop qui en fait n’en sont pas, vu qu’ils codent pour une sélénocystéine, parce que le code génétique universel, ha ha ha, laissez moi rire.
Il est temps de s'intéresser à la colonne gene_type (types de gènes), en restreignant l’analyse aux lignes du tableau de type gene uniquement.
|
1 2 3 4 5 6 7 8 9 10 11 |
# on ne guarde que les lignes du tableau dont le type est "gene" fig2 <- filter(gencode, type == "gene") %>% # trions les types de gènes en fonction de leur effectif ggplot(aes(x = factor( gene_type, levels = names(sort(table(gene_type))) ))) + geom_bar(stat = "count", fill = "darkorange") + geom_text(aes(label = ..count..), y = 10, hjust = 0, stat = "count") + labs(x = "", y = "", title = "Le génome humain") + coord_flip() |
Figure 2 : Nombre de gènes en fonction de leur types. Où l’on découvre que ça n'embête personne chez GENCODE de créer des nouveaux types de gènes pour UN SEUL gène.
GENCODE recense 19 950 gènes protéiques, un nombre plus raisonnable (même si il est un peu faible pour notre égo d'espèce ;-)). Pas mal de pseudogènes trainent un peu partout, ce qui ne fait pas très propre pour un génome d'espèce dominante. Notez aussi les nombreuses catégories de gènes en lien avec le système de recombinaison VDJ, parce que les immunologues se complaisent à noyer sous une nomenclature imbitable absolument tout ce qu’ils touchent. Détaillons quelques petites catégories mystérieuses :
- antisense : “Has transcripts that overlap the genomic span (i.e. exon or introns) of a protein-coding locus on the opposite strand.” (Gène dont les transcrits chevauchent les introns et ou exons d’un gène protéique codé dans le brin opposé).
- TEC : To be Experimentaly Confirmed: gène en attente de confirmation expérimentale. Pauvres 1048 gènes putatifs auxquels personne ne semble vouloir s’intéresser. 🙁
- scaRNA, sRNA, vaultRNA, scRNA, macro_lnc_RNA : parce que bon, il ne faudrait pas non plus mettre tous les gènes ARN dans le même sac. Il y a des différences importantes entre les différentes catégories de gènes ARN, donc il faut bien tout ranger soigneusement dans autant de petites cases que nécessaires. La classification fine des ARN, c’est du sérieux, du travail d’expert. D'orfèvre, dirais-je même. Ne mélangeons pas tout, un peu de rigueur. Quoi la catégorie misc_RNA (de miscellaneous: divers en Anglais) de 2213 membres ? J’entends mal. Mettre le vaultRNA, le scRNA, et le macro_lcnRNA dans la catégorie misc_RNA ? Nan.
Bref, c’est un peu le bazar cette classification, alors pour la suite de l’article, faisons nous notre petite classification grossière (libre à vous de l'adapter pour mettre en évidence votre type de gène préféré !):
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
gencode$simple_gene_type <- gencode$gene_type # si il y a pseudogene dans le nom, c'est un pseudogène gencode$simple_gene_type[grepl("pseudogene", gencode$simple_gene_type)] <- "pseudogène" # si il y a ARN dans le nom, c'est un gène ARN gencode$simple_gene_type[grepl("RNA", gencode$simple_gene_type)] <- "gène ARN" # si c'est un autre type de gène, c'est un autre type de gène gencode$simple_gene_type[!(gencode$simple_gene_type %in% c("protein_coding", "pseudogène", "gène ARN"))] <- "autre" # en Français, c'est plus jolie gencode$simple_gene_type[gencode$simple_gene_type == "protein_coding"] <- "gène protéique" # trions par effectif, ça simplifiera le code des figures gencode <- mutate(gencode, simple_gene_type = factor( simple_gene_type, levels = names(sort(table(simple_gene_type), decreasing = TRUE)) )) fig3 <- filter(gencode, type == "gene") %>% ggplot(aes(x = factor(simple_gene_type, levels = rev(levels(simple_gene_type))), fill = simple_gene_type)) + geom_bar(stat = "count") + geom_text(aes(label = ..count..), y = 100, hjust = 0, stat = "count") + labs(x = "", y = "", title = "Le génome humain") + coord_flip() + theme(legend.position="none") |
Figure 3 : Nombres de gènes en fonction de leur type, mais des types plus simples, parce que sinon c’est trop touffu. Il y a autant de gènes ARN que de pseudo-gènes… hum, je me demande si… heu, non, rien.
Et le premier qui vient me dire que certaines catégories que j’ai mis dans “autre” seraient mieux dans “gène ARN” peut aller se retaper cette analyse en Rust.
dplyr et l'évaluation non standard
dplyr, comme ggplot2, est un adepte de l'évaluation non standard des paramètres. Il faut donc indiquer le nom des colonnes sur lequel on travaille sans apostrophe, R cherchant alors cet objet d'abord au sein de l’environnement du tableau de données, avant de chercher dans l’environnement global. Un des inconvénients de ce comportement (oui, je sais.), c'est que si le nom de la colonne en question est dans une variable, la fonction risque de ne plus fonctionner. Il existe pour cela des versions à évaluation standard des paramètres pour chaque fonction de dplyr. Ces fonctions ont le même nom mais se finissent par "_". En résumé :
|
1 2 3 4 5 6 7 8 9 |
myVar <- "seqnames" select(gencode, seqnames ) # fonctionne select(gencode, "seqnames") # erreur select(gencode, myVar ) # erreur select_(gencode, seqnames ) # erreur select_(gencode, "seqnames") # fonctionne select_(gencode, myVar ) # fonctionne |
L'origine des annotations
Mais d'où viennent toutes ces annotations de gènes ?
|
1 2 3 4 5 6 |
fig4 <- filter(gencode, type == "gene") %>% # la colonne source contient la source des annotations ggplot(aes(x = source, fill = simple_gene_type)) + geom_bar(stat = "count") + labs(x = "", y = "", title = "Origine des annotations", fill = "") + theme(axis.text.x = element_text(angle = 45, hjust = 1)) |
Figure 4: Source des annotations en fonction du type de gène. HAVANA: consortium d’annotation manuelle. Ensembl: annotation algorithmique.
Je dois dire que je suis assez impressionné. L'annotation de presque tous les pseudogènes et “autre” gènes a été validée plus ou moins manuellement. Il ne reste qu’une moitié de gènes ARN prédite algorithmiquement, mais non validée manuellement. Du beau travail !
Répartition des gènes par chromosomes
Jusque là, dplyr ne nous a pas beaucoup servi, à peine à filtrer un tableau selon la valeur contenue dans une de ses colonnes. Nous allons maintenant monter un peu en complexité.
Par exemple, étudions la répartition de nos gènes sur nos chromosomes. Nous comparerons aussi la densité en gène de chaque chromosome. Pour cela nous récupérerons la longueur de chaque chromosome en lisant un petit fichier texte sur les serveurs de l'UCSC.
|
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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
# création d'une table contenant le nombre de gènes par chromosome, pour chaque type simplifié de gènes genes_by_chr <- filter(gencode, type == "gene") %>% # table() compte les effectifs pour chaque combinaison with(., table(simple_gene_type, seqnames)) %>% # tidy transforme l’affreux objet retourné par table() en un joli tableau tidy fig5A <- ggplot(genes_by_chr, aes(x = seqnames, y = Freq, fill = simple_gene_type)) + # pour une fois, on ne veut pas stat = "count" geom_bar(stat = "identity") + labs(x = "", y = "", title = "Gènes par chromosomes") + theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5)) # lecture du fichier texte directement depuis le serveur chr_length <- read_tsv("http://hgdownload-test.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes", col_names = FALSE) %>% # on nome les colonnes dplyr::rename(seqnames = X1, length = X2) %>% # on ne guarde que les chromosomes standards filter(seqnames %in% c(paste0("chr", 1:22), "chrX", "chrY", "chrM")) %>% # que l'on ordonne mutate(seqnames = factor(seqnames, levels = c(paste0("chr", 1:22), "chrX", "chrY", "chrM"))) fig5B <- ggplot(chr_length, aes(x = seqnames, y = length)) + geom_bar(stat = "identity", fill = "darkorange") + labs(x = "", y = "", title = "Taille des chromosomes (pb)") + theme(axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5)) # fusion des deux tableau grâce a dplyr::left_join genes_by_chr <- left_join(genes_by_chr, chr_length, by = "seqnames") %>% # création d'une nouvelle colonne: le nombre de gènes par mégabase mutate(perMb = Freq * 1E6 / length) fig5C <- ggplot(genes_by_chr, aes(x = seqnames, y = perMb, fill = simple_gene_type)) + geom_bar(stat = "identity") + coord_cartesian(ylim = c(0, 25)) + # une facette par type de gènes facet_wrap(~simple_gene_type, ncol = 1) + labs(x = "", y = "", title = "Gènes par mégabase par chromosomes") + theme(legend.position = "none", axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5)) # assemblage des panneaux avec cowplot # la zone fait 1 x 1, avec le point (0,0) en bas à gauche fig5 <- ggdraw() + draw_plot(fig5A, x = 0 , y = 0.5, w = 0.5, h = 0.5) + draw_plot(fig5B, x = 0 , y = 0 , w = 0.5, h = 0.5) + draw_plot(fig5C, x = 0.5, y = 0 , w = 0.5, h = 1 ) + draw_plot_label(c("A", "B", "C"), c(0, 0, 0.5), c(1, 0.5, 1), size = 20) |
Figure 5: Répartition des gènes entre les différents chromosomes. A. Nombre de gènes par chromosomes. B. Taille de chaque chromosome. C. Nombre de gènes par mégabase par chromosome. Le chromosome M est hors catégories, avec 785 gènes protéiques par mégabase et 1449 gènes ARN par mégabase (il ne fait que 16.5 kilobases).
Remarquons que la densité en gènes protéiques est assez variable d’un chromosome à l’autre, le chromosome 13 est assez pauvre, le chromosome 19 bien plus riche. La densité en gènes ARN et en pseudogènes semble plus homogène. Il est amusant de comparer le destin du chromosome Y, à transmission uniquement paternelle, avec celui du chromosome M, à transmission uniquement maternelle. Le chromosome Y a perdu presque tous ses gènes, il ne reste en grande partie que des pseudogènes (mais sa densité en pseudogènes n’a rien de scandaleusement élevée par rapport aux autosomes), et des régions intergéniques, lui conférant une faible densité en gène. Le chromosome M, qui a également subi une perte de gènes au cours de son évolution, est cependant hors classe en terme de densité génique avec plus de 500 gènes protéiques par mégabase ! Il faut dire qu'il a commencé son histoire avec très peu de régions inter-géniques.
Saluons aussi le bel effort de nomenclature: les autosomes humains étaient censés être numérotés par ordre décroissant de leur taille. Malheureusement, l’estimation a eu quelques légers ratés : le chromosome 11 est plus long que le 10, le 20 plus long que le 19, et le 22 plus long que le 21. Rha, et c'est trop tard pour les renommer, maintenant !
Taille des gènes
Nous l’avons vu plus haut, il y a en moyenne 3.4 transcrits par gène, et 5.9 exons par transcrits. Rentrons dans le détail de ces distributions. Nous allons aussi nous intéresser au nombre de fois que l'annotation d'un gène a été modifié. Cette information est contenue dans la dernière partie du code Ensembl du gène. Par exemple, mon gène préféré, MBD2 (ENSG00000134046.11) a vu son annotation modifiée 11 fois.
|
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 34 35 36 37 38 39 40 41 42 43 44 45 |
fig6A <- filter(gencode, type == "transcript") %>% # group_by groupe le tableau de donnée en fonction de colonnes particuliéres # les fonction dplyr appelées aprés un group_by s'appliquent individuellement pour chaque groupe group_by(gene_id, simple_gene_type) %>% # en l'occurence, pour chaque gène, on demande le nombre de transcrit via la fonction n() summarise(n_transcript = n()) %>% # pour rendre le plot plus lisible, la valeur maximale sera 20 transcrits mutate(n_transcript = if_else(n_transcript > 20, 20L, n_transcript)) %>% ggplot(aes(x = n_transcript, fill = simple_gene_type)) + geom_bar() + scale_x_continuous(breaks = c(1, 10, 20), labels = c(1, 10, "\u2265 20")) + facet_wrap(~simple_gene_type, scales = "free_y", ncol = 1) + theme(legend.position = "none") + labs(x = "Nombre de\ntranscrits", y = "Nombre de gènes", title = "Transcrits\npar gènes") fig6B <- filter(gencode, type == "exon") %>% mutate(exon_number = as.integer(exon_number)) %>% group_by(gene_id, simple_gene_type) %>% # cette fois, le nombre d'exon d'un gènes et définie comme le maximum du numéro d'exon de ce gène summarise(n_exon = max(exon_number)) %>% mutate(n_exon = if_else(n_exon > 20, 20L, n_exon)) %>% ggplot(aes(x = n_exon, fill = simple_gene_type)) + geom_bar() + scale_x_continuous(breaks = c(1, 10, 20), labels = c(1, 10, "\u2265 20")) + facet_wrap(~simple_gene_type, scales = "free_y", ncol = 1) + theme(legend.position = "none") + labs(x = "Nombre\nd'exons", y = "Nombre de gènes", title = "Exons\npar gènes") fig6C <- filter(gencode, type == "gene") %>% mutate( # on extrait le numéro de version de l'annotation de chaque gène gene_version = strsplit(gene_id, ".", fixed = TRUE) %>% sapply(last) %>% as.integer() ) %>% mutate(gene_version = if_else(gene_version > 20, 20L, gene_version)) %>% ggplot(aes(x = gene_version, fill = simple_gene_type)) + geom_bar() + scale_x_continuous(breaks = c(1, 10, 20), labels = c(1, 10, "\u2265 20")) + facet_wrap(~simple_gene_type, scales = "free_y", ncol = 1) + theme(legend.position = "none") + labs(x = "Nombre de\nrévisions", y = "Nombre de gènes", title = "Révisions\npar gènes") # plot_grid est une version simplifié de ggdraw(), dans lequel les plots sont rangés dans une grille fig6 <- plot_grid(fig6A, fig6B, fig6C, labels = LETTERS[1:3], ncol = 3, label_size = 20, hjust = -3) |
Figure 6: A. Nombre de gènes en fonction de leur nombre de transcrits. B. Nombre de gènes en fonction de leur nombre d’exons. C. Nombre de gènes en fonction du nombre de fois où leur annotation a été corrigée.
Les graphes sont éloquents : les gènes protéiques ont souvent plusieurs transcrits différents annotés (plus de 1000 gènes protéiques ont plus de 20 transcrits différents annotés !), tandis que les autres catégories se contentent dans leur plus grande majorité d’un seul transcrit. Il en va de même pour le nombre d'exons.
Là où c’est un peu plus inquiétant, c’est quand on regarde le nombre de fois que l’annotation d’un gène a été modifiée par GENCODE: la majorité des gènes protéiques ont vu leur annotation modifiée plus de dix fois, tandis les gènes ARN, les pseudogènes et les “autres” gènes sont en général inscrits une fois dans GENCODE et plus personne n’y touche. Alors on peut se dire que :
- GENCODE est exceptionnellement doué pour annoter des gènes ARN, pseudogènes et autre gènes, et du coup l’annotation est parfaite du premier coup.
- GENCODE est exceptionnellement mauvais pour annoter des gènes protéiques, et du coup doit sans arrêt corriger toute les bêtises. Notons qu’ils ne sont pas très doués pour corriger non plus, vu que les corrections s’empilent les unes sur les autres. Notez aussi que nous ne somme qu'à la 25éme version de GENCODE, et que le nombre maximum de correction pour un gènes est donc 25…
- Mais en fait, GENCODE se base sur le travail de la communauté. Oui, sur votre travail à vous tous. Et il faut le dire, autant pour travailler sur des gènes protéiques, il y a du monde, autant lorsqu’il faut travailler sur des gènes ARN, des pseudogènes et autres, il n’y a plus personne ! Peut-être devrions nous faire un effort.
Quelle est la taille moyenne de nos gènes ? Et bien, cela dépend de comment on la mesure. On peut définir la longueur totale du gène entre son site d'initiation de la transcription et son site de terminaison de la transcription, ou bien ne guarder que la somme de la longueur de ses exons (c'est à dire sans les introns). Nous allons faire les deux. Nous allons même calculer la distribution de la longueur des introns. C'est un peu plus complexe car il n'y a pas de ligne correspondant aux introns dans notre fichier d'annotation ! Mais vous allez voir qu'avec quelques fonctions dplyr, nous allons nous en sortir.
|
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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
# la longueur exons + introns est triviale à obtenir grâce à la colonne width fig7A <- filter(gencode, type == "gene") %>% ggplot(aes(x = width, fill = simple_gene_type, color = simple_gene_type)) + # geom_denisty plot la distribution des longueurs geom_density(alpha = 0.2) + # on fixe la même échelle pour tout les plots scale_x_log10(limits = c(10, 10000000)) + annotation_logticks(sides = "b") + labs(x = "Paire de bases", y = "Densité", title = "Longueur des gènes (exons + introns)", fill = "", color = "") # pour la longueur des transcrits, on somme tout les exons de chaque transcrits fig7B <- filter(gencode, type == "exon") %>% # on groupe par transcrits group_by(gene_id, transcript_id, simple_gene_type) %>% summarise(transcript_length = sum(width)) %>% # puis par gènes group_by(gene_id, simple_gene_type) %>% # et on prend la longueur du transcrit médian summarise(gene_length = median(transcript_length)) %>% ggplot(aes(x = gene_length, fill = simple_gene_type, color = simple_gene_type)) + geom_density(alpha = 0.2) + scale_x_log10(limits = c(10, 10000000)) + annotation_logticks(sides = "b") + labs(x = "Paire de bases", y = "Densité", title = "Longueur des gènes (exons seuls)", fill = "", color = "") # la longueur des exons est triviale à obtenir fig7C <- filter(gencode, type == "exon") %>% ggplot(aes(x = width, fill = simple_gene_type, color = simple_gene_type)) + geom_density(alpha = 0.2) + scale_x_log10(limits = c(10, 10000000)) + annotation_logticks(sides = "b") + labs(x = "Paire de bases", y = "Densité", title = "Longueur des exons", fill = "", color = "") # pour la longueur des introns, il faut ruser ! fig7D <- filter(gencode, type == "exon") %>% # on groupe par transcrit group_by(gene_id, transcript_id, simple_gene_type) %>% # pour chaque transcrit, on trie bien chaque exons dans l'ordre arrange(start) %>% # nouvelle colonne qui contient le start de la ligne du dessous, cad de l'exon suivant, # via la fonction lead(). # puis, dans le même mutate, nouvelle colonne qui contient la différence entre la fin # d'un exon, et le début de l'exon suivant. mutate(next_exon_start = lead(start), intron_length = next_exon_start - end) %>% ggplot(aes(x = intron_length, fill = simple_gene_type, color = simple_gene_type)) + geom_density(alpha = 0.2) + scale_x_log10(limits = c(10, 10000000)) + annotation_logticks(sides = "b") + labs(x = "Paire de bases", y = "Densité", title = "Longueur des introns", fill = "", color = "") fig7 <- plot_grid(fig7A, fig7B, fig7C, fig7D, labels = LETTERS[1:4], ncol = 1, label_size = 20) |
Figure 7: distribution de la longueur des gènes (A), des transcrits après épissage (B), des exons (C) et des introns (D).
La taille totale des gènes (exons + introns) est bien différente en fonction du type de gènes : les gènes protéiques sont les plus gros, suivit des "autres", des pseudogènes, et des gènes ARN et leur jolie distribution bimodale (quoi j'aurais du faire deux groupes de gènes ARN !?). La longueur des transcrits (soit la somme de la longueur des exons d'un gène) présente des différences bien moins marquées, même si on y voit toujours une belle distribution bimodale pour les gènes ARN. Pour ce qui est de la longueur des exons et des introns, les distributions sont quasi identiques pour chaque types de gènes. Les gènes protéiques sont donc plus gros car ils ont plus d'introns et plus d'exons que les autres (cf figure 6B), mais leurs exons et leurs introns ont une taille tout à fait standard.
Le modulo 3 des CDS
Chaque CDS (coding DNA sequence) possède sa petite ligne dans le fichier fournit par GENCODE. Les codons mesurant 3 nucléotides, il est peut-être intéressant de regarder le modulo 3 de la longueur des CDS ? Nous allons comparer le modulo 3 des CDS contenu dans chaque exons, mais aussi le modulo 3 des CDS pour chaque gène.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
# on ne garde que les lignes de type "CDS" fig8A <- filter(gencode, type == "CDS") %>% # nouvelle colonne contenant le modulo 3 mutate(modulo_3 = width %% 3) %>% ggplot(aes(x = modulo_3)) + geom_bar(fill = "darkorange") + labs(x = "", y = "effectif", title = "Modulo 3 des\nCDS (exons)") fig8B <- filter(gencode, type == "CDS") %>% # on groupe par transcrit group_by(gene_id, transcript_id, simple_gene_type) %>% # on somme la taille de chaque CDS pour chaque transcrit summarise(total_CDS = sum(width)) %>% mutate(modulo_3 = total_CDS %% 3) %>% group_by(gene_id, simple_gene_type) %>% # pour chaque gène, on prends la mediane des modulos 3 des transcrits # qu'on arrondit au plus bas en cas de chiffre non entier summarise(modulo_3 = floor(median(modulo_3))) %>% ggplot(aes(x = modulo_3)) + geom_bar(fill = "darkorange") + labs(x = "", y = "effectif", title = "Modulo 3 des\nCDS (gènes)") fig8 <- plot_grid(fig8A, fig8B, labels = LETTERS[1:2], ncol = 2, label_size = 20) |
Figure 8: modulo 3 des séquences codantes. A. au niveau des exons. B. au niveau des gènes.
La grande majorité des séquences codantes des gènes est divisible par trois. Ouf ! Tout va bien. Il reste quelques gènes dont la séquences codante est non divisible par trois. Il faudrait regarder plus en détail : pseudogènes ? Annotations d'isoformes non sens dégradées ? Ou bien erreurs d'annotations ?
Les séquences codantes des exons sont réparties presque équitablement entre les trois possibilité de modulo 3. Donc si un exon codant d'un gène est inséré dans l'intron d'un autre gène, il y a grossièrement deux chance sur trois de casser la phase du gène d'arrivé. Ce qui rend le ré-assemblage évolutif d'exons un poil compliqué. Le léger enrichissement en exons codant divisible par trois est sans doute lié à la présence de gènes dont la séquence codante est contenue dans un seul exons.
Les plus grandes (pseudo) familles de gènes
Dernière analyse de ce billet, intéressons nous aux pseudo-familles des gènes. Je ne vais pas faire ici de comparaison de séquences ou de domaines protéiques, mais une analyse beaucoup plus simple : regarder les gènes dont le nom est le même à l’exception du numéro final (par exemple MBD1, MBD2, MBD3, etc.). C'est bien entendu une analyse grossière, relevant bien plus de nos biais de nommage que d'une éventuelle relation paralogue.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
fig9 <- filter(gencode, type == "gene") %>% # on enléve les chiffres à la fin mutate(family_name = sub("[0-9]*$", "", gene_name)) %>% # on compte count(family_name, gene_type = gene_type) %>% # on ungroup avant de trier ungroup %>% arrange(desc(n)) %>% # on ne guarde que les familles de plus de 100 membres filter(n >= 100) %>% ggplot(aes(x = factor(family_name, levels = rev(family_name)), y = n, fill = gene_type)) + geom_bar(stat = "identity") + geom_text(aes(label = gene_type, y = 10), hjust = 0, color = "grey40", size = 3) + coord_flip() + labs(x = "", y = "Nombre de membres", title = "Pseudo-familles de gènes") + theme(legend.position = "none") |
Figure 9: Noms de gènes les plus fréquents.
Nous découvrons donc que un gros tas de miRNA s'appellent MIR, des lincRNA qui s'appelent LINC, des snoRNA qui s'appellent SNORA et SNORD, quelle surprise ! Vous m'étonnez qu'avec une nomenclature pareille personne ne veuille travailler dessus... Les ARN Y sont apparemment assez nombreux, je ne les connaissais pas, mais Wikipedia si. C'est étonnement intéressant. Même histoire pour les ARN 7SK, vous connaissiez, vous ? Autres ARN dont j'ignorai l’existence et l'importance: les ARN SRP. Comme quoi, on ne nous dit pas tout.
Trois grands groupes de gènes protéiques ont plus de 100 membres ayant le même nom : les protéines à doigt de zinc (ZNF - zinc finger proteins), ou potentiellement autant de facteurs de transcription auquel personne ne s'est assez intéressé pour leur donner un nom moins générique, les protéines transmembranaires (TMEM - transmembrane protein), et enfin les protéines contenant une superhélice (CCDC - coiled coil domain containing protein).
Je le répète : ce classement par popularité du nom des gènes ne reflète pas l’effectif des familles de gènes. Il manque par exemple la famille des récepteurs olfactifs, dont les membres ont un nom différent en fonction de leur sous-famille.
Exercices
dplyr vous intrigue, et vous avez envie d'approfondir le sujet ? Voici quelques questions que vous pourriez résoudre. N'hésitez pas à partager votre réponse en commentaire !
- Sur quels chromosomes sont situés les 117 codons stop qui codent pour une sélénocystéine ?
- Il y a t'il une différence notable entre la distribution des tailles des 5'UTR et celle des 3'UTR ? Quelle est la proportion de 5'UTR et de 3'UTR qui s'étendent sur plus d'un exons ?
- Pour chaque chromosome, quelle est la répartition des gènes codés sur le brin + et sur le brin - ? Notez-vous quelques cas extrêmes ? Il y a-t-il une différence entre types de gènes ?
- Quels sont donc ces gènes protéiques dont la longueur de la séquence codante n'est pas divisible par trois ?
Conclusion
J’espère que cet article vous aura appris quelques informations pertinentes sur l'état actuel de l'annotation des gènes du génome humain. Ce fichier d'annotation est une merveilleuse réussite collective, mais des zones d'ombres restent à être éclaircies par la communauté. Comment, en tant que civilisation, pouvons nous accepter que plusieurs milliers de gènes protéiques du génome humain restent non étudiés, alors que nous n'en possédons que 20 000 aux dernières nouvelles ?
Certains auront peut-être réalisé le potentiel de la combinaison dplyr + ggplot2 pour travailler sur des données tabulaires. En s'y mettant petit à petit, avec les feuilles d'astuces sous les yeux, c'est vraiment un plaisir. Le package data.table est une alternative populaire à dplyr. data.table est généralement plus rapide, avec une philosophie et une syntaxe différente. Le super auteur et relecteur Lelouar a reproduit cette analyse en utilisant data.table et en mesurant les temps d'execution, vous pouvez donc vous faire une idée en regardant ces deux notebook Jupyter.
Lors de la rédaction de cet article, je suis tombé sur ce billet de blog en anglais, reproduisant certaines de ces analyses avec Python et Panda. Un de nos rédacteurs préférés a aussi extrait quelques chiffres similaires en utilisant bash et quelques utilitaires.
Milles merci aux relecteurs et reléctrice Lelouar, Akira et Bebatut ! <3 <3 <3
Devailly Elisabeth
janvier 26, 2017 à 9:08
Article bien écrit et accessible à tous ..... Surtout par l'approche rigolote.
Mais bon dans une famille pas très utile de charger ce fichier un peu lourd en mémoire vive. !
ZaZo0o
août 16, 2017 à 2:25
Merci Guillaume !
J'ai adapté à mes données et maintenant j'ai une jolie figure pour mon article \o/
[Remarque] j'ai utilisé un fichier GTF et non GFF, du coup je l'ai chargé comme ça:
gencode <- import("gencode.annotation.gtf.gz", format = "GTF")
Bastien
octobre 17, 2017 à 2:20
Merci pour le post qui est très agréable à suivre !
Pour créer la figure 6B, la fonction summarise me génère un n_exon de type double alors que la fonction mutate juste après fait des comparaisons d'entier.
J'ai modifié la ligne
"summarise(n_exon = max(exon_number))"
par
"summarise(n_exon = as.integer(max(exon_number)))"
Si jamais c'est utile à quelqu'un 🙂