Suivez l'guide :
dplyr et le génome humain

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 :

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 !

%>%: 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 :

Le raccourcis . permet de définir la position de l’objet dans le cas où il n’est pas le premier argument :

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é.

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.

Figure 1: décomposition du tableau d'annotation GENCODE. Beaucoup de lignes, beaucoup d’éléments.

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.

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.

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é !):

Figure 3 : Nombres de gènes en fonction de leur type, mais des types plus simples, parce que sinon c’est trop touffu. Autant de gènes ARN que de pseudo-gènes… hum, je me demande si… heu, non, rien.

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é :

L'origine des annotations

Mais d'où viennent toutes ces annotations de gènes ?

Figure 4: Source des annotations en fonction du type de gène. HAVANA: consortium d’annotation manuelle. Ensembl: annotation algorithmique.

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.

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 et 1449 gènes ARN par mégabase (il ne fait que 16.5 kilobases).

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.

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’exon. C. Nombre de gène en fonction du nombre de fois où leur annotation a été corrigé.

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.

Figure 7: distribution de la longueur des gènes (A), des transcrits après épissage (B), des exons (C) et des introns (D).

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.

Figure 8: modulo 3 des séquences codantes

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.

Figure 9: Noms de gènes les plus fréquents.

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 !

  1. Sur quels chromosomes sont situés les 117 codons stop qui  codent pour une sélénocystéine ?
  2. 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 ?
  3. 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 ?
  4. 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

 

  • À 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

3 commentaires sur “dplyr et le génome humain

  1. 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. !

  2. 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")

  3. 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 🙂

Laisser un commentaire