dplyr et le génome humain

Introduction

Non, ne fuyez pas tout de suite, chers lec­teurs, tout va s'éclaircir : dplyr, c’est plyr pour les data.frame (les tableaux de don­nées). Atten­dez, j’y viens, plyr, c’est un package R pour appli­quer (apply) des fonc­tions. Donc, dplyr (pro­non­cez “diplir”), c’est un package R, pour appli­quer des fonc­tions à un tableau de don­nées.

Et ça, mes amis, c'est super cool.

Ren­dez-vous ser­vice : impri­mez immé­dia­te­ment cette feuille d’astuces le concer­nant (pdf, 0.5Mo). En cou­leur. En fait impri­mez en 10, et don­nez en 8 à vos col­lègues, et met­tez en une chez vous, et une sur votre lieu de tra­vail. dplyr est cen­tré autour du concept de don­nées propres (ou tidy) : chaque ligne est une obser­va­tion, chaque colonne un para­mètre obser­vé. Finis­sons notre petite digres­sion en évo­quant broom, un package qui per­met de rendre propres (tidy) cer­tains retour de fonc­tions sales (mes­sy) de R, via la fonc­tion

tidy()

 .

Pour vous ini­tier à dplyr, je vous pro­pose une petite série d'exercices appli­qués à notre génome à tous. Ici on ne par­le­ra pas tant du génome dans sa défi­ni­tion moderne (ensemble des séquences d’ADN d’un orga­nisme), mais dans une ver­sion plus éty­mo­lo­gique, et plus désuète : ensemble des gènes d’un orga­nisme (ici, l’humain).

Si vous n’en avez rien à faire de dplyr, ou a for­tio­ri de R, ne fuyez pas pour autant : vous per­driez une occa­sion de mieux connaître votre génome !

GENCODE et l’annotation du génome humain

Nous allons uti­li­ser les don­nées d’annotations de GENCODE, un consor­tium pour l’annotation du génome humain et du génome murin, notam­ment uti­li­sé par Ensem­bl.

Pour télé­char­ger le der­nier fichier d’annotation dis­po­nible, il suf­fit de se rendre sur leur site inter­net, et de sélec­tion­ner Data > Human > Cur­rent release. Lors de l’écriture de cet article, la der­nière ver­sion est la 25, cor­res­pon­dant aux coor­don­nées géno­miques GRCh38 (aus­si appe­lé hg38). Plu­sieurs fichiers sont dis­po­nibles mais un seul suf­fi­ra aujourd’hui, le pre­mier : “Com­pre­hen­sive gene anno­ta­tion | CHR”, en ver­sion GFF3. Le lec­teur curieux se deman­de­ra quelle est la dif­fé­rence entre GTF et GFF3 ? C’est com­pli­qué. Le GTF est une variante mineure du GFF2. Le GFF3 est plus récent et mieux spé­ci­fié, nous allons donc uti­li­ser cette ver­sion. Vous devez télé­char­ger ce fichier (45Mo) si vous sou­hai­tez repro­duire chez vous l’analyse ci des­sous (ce que je vous conseille !). Le fichier GFF3 est un fichier texte, déli­mi­té par des tabu­la­tions, dont les colonnes suivent un for­ma­tage défi­ni.

Lec­teur du futur, peut-être aurez-vous la chance de dis­po­ser d’un fichier d’annotation plus récent ? Deux choix s’offrent alors à vous : refaire l’analyse avec les jolies anno­ta­tions du futur, vous pour­rez alors com­pa­rer avec l’état des connais­sances de 2017 et bien vous mar­rer ; ou bien télé­char­ger la ver­sion que nous uti­li­se­rons ici (le lien devrait res­ter valide — 45Mo).

Pré-requis

Une connais­sance basique de R, et notam­ment quelques notions de ggplot2, vous aide­ra à com­prendre le code pro­duit ici (pour s’initier à R, pour­quoi ne pas com­men­cer par ici, et pour ggplot2 par là ?). Je vais cepen­dant com­men­ter le code dans l'espoir de le rendre plus intel­li­gible. Si vous sou­hai­tez repro­duire l’analyse en paral­lèle de votre lec­ture, vous aurez besoin d’une ver­sion de R un peu récente (j'utilise la  3.3.1 “Bug in your hair”), et de quelques packages. L'installation de ces der­niers ne devrait pas être trop longue :

Il vous faut aus­si, mal­heu­reu­se­ment, pas mal de mémoire vive, je dirai au moins 3 ou 4Go pour R (j’ai écrit cet article sur un ordi­na­teur ayant 8Go). Un des défaut majeur de R étant qu’il tra­vaille essen­tiel­le­ment sur la mémoire vive… Oui, je sais, ça fait pas mal de défauts.

Analyse

Parsons le fichier GENCODE

Pour char­ger le fichier dans R, plu­tôt que d’écrire notre propre par­seur, nous allons uti­li­ser celui du package rtra­ck­layer en uti­li­sant la fonc­tion

import()

 . Cette fonc­tion retourne un objet du type GRanges, que nous allons conver­tir en data_​frame via la fonc­tion

as_​data_​frame()

  du package dplyr. Les data_​frame, (aus­si appe­lés tibble) (ou encore tbl_​df) (oui, je sais) de dplyr sont légè­re­ment dif­fé­rents des data.frame de R tout court : la méthode d’affichage (print) par défaut est plus sym­pa, et les row.names sont inter­dits !

%>%: le pipe R

Petite digres­sion à pro­pos du pipe :

%>%

 , inclus dans le package dplyr (et issus de magr­tit­tr). Il s’agit de l’équivalent R du pipe

|

  de shell/​bash. Il passe le résul­tat de la fonc­tion de gauche comme pre­mier argu­ment (par défaut) de la fonc­tion de droite :

Le rac­cour­cis

.

per­met de défi­nir la posi­tion de l’objet dans le cas où il n’est pas le pre­mier argu­ment :

Mes exemples bidons ne sont pas trop flat­teurs, mais ce petit su-sucre syn­taxique devient vite indis­pen­sable dès qu’on l'essaye. Il per­met d’accomplir des tâches com­plexes sans mul­ti­plier les variables inter­mé­diaires ou les fonc­tions nichées, et comme nous allons le voir, il se marie admi­ra­ble­ment bien avec les fonc­tions du package dplyr, qui prennent tou­jours comme pre­mier para­mètre un tableau de don­nées, et retournent tou­jours un tableau de don­nées.

Un bien beau jeu de données

Exa­mi­nons un peu le jeu de don­nées tout juste char­gé.

seq­names 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 ccd­sid
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 trans­cript 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 trans­cript 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 mil­lions de lignes !

29 colonnes !

Voi­là un bien joli jeu de don­nées. Presque un gros jeu de don­nées, même !

Les types de gènes

Com­bien de gènes dans le génome humain : un par ligne ? 2.5 mil­lions de gènes ? Sans doute pas : regar­dons d’un peu plus près la colonne type. Pour cela, nous allons pro­duire un petit his­to­gramme comp­tant l'occurrence de chaque valeur pos­sible de type dans le tableau. Et là, sur­prise, 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écom­po­si­tion du tableau d'annotation GENCODE. Beau­coup de lignes, beau­coup d’éléments.

GENCODE compte donc 58 037 gènes. C’est beau­coup, on entend plu­tôt par­ler de 20 ou 25 000 gènes d’habitude. Une petite divi­sion nous don­ne­ra un nombre moyen de 3,4 trans­crits dif­fé­rents par gènes, et de 5,9 exons par trans­crits. 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 tra­duites en 5’ qu’en 3’. Nous trou­vons aus­si 117 codons stop qui en fait n’en sont pas, vu qu’ils codent pour une sélé­no­cys­téine, parce que le code géné­tique uni­ver­sel, ha ha ha, lais­sez moi rire.

Il est temps de s'intéresser à la colonne gene_​type (types de gènes), en restrei­gnant l’analyse aux lignes du tableau de type gene uni­que­ment.

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 fonc­tion de leur types. Où l’on découvre que ça n'embête per­sonne chez GENCODE de créer des nou­veaux types de gènes pour UN SEUL gène.

GENCODE recense 19 950 gènes pro­téiques, un nombre plus rai­son­nable (même si il est un peu faible pour notre égo d'espèce ;-)). Pas mal de pseu­do­gènes trainent un peu par­tout, ce qui ne fait pas très propre pour un génome d'espèce domi­nante. Notez aus­si les nom­breuses caté­go­ries de gènes en lien avec le sys­tème de recom­bi­nai­son VDJ, parce que les immu­no­logues se com­plaisent à noyer sous une nomen­cla­ture imbi­table abso­lu­ment tout ce qu’ils touchent. Détaillons quelques petites caté­go­ries mys­té­rieuses :

  • anti­sense : “Has trans­cripts that over­lap the geno­mic span (i.e. exon or introns) of a pro­tein-coding locus on the oppo­site strand.” (Gène dont les trans­crits che­vauchent les introns et ou exons d’un gène pro­téique codé dans le brin oppo­sé).
  • TEC : To be Expe­ri­men­ta­ly Confir­med : gène en attente de confir­ma­tion expé­ri­men­tale. Pauvres 1048 gènes puta­tifs aux­quels per­sonne ne semble vou­loir s’intéresser. 🙁
  • scaR­NA, sRNA, vaul­tR­NA, scR­NA, macro_​lnc_​RNA : parce que bon, il ne fau­drait pas non plus mettre tous les gènes ARN dans le même sac. Il y a des dif­fé­rences impor­tantes entre les dif­fé­rentes caté­go­ries de gènes ARN, donc il faut bien tout ran­ger soi­gneu­se­ment dans autant de petites cases que néces­saires. La clas­si­fi­ca­tion fine des ARN, c’est du sérieux, du tra­vail d’expert. D'orfèvre, dirais-je même. Ne mélan­geons pas tout, un peu de rigueur. Quoi la caté­go­rie misc_​RNA (de mis­cel­la­neous : divers en Anglais) de 2213 membres ? J’entends mal. Mettre le vaul­tR­NA, le scR­NA, et le macro_​lcnRNA dans la caté­go­rie misc_​RNA ? Nan.

Bref, c’est un peu le bazar cette clas­si­fi­ca­tion, alors pour la suite de l’article,  fai­sons nous notre petite clas­si­fi­ca­tion gros­sière (libre à vous de l'adapter pour mettre en évi­dence 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 fonc­tion de leur type, mais des types plus simples, parce que sinon c’est trop touf­fu. Il y a autant de gènes ARN que de pseu­do-gènes… hum, je me demande si… heu, non, rien.

Et le pre­mier qui vient me dire que cer­taines caté­go­ries que j’ai mis dans “autre” seraient mieux dans “gène ARN” peut aller se reta­per cette ana­lyse en Rust.

dplyr et l'évaluation non standard

dplyr, comme ggplot2, est un adepte de l'éva­lua­tion non stan­dard des para­mètres. Il faut donc indi­quer le nom des colonnes sur lequel on tra­vaille sans apos­trophe, R cher­chant alors cet objet d'abord au sein de l’environnement du tableau de don­nées, avant de cher­cher dans l’environnement glo­bal. Un des incon­vé­nients de ce com­por­te­ment (oui, je sais.), c'est que si le nom de la colonne en ques­tion est dans une variable, la fonc­tion risque de ne plus fonc­tion­ner. Il existe pour cela des ver­sions à éva­lua­tion stan­dard des para­mètres pour chaque fonc­tion de dplyr. Ces fonc­tions ont le même nom mais se finissent par "_​". En résu­mé :

L'origine des annotations

Mais d'où viennent toutes ces anno­ta­tions 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 anno­ta­tions en fonc­tion du type de gène. HAVANA : consor­tium d’annotation manuelle. Ensem­bl : anno­ta­tion algo­rith­mique.

Je dois dire que je suis assez impres­sion­né. L'annotation de presque tous les pseu­do­gènes et “autre” gènes a été vali­dée plus ou moins manuel­le­ment. Il ne reste qu’une moi­tié de gènes ARN pré­dite algo­rith­mi­que­ment, mais non vali­dée manuel­le­ment. Du beau tra­vail !

Répartition des gènes par chromosomes

Jusque là, dplyr ne nous a pas beau­coup ser­vi, à peine à fil­trer un tableau selon la valeur conte­nue dans une de ses colonnes. Nous allons main­te­nant mon­ter un peu en com­plexi­té.

Par exemple, étu­dions la répar­ti­tion de nos gènes sur nos chro­mo­somes. Nous com­pa­re­rons aus­si la den­si­té en gène de chaque chro­mo­some. Pour cela nous récu­pé­re­rons la lon­gueur de chaque chro­mo­some en lisant un petit fichier texte sur les ser­veurs 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épar­ti­tion des gènes entre les dif­fé­rents chro­mo­somes. A. Nombre de gènes par chro­mo­somes. B. Taille de chaque chro­mo­some. C. Nombre de gènes par méga­base par chro­mo­some. Le chro­mo­some M est hors caté­go­ries, avec 785 gènes pro­téiques par méga­base et 1449 gènes ARN par méga­base (il ne fait que 16.5 kilo­bases).

Remar­quons que la den­si­té en gènes pro­téiques est assez variable d’un chro­mo­some à l’autre, le chro­mo­some 13 est assez pauvre, le chro­mo­some 19 bien plus riche. La den­si­té en gènes ARN et en pseu­do­gènes semble plus homo­gène. Il est amu­sant de com­pa­rer le des­tin du chro­mo­some Y, à trans­mis­sion uni­que­ment pater­nelle, avec celui du chro­mo­some M, à trans­mis­sion uni­que­ment mater­nelle. Le chro­mo­some Y a per­du presque tous ses gènes, il ne reste en grande par­tie que des pseu­do­gènes (mais sa den­si­té en pseu­do­gènes n’a rien de scan­da­leu­se­ment éle­vée par rap­port aux auto­somes), et des régions inter­gé­niques, lui confé­rant une faible den­si­té en gène. Le chro­mo­some M, qui a éga­le­ment subi une perte de gènes au cours de son évo­lu­tion, est cepen­dant hors classe en terme de den­si­té génique avec plus de 500 gènes pro­téiques par méga­base ! Il faut dire qu'il a com­men­cé son his­toire avec très peu de régions inter-géniques.

Saluons aus­si le bel effort de nomen­cla­ture : les auto­somes humains étaient cen­sés être numé­ro­tés par ordre décrois­sant de leur taille. Mal­heu­reu­se­ment, l’estimation a eu quelques légers ratés : le chro­mo­some 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 renom­mer, main­te­nant !

Taille des gènes

Nous l’avons vu plus haut, il y a en moyenne 3.4 trans­crits par gène, et 5.9 exons par trans­crits. Ren­trons dans le détail de ces dis­tri­bu­tions. Nous allons aus­si nous inté­res­ser au nombre de fois que l'annotation d'un gène a été modi­fié. Cette infor­ma­tion est conte­nue dans la der­nière par­tie du code Ensem­bl du gène. Par exemple, mon gène pré­fé­ré, MBD2 (ENSG00000134046.11) a vu son anno­ta­tion modi­fié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 fonc­tion de leur nombre de trans­crits. B. Nombre de gènes en fonc­tion de leur nombre d’exons. C. Nombre de gènes en fonc­tion du nombre de fois où leur anno­ta­tion a été cor­ri­gée.

Les graphes sont élo­quents : les gènes pro­téiques ont sou­vent plu­sieurs trans­crits dif­fé­rents anno­tés (plus de 1000 gènes pro­téiques ont plus de 20 trans­crits dif­fé­rents anno­tés !), tan­dis que les autres caté­go­ries se contentent dans leur plus grande majo­ri­té d’un seul trans­crit. 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é modi­fiée par GENCODE : la majo­ri­té des gènes pro­téiques ont vu leur anno­ta­tion modi­fiée plus de dix fois, tan­dis les gènes ARN, les pseu­do­gènes et les “autres” gènes sont en géné­ral ins­crits une fois dans GENCODE et plus per­sonne n’y touche. Alors on peut se dire que :

  • GENCODE est excep­tion­nel­le­ment doué pour anno­ter des gènes ARN, pseu­do­gènes et autre gènes, et du coup l’annotation est par­faite du pre­mier coup.
  • GENCODE est excep­tion­nel­le­ment mau­vais pour anno­ter des gènes pro­téiques, et du coup doit sans arrêt cor­ri­ger toute les bêtises. Notons qu’ils ne sont pas très doués pour cor­ri­ger non plus, vu que les cor­rec­tions s’empilent les unes sur les autres. Notez aus­si que nous ne somme qu'à la 25éme ver­sion de GENCODE, et que le nombre maxi­mum de cor­rec­tion pour un gènes est donc 25…
  • Mais en fait, GENCODE se base sur le tra­vail de la com­mu­nau­té. Oui, sur votre tra­vail à vous tous. Et il faut le dire, autant pour tra­vailler sur des gènes pro­téiques, il y a du monde, autant lorsqu’il faut tra­vailler sur des gènes ARN, des pseu­do­gènes et autres, il n’y a plus per­sonne ! Peut-être devrions nous faire un effort.

Quelle est la taille moyenne de nos gènes ? Et bien, cela dépend de com­ment on la mesure. On peut défi­nir la lon­gueur totale du gène entre son site d'initiation de la trans­crip­tion et son site de ter­mi­nai­son de la trans­crip­tion, ou bien ne guar­der que la somme de la lon­gueur de ses exons (c'est à dire sans les introns). Nous allons faire les deux. Nous allons même cal­cu­ler la dis­tri­bu­tion de la lon­gueur des introns. C'est un peu plus com­plexe car il n'y a pas de ligne cor­res­pon­dant aux introns dans notre fichier d'annotation ! Mais vous allez voir qu'avec quelques fonc­tions dplyr, nous allons nous en sor­tir.

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 : dis­tri­bu­tion de la lon­gueur des gènes (A), des trans­crits après épis­sage (B), des exons (C) et des introns (D).

La taille totale des gènes (exons + introns) est bien dif­fé­rente en fonc­tion du type de gènes : les gènes pro­téiques sont les plus gros, sui­vit des "autres", des pseu­do­gènes, et des gènes ARN et leur jolie dis­tri­bu­tion bimo­dale (quoi j'aurais du faire deux groupes de gènes ARN !?). La lon­gueur des trans­crits (soit la somme de la lon­gueur des exons d'un gène) pré­sente des dif­fé­rences bien moins mar­quées, même si on y voit tou­jours une belle dis­tri­bu­tion bimo­dale pour les gènes ARN. Pour ce qui est de la lon­gueur des exons et des introns, les dis­tri­bu­tions sont qua­si iden­tiques pour chaque types de gènes. Les gènes pro­té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 stan­dard.

Le modulo 3 des CDS

Chaque CDS (coding DNA sequence) pos­sède sa petite ligne dans le fichier four­nit par GENCODE. Les codons mesu­rant 3 nucléo­tides, il est peut-être inté­res­sant de regar­der le modu­lo 3 de la lon­gueur des CDS ? Nous allons com­pa­rer le modu­lo 3 des CDS conte­nu dans chaque exons, mais aus­si le modu­lo 3 des CDS pour chaque gène.

Figure 8: modulo 3 des séquences codantes
Figure 8 : modu­lo 3 des séquences codantes. A. au niveau des exons. B. au niveau des gènes.

La grande majo­ri­té des séquences codantes des gènes est divi­sible par trois. Ouf ! Tout va bien. Il reste quelques gènes dont la séquences codante est non divi­sible par trois. Il fau­drait regar­der plus en détail : pseu­do­gènes ? Anno­ta­tions d'isoformes non sens dégra­dées ? Ou bien erreurs d'annotations ?

Les séquences codantes des exons sont répar­ties presque équi­ta­ble­ment entre les trois pos­si­bi­li­té de modu­lo 3. Donc si un exon codant d'un gène est insé­ré dans l'intron d'un autre gène, il y a gros­siè­re­ment deux chance sur trois de cas­ser la phase du gène d'arrivé. Ce qui rend le ré-assem­blage évo­lu­tif d'exons un poil com­pli­qué. Le léger enri­chis­se­ment en exons codant divi­sible par trois est sans doute lié à la pré­sence de gènes dont la séquence codante est conte­nue dans un seul exons.

Les plus grandes (pseudo) familles de gènes

Der­nière ana­lyse de ce billet, inté­res­sons nous aux pseu­do-familles des gènes. Je ne vais pas faire ici de com­pa­rai­son de séquences ou de domaines pro­téiques, mais une ana­lyse beau­coup plus simple : regar­der 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 enten­du une ana­lyse gros­sière, rele­vant bien plus de nos biais de nom­mage que d'une éven­tuelle rela­tion para­logue.

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

Nous décou­vrons donc que un gros tas de miR­NA s'appellent MIR, des lin­cR­NA qui s'appelent LINC, des snoR­NA qui s'appellent SNORA et SNORD, quelle sur­prise ! Vous m'étonnez qu'avec une nomen­cla­ture pareille per­sonne ne veuille tra­vailler des­sus… Les ARN Y sont appa­rem­ment assez nom­breux, je ne les connais­sais pas, mais Wiki­pe­dia si. C'est éton­ne­ment inté­res­sant. Même his­toire pour les ARN 7SK, vous connais­siez, 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 pro­téiques ont plus de 100 membres ayant le même nom : les pro­téines à doigt de zinc (ZNF — zinc fin­ger pro­teins), ou poten­tiel­le­ment autant de fac­teurs de trans­crip­tion auquel per­sonne ne s'est assez inté­res­sé pour leur don­ner un nom moins géné­rique, les pro­téines trans­mem­bra­naires (TMEM — transmembrane pro­tein), et enfin les pro­téines conte­nant une super­hé­lice (CCDC — coiled coil domain contai­ning pro­tein).

Je le répète : ce clas­se­ment par popu­la­ri­té du nom des gènes ne reflète pas l’effectif des familles de gènes. Il manque par exemple la famille des récep­teurs olfac­tifs, dont les membres ont un nom dif­fé­rent en fonc­tion de leur sous-famille.

Exercices

dplyr vous intrigue, et vous avez envie d'approfondir le sujet ? Voi­ci quelques ques­tions que vous pour­riez résoudre. N'hésitez pas à par­ta­ger votre réponse en com­men­taire !

  1. Sur quels chro­mo­somes sont situés les 117 codons stop qui  codent pour une sélé­no­cys­téine ?
  2. Il y a t'il une dif­fé­rence notable entre la dis­tri­bu­tion des tailles des 5'UTR et celle des 3'UTR ? Quelle est la pro­por­tion de 5'UTR et de 3'UTR qui s'étendent sur plus d'un exons ?
  3. Pour chaque chro­mo­some, quelle est la répar­ti­tion des gènes codés sur le brin + et sur le brin — ? Notez-vous quelques cas extrêmes ? Il y a‑t-il une dif­fé­rence entre types de gènes ?
  4. Quels sont donc ces gènes pro­téiques dont la lon­gueur de la séquence codante n'est pas divi­sible par trois ?

Conclusion

J’espère que cet article vous aura appris quelques infor­ma­tions per­ti­nentes sur l'état actuel de l'annotation des gènes du génome humain. Ce fichier d'annotation est une mer­veilleuse réus­site col­lec­tive, mais des zones d'ombres res­tent à être éclair­cies par la com­mu­nau­té. Com­ment, en tant que civi­li­sa­tion, pou­vons nous accep­ter que plu­sieurs mil­liers de gènes pro­téiques du génome humain res­tent non étu­diés, alors que nous n'en pos­sé­dons que 20 000 aux der­nières nou­velles ?

Cer­tains auront peut-être réa­li­sé le poten­tiel de la com­bi­nai­son dplyr + ggplot2 pour tra­vailler sur des don­nées tabu­laires. En s'y met­tant petit à petit, avec les feuilles d'astuces sous les yeux, c'est vrai­ment un plai­sir. Le package data.table est une alter­na­tive popu­laire à dplyr. data.table est géné­ra­le­ment plus rapide, avec une phi­lo­so­phie et une syn­taxe dif­fé­rente. Le super auteur et relec­teur Lelouar a repro­duit cette ana­lyse en uti­li­sant data.table et en mesu­rant les temps d'execution, vous pou­vez donc vous faire une idée en regar­dant ces deux note­book Jupy­ter.

Lors de la rédac­tion de cet article, je suis tom­bé sur ce billet de blog en anglais, repro­dui­sant cer­taines de ces ana­lyses avec Python et Pan­da. Un de nos rédac­teurs pré­fé­rés a aus­si extrait quelques chiffres simi­laires en uti­li­sant bash et quelques uti­li­taires.

Milles mer­ci aux relec­teurs et reléc­trice Lelouar, Aki­ra et Beba­tut ! <3 <3 <3

 



Pour continuer la lecture :


Commentaires

3 réponses à “dplyr et le génome humain”

  1. Avatar de Devailly Elisabeth
    Devailly Elisabeth

    Article bien écrit et acces­sible à tous .…. Sur­tout par l'approche rigo­lote.
    Mais bon dans une famille pas très utile de char­ger ce fichier un peu lourd en mémoire vive. !

  2. Mer­ci Guillaume !
    J'ai adap­té à mes don­nées et main­te­nant j'ai une jolie figure pour mon article \o/​

    [Remarque] j'ai uti­li­sé un fichier GTF et non GFF, du coup je l'ai char­gé comme ça :
    gen­code <- import("gencode.annotation.gtf.gz", for­mat = "GTF")

  3. Avatar de Bastien

    Mer­ci pour le post qui est très agréable à suivre !
    Pour créer la figure 6B, la fonc­tion sum­ma­rise me génère un n_​exon de type double alors que la fonc­tion mutate juste après fait des com­pa­rai­sons d'entier.
    J'ai modi­fié 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