L'annotation de régions génomiques et les analyses d’enrichissement

Les anno­ta­tions sont essen­tielles lors d'analyses fonc­tion­nelles à large échelle sur le génome. 

Lorsque l’on pra­tique des ana­lyses en géno­mique, basées sur des tech­niques comme le RNA-seq ou le ChIP-seq, on se retrouve avec res­pec­ti­ve­ment une liste de trans­crits ou de pics (régions géno­miques). Dans le cas des ana­lyses de ChIP-seq, on sou­haite carac­té­ri­ser les gènes cibles du fac­teur de trans­crip­tion étu­dié sur tout le génome (genome-wide), pour com­prendre la fonc­tion bio­lo­gique de ce fac­teur. Dans le cas du RNA-seq, on obtient une liste de trans­crits dif­fé­ren­tiel­le­ment expri­més dont on sou­haite carac­té­ri­ser la fonc­tion.

Dans cet article nous allons uti­li­ser plu­sieurs librai­ries R pour auto­ma­ti­ser cette ana­lyse, à par­tir d’un fichier .bed (voir article Bed & Fas­ta), ou d’une liste de trans­cripts Ensem­bl.

Pour com­men­cer, nous allons uti­li­ser des sites de ChIP-seq du récep­teur aux glu­co­cor­ti­coïdes GR (Grønt­ved L, John S, Baek S, Liu Y et al. , EMBO J 2013, GSE46047). GR est un récep­teur nucléaire impor­tant impli­qué dans la glu­co­néo­ge­nèse, la gly­co­lyse, le méta­bo­lisme des acides gras et la réponse immu­ni­taire et inflam­ma­toire. Nous allons anno­ter les pics de GR avec les gènes les plus proches sur le génome de la sou­ris (mm9) à l’aide de la librai­rie ChIP­pea­kan­no.

Ensuite, nous allons voir com­ment conver­tir des iden­ti­fiants Ensem­bl en sym­boles de gènes à l’aide de la librai­rie Bio­maRt. Enfin, nous allons faire une ana­lyse d’enrichissement d’annotations à l’aide de la librai­rie RDA­VID­Web­Ser­vice. DAVID est un très bon site d’analyse d’annotations qui per­met de tra­vailler avec dif­fé­rentes sources comme les onto­lo­gies de gènes (GO terms), que nous avions intro­duites dans un article pré­cé­dent et les voies de signa­li­sa­tions entre autres. La base de don­nées de l’outil DAVID per­met de faire des requêtes sur 82 sources, dont notam­ment REACTOME, KEGG et PANTHER, qui sont main­te­nues par des bio­cu­ra­teurs.

DAVID vous per­met d’utiliser son inter­face web via son site ou des ser­vices web (accès pro­gram­ma­tique). D’autres appli­ca­tions web per­mettent de tra­vailler direc­te­ment avec des fichiers .bed, comme l’excellent outil GREAT du labo­ra­toire Beje­ra­no de Stan­ford.

Un peu de code R : mise en place de l'environnement de travail 

Vous pou­vez télé­char­ger les don­nées et le script ici.

Dans un pre­mier temps, il faut ins­tal­ler les librai­ries néces­saires pour l’analyse avec les com­mandes R sui­vantes :

Dans un deuxième temps, nous devons impor­ter les don­nées des pics de GR au for­mat .bed

Nous pou­vons uti­li­ser bio­maRt pour obte­nir les sites d’initiation de la trans­crip­tion (TSS). Nous uti­li­sons les don­nées de Ensem­bl NCBIM37 pour l’assemblage du génome “mm9”.

Grâce aux com­mandes ci-des­sus, nous pou­vons obte­nir une table ayant la struc­ture sui­vante :

Le champ « seq­names » repré­sente le chro­mo­some, les champs « start » et « end » repré­sentent les coor­don­nées géno­miques des pics de GR. Ces pics sont anno­tés avec le TSS du gène Ensem­bl le plus proche en indi­quant la dis­tance et le che­vau­che­ment (over­lap).

Com­ment faire des requêtes sur Bio­mart pour conver­tir des iden­ti­fiants Ensem­bl en sym­boles de gènes ?

En regar­dant à nou­veau la table, on voit que les sym­boles des gènes sont à pré­sent affi­chés sur la der­nière colonne. Pour tra­vailler avec les iden­ti­fiants des trans­crits et non des gènes, il faut sim­ple­ment uti­li­ser "ensembl_​transcript_​id".

Comment faire une analyse d’enrichissement grâce à DAVID ? Et qu'est-ce qu’une analyse d’enrichissement ?

Les ana­lyses d’enrichissement d’annotations ont pour but de cal­cu­ler une pro­ba­bi­li­té : sachant qu’un groupe de gènes est anno­té avec un terme spé­ci­fique de KEGG_​PATHWAY, quelle est la pro­ba­bi­li­té que la tota­li­té ou une frac­tion de ces gènes soit dans le groupe des cibles de GR (dans ce cas). David uti­lise le test hyper géo­mé­trique pour cal­cu­ler cette valeur p.

On peut faire une ana­lyse d'enrichissement à par­tir de n'importe quelle liste de trans­crits pro­ve­nant soit de RNA-seq soit de pics de ChIP anno­tés avec les gènes les plus proches (comme ci-des­sus).

Vous devez d’abord vous ins­crire sur le site de DAVID pour uti­li­ser le ser­vice web.
Puis vous pou­vez vous connec­ter :

Vous pou­vez accé­der à la liste des anno­ta­tions dis­po­nibles avec la com­mande :

Lorsque vous avez choi­si les anno­ta­tions utiles (une ou plu­sieurs) pour votre ana­lyse (dans ce cas KEGG_​PATHWAY), vous pou­vez uti­li­ser le code sui­vant :

Vous pou­vez éga­le­ment ajou­ter une liste de gènes en arrière-plan (avec listType="Background") pour amé­lio­rer votre ana­lyse, en ne consi­dé­rant que les gènes expri­més dans un tis­su pré­cis (par exemple).

Le ser­vice web de DAVID va vous envoyer une table conte­nant dans chaque ligne les anno­ta­tions triées par la P-valeur, avec dif­fé­rentes métriques (FDR, Bon­fer­ro­ni) et les iden­ti­fiants Ensem­bl.

On peut enfin visua­li­ser l’analyse à l’aide d’un bar plot.

On constate sur ce gra­phique que, par­mi les anno­ta­tions enri­chies, on retrouve une majo­ri­té de fonc­tions, décrites dans la lit­té­ra­ture, des récep­teurs aux glu­co­cor­ti­coïdes (GR). Nous voyons par exemple que les adi­po­cy­to­kines (inflam­ma­tion) et les récep­teurs des cel­lules T (réponse immu­ni­taire) sont par­mi les cibles les plus impor­tantes de GR.

En conclu­sion, nous avons un joli script qui per­met d'automatiser des requêtes sur DAVID à par­tir d’un fichier de régions géno­miques au for­mat .bed. D’autres packages R offrent la pos­si­bi­li­té de faire des ana­lyses simi­laires comme top­GO ou encore GAGE. N’hésitez pas à les tes­ter. Don­nez-nous votre avis sur les outils que vous connais­sez !

Mer­ci aux relec­teurs : Yoann M, ook4mi, NiGo­Pol, mura­veill, Zazo0o et Estel