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

1920px-Gas_centrifuge_cascade
Non il ne s'agit pas d'enrichissement d'uranium ! (U.S. Depart­ment of Ener­gy, Domaine Public)

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 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­viceDAVID 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 REACTOMEKEGG 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 tra­vail 

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 :

source("http://bioconductor.org/biocLite.R")
biocLite("ChIPpeakAnno")
biocLite("rtracklayer")
biocLite("biomaRt")
biocLite("RDAVIDWebService")
library("rtracklayer")
library("biomaRt")
library("ChIPpeakAnno")
library("RDAVIDWebService")
Dans un deuxième temps, nous devons importer les données des pics de GR au format .bed
GR_bed=import.bed("GSE46047_GR_peaks_mouse_liver_mm9.txt", genome = "mm9")

Nous pouvons utiliser biomaRt pour obtenir les sites d’initiation de la transcription (TSS). Nous utilisons les données de Ensembl NCBIM37 pour l’assemblage du génome “mm9”.

ensembl = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
TSS.mouse.NCBIM37 = getAnnotation(mart=ensembl, featureType="TSS")
Comment annoter des pics de GR sur le génome avec les gènes les plus proches ?
annotatedPeak = annotatePeakInBatch(GR_bed, AnnotationData = TSS.mouse.NCBIM37)
Annotatedbed = as.data.frame(annotatedPeak)
Grâce aux commandes ci-dessus nous pouvons obtenir une table ayant la structure suivante :
head(Annotatedbed)
seqnames start end width strand peak feature start_position end_position insideFeature distancetoFeature
00002 ENSMUSG00000103922 1 4759134 4759283 150 + 00002 ENSMUSG00000103922 4771131 4772199 upstream -11997
00006 ENSMUSG00000051285 1 7079940 7080089 150 + 00006 ENSMUSG00000051285 7088920 7173628 upstream -8980
00009 ENSMUSG00000104504 1 8862999 8863148 150 + 00009 ENSMUSG00000104504 8856763 8857102 downstream 6236
00010 ENSMUSG00000025911 1 9554202 9554351 150 + 00010 ENSMUSG00000025911 9547948 9580673 inside 6254
00011 ENSMUSG00000081441 1 9574861 9575010 150 + 00011 ENSMUSG00000081441 9574548 9575649 inside 313
00012 ENSMUSG00000067879 1 9611749 9611898 150 + 00012 ENSMUSG00000067879 9601199 9627143 inside 10550
shortestDistance fromOverlappingOrNearest
00002 ENSMUSG00000103922 11848 NearestStart
00006 ENSMUSG00000051285 8831 NearestStart
00009 ENSMUSG00000104504 5897 NearestStart
00010 ENSMUSG00000025911 6254 NearestStart
00011 ENSMUSG00000081441 313 NearestStart
00012 ENSMUSG00000067879 10550 NearestStart

Le champ « seqnames » représente le chromosome , les champs « start » et « end » représentent les coordonnées génomiques des pics de GR. Ces pics sont annotés avec le TSS du gène Ensembl le plus proche en indiquant la distance et le chevauchement (overlap).

Comment faire des requêtes sur Biomart pour convertir des identifiants Ensembl en symboles de gènes ?

maptable = getBM(attributes = c("mgi_symbol", "ensembl_gene_id"), filters = "ensembl_gene_id", values = Annotatedbed$feature, mart = ensembl)
ii = match(Annotatedbed$feature,maptable[,2])
Annotatedbed$feature_symbol = maptable[ii,1]
En regardant à nouveau la table, on voit que les symboles des gènes sont à présent affichés sur la dernière colonne. Pour travailler avec les identifiants des transcrits et non des gènes, il faut simplement utiliser "ensembl_transcript_id".
head(Annotatedbed)
Comment faire une analyse d’enrichissement grâce à DAVID ? Et qu'est-ce qu’une analyse d’enrichissement ?

Les analyses d’enrichissement d’annotations ont pour but de calculer une probabilité : sachant qu’un groupe de gènes est annoté avec un terme spécifique de KEGG_PATHWAY, quelle est la probabilité que la totalité ou une fraction de ces gènes soit dans le groupe des cibles de GR (dans ce cas). David utilise le test hyper géométrique pour calculer cette valeur p.

On peut faire une analyse d'enrichissement à partir de n'importe quelle liste de transcrits provenant soit de RNA-seq soit de pics de ChIP annotés avec les gènes les plus proches (comme ci-dessus).

Vous devez d’abord vous inscrire sur le site de DAVID pour utiliser le service web.
Puis vous pouvez vous connecter :

david = DAVIDWebService$new(email="votre_mail@institution.fr")

Vous pouvez accéder à la liste des annotations disponibles avec la commande:

getAllAnnotationCategoryNames (david)
Lorsque vous avez choisi les annotations utiles (une ou plusieurs) pour votre analyse (dans ce cas KEGG_PATHWAY), vous pouvez utiliser le code suivant :
setAnnotationCategories(david, c("KEGG_PATHWAY"))
result1 = addList(david, Annotatedbed$feature,idType="ENSEMBL_GENE_ID", listName="GR_targets", listType="Gene")
t_p1=getFunctionalAnnotationChart(david)

Vous pouvez également ajouter une liste de gènes en arrière plan (avec listType="Background") pour améliorer votre analyse, en ne considérant que les gènes exprimés dans un tissus précis (par exemple).

Le service web de DAVID va vous envoyer une table contenant dans chaque ligne les annotations triées par la P-valeur, avec différentes métriques (FDR, Bonferroni) et les identifiants Ensembl.

head(t_p1)
On peut enfin visualiser l’analyse à l’aide d’un bar plot.
par(mar=c(5, 20, 4, 2) + 0.1)
barplot(rev(-log10(t_p1$PValue)), names.arg=rev(t_p1$Term),col="darkgreen",horiz=T,las=2,cex.lab=1,cex.names=0.8,xlab="-log10 P valeur")
title("analyse d'enrichissement d'annotations")
Screenshot 2015-04-07 20.51.41

On constate sur ce graphique que, parmi les annotations enrichies, on retrouve une majorité de fonctions, décrites dans la littérature, des récepteurs aux glucocorticoïdes (GR). Nous voyons par exemple que les adipocytokines (inflammation) et les récepteurs des cellules T (réponse immunitaire) sont parmi les cibles les plus importantes de GR.

En conclusion, nous avons un joli script qui permet d' automatiser des requêtes sur DAVID à partir d’un fichier de régions génomiques au format .bed. D’autres packages R offrent la possibilité de faire des analyses similaires comme topGO ou encore GAGE. N’hésitez pas à les tester. Donnez-nous votre avis sur les outils que vous connaissez !

Merci aux relecteurs: Yoann M, ook4mi, NiGoPol, muraveill, Zazo0o et Estel



Pour continuer la lecture :


Commentaires

Laisser un commentaire