- Le blog participatif de bioinformatique francophone depuis 2012 -

R : représenter des genomic tracks avec Gviz

Si vous ana­ly­sez des don­nées d'épigénomique telles que de l'ATAC-seq ou des ChIP-seq, vous sou­hai­tez sûre­ment pou­voir repré­sen­ter des exemples de pics sous forme de geno­mic tracks comme on en voit sou­vent dans les publi­ca­tions.

Exemple de geno­mic tracks tiré d'un article scien­ti­fique (CC-BY Davie et al. 2015)

Lorsque l'on ins­pecte ses don­nées de cove­rage (ou cou­ver­ture), on charge géné­ra­le­ment le fichier Bam ou le fichier Big­Wig dans notre navi­ga­teur de génome pré­fé­ré (pour ma part IGV, mais aus­si IGB ou UCSC). Après quelques para­mé­trages, on peut visua­li­ser et com­pa­rer plu­sieurs échan­tillons ensemble comme dans l'exemple ci-des­sous :

Screen­shot d'IGV avec des don­nées d'épigénomique dont les répli­cats ont été sup­per­po­sés ("Over­lay").

IGV, c'est super pra­tique pour ins­pec­ter ses don­nées, se bala­der le long du génome et regar­der ses signaux, mais c'est moins pra­tique si l'on veut en faire une figure pour un papier. Il y a bien la pos­si­bi­li­té de sau­ve­gar­der sa vue en PNG ou en SVG, mais ensuite il faut retou­cher le résul­tat pour pou­voir le mon­ter en figure qua­li­té publi­ca­tion, et ça quand on est bio­in­for­ma­ti­cien, on aime pas trop…

J'ai long­temps cher­ché un moyen d'automatiser la géné­ra­tion de geno­mic tracks pas trop moches avec R. Il existe bien des paquets pour faire cela, mais cha­cun a ses propres limi­ta­tions et/​ou son esthé­tique bien à lui. Je ne vais pas vous en faire une liste exhaus­tive, mais sachez que mal­gré mes recherches, je n'ai tou­jours pas trou­vé l'outil par­fait qui répond à tous mes besoins sans devoir faire des com­pro­mis.

Voi­ci une liste des paquets que j'ai regar­dé (sans for­cé­ment les tes­ter) :

Vous l'aurez com­pris, mon choix s'est por­té sur ce der­nier, Gviz, pour construire mes figures, et je vais vous par­ta­ger mes para­mètres pour arri­ver à ce résul­tat :

Le but de cet article n'étant pas de vous faire un tour exhaus­tif de Gviz, je vous conseille for­te­ment de lire, voire même de tes­ter l'excellente vignette du paquet qui est dis­po­nible depuis bio­Con­duc­tor : https://​bio​con​duc​tor​.org/​p​a​c​k​a​g​e​s​/​r​e​l​e​a​s​e​/​b​i​o​c​/​v​i​g​n​e​t​t​e​s​/​G​v​i​z​/​i​n​s​t​/​d​o​c​/​G​v​i​z​.​h​tml

Le génome de référence

Tout d'abord, on va com­men­cer par pré­pa­rer la track d'annotation du génome de réfé­rence, celle où l'on voit les gènes.

Pour cela, vous pou­vez soit uti­li­ser les génomes déjà tous prêts pré­sent dans bio­Con­duc­tor, ou bien char­ger le même génome que vous avez uti­li­sé pour votre ana­lyse (option que je pré­fère, pour une ques­tion d'exactitude de la pré­sen­ta­tion de résul­tats).

Comme je cherche à vous par­ta­ger les moments où j'ai eu des dif­fi­cul­tés afin de vous sim­pli­fier la vie (et aus­si pour mon moi du futur), je vais vous mon­trer un cas un peu com­pli­qué avec un GTF de chez Gen­code.

Pour des­si­ner la track du gène Tram1 chez la sou­ris (mm10), voi­ci le code mini­mal par défaut :

Voi­ci ce que l'on obtient :

Région du gène Tram1 à par­tir d'un fichier GTF Gen­code

Pas très pro­bant… Le gène est repré­sen­té par de gros rec­tangles, alors qu'on s'attend à voir des exons et des introns.

Les fichiers GTF de chez Gen­code contiennent des entrées pour les gènes (et pas que exon et CDS par exemple), ce qui génère de gros rec­tangles qui masquent la struc­ture exon/​intron des gènes.

Voi­ci une des lignes en ques­tion dans mon fichier GTF :

Pour s'en débar­ras­ser, l'astuce consiste à char­ger le génome comme un objet TxDb, lisible par Gviz. Cela va enle­ver les gros rec­tangles qui masquent la struc­ture des gènes :

Région du gène Tram1 à par­tir d'un fichier GTF Gen­code trans­for­mé en TxDb

C'est beau­coup mieux, mais il manque encore quelque chose, le nom des gènes ! En temps nor­mal, voi­ci à quoi devrait res­sem­bler le code pour affi­cher les noms des gènes :

Mais voi­là ce que ça donne avec un objet TxDb, on se retrouve avec les ID Ensem­bl à la place des noms.

Région du gène Tram1 avec les ID des genes à la place des noms des gènes

Là encore, il va fal­loir avoir recours à une astuce. Une fois trans­for­mé en TxDb, nous avons per­du l'information du nom des gènes car les sym­boles ne font pas par­ti du for­mat TxDb. Il va alors fal­loir aller les récu­pé­rer depuis le fichier GTF :

Nous obte­nons alors une table de cor­res­pon­dance des ID des gènes et de leurs noms res­pec­tifs :

Pour affi­cher les noms des gènes sur la track, il va fal­loir rem­pla­cer les ID par les sym­boles, comme ceci :

Région du gène Tram1 avec les noms des gènes

Enfin, ici nous voyons les trans­crits alter­na­tifs des gènes. Ce n'est pas trop gênant quand il y en a peu, mais ça le devient quand il y en a vrai­ment plein. Nous allons alors réduire l'affichage pour ne mon­trer qu'un seul modèle de gène avec tous les exons alter­na­tifs :

Région du gène Tram1 dont les trans­crits alter­na­tifs sont réduit en un modèle de gène

Vous note­rez au pas­sage que l'on a per­du les boîtes plus fines qui montrent les UTR des gènes. Comme il existe plu­sieurs UTR pos­sibles, Gviz fait le choix de ne pas les dis­tin­guer et de des­si­ner les UTR aus­si épais que des exons.

Tra­vaillons main­te­nant un peu l'esthétique de la track, parce que le diable se cache dans les détails. Pre­miè­re­ment, je n'aime pas le jaune pour les boîtes des gènes, je n'aime pas non plus qu'il y ait un contours autours des boîtes, parce qu'ils font appa­raître les exons plus gros qu'ils ne le sont réel­le­ment, je n'aime pas le fond gri­sé du titre de la track qui manque cruel­le­ment de contraste, et enfin chez la sou­ris les gènes doivent être notés en ita­lique. Voi­là à quoi res­semble ma ver­sion "pim­pée" :

Région du gène Tram1 avec mes para­mètres esthé­tiques

Le résul­tat est plus sobre, et c'est ce que l'on veut pour une figure de papier, parce que l'on veut que l'attention du lec­teur se fixe sur les don­nées que l'on va pré­sen­ter avec.

Les fichiers Bed

Les ana­lyses d'ATAC ou ChIP-seq consistent à appe­ler des pics, c'est à dire à détec­ter les signaux enri­chis d'ouverture de la chro­ma­tine ou de site de fixa­tion de pro­téines sur l'ADN. Le résul­tats de cette ana­lyse se concré­tise en un ficher Bed conte­nant les coor­don­nées des pics. Nous pou­vons repré­sen­ter ces régions géno­mique autours de nos gènes pré­fé­ré dans Gviz comme ceci :

Pics dans la région autours du gène Tram1

Les bigWig

Charger les bigWig dans R

Les fichiers big­Wig sont sou­vent lourds (autours de 100 Mb, voire beau­coup plus selon le type d'expérience), alors char­ger ces fichiers en entier dans R peut vite deve­nir limi­tant en terme de RAM. Pour ne pas explo­ser les res­sources de votre machine, il est pos­sible de n'importer que les don­nées d'un locus pré­cis grâce à rtracklayer. Voyons com­ment char­ger de façon rela­ti­ve­ment auto­ma­ti­sée des fichiers big­Wig.

Dans un pre­mier temps, on va lis­ter les fichiers à notre dis­po­si­tion qui sont ran­gés ensemble dans un seul dos­sier. Les fichiers big­Wig sont tous nom­més comme suit : "Condition_Replicat.bw" afin de faci­li­ter l'extraction des infor­ma­tions depuis leur nom (c'est une habi­tude que j'ai prise, vous faites peut-être dif­fé­rem­ment) :

Nous allons ensuite extraire les condi­tions des noms de fichiers, défi­nir une cou­leur pour cha­cune d'entre elles, et impor­ter les don­nées rela­tive à notre gène d'intérêt, ici Tram1 :

Visualiser les bigWig avec Gviz

Nous avons main­te­nant une liste d'objets GRanges avec les don­nées de cou­ver­tures autours de notre gène pré­fé­ré. Nous pou­vons main­te­nant les repré­sen­ter avec Gviz :

Et voi­là le résul­tat :

Cou­ver­ture dans la région autours du gène Tram1

Superposer les réplicats

Parce que je suis tou­jours plus exi­geante, j'aimerais super­po­ser mes répli­cats, pour n'avoir que 2 tracks à mon­trer, ce qui fera gagner beau­coup de place. Je vais repro­duire la fonc­tion "Over­lay" de IGV, avec mes répli­cats qui consiste à super­po­ser les tracks, et je vais appli­quer une légère trans­pa­rence pour voir les dif­fé­rences entre les répli­cats :

Cou­ver­ture dans la région autours du gène Tram1 avec les répli­cat sup­per­po­sés

Normaliser les valeurs des tracks

Enfin, pour aller encore plus dans les détails, vous note­rez que la hau­teur des tracks est réglée sur les valeurs maxi­males des don­nées de cou­ver­tures, mais que ces valeurs ne sont pas les mêmes d'une condi­tion à l'autre. Pour faci­li­ter la com­pa­rai­son des deux condi­tions, on va fixer la hau­teur des tracks à la valeur maxi­male toutes condi­tions confon­dues, pour pou­voir appré­cier visuel­le­ment les dif­fé­rences entre les deux :

On va main­te­nant pou­voir pas­ser cette valeur à notre fonc­tion bigwig_tracks_overlay :

Cou­ver­ture dans la région autours du gène Tram1 avec les répli­cat sup­per­po­sés, cette fois avec les valeurs en Y nor­ma­li­sées entre les deux condi­tions

Mettre en évidence des régions ou des pics

On y est presque, réglons les der­niers détails. Pour mettre en valeur les régions où il y a des pics, on peut des­si­ner des rec­tangles en arrière plan qui tra­versent les tracks afin de mieux voir leur déli­mi­ta­tions. Pour cela, on uti­lise la fonc­tion HighlightTrack qui prend une liste de posi­tions start et end pour déli­mi­ter les régions à mettre en évi­dence, et qui englobe les tracks sous les­quelles on veut faire pas­ser les rec­tangles :

Et voi­là ce que ça donne :

Pics dans la région autours du gène Tram1 mis en évi­dence par des rec­tangles gri­sés

On voit main­te­nant bien où sont les pics par rap­port aux signaux des big­Wigs et aus­si par rap­port au gène. Notez tou­te­fois qu'il y a une limi­ta­tion à l'usage des high­lights si vous expor­tez vos figures en PNG avec des uni­tés autres que des pixels. Voi­ci ce que ça donne quand on l'exporte en don­nant des tailles en cm :

Export de la figure en PNG avec des uni­tées en cm. Les high­lights ne sont plus ali­gnées cor­rec­te­ment…

Les rec­tangles gris vers la gauche ne sont plus ali­gnés cor­rec­te­ment avec les pics ! Je ne com­pred pas d'où vient ce pro­blème, alors il faut trou­ver un autre moyen d'exporter ses figures. Je n'ai pas tes­té toutes les options, mais l’export en PDF fonc­tionne bien.

Limitations et conclusion

Comme vous l'avez vu, pour géné­rer des figures de geno­mic tracks de bonne qua­li­té visuelle, il faut en suer un petit peu, et c'est pour ça que j'ai vou­lu par­ta­ger mon code avec vous, afin de vous épar­gner du temps et des tâton­ne­ments. Bien que très modu­lable, Gviz réclame pas mal de petits bidouillages à gauche à droite pour obte­nir le résul­tat que l'on sou­haite, sur­tout quand on est un poil maniaque des gra­phiques comme moi… Je n'ai pas pré­sen­té toutes mes astuces ici mais si vous vou­lez par exemple avoir les titres des tracks à l'horizontale, sachez que l'option rotation.title=0 existe bien, mais n'est que par­tiel­le­ment implé­men­tée et il fau­dra bidouiller avec la lar­geur de la colonne de titre et rajou­ter des espaces dans les noms des titres pour que ça passe…

Néan­moins, je trouve quand même que l'effort en vaut la peine, le résul­tat est visuel­le­ment satis­fai­sant, et je ne pense pas que les autres paquets R puissent faire mieux à l'heure de l'écriture de cet article.

Mer­ci à Léo­pold, Ista, et Aze­rin pour la relec­ture.

Vous avez aimé ? Dites-le nous !

Moyenne : 4.8 /​ 5. Nb de votes : 4

Pas encore de vote pour cet article.

Partagez cet article :




Commentaires

2 réponses à “R : représenter des genomic tracks avec Gviz”

    1. Des para­mètres, plein ! Mer­ci, je vais y jeter un œil.

Laisser un commentaire