Créer des Heatmaps à partir de grosses matrices en R

En géno­mique, et sans doute dans tout un tas d'autres domaines omiques ou big data, nous essayons sou­vent de tra­cer des grosses matrices sous forme d'heat­map. Par grosse matrice, j'entends une matrice dont le nombre de lignes et/​ou de colonnes est plus grand que le nombre de pixels sur l'écran que vous uti­li­sez. Par exemples, si vous avez une matrice de 50 colonnes et de 20 000 lignes (cas assez fré­quent quand il y a une ligne par gène), il y a de forte chances que cette matrice aura plus de lignes qu'il n'y a de pixels sur votre écran — 1080 pixels ver­ti­caux sur un écran HD (à moins bien sûr que vous lisiez ceci dans un futur loin­tain d'hyper haute défi­ni­tion).

Le pro­blème lorsqu'on affiche des matrices qui ont plus de lignes que de pixel à l'écran, c'est jus­te­ment que chaque pixel va devoir repré­sen­ter plu­sieurs cel­lules de la matrice, et que le com­por­te­ment par défaut de R sur ce point-là n'est pas for­cé­ment opti­mal.

Un exemple de données numériques

Com­men­çons par géné­rer un faux jeu de don­nées, imi­tant ce qu'on peut obte­nir en épi­gé­no­mique. J'essaye de pro­duire un signal cen­tré au milieu des rangs, de plus en plus fort, tout en gar­dant une part d'aléatoire. Je laisse le code pour que vous puis­siez jouer chez vous à repro­duire les figures de cet article, mais vous n'avez pas besoin de com­prendre cette sec­tion pour com­prendre la suite.

Il existe de nom­breuses fonc­tions en R pour affi­cher cette matrice, par exemple heatmap(), gplots::heatmap.2(), ggplot2::geom_raster(), ou ComplexHeatmap::Heatmap(). La plu­part de ces fonc­tions font appel à la fonc­tion image() de plus bas niveau, qui est celle qui trace la matrice colo­rée. C'est cette fonc­tion que nous allons uti­li­ser dans ce billet :

Figure 1 : Une bien belle heat­map ?

On pour­rait pen­ser que tout va bien ici. On arrive bien à voir un signal, et on en tire la conclu­sion qu'il est au centre, plus fort en haut qu'en bas. Le sou­ci c'est qu'avec 20 000 lignes dans notre matrice on devrait avoir une image beau­coup moins brui­tée. Comme il n'y a que quelques cen­taines de pixels de hau­teur dans le png (par défaut), R doit déci­der d'une manière ou d'une autre com­ment résu­mer l'information de plu­sieurs cel­lules en un seul pixel. Il semble que R choi­sisse plus ou moins au hasard une seule cel­lule à affi­cher par pixel (pro­ba­ble­ment la pre­mière ou la der­nière dans la pile). Il y a donc un sous-échan­tillon­nage impor­tant.

On pour­rait ima­gi­ner géné­rer un png de plus de 20 000 pixels de haut pour com­pen­ser, mais ça fait des fichiers lourds à mani­pu­ler, et il faut pen­ser à agran­dir d'autant la taille du texte et l’épaisseur des traits pour un résul­tat potable.

Autre idée, cer­taines devices gra­phiques (par exemple pdf(), mais pas png()) per­mettent jouer avec le para­mètre useRaster = TRUE de la fonc­tion image(), ce qui peut aider dans quelques situa­tions. La ras­te­ri­sa­tion, d'après wiki­pe­dia, "est un pro­cé­dé qui consiste à conver­tir une image vec­to­rielle en une image matri­cielle". L’algorithme de ras­te­ri­sa­tion va donc essayer de conver­tir plu­sieurs lignes de don­nées en un seul pixel.

Le fichier pdf géné­ré est dis­po­nible ici. Mais les dif­fé­rents lec­teurs pdf n'affichent pas le même ren­du des plots en ques­tions :

Figure 2 : use­Ras­ter = TRUE, une solu­tion loin d'être idéale

Acro­bat, Edge et Oku­lar donnent le ren­du atten­du : une repré­sen­ta­tion bien plus fine des don­nées ori­gi­nales lorsque la ras­te­ri­sa­tion est acti­vée. Evince et Suma­traPDF inversent les ren­dus, et voilent la ver­sion "non ras­te­ri­sée" ! Le lec­teur de pdf de Fire­fox aban­donne car­ré­ment (en tout cas sous Win­dows 10, sous GNU/​Linux il affiche le même résul­tat qu'Acrobat, Edge et Oku­lar). Si votre lec­teur de pdf pré­fé­ré n'est pas par­mi ceux que j'ai tes­té, je serai curieux d'avoir le résul­tat que vous obte­nez en com­men­taire.

Pour info, alors que le pdf fait 5 Mo, le même code expor­tant du svg génère un fichier de 200 Mo ! Je n'ai lâche­ment pas eu le cou­rage de l'ouvrir pour voir le ren­du obte­nu…

Au final, la ras­te­ri­sa­tion essaie de résu­mer les infor­ma­tions conte­nues der­rières chaque pixel en en fai­sant une moyenne. Mais c'est un pro­ces­sus qu'on peut essayer de faire nous même, ce qui a deux avan­tages : on s'affranchit des dif­fé­rences de ren­dus entre lec­teurs de pdf, et ça mar­che­ra même sur les devices non vec­to­riels, du genre png, ce qui évite de géné­rer des images trop lourdes.

L'idée est donc de redi­men­sion­ner la matrice avant le plot, en la ren­dant plus petite et en appli­quant une fonc­tion qui "résu­me­ra" les cel­lules cor­res­pon­dantes à chaque pixel (par exemple la fonc­tion mean()). Je vous pro­pose cette petite fonc­tion (aus­si dis­po­nible sur canS­nip­pet) :

Com­pa­rons le ren­du Avant /​ Après :

Figure 3 : Un ren­du bien plus fin lorsqu'on réduit la taille de la matrice nous même.

Para­doxa­le­ment, on "dis­cerne" bien mieux les détails des 20 000 lignes de la matrice en rédui­sant la taille de la matrice nous même, plu­tôt qu'en lais­sant R (mal) affi­cher les 20 000 lignes.

Matrices creuses (sparse)

Dans cer­taines situa­tions, faire la moyenne des cel­lules par pixel n'est pas la manière la plus maline de résu­mer les don­nées. Par exemple, dans le cas de matrices creuses, on ne sou­haite pas moyen­ner nos quelques valeurs iso­lées par tous les zéros les entou­rant. Dans ces cas-là, prendre la valeur maxi­male cor­res­pon­dra mieux à ce qu'on cherche à mon­trer.

J'ai ren­con­tré ce cas dans une étude d'eQTL (ana­lyse QTL en uti­li­sant les niveaux d'expressions des gènes comme phé­no­types). L'idée est d'identifier des Single Nucleo­tide Poly­mor­phisms (SNP, des variants /​ muta­tions) qui sont asso­ciés à des chan­ge­ments d'expression des gènes. Pour cela, on fait un test sta­tis­tique d'association entre chaque SNP et chaque niveau d'expression de gènes, ce qui nous donne autant de p‑valeurs.

Nous avions l'expression d'environ 20 000 gènes, et envi­ron 45 000 SNP, résul­tant en une matrice de 20 000 x 45 000 p‑valeurs. La plu­part des p‑valeurs sont non signi­fi­ca­tives, et seule une mino­ri­té était très petite (ou très grande après trans­for­ma­tion en -log10(p‑valeur)). Or, ce qu'on sou­haite c'est affi­cher les p‑valeurs des SNP prin­ci­paux (lead SNP). On va donc plu­tôt prendre le maxi­mum des -log10(p‑valeur) plu­tôt que leur moyenne :

Figure 4 : À gauche : On résume la grosse matrice en cal­cu­lant la moyenne des p‑valeurs par pixel. À droite : On prend la plus petite p‑valeur (la plus grande après trans­for­ma­tion ‑log10) pour mieux voir la signi­fi­ca­ti­vi­té des lead SNP. Les p‑valeurs réelles sont bien mieux repré­sen­tés.

Données catégorielles

Dans le cas de don­nées caté­go­rielles, on ne peut pas vrai­ment prendre une moyenne des valeurs. Il faut plu­tôt faire les moyennes des cou­leurs asso­ciées à chaque caté­go­rie (Ici, je le fais dans l'espace colo­ri­mé­trique RGB, mais ça fonc­tionne peut-être encore mieux si la moyenne est faite en espace HCL ?). Pour affi­cher une matrice de cou­leurs, il faut uti­li­ser rasterImage() au lieu d'image().

Figure 5 : Réduire la taille d'une matrice caté­go­rielle avant le plot, en fai­sant la moyenne des cou­leurs par pixel, per­met de repré­sen­ter plus fidè­le­ment les don­nées.

ggplot2

D’après mes tests, ggplot2 est aus­si affec­té par ce sou­ci d'overplotting, que ce soit geom_tile() ou geom_raster() (qui est une ver­sion opti­mi­sée de geom_​tile() quand les cases sont régu­lières).

Figure 6 : ggplot2 vic­time de l'overplotting. Notez de sub­tiles dif­fé­rences entre geom_​tile() et geom_​raster().

ComplexHeatmap

Le package Bio­con­duc­tor Com­plex­Heat­map est vrai­ment top pour géné­rer des Heat­maps un peu com­plexe, avec des anno­ta­tions dans tous les sens.

Cela dit, mes quelques tests sug­gèrent qu'il souffre du même pro­blème d'overplotting que les autres fonc­tions. Il réa­lise un sous-échan­tillon­nage des cel­lules à affi­cher, au lieu de moyen­ner les don­nées par pixel :

Figure 7 : Com­plex­Heat­map nous déçoit.

La fonc­tion Heatmap() a bien des para­mètres qui per­mettent une ras­te­ri­za­tion dans le cas de grosses matrices, mais ils semblent plus utiles pour réduire le poids des fichiers vec­to­riels que pour résoudre le pro­blème de sous-échan­tillon­nage :

Figure 8 : Com­plex­Heat­map nous déçoit même avec use_​raster = TRUE

Conclusion

En R, rédui­sez vos grosses matrices avant de les affi­cher, vous ver­rez mieux les petits détails. Sinon vous obtien­drez des heat­maps un peu approxi­ma­tives.

Le prin­ci­pal sou­ci de cette solu­tion, c'est qu'il faut faire les déco­ra­tions (axes, barres de cou­leurs sur les côtés, den­dro­grammes, etc.) à la main, ce qui est un peu labo­rieux.

Mer­ci a mes talen­tueux relec­teurs Mathu­rin, Gwe­naëlle, et lhtd. Une ver­sion de cet article tra­duite en anglais sera publiée pro­chai­ne­ment sur le blog de l'auteur.



Pour continuer la lecture :


Commentaires

Laisser un commentaire