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

En génomique, et sans doute dans tout un tas d'autres domaines omiques ou big data, nous essayons souvent de tracer des grosses matrices sous forme d'heatmap. 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 utilisez. 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 verticaux sur un écran HD (à moins bien sûr que vous lisiez ceci dans un futur lointain d'hyper haute définition).

Le problème lorsqu'on affiche des matrices qui ont plus de lignes que de pixel à l'écran, c'est justement que chaque pixel va devoir représenter plusieurs cellules de la matrice, et que le comportement par défaut de R sur ce point-là n'est pas forcément optimal.

Un exemple de données numériques

Commençons par générer un faux jeu de données, imitant ce qu'on peut obtenir en épigénomique. J'essaye de produire un signal centré au milieu des rangs, de plus en plus fort, tout en gardant une part d'aléatoire. Je laisse le code pour que vous puissiez jouer chez vous à reproduire les figures de cet article, mais vous n'avez pas besoin de comprendre cette section pour comprendre la suite.

Il existe de nombreuses fonctions en R pour afficher cette matrice, par exemple heatmap(), gplots::heatmap.2(), ggplot2::geom_raster(), ou ComplexHeatmap::Heatmap(). La plupart de ces fonctions font appel à la fonction image() de plus bas niveau, qui est celle qui trace la matrice colorée. C'est cette fonction que nous allons utiliser dans ce billet:

Figure 1: Une bien belle heatmap ?

On pourrait penser que tout va bien ici. On arrive bien à voir un signal, et on en tire la conclusion qu'il est au centre, plus fort en haut qu'en bas. Le souci c'est qu'avec 20 000 lignes dans notre matrice on devrait avoir une image beaucoup moins bruitée. Comme il n'y a que quelques centaines de pixels de hauteur dans le png (par défaut), R doit décider d'une manière ou d'une autre comment résumer l'information de plusieurs cellules en un seul pixel. Il semble que R choisisse plus ou moins au hasard une seule cellule à afficher par pixel (probablement la première ou la dernière dans la pile). Il y a donc un sous-échantillonnage important.

On pourrait imaginer générer un png de plus de 20 000 pixels de haut pour compenser, mais ça fait des fichiers lourds à manipuler, et il faut penser à agrandir d'autant la taille du texte et l’épaisseur des traits pour un résultat potable.

Autre idée, certaines devices graphiques (par exemple pdf(), mais pas png()) permettent jouer avec le paramètre useRaster = TRUE de la fonction image(), ce qui peut aider dans quelques situations. La rasterisation, d'après wikipedia, "est un procédé qui consiste à convertir une image vectorielle en une image matricielle". L’algorithme de rasterisation va donc essayer de convertir plusieurs lignes de données en un seul pixel.

Le fichier pdf généré est disponible ici. Mais les différents lecteurs pdf n'affichent pas le même rendu des plots en questions :

Figure 2 : useRaster = TRUE, une solution loin d'être idéale

Acrobat, Edge et Okular donnent le rendu attendu: une représentation bien plus fine des données originales lorsque la rasterisation est activée. Evince et SumatraPDF inversent les rendus, et voilent la version "non rasterisée" ! Le lecteur de pdf de Firefox abandonne carrément (en tout cas sous Windows 10, sous GNU/Linux il affiche le même résultat qu'Acrobat, Edge et Okular). Si votre lecteur de pdf préféré n'est pas parmi ceux que j'ai testé, je serai curieux d'avoir le résultat que vous obtenez en commentaire.

Pour info, alors que le pdf fait 5 Mo, le même code exportant du svg génère un fichier de 200 Mo ! Je n'ai lâchement pas eu le courage de l'ouvrir pour voir le rendu obtenu...

Au final, la rasterisation essaie de résumer les informations contenues derrières chaque pixel en en faisant une moyenne. Mais c'est un processus qu'on peut essayer de faire nous même, ce qui a deux avantages : on s'affranchit des différences de rendus entre lecteurs de pdf, et ça marchera même sur les devices non vectoriels, du genre png, ce qui évite de générer des images trop lourdes.

L'idée est donc de redimensionner la matrice avant le plot, en la rendant plus petite et en appliquant une fonction qui "résumera" les cellules correspondantes à chaque pixel (par exemple la fonction mean()). Je vous propose cette petite fonction (aussi disponible sur canSnippet) :

Comparons le rendu Avant / Après :

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

Paradoxalement, on "discerne" bien mieux les détails des 20 000 lignes de la matrice en réduisant la taille de la matrice nous même, plutôt qu'en laissant R (mal) afficher les 20 000 lignes.

Matrices creuses (sparse)

Dans certaines situations, faire la moyenne des cellules par pixel n'est pas la manière la plus maline de résumer les données. Par exemple, dans le cas de matrices creuses, on ne souhaite pas moyenner nos quelques valeurs isolées par tous les zéros les entourant. Dans ces cas-là, prendre la valeur maximale correspondra mieux à ce qu'on cherche à montrer.

J'ai rencontré ce cas dans une étude d'eQTL (analyse QTL en utilisant les niveaux d'expressions des gènes comme phénotypes). L'idée est d'identifier des Single Nucleotide Polymorphisms (SNP, des variants / mutations) qui sont associés à des changements d'expression des gènes. Pour cela, on fait un test statistique 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 environ 45 000 SNP, résultant en une matrice de 20 000 x 45 000 p-valeurs. La plupart des p-valeurs sont non significatives, et seule une minorité était très petite (ou très grande après transformation en -log10(p-valeur)). Or, ce qu'on souhaite c'est afficher les p-valeurs des SNP principaux (lead SNP). On va donc plutôt prendre le maximum des -log10(p-valeur) plutôt que leur moyenne :

Figure 4: À gauche : On résume la grosse matrice en calculant la moyenne des p-valeurs par pixel. À droite : On prend la plus petite p-valeur (la plus grande après transformation -log10) pour mieux voir la significativité des lead SNP. Les p-valeurs réelles sont bien mieux représentés.

Données catégorielles

Dans le cas de données catégorielles, on ne peut pas vraiment prendre une moyenne des valeurs. Il faut plutôt faire les moyennes des couleurs associées à chaque catégorie (Ici, je le fais dans l'espace colorimétrique RGB, mais ça fonctionne peut-être encore mieux si la moyenne est faite en espace HCL ?). Pour afficher une matrice de couleurs, il faut utiliser rasterImage() au lieu d'image().

Figure 5: Réduire la taille d'une matrice catégorielle avant le plot, en faisant la moyenne des couleurs par pixel, permet de représenter plus fidèlement les données.

ggplot2

D’après mes tests, ggplot2 est aussi affecté par ce souci d'overplotting, que ce soit geom_tile() ou geom_raster() (qui est une version optimisée de geom_tile() quand les cases sont régulières).

Figure 6 : ggplot2 victime de l'overplotting. Notez de subtiles différences entre geom_tile() et geom_raster().

ComplexHeatmap

Le package Bioconductor ComplexHeatmap est vraiment top pour générer des Heatmaps un peu complexe, avec des annotations dans tous les sens.

Cela dit, mes quelques tests suggèrent qu'il souffre du même problème d'overplotting que les autres fonctions. Il réalise un sous-échantillonnage des cellules à afficher, au lieu de moyenner les données par pixel :

Figure 7 : ComplexHeatmap nous déçoit.

La fonction Heatmap() a bien des paramètres qui permettent une rasterization dans le cas de grosses matrices, mais ils semblent plus utiles pour réduire le poids des fichiers vectoriels que pour résoudre le problème de sous-échantillonnage :

Figure 8: ComplexHeatmap nous déçoit même avec use_raster = TRUE

Conclusion

En R, réduisez vos grosses matrices avant de les afficher, vous verrez mieux les petits détails. Sinon vous obtiendrez des heatmaps un peu approximatives.

Le principal souci de cette solution, c'est qu'il faut faire les décorations (axes, barres de couleurs sur les côtés, dendrogrammes, etc.) à la main, ce qui est un peu laborieux.

Merci a mes talentueux relecteurs Mathurin, Gwenaëlle, et lhtd. Une version de cet article traduite en anglais sera publiée prochainement sur le blog de l'auteur.

  • À propos de
  • Après une thèse en cancérologie à Lyon et un postdoc en bioinformatique à Édimbourg, je suis chercheur à l'INRA Toulouse depuis fin 2017. Régulation transcriptionnelle et épigénétique. Twitter: @G_Devailly

Catégorie: Astuce | Tags: , , , ,

Laisser un commentaire