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.

library(dplyr)
library(purrr)
Ncol <- 50
genmat <- map(
    1:20000,
    function(i) {
        runif(Ncol) + c(sort(abs(rnorm(50)))[1:25], rev(sort(abs(rnorm(50)))[1:25]))  * i/5000
    }
) %>% do.call(rbind, .)
genmat[1:5, 1:5]
##            [,1]      [,2]       [,3]      [,4]      [,5]
## [1,] 0.24750229 0.8144309 0.31405005 0.2787540 0.8435071
## [2,] 0.44266149 0.1147394 0.28464511 0.6437944 0.7597911
## [3,] 0.11495737 0.6750608 0.04393633 0.5712240 0.2088942
## [4,] 0.16660166 0.5508895 0.75274403 0.7340737 0.9325773
## [5,] 0.07285492 0.4573314 0.09437322 0.1534962 0.4939674
dim(genmat)
## [1] 20000    50

Il existe de nom­breuses fonc­tions en R pour affi­cher cette matrice, par exemple

,

,

, ou

. La plu­part de ces fonc­tions font appel à la fonc­tion

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 :

oldpar <- par(mar = rep(0.2, 4)) # reducing plot margins
image(
    t(genmat), # image() has some weird opinions about how your matrix will be plotted
    axes = FALSE,
    col = colorRampPalette(c("white", "darkorange", "black"))(30), # our colour palette
    breaks = c(seq(0, 3, length.out = 30), 100) # colour-to-value mapping
)
box() # adding a box arround the heatmap

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

, mais pas

) per­mettent jouer avec le para­mètre

de la fonc­tion

, 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.

pdf("big_hm_1.pdf")
layout(matrix(c(1, 2), nrow = 1)) # side by side plot
# Left plot, no rasterisation
image(
    t(genmat),
    axes = FALSE,
    col = colorRampPalette(c("white", "darkorange", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Original matrix"
)
# Right plot, with rasterisation
image(
    t(genmat),
    axes = FALSE,
    col = colorRampPalette(c("white", "darkorange", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    useRaster = TRUE,
    main = "With rasterization"
)
dev.off()

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 :

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

). Je vous pro­pose cette petite fonc­tion (aus­si dis­po­nible sur canS­nip­pet) :

# reduce matrix size, using a summarizing function (default, mean)
redim_matrix <- function(
    mat,
    target_height = 100,
    target_width = 100,
    summary_func = function(x) mean(x, na.rm = TRUE),
    output_type = 0.0, #vapply style
    n_core = 1 # parallel processing
    ) {
    if(target_height > nrow(mat) | target_width > ncol(mat)) {
        stop("Input matrix must be bigger than target width and height.")
    }
    seq_height <- round(seq(1, nrow(mat), length.out = target_height + 1))
    seq_width  <- round(seq(1, ncol(mat), length.out = target_width  + 1))
    # complicated way to write a double for loop
    do.call(rbind, parallel::mclapply(seq_len(target_height), function(i) { # i is row
        vapply(seq_len(target_width), function(j) { # j is column
            summary_func(
                mat[
                    seq(seq_height[i], seq_height[i + 1]),
                    seq(seq_width[j] , seq_width[j + 1] )
                    ]
            )
        }, output_type)
    }, mc.cores = n_core))
}
genmatred <- redim_matrix(genmat, target_height = 600, target_width = 50) # 600 is very rougthly the pixel height of the image.
genmatred[1:5, 1:5]
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.4530226 0.4911097 0.4927123 0.5302643 0.5561331
## [2,] 0.5263392 0.5138786 0.5324716 0.5354325 0.5050932
## [3,] 0.4196155 0.4887105 0.5238630 0.5183627 0.5296764
## [4,] 0.5024431 0.5015508 0.5155568 0.5537814 0.5318501
## [5,] 0.5121447 0.5533040 0.4882006 0.4877140 0.5222805
dim(genmatred)
## [1] 600  50

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

layout(matrix(c(1, 2), nrow = 1))

# left plot, original matrix
image(
    t(genmat),
    axes = FALSE,
    col = colorRampPalette(c("white", "darkorange", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Original matrix"
)
box()

# Right plot, reduced matrix
image(
    t(genmatred),
    axes = FALSE,
    col = colorRampPalette(c("white", "darkorange", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Reduced matrix"
)
box()

par(oldpar) # restoring margin size to default values

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 :

# left matrix, we take the mean of the -log10 of the p-values
redim_matrix(
  eqtls,
  target_height = 600, target_width = 600,
  summary_func = function(x) mean(x, na.rm = TRUE),
  n_core = 14
)
# right matrix we take the maximum of the -log10 of the p-values
redim_matrix(
  eqtls,
  target_height = 600, target_width = 600,
  summary_func = function(x) max(x, na.rm = TRUE),
  n_core = 14
)

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

au lieu d'

.

# some fake data
mycolors <- matrix(c(
  sample(c("#0000FFFF", "#FFFFFFFF", "#FF0000FF"), size = 5000, replace = TRUE, prob = c(2, 1, 1)),
  sample(c("#0000FFFF", "#FFFFFFFF", "#FF0000FF"), size = 5000, replace = TRUE, prob = c(1, 2, 1)),
  sample(c("#0000FFFF", "#FFFFFFFF", "#FF0000FF"), size = 5000, replace = TRUE, prob = c(1, 1, 2))
), ncol = 1)
color_mat <- t(as.matrix(mycolors))

# custom funtion to average HTML colors
mean_color <- function(mycolors) {
    R     <- strtoi(x = substr(mycolors,2,3), base = 16)
    G     <- strtoi(x = substr(mycolors,4,5), base = 16)
    B     <- strtoi(x = substr(mycolors,6,7), base = 16)
    alpha <- strtoi(x = substr(mycolors,8,9), base = 16)

    return(
        rgb(
            red   = round(mean(R)),
            green = round(mean(G)),
            blue  = round(mean(B)),
            alpha = round(mean(alpha)),
            maxColorValue = 255
        )
    )
}

# Let's apply the redim_matrix() function using ou newly defined mean_color() function:
color_mat_red <- redim_matrix(
  color_mat,
  target_height = 1,
  target_width = 500,
  summary_func = mean_color,
  output_type = "string"
)

# And do the plotting
layout(matrix(c(1, 2), nrow = 2))

# left plot, original matrix
plot(c(0,1), c(0,1), axes = FALSE, type = "n",  xlab = "", ylab = "", xlim = c(0, 1), ylim = c(0,1), xaxs="i", yaxs="i", main = "Full matrix")
rasterImage(
    color_mat,
    xleft   = 0,
    xright  = 1,
    ybottom = 0,
    ytop    = 1
)
box()

# right plot, summarised matrix
plot(c(0,1), c(0,1), axes = FALSE, type = "n",  xlab = "", ylab = "", xlim = c(0, 1), ylim = c(0,1), xaxs="i", yaxs="i", main = "Reduced matrix")
rasterImage(
    color_mat_red,
    xleft   = 0,
    xright  = 1,
    ybottom = 0,
    ytop    = 1
)
box()

ggplot2

D’après mes tests,

est aus­si affec­té par ce sou­ci d'overplotting, que ce soit

ou

(qui est une ver­sion opti­mi­sée de geom_​tile() quand les cases sont régu­lières).

library(ggplot2)
library(patchwork)
# Wide to long transformation
data_for_ggplot <- as.data.frame(genmat) %>% 
    mutate(row = rownames(.)) %>% 
    tidyr::pivot_longer(-row, names_to = "col") %>%
    mutate(row = as.numeric(row), col = readr::parse_number(col))
    
# with geom_tile()
p1 <- ggplot(data_for_ggplot, aes(x = col, y = row, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(
      low = "white", mid = "darkorange", high = "black",
      limits = c(0, 3), midpoint = 1.5, oob = scales::squish
  ) +
  labs(title = "geom_tile") +
  theme_void() +
  theme(legend.position = "none")
# with geom_raster()
p2 <- ggplot(data_for_ggplot, aes(x = col, y = row, fill = value)) +
  geom_raster() +
  scale_fill_gradient2(
      low = "white", mid = "darkorange", high = "black",
      limits = c(0, 3), midpoint = 1.5, oob = scales::squish
  ) +
  labs(title = "geom_raster") +
  theme_void()
  theme(legend.position = "none")
p1 + p2

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 :

library(ComplexHeatmap)
Heatmap(
  genmat[nrow(genmat):1, ], # putting the top on top
  col = circlize::colorRamp2(c(0, 1.5, 3), c("white", "darkorange", "black")),
  cluster_rows = FALSE, cluster_columns = FALSE,
  show_heatmap_legend = FALSE,
  column_title = "No rasterisation"
)

La fonc­tion

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 :

Heatmap(
  genmat[nrow(genmat):1, ],
  col = circlize::colorRamp2(c(0, 1.5, 3), c("white", "darkorange", "black")),
  cluster_rows = FALSE, cluster_columns = FALSE,
  show_heatmap_legend = FALSE,
  use_raster = TRUE,
  raster_resize = TRUE, raster_device = "png",
  column_title = "With rasterisation"
)

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.