Accessibility Tools

Web Accessibility plugin by DJ-Extensions.com

- Le blog participatif de bioinformatique francophone depuis 2012 -

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

-

par

Guillaume Devailly
dans

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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
[crayon-67e51d0157ac4833411460  ]libra­ry(dplyr)
libra­ry(purrr)
Ncol <- 50
gen­mat <- map(
    1 :20000,
    func­tion(i) {
        runif(Ncol) + c(sort(abs(rnorm(50)))[1 :25], rev(sort(abs(rnorm(50)))[1 :25]))  * i/​5000
    }
) %>% do.call(rbind, .)
gen­mat[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(gen­mat)
## [1] 20000    50
[/​crayon]

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

1
heat­map()

,

1
gplots::heat­map.2()

,

1
ggplot2::geom_​raster()

, ou

1
Com­plex­Heat­map::Heat­map()

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

1
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 :

1
2
3
4
5
6
7
8
[crayon-67e51d0157ace497381598  ]old­par <- par(mar = rep(0.2, 4)) # redu­cing plot mar­gins
image(
    t(gen­mat), # image() has some weird opi­nions about how your matrix will be plot­ted
    axes = FALSE,
    col = color­Ramp­Pa­lette(c("white", "dar­ko­range", "black"))(30), # our colour palette
    breaks = c(seq(0, 3, length.out = 30), 100) # colour-to-value map­ping
)
box() # adding a box arround the heat­map
[/​crayon]
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

1
pdf()

, mais pas

1
png()

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

1
use­Ras­ter = TRUE

de la fonc­tion

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
[crayon-67e51d0157af0590115189  ]pdf("big_hm_1.pdf")
layout(matrix(c(1, 2), nrow = 1)) # side by side plot
# Left plot, no ras­te­ri­sa­tion
image(
    t(gen­mat),
    axes = FALSE,
    col = color­Ramp­Pa­lette(c("white", "dar­ko­range", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Ori­gi­nal matrix"
)
# Right plot, with ras­te­ri­sa­tion
image(
    t(gen­mat),
    axes = FALSE,
    col = color­Ramp­Pa­lette(c("white", "dar­ko­range", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    use­Ras­ter = TRUE,
    main = "With ras­te­ri­za­tion"
)
dev.off()
[/​crayon]

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

1
mean()

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
[crayon-67e51d0157af3082743101  ]# reduce matrix size, using a sum­ma­ri­zing func­tion (default, mean)
redim_​matrix <- func­tion(
    mat,
    target_​height = 100,
    target_​width = 100,
    summary_​func = func­tion(x) mean(x, na.rm = TRUE),
    output_​type = 0.0, #vap­ply style
    n_​core = 1 # paral­lel pro­ces­sing
    ) {
    if(target_​height > nrow(mat) | target_​width > ncol(mat)) {
        stop("Input matrix must be big­ger than tar­get 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))
    # com­pli­ca­ted way to write a double for loop
    do.call(rbind, paral­lel::mclap­ply(seq_​len(target_​height), func­tion(i) { # i is row
        vap­ply(seq_​len(target_​width), func­tion(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))
}
gen­ma­tred <- redim_​matrix(gen­mat, target_​height = 600, target_​width = 50) # 600 is very rougth­ly the pixel height of the image.
gen­ma­tred[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(gen­ma­tred)
## [1] 600  50
[/​crayon]

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
[crayon-67e51d0157af6303114101  ]layout(matrix(c(1, 2), nrow = 1))
 
# left plot, ori­gi­nal matrix
image(
    t(gen­mat),
    axes = FALSE,
    col = color­Ramp­Pa­lette(c("white", "dar­ko­range", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Ori­gi­nal matrix"
)
box()
 
# Right plot, redu­ced matrix
image(
    t(gen­ma­tred),
    axes = FALSE,
    col = color­Ramp­Pa­lette(c("white", "dar­ko­range", "black"))(30),
    breaks = c(seq(0, 3, length.out = 30), 100),
    main = "Redu­ced matrix"
)
box()
 
par(old­par) # res­to­ring mar­gin size to default values
[/​crayon]
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 :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
[crayon-67e51d0157af9125182088  ]# left matrix, we take the mean of the -log10 of the p-values
redim_​matrix(
  eqtls,
  target_​height = 600, target_​width = 600,
  summary_​func = func­tion(x) mean(x, na.rm = TRUE),
  n_​core = 14
)
# right matrix we take the maxi­mum of the -log10 of the p-values
redim_​matrix(
  eqtls,
  target_​height = 600, target_​width = 600,
  summary_​func = func­tion(x) max(x, na.rm = TRUE),
  n_​core = 14
)
[/​crayon]
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

1
ras­te­rI­mage()

au lieu d'

1
image()

.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
[crayon-67e51d0157afe623413478  ]# some fake data
myco­lors <- 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(myco­lors))
 
# cus­tom fun­tion to ave­rage HTML colors
mean_​color <- func­tion(myco­lors) {
    R     <- str­toi(x = sub­str(myco­lors,2,3), base = 16)
    G     <- str­toi(x = sub­str(myco­lors,4,5), base = 16)
    B     <- str­toi(x = sub­str(myco­lors,6,7), base = 16)
    alpha <- str­toi(x = sub­str(myco­lors,8,9), base = 16)
 
    return(
        rgb(
            red   = round(mean(R)),
            green = round(mean(G)),
            blue  = round(mean(B)),
            alpha = round(mean(alpha)),
            max­Co­lor­Va­lue = 255
        )
    )
}
 
# Let's apply the redim_​matrix() func­tion using ou new­ly defi­ned mean_​color() func­tion :
color_​mat_​red <- redim_​matrix(
  color_​mat,
  target_​height = 1,
  target_​width = 500,
  summary_​func = mean_​color,
  output_​type = "string"
)
 
# And do the plot­ting
layout(matrix(c(1, 2), nrow = 2))
 
# left plot, ori­gi­nal 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")
ras­te­rI­mage(
    color_​mat,
    xleft   = 0,
    xright  = 1,
    ybot­tom = 0,
    ytop    = 1
)
box()
 
# right plot, sum­ma­ri­sed 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 = "Redu­ced matrix")
ras­te­rI­mage(
    color_​mat_​red,
    xleft   = 0,
    xright  = 1,
    ybot­tom = 0,
    ytop    = 1
)
box()
[/​crayon]
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,

1
ggplot2

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

1
geom_​tile()

ou

1
geom_​raster()

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
[crayon-67e51d0157b0a567413622  ]libra­ry(ggplot2)
libra­ry(patch­work)
# Wide to long trans­for­ma­tion
data_​for_​ggplot <- as.data.frame(gen­mat) %>%
    mutate(row = row­names(.)) %>%
    tidyr::pivot_​longer(-row, names_​to = "col") %>%
    mutate(row = as.nume­ric(row), col = rea­dr::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 = "dar­ko­range", high = "black",
      limits = c(0, 3), mid­point = 1.5, oob = scales::squish
  ) +
  labs(title = "geom_​tile") +
  theme_​void() +
  theme(legend.posi­tion = "none")
# with geom_​raster()
p2 <- ggplot(data_​for_​ggplot, aes(x = col, y = row, fill = value)) +
  geom_​raster() +
  scale_​fill_​gradient2(
      low = "white", mid = "dar­ko­range", high = "black",
      limits = c(0, 3), mid­point = 1.5, oob = scales::squish
  ) +
  labs(title = "geom_​raster") +
  theme_​void()
  theme(legend.posi­tion = "none")
p1 + p2
[/​crayon]
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 :

1
2
3
4
5
6
7
8
[crayon-67e51d0157b0c478134200  ]libra­ry(Com­plex­Heat­map)
Heat­map(
  gen­mat[nrow(gen­mat):1, ], # put­ting the top on top
  col = cir­clize::colorRamp2(c(0, 1.5, 3), c("white", "dar­ko­range", "black")),
  cluster_​rows = FALSE, cluster_​columns = FALSE,
  show_​heatmap_​legend = FALSE,
  column_​title = "No ras­te­ri­sa­tion"
)
[/​crayon]
Figure 7 : Com­plex­Heat­map nous déçoit.

La fonc­tion

1
Heat­map()

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 :

1
2
3
4
5
6
7
8
9
[crayon-67e51d0157b12418395729  ]Heat­map(
  gen­mat[nrow(gen­mat):1, ],
  col = cir­clize::colorRamp2(c(0, 1.5, 3), c("white", "dar­ko­range", "black")),
  cluster_​rows = FALSE, cluster_​columns = FALSE,
  show_​heatmap_​legend = FALSE,
  use_​raster = TRUE,
  raster_​resize = TRUE, raster_​device = "png",
  column_​title = "With ras­te­ri­sa­tion"
)
[/​crayon]
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.

Vous avez aimé ? Dites-le nous !

Moyenne : 0 /​ 5. Nb de votes : 0

Pas encore de vote pour cet article.

We are sor­ry that this post was not use­ful for you !

Let us improve this post !

Tell us how we can improve this post ?

Partagez cet article

Guillaume Devailly

Guillaume Devailly

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

Pour continuer la lecture :

Commentaires

Laisser un commentaire

Pour insérer du code dans vos commentaires, utilisez les balises <code> et <\code>.

fr fr