歡迎關注天下博客：http://blog.genesino.com/2017/06/heatmap-beautiful/
Jump to…
熱圖美化
對數轉換
Z-score轉換
抹去異常值
非線性顏色
調整行或列的順序
生信寶典，一起換個角度學生信
熱圖美化
上一期的繪圖命令中，最后一行的操作抹去了之前設定的橫軸標記的旋轉，最后出來的圖比較難看。

上次我們是這么寫的

p <- p + xlab(“samples”) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank())
為了使橫軸旋轉45度，需要把這句話theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1))放在theme_bw()的后面。

p <- p + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1))
最后的圖應該是下邊樣子的。

上次的測試數據，數值的分布比較均一，相差不是太大，但是Gene_4和Gene_5由于整體的值低于其它的基因，從顏色上看，不仔細看，看不出差別。

data <- c(rnorm(5,mean=5), rnorm(5,mean=20), rnorm(5, mean=100), c(600,700,800,900,10000))
data <- matrix(data, ncol=5, byrow=T)
data <- as.data.frame(data)
rownames(data) <- letters[1:4]
colnames(data) <- paste(“Grp”, 1:5, sep=”_”)
data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 5.958073 5.843652 3.225465 4.886184 3.411362
b 19.630582 20.376791 20.744580 18.534027 20.638288
c 100.351299 99.849900 102.197343 98.583629 99.540488
d 600.000000 700.000000 800.000000 900.000000 10000.000000
data$ID <- rownames(data)
data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5 ID
a 5.958073 5.843652 3.225465 4.886184 3.411362 a
b 19.630582 20.376791 20.744580 18.534027 20.638288 b
c 100.351299 99.849900 102.197343 98.583629 99.540488 c
d 600.000000 700.000000 800.000000 900.000000 10000.000000 d
data_m <- melt(data, id.vars=c(“ID”))
head(data_m)
ID variable value
1 a Grp_1 5.958073
2 b Grp_1 19.630582
3 c Grp_1 100.351299
4 d Grp_1 600.000000
5 a Grp_2 5.843652
6 b Grp_2 20.376791
p <- ggplot(data_m, aes(x=variable,y=ID)) + xlab(“samples”) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + theme(legend.position=”top”) + geom_tile(aes(fill=value)) + scale_fill_gradient(low = “white”, high = “red”)
p
dev.off()
輸出的結果是這個樣子的

圖中只有右上角可以看到紅色，其他地方就沒了顏色的差異。這通常不是我們想要的。為了更好的可視化效果，需要對數據做些預處理，主要有 對數轉換，Z-score轉換，抹去異常值，非線性顏色等方式。

對數轉換

假設下面的數據是基因表達數據，4個基因 (a, b, c, d)和5個樣品 (Grp_1, Grp_2, Grp_3, Grp_4)，矩陣中的值代表基因表達FPKM值。

data <- c(rnorm(5,mean=5), rnorm(5,mean=20), rnorm(5, mean=100), c(600,700,800,900,10000))
data <- matrix(data, ncol=5, byrow=T)
data <- as.data.frame(data)
rownames(data) <- letters[1:4]
colnames(data) <- paste(“Grp”, 1:5, sep=”_”)
data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 6.61047 20.946720 100.133106 600.000000 5.267921
b 20.80792 99.865962 700.000000 3.737228 19.289715
c 100.06930 800.000000 6.252753 21.464081 98.607518
d 900.00000 3.362886 20.334078 101.117728 10000.000000
data_log <- log2(data+1)
data_log
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 2.927986 4.455933 6.660112 9.231221 2.647987
b 4.446780 6.656296 9.453271 2.244043 4.342677
c 6.659201 9.645658 2.858529 4.489548 6.638183
d 9.815383 2.125283 4.415088 6.674090 13.287857
data_log$ID = rownames(data_log)
data_log_m = melt(data_log, id.vars=c(“ID”))

p <- ggplot(data_log_m, aes(x=variable,y=ID)) + xlab(“samples”) + ylab(NULL) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + theme(legend.position=”top”) + geom_tile(aes(fill=value)) + scale_fill_gradient(low = “white”, high = “red”)
ggsave(p, filename=”heatmap_log.pdf”, width=8, height=12, units=c(“cm”),colormodel=”srgb”)
對數轉換后的數據，看起來就清晰的多了。而且對數轉換后，數據還保留著之前的變化趨勢，不只是基因在不同樣品之間的表達可比 (同一行的不同列)，不同基因在同一樣品的值也可比 (同一列的不同行) (不同基因之間比較表達值存在理論上的問題，即便是按照長度標準化之后的FPKM也不代表基因之間是完全可比的)。

Z-score轉換

Z-score又稱為標準分數，是一組數中的每個數減去這一組數的平均值再除以這一組數的標準差，代表的是原始分數距離原始平均值的距離，以標準差為單位。可以對不同分布的各原始分數進行比較，用來反映數據的相對變化趨勢，而非絕對變化量。

data_ori <- “Grp_1;Grp_2;Grp_3;Grp_4;Grp_5
a;6.6;20.9;100.1;600.0;5.2
b;20.8;99.8;700.0;3.7;19.2
c;100.0;800.0;6.2;21.4;98.6
d;900;3.3;20.3;101.1;10000”

data <- read.table(text=data_ori, header=T, row.names=1, sep=”;”, quote=”“)

去掉方差為0的行，也就是值全都一致的行

data <- data[apply(data,1,var)!=0,]

data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 6.6 20.9 100.1 600.0 5.2
b 20.8 99.8 700.0 3.7 19.2
c 100.0 800.0 6.2 21.4 98.6
d 900.0 3.3 20.3 101.1 10000.0

標準化數據，并轉換為data.frame

data_scale <- as.data.frame(t(apply(data,1,scale)))

重命名列

colnames(data_scale) <- colnames(data)
data_scale
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a -0.5456953 -0.4899405 -0.1811446 1.7679341 -0.5511538
b -0.4940465 -0.2301542 1.7747592 -0.5511674 -0.4993911
c -0.3139042 1.7740182 -0.5936858 -0.5483481 -0.3180801
d -0.2983707 -0.5033986 -0.4995116 -0.4810369 1.7823177
data_scale$ID = rownames(data_scale)
data_scale_m = melt(data_scale, id.vars=c(“ID”))

p <- ggplot(data_scale_m, aes(x=variable,y=ID)) + xlab(“samples”) + ylab(NULL) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + geom_tile(aes(fill=value)) + scale_fill_gradient(low = “white”, high = “red”)
ggsave(p, filename=”heatmap_scale.pdf”, width=8, height=12, units=c(“cm”),colormodel=”srgb”)
Z-score轉換后，顏色分布也相對均一了，每個基因在不同樣品之間的表達的高低一目了然。但是不同基因之間就完全不可比了。

抹去異常值

粗暴一點，假設檢測飽和度為100，大于100的值都視為100對待。

data_ori <- “Grp_1;Grp_2;Grp_3;Grp_4;Grp_5
a;6.6;20.9;100.1;600.0;5.2
b;20.8;99.8;700.0;3.7;19.2
c;100.0;800.0;6.2;21.4;98.6
d;900;3.3;20.3;101.1;10000”

data <- read.table(text=data_ori, header=T, row.names=1, sep=”;”, quote=”“)

data[data>100] <- 100
data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 6.6 20.9 100.0 100.0 5.2
b 20.8 99.8 100.0 3.7 19.2
c 100.0 100.0 6.2 21.4 98.6
d 100.0 3.3 20.3 100.0 100.0
data$ID = rownames(data)
data_m = melt(data, id.vars=c(“ID”))

p <- ggplot(data_m, aes(x=variable,y=ID)) + xlab(“samples”) + ylab(NULL) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + geom_tile(aes(fill=value)) + scale_fill_gradient(low = “white”, high = “red”)
ggsave(p, filename=”heatmap_nooutlier.pdf”, width=8, height=12, units=c(“cm”),colormodel=”srgb”)
雖然損失了一部分信息，但整體模式還是出來了。只是在選擇異常值標準時需要根據實際確認。

非線性顏色

正常來講，顏色的賦予在最小值到最大值之間是均勻分布的。非線性顏色則是對數據比較小但密集的地方賦予更多顏色，數據大但分布散的地方賦予更少顏色，這樣既能加大區分度，又最小的影響原始數值。通常可以根據數據模式，手動設置顏色區間。為了方便自動化處理，我一般選擇用四分位數的方式設置顏色區間。

data_ori <- “Grp_1;Grp_2;Grp_3;Grp_4;Grp_5
a;6.6;20.9;100.1;600.0;5.2
b;20.8;99.8;700.0;3.7;19.2
c;100.0;800.0;6.2;21.4;98.6
d;900;3.3;20.3;101.1;10000”

data <- read.table(text=data_ori, header=T, row.names=1, sep=”;”, quote=”“)

data
Grp_1 Grp_2 Grp_3 Grp_4 Grp_5
a 6.6 20.9 100.1 600.0 5.2
b 20.8 99.8 700.0 3.7 19.2
c 100.0 800.0 6.2 21.4 98.6
d 900.0 3.3 20.3 101.1 10000.0
data$ID = rownames(data)
data_m = melt(data, id.vars=c(“ID”))

獲取數據的最大、最小、第一四分位數、中位數、第三四分位數

summary_v <- summary(data_m$value)
summary_v
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.30 16.05 60.00 681.40 225.80 10000.00

在最小值和第一四分位數之間劃出6個區間，第一四分位數和中位數之間劃出6個區間，中位數和第三四分位數之間劃出5個區間，最后的數劃出5個區間

break_v <- unique(c(seq(summary_v[1]*0.95,summary_v[2],length=6),seq(summary_v[2],summary_v[3],length=6),seq(summary_v[3],summary_v[5],length=5),seq(summary_v[5],summary_v[6]*1.05,length=5)))
break_v
[1] 3.135 5.718 8.301 10.884 13.467 16.050 24.840
[8] 33.630 42.420 51.210 60.000 101.450 142.900 184.350
[15] 225.800 2794.350 5362.900 7931.450 10500.000

安照設定的區間分割數據

原始數據替換為了其所在的區間的數值

data_mvalue<?cut(datamvalue<?cut(datamvalue, breaks=break_v,labels=break_v[2:length(break_v)])
break_v=unique(data_m$value)

data_m
ID variable value
1 a Grp_1 8.301
2 b Grp_1 24.84
3 c Grp_1 101.45
4 d Grp_1 2794.35
5 a Grp_2 24.84
6 b Grp_2 101.45
7 c Grp_2 2794.35
8 d Grp_2 5.718
9 a Grp_3 101.45
10 b Grp_3 2794.35
11 c Grp_3 8.301
12 d Grp_3 24.84
13 a Grp_4 2794.35
14 b Grp_4 5.718
15 c Grp_4 24.84
16 d Grp_4 101.45
17 a Grp_5 5.718
18 b Grp_5 24.84
19 c Grp_5 101.45
20 d Grp_5 10500

雖然看上去還是數值，但已經不是數字類型了

而是不同的因子了，這樣就可以對不同的因子賦予不同的顏色了

is.numeric(data_mvalue)[1]FALSEis.factor(datamvalue)[1]FALSEis.factor(datamvalue)
[1] TRUE
break_v
[1] 8.301 24.84 101.45 2794.35 5.718 10500
18 Levels: 5.718 8.301 10.884 13.467 16.05 24.84 33.63 42.42 51.21 … 10500

產生對應數目的顏色

gradientC=c(‘green’,’yellow’,’red’)
col <- colorRampPalette(gradientC)(length(break_v))
col
[1] “#00FF00” “#66FF00” “#CCFF00” “#FFCB00” “#FF6500” “#FF0000”
p <- ggplot(data_m, aes(x=variable,y=ID)) + xlab(“samples”) + ylab(NULL) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + geom_tile(aes(fill=value))

與上面不同的地方，使用的是scale_fill_manual逐個賦值

p <- p + scale_fill_manual(values=col)
ggsave(p, filename=”heatmap_nonlinear.pdf”, width=8, height=12, units=c(“cm”),colormodel=”srgb”)

調整行或列的順序

如果想保持圖中每一行的順序與輸入的數據框一致，需要設置因子的水平。這也是ggplot2中調整圖例或橫縱軸字符順序的常用方式。

data_rowname <- rownames(data)
data_rowname <- as.vector(rownames(data))
data_rownames <- rev(data_rowname)
data_log_mID<?factor(datalogmID<?factor(datalogmID, levels=data_rownames, ordered=T)
p <- ggplot(data_log_m, aes(x=variable,y=ID)) + xlab(NULL) + ylab(NULL) + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) + theme(legend.position=”top”) + geom_tile(aes(fill=value)) + scale_fill_gradient(low = “white”, high = “red”)
ggsave(p, filename=”heatmap_log.pdf”, width=8, height=12, units=c(“cm”),colormodel=”srgb”)

基于ggplot2的heatmap繪制到現在就差不多了，但總是這么畫下去也會覺得有點累，有沒有辦法更簡化呢？ 且聽下回分解。

生信寶典，一起換個角度學生信

R語言學習 - 熱圖繪制 (heatmap)
R語言學習 - 基礎概念和矩陣操作
R語言學習 - 熱圖美化
R語言學習 - 熱圖簡化
RBIOINFO
CHENTONG
版權聲明：本文為博主原創文章，轉載請注明出處。
alipay.png WeChatPay.png

總結

以上是生活随笔為你收集整理的R统计绘图 - 热图美化的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。