专栏名称: 生信宝典
生物信息分析入门、晋级和经验分享。Linux、R、Python学习教程;高通量测序数据分析学习教程;生信软件安装教程。所有内容均为原创分享,致力于从基础学习到提高整个过程。
目录
相关文章推荐
生信宝典  ·  微生物组-扩增子16S分析和可视化(2025.4) ·  昨天  
生物制品圈  ·  刚刚!科伦博泰ADC新药获批肺癌新适应症 ·  2 天前  
生信人  ·  微生物助力中药研究,拿下一区顶刊 ·  4 天前  
51好读  ›  专栏  ›  生信宝典

R语言学习 - 热图美化

生信宝典  · 公众号  · 生物  · 2017-07-03 07:51

正文

热图美化

上一期的绘图命令中,最后一行的操作抹去了之前设定的横轴标记的旋转,最后出来的图比较难看。 上次我们是这么写的

p "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 45,hjust=1, vjust=1))

最后的图应该是下边样子的。

上图的测试数据,数值的分布比较均一,相差不是太大,但是Gene_4和Gene_5由于整体的值低于其它的基因,从颜色上看,不仔细看,看不出差别。

实际应用中,异常值的出现会毁掉一张热图,如下是一个例子。

data 5,mean=5), rnorm(5,mean=20), rnorm(5, mean=100), c(600,700,800,900,10000))
data 5, byrow=T)
data rownames(data) 1:4]
colnames(data) "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 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 "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 "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 5,mean=5), rnorm(5,mean=20), rnorm(5, mean=100), c(600,700,800,900,10000))
data 5, byrow=T)
data rownames(data) 1:4]
colnames(data) "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
# 对数转换
# +1是为了防止对0取对数;是加1还是加个更小的值取决于数据的分布。
# 加的值一般认为是检测的低阈值,低于这个值的数字之间的差异可以忽略。
data_log 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 "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 T, row.names=1, sep=";", quote="")

# 去掉方差为0的行,也就是值全都一致的行
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
# 标准化数据,获得Z-score,并转换为data.frame
data_scale 1,scale)))

# 重命名列
colnames(data_scale) 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"






请到「今天看啥」查看全文