百度360必应搜狗淘宝本站头条
mqtt服务器xftp6mybatiscollectionkeyerrorc#var
当前位置：网站首页 > 热门文章 > 正文

R语言学习 - 热图美化

bigegpt 2024-08-05 11:51 2 浏览

热图美化

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

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
# 对数转换
# +1是为了防止对0取对数；是加1还是加个更小的值取决于数据的分布。
# 加的值一般认为是检测的低阈值，低于这个值的数字之间的差异可以忽略。
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
# 标准化数据，获得Z-score，并转换为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_m$value <- cut(data_m$value, 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_m$value)
[1] FALSE
> is.factor(data_m$value)
[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_m$ID <- factor(data_log_m$ID, 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绘制到现在就差不多了，但总是这么画下去也会觉得有点累，有没有办法更简化呢？ 且听下回分解。

相关推荐

程序员请收好:10个非常有用的 Visual Studio Code 插件

一个插件列表，可以让你的程序员生活变得轻松许多。作者|Daan译者|Elle出品|CSDN（ID：CSDNnews）以下为译文：无论你是经验丰富的开发人员还是刚刚开始第一份工作的初级开发人...

PADS在WIN10系统中菜单显示不全的解决方法

决定由AD转PADS，打开发现菜单显示不正常，如下图所示：这个是由于系统的默认字体不合适导致，修改一下系统默认字体即可，修改方法如下：打开开始菜单-->所有程序-->Windows系统--...

一文讲解Web前端开发基础环境配置

先从基本的HTML语言开始学习。一个网页的所有内容都是基于HTML，为了学好HTML，不使用任何集成工具，而用一个文本编辑器，直接从最简单的HTML开始编写HTML。先在网上下载notepad++文...

TCP/IP协议栈在Linux内核中的运行时序分析

本文主要是讲解TCP/IP协议栈在Linux内核中的运行时序，文章较长，里面有配套的视频讲解，建议收藏观看。1Linux概述　　1.1Linux操作系统架构简介Linux操作系统总体上由Linux...

从 Angular Route 中提前获取数据

#头条创作挑战赛#介绍提前获取意味着在数据呈现在屏幕之前获取到数据。本文中，你将学到，在路由更改前怎么获取到数据。通过本文，你将学会使用resolver，在AngularApp中应用re...

边做游戏边划水: 基于浅水方程的水面交互、河道交互模拟方法

以下文章来源于腾讯游戏学堂，作者Byreave篇一：基于浅水方程的水面交互本文主要介绍一种基于浅水方程的水体交互算法，在基本保持水体交互效果的前提下，实现了一种极简的水面模拟和物体交互方法。真实感的...

Nacos介绍及使用

一、Nacos介绍Nacos是SpringCloudAlibaba架构中最重要的组件。Nacos是一个更易于帮助构建云原生应用的动态服务发现、配置和服务管理平台，提供注册中心、配置中心和动态DNS...

Spring 中@Autowired，@Resource，@Inject 注解实现原理

使用案例前置条件：现在有一个Vehicle接口，它有两个实现类Bus和Car，现在还有一个类VehicleService需要注入一个Vehicle类型的Bean：publicinte...

一文带你搞懂Vue3 底层源码

作者：妹红大大转发链接：https://mp.weixin.qq.com/s/D_PRIMAD6i225Pn-a_lzPA前言vue3出来有一段时间了。今天正式开始记录一下梗vue3.0.0-be...

一线开发大牛带你深度解析探讨模板解释器，解释器的生成

解释器生成解释器的机器代码片段都是在TemplateInterpreterGenerator::generate_all()中生成的，下面将分小节详细展示该函数的具体细节，以及解释器某个组件的机器代码...

Nacos源码—9.Nacos升级gRPC分析五

大纲10.gRPC客户端初始化分析11.gRPC客户端的心跳机制(健康检查)12.gRPC服务端如何处理客户端的建立连接请求13.gRPC服务端如何映射各种请求与对应的Handler处理类14.gRP...

聊聊Spring AI的Tool Calling

序本文主要研究一下SpringAI的ToolCallingToolCallbackorg/springframework/ai/tool/ToolCallback.javapublicinter...

「云原生」Containerd ctr，crictl 和 nerdctl 命令介绍与实战操作

一、概述作为接替Docker运行时的Containerd在早在Kubernetes1.7时就能直接与Kubelet集成使用，只是大部分时候我们因熟悉Docker，在部署集群时采用了默认的dockers...

在MySQL登录时出现Access denied for user ~~ (using password: YES)

Windows~~~在MySQL登录时出现Accessdeniedforuser‘root‘@‘localhost‘(usingpassword:YES)，并修改MySQL密码目录适用...

mysql 8.0多实例批量部署script

背景最近一个项目上，客户需要把阿里云的rdsformysql数据库同步至线下，用作数据的灾备，需要在线下的服务器上部署mysql8.0多实例，为了加快部署的速度，写了一个脚本。解决方案#!/bi...

一周热门
最近发表
标签列表