热图在单细胞数据分析中的应用
分享是一种态度
作者 | 周运来
男,
一个长大了才会遇到的帅哥,
稳健,潇洒,大方,靠谱。
一段生信缘,一棵技能树,
一枚大型测序工厂的螺丝钉,
一个随机森林中提灯觅食的津门旅客。
什么是热图?
热图是一个以颜色变化来显示数据的可视化矩阵,Toussaint Loua在1873年就曾使用过热图来绘制对巴黎各区的社会学统计。我们就拿这张简单朴素的热图来讲一下热图怎么看。
首先映入我们眼帘的是有的地方是黑的,有的地方是白的(颜色),每一块颜色都有对应的XY轴。言下之意,对象X的属性Y的值是用颜色表征的。颜色的聚集代表相应对象X的属性Y具有相似性(模式,pattern)。本质上它是表现一个数值矩阵,图上每一个小方格都是一个数值,按一条预设好的色彩变化尺(称为色键,Color Key),给每个数值分配颜色。
有时候我们还能看到对象X或者属性Y的聚类结果也绘制在热图的旁边,但是这就不属于热图的部分了,因为他已经不热了(热,就是有的地方冷,有的地方热)。
热图能说明哪些问题
表达量
广泛的应用就是用热图来可视化表达量。我们想象一下一个9个样本50个基因的表达谱,人类一眼看过去就是一堆数字,而表达量数值大小映射到颜色的深浅上,看起来就很清楚了。
很多时候,为了同一个基因在不同样本中的表达量有可比性,需要对表达量取对数,或取Z-score,把数据标准化到一个水平上。
相关性
计算两个矩阵的相关性,可以得到两两的相关性,这时,用热图的颜色来表示相关性可以看出哪些配对相关性较高。
在单细胞中的应用
表达量
seurat主题
这是一张典型的seurat做的热图,可以清楚地看出不同分群有着不同的表达模式。这里的每一个色块是一个细胞某基因的表达量。cluster可以看做是细胞的聚类,Y轴的基因,我们看到也是聚类了的(很可能是手动的,每一类基因作者都给出了注释)。所以这张热图的关键是什么?细胞和基因及其顺序。选择合适的细胞和基因(一般是每个群的差异高表达基因)后,为什么我做的图是一团黑?很可能是因为:
- 基因的顺序没排好。
- 没有标准化
人们经常需要根据差异分析的结果来探索基因列表的排序,如SC3的策略。差异基因的计算采用非参数Kruskal-Wallis检验。SC3提供了调整p值< 0.01的所有差异表达基因的列表,并绘制了p值最低的50个基因的基因表达谱。值得注意的是,聚类后的差异表达计算可能会在p值的分布中引入偏差,因此我们建议仅使用p值对基因进行排序。
SC3主题
这类图无疑反映了某geneList在某cluster的表达情况。如果巧了,这个geneList是某个细胞类型的marker基因,或者是某个功能的主要集合,热图有助于细胞群功能和类型的鉴定。热图很好地将对象(X,一般是我们的细胞)与它的属性(Y,一般是我们的基因)联系起来。
scanpy主题
在monocle2 中我们还看到一种热图将基因的表达情况与细胞发育轨迹结合到一起。可视化所有明显依赖于分支的基因的变化(如果愿意也可以自己定义geneList)。这张热图同时显示了两种命运的变化,它还要求选择分支点(branch_point )。列是伪时间中的点,行是基因,伪时间的开始在热图的中间。当你从热图的中间读到右边的时候,你正在跟随一个伪时间谱系。当你读到左边时,另一个。这些基因是分层聚类的,因此您可以可视化具有类似的依赖于序列的基因模块。
monocle2主题
logFC
有时候为了突出展示群与群之间的差异基因,可以把差异分析的指标反映在热图上,如10X的loup软件那样:
Question 2: Which genes are displayed in the Loupe Cell Browser heatmap? How can I export the list of genes displayed?
Answer: In Loupe Cell Browser (version 2.0.0), the heatmap is a compact display of a subset of differentially expressed genes per cluster. Specifically, the gene list is the union of the top 120/N upregulated genes for each cluster ranked by log2 fold-change (N=total number of clusters).The gene names are on the plot when you export the heatmap. However, as of version 2.0.0, there is no 1-click function to export the associated information for the subset of heatmap genes.
cloup主题
相关性
提到相关性,我们很容易注意到WGCNA(weighted correlation network analysis,加权基因共表达网络分析), 用于提取与性状或临床特征相关的基因模块,解析与表达量相关生物学过程。这是除了富集分析之外另一个寻找好的geneList的方法。这里的颜色不再是表达量的度量而是相似性的度量。
WGCNA主题
ComplexHeatmap在单细胞数据可视化中的应用
人们针对单细胞发展了相应的数据结构如seurat的S4类,monocle的CDS,SingleCellExperiment的sce,scanpy的anndata等,可见单细胞的故事远大于一张二维的表达谱。那么一张热图往往也不能完全的说明问题,于是我们希望能够灵活地操纵热图来讲更多的故事。于是,我们发现ComplexHeatmap这个R包真的是热图神器。
ComplexHeatmap 主题
library(gdata)
library(ComplexHeatmap)
library(GetoptLong)
library(circlize)
getwd()
#?read.xls
expr <-readxl::read_excel("ComplexHeatmap//data//41587_2015_BFnbt3102_MOESM18_ESM.xlsx",sheet=2)
#View(expr)
#expr <- as.matrix(expr)
expr = expr[!duplicated(expr[[1]]), ]
rownames(expr) = expr[[1]]
#expr = expr[-1]
colnames(expr<-as.data.frame(expr))
expr = expr[,-1]
expr = t(expr)
#expr = expr[-1,]
expr[1:3,1:3];dim(expr)
Cell 1 Cell 2 Cell 3
Gnai3 3.23220 1.98320 2.2482
Cdc45 3.19810 1.17300 3.1705
Narf 0.29411 0.49389 1.6279
[1] 7073 81
expr = expr[apply(expr, 1, function(x) sum(x > 0)/length(x) > 0.5), , drop = FALSE]
get_correlated_variable_genes = function(mat, n = nrow(mat), cor_cutoff = 0, n_cutoff = 0) {
ind = order(apply(mat, 1, function(x) {
q = quantile(x, c(0.1, 0.9))
x = x[x < q[1] & x > q[2]]
var(x)/mean(x)
}), decreasing = TRUE)[1:n]
mat2 = mat[ind, , drop = FALSE]
dt = cor(t(mat2), method = "spearman")
diag(dt) = 0
dt[abs(dt) < cor_cutoff] = 0
dt[dt < 0] = -1
dt[dt > 0] = 1
i = colSums(abs(dt)) > n_cutoff
mat3 = mat2[i, ,drop = FALSE]
return(mat3)
}
mat = get_correlated_variable_genes(expr, cor_cutoff = 0.5, n_cutoff = 20)
mat[1:4,1:4];dim(mat)
Cell 1 Cell 2 Cell 3 Cell 4
Sdhd 3.7634 2.7851 3.0983 3.9585
Ccnd2 3.5098 4.2493 3.7505 4.5186
Dazap2 2.9348 2.0173 3.0520 4.7106
Lck 3.8987 1.8177 3.6390 4.2300
[1] 722 81
mat2 = t(apply(mat, 1, function(x) {
q10 = quantile(x, 0.1)
q90 = quantile(x, 0.9)
x[x < q10] = q10
x[x > q90] = q90
scale(x)
}))
colnames(mat2) = colnames(mat)
mat2[1:4,1:4]
Cell 1 Cell 2 Cell 3 Cell 4
Sdhd 1.4222656 -0.1033475 0.38507327 1.458445
Ccnd2 -1.1291769 0.8112796 -0.49757737 1.517926
Dazap2 -0.2861832 -1.5546649 -0.05732535 1.585296
Lck 0.5384592 -1.5226640 -0.13645403 1.399448
base_mean = rowMeans(mat)
cc = readRDS("ComplexHeatmap//data//mouse_cell_cycle_gene.rds")
head(cc)
ENSMUSG00000019256 ENSMUSG00000021866 ENSMUSG00000017716 ENSMUSG00000007815 ENSMUSG00000034218 ENSMUSG00000042750
"Ahr" "Anxa11" "Birc5" "Rhoa" "Atm" "Bex2"
ccl = rownames(mat) %in% cc
table(cc1)
cc_gene = rownames(mat)[ccl]
rp = readRDS("ComplexHeatmap//data//mouse_ribonucleoprotein.rds")
head(rp)
ENSMUSG00000029580 ENSMUSG00000028427 ENSMUSG00000067274 ENSMUSG00000000568 ENSMUSG00000022557 ENSMUSG00000017146
"Actb" "Aqp7" "Rplp0" "Hnrnpd" "Bop1" "Brca1"
rpl = rownames(mat) %in% rp
ht_list = Heatmap(mat2, col = colorRamp2(c(-1.5, 0, 1.5), c("blue", "white", "red")),
name = "scaled_expr", column_title = qq("relative expression for @{nrow(mat)} genes"),
show_column_names = FALSE, width = unit(8, "cm"),
heatmap_legend_param = list(title = "Scaled expr")) +
Heatmap(base_mean, name = "base_expr", width = unit(5, "mm"),
heatmap_legend_param = list(title = "Base expr")) +
Heatmap(rpl + 0, name = "ribonucleoprotein", col = c("0" = "white", "1" = "purple"),
show_heatmap_legend = FALSE, width = unit(5, "mm")) +
Heatmap(ccl + 0, name = "cell_cycle", col = c("0" = "white", "1" = "red"),
show_heatmap_legend = FALSE, width = unit(5, "mm")) +
rowAnnotation(link = anno_mark(at = which(ccl & base_mean > quantile(base_mean, 0.25)),
labels = rownames(mat)[ccl & base_mean > quantile(base_mean, 0.25)],
labels_gp = gpar(fontsize = 10), padding = unit(1, "mm"))) +
Heatmap(cor(t(mat2)), name = "cor",
col = colorRamp2(c(-1, 0, 1), c("green", "white", "red")),
show_row_names = FALSE, show_column_names = FALSE, row_dend_side = "right",
show_column_dend = FALSE, column_title = "pairwise correlation between genes",
heatmap_legend_param = list(title = "Correlation"))
ht_list = draw(ht_list, main_heatmap = "cor")
decorate_column_dend("scaled_expr", {
tree = column_dend(ht_list)$scaled_expr
ind = cutree(as.hclust(tree), k = 2)[order.dendrogram(tree)]
first_index = function(l) which(l)[1]
last_index = function(l) { x = which(l); x[length(x)] }
x1 = c(first_index(ind == 1), first_index(ind == 2)) - 1
x2 = c(last_index(ind == 1), last_index(ind == 2))
grid.rect(x = x1/length(ind), width = (x2 - x1)/length(ind), just = "left",
default.units = "npc", gp = gpar(fill = c("#FF000040", "#00FF0040"), col = NA))
})
dev.off()
数据可视化的过程就是一段探索意义的旅程,给每一种颜色、每一种形状、每一种聚集和离散找到一种生物学意义。这让我想起海子的《面朝大海,春暖花开》:
给每一条河每一座山取一个温暖的名字
陌生人,我也为你祝福
愿你有一个灿烂的前程
愿你有情人终成眷属
愿你在尘世获得幸福
我只愿面朝大海,春暖花开
References
[1]
ComplexHeatmap: https://links.jianshu.com/go?to=https%3A%2F%2Fjokergoo.github.io%2FComplexHeatmap-reference%2Fbook%2F
[2]
R数据可视化3:热图: https://links.jianshu.com/go?to=https%3A%2F%2Fdb.cngb.org%2Fscience%2FARTtoo54%2F
[3]
如何画热图: https://links.jianshu.com/go?to=https%3A%2F%2Fzhuanlan.zhihu.com%2Fp%2F21258579
[4]
10秒钟-完美掌握-热图(heatmap)绘制 - 所有人都可以!: https://www.jianshu.com/p/fb1170d03b97
推荐阅读
-
手绘 | 深入分析机器学习在 8 个风控场景中的应用
-
ZH-HPLC/RS485 传输模块在光伏逆变器发电数据采集器中的应用 + 宽带电力线载波 HPLC + 站式智能汇聚终端在分布式光伏并网数据采集与编程中的应用特点
-
35 岁实现财务*,腾讯程序员手握2300万提前退休?-1000万房产、1000万腾讯股票、加上300万的现金,一共2300万的财产。有网友算了一笔账,假设1000万的房产用于自住,剩下1300万资产按照平均税后20-50万不等进行计算,大约花上26-60年左右的时间才能赚到这笔钱。也就是说,普通人可能奋斗一辈子,才能赚到这笔钱。在很多人还在为中年危机而惶惶不可终日的时候,有的人的35岁,就已经安全着陆,试问哪个打工人不羡慕?但问题是有这样财富积累必然有像样的实力做靠山。没有人可以不劳而获。 看到这里,肯定有人说,那么对于普通人来说,卷可能真就成了唯一的出路。但是卷也有轻松的卷,“偷懒”的卷法,对于程序员而言,刨除掉一时无法改掉的开会传统占用的大部分时间,如何把有限的时间和精力放在真正重要的架构设计、需求设计上,而不是重复的造*,编码、改bug、手动测试。因此在科技改变生活的今天,学会使用AI工具成为程序员们的必备技能。 以全栈式全自动的软件开发工具飞算SoFlu软件机器人为例,作为全球首款面向微服务架构设计和最佳实践的软件机器人,SoFlu软件机器人改变了原来手工编码的作业模式,通过可视化拖拽方式以及参数配置就能实现等同于编写复杂代码的业务逻辑,在设计业务逻辑时就完成了微服务应用开发,做到“业务即图,图即代码”,实现“软件开发,十倍提效”。 而作为飞算SoFlu软件机器人的一个重要组成部分,FuncGPT(慧函数)支持所有类型函数创建。通过自然语言描述Java函数需求,实时生成高质量、高可读性的Java函数代码。生成代码可直接复制到IDEA,或一键导入Java全自动开发工具函数库。 FuncGPT(慧函数)具有五大核心能力: ● 自然语言:通过自然语言即可生成函数,降低软件开发门槛。 ● 秒级函数开发:全自动秒级开发函数,秒级完成,效率千倍提升。 ● 一站式函数开发应用:函数生成即可直接应用,效率、安全有保障。 ● 质量可靠:生成的函数代码符合业界规范、具备良好的可读性,符合最优解。 ● 灵活生成:可以根据开发人员的具体需求生成并支持调整和修改。 以“将文件流转换为File对象”这个需求为例,以下是FuncGPT(慧函数)秒级生成的代码: / 类名:FileConverter // 函数名:convertStreamToFile // 函数功能:将文件流转换为File对象 // POM依赖包:无 import java.io.File; import java.io.FileOutputStream; import java.io.IOException; import java.io.InputStream; public class FileConverter { /** * 将文件流转换为File对象 * @param inputStream 文件流 * @param filePath 文件路径 * @param fileName 文件名 * @return 转换后的File对象 * @throws IOException 如果发生I/O异常 */ public static File convertStreamToFile(InputStream inputStream, String filePath, String fileName) throws IOException { File file = new File(filePath + File.separator + fileName); // 创建File对象 try (FileOutputStream outputStream = new FileOutputStream(file)) { // 创建文件输出流 byte buffer = new byte[1024]; int bytesRead; while ((bytesRead = inputStream.read(buffer)) != -1) { // 从文件流读取数据并写入文件 outputStream.write(buffer, 0, bytesRead); } } return file; // 返回转换后的File对象 } } // 函数示例 // 将文件流转换为File对象示例 // 入参:inputStream,文件流 // 入参:filePath,文件路径 // 入参:fileName,文件名 // 出参:file,转换后的File对象 // 调用示例: // InputStream inputStream = new FileInputStream("example.txt"); // String filePath = "C:\\Users\\User\\Documents"; // String fileName = "example.txt"; // File file = FileConverter.convertStreamToFile(inputStream, filePath, fileName); // System.out.println(file.getAbsolutePath); // 输出结果:例如,将文件流转换为File对象后,文件的绝对路径为:C:\Users\User\Documents\example.txt // 则输出结果为:C:\Users\User\Documents\example.txt 通过分析,不难发现以上代码:
-
Apache Spark 的基本概念和在大数据分析中的应用
-
python 数据分析在电力行业中的应用 python 电影数据分析报告
-
正负偏差变量 即 d2+、d2- 分别表示决策值中超出和未达到目标值的部分。而 di+、di- 均大于 0 刚性约束和目标约束(柔性目标约束有偏差) 在多目标规划中,>=/<= 在刚性约束中保持不变。当需要将约束条件转换为柔性约束条件时,需要将 >=/<= 更改为 =(因为已经有 d2+、d2- 用来表示正负偏差),并附加上 (+dii-di+) 注意这里是 +di、-di+!之所以是 +di,-di+,是因为需要将目标还原为最接近的原始刚性约束条件 优先级因素和权重因素 对多个目标进行优先排序和优先排序 目标规划的目标函数 是所有偏差变量的加权和。值得注意的是,这个加权和都取最小值。而 di+ 和 dii- 并不一定要出现在每个不同的需求层次中。具体分析需要具体问题具体分析 下面是一个例子: 题目中说设备 B 既要求充分利用,又要求尽可能不加班,那么列出的时间计量表达式即为:min z = P3 (d3- + d3 +) 使用 + 而不是 -d3 + 的原因是:正负偏差不可能同时存在,必须有 di+di=0 (因为判定值不可能同时大于目标值和小于目标值),而前面是 min,所以只要取 + 并让 di+ 和 dii- 都为正值即可。因此,得出以下规则: 最后,给出示例和相应的解法: 问题:某企业生产 A 和 B 两种产品,需要使用 A、B、C 三种设备。下表显示了与工时和设备使用限制有关的产品利润率。问该企业应如何组织生产以实现下列目标? (1) 力争利润目标不低于 1 500 美元; (2) 考虑到市场需求,A、B 两种产品的生产比例应尽量保持在 1:2; (3)设备 A 是贵重设备,严禁超时使用; (4)设备 C 可以适当加班,但要控制;设备 B 要求充分利用,但尽量不加班。 从重要性来看,设备 B 的重要性是设备 C 的三倍。 建立相应的目标规划模型并求解。 解:设企业生产 A、B 两种产品的件数分别为 x1、x2,并建立相应的目标计划模型: 以下为顺序求解法,利用 LINGO 求解: 1 级目标: 模型。 设置。 variable/1..2/:x;! s_con_num/1...4/:g,dplus,dminus;!所需软约束数量(g=dplus=dminus 数量)及相关参数; s_con(s_con_num);! s_con(s_con_num,variable):c;!软约束系数; 结束集 数据。 g=1500 0 16 15. c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(1);!第一个目标函数;!对应于 min=z 的第一小部分;! 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); !使用设置完成的数据构建软约束表达式; ! !软约束表达式 @for(variable:@gin(x)); !将变量约束为整数; ! 结束 此时,第一级目标的最优值为 0,第一级偏差为 0: 第二级目标: !求 dminus(1)=0,然后求解第二级目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(2)+dplus(2);!第二个目标函数 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标结果 @for(variable:@gin(x)); ! 结束 此时,第二个目标的最优值为 0,偏差为 0: 第三目标 !求 dminus(2)=0,然后求解第三个目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=3*dminus(3)+3*dplus(3)+dminus(4);!第三个目标函数。 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标约束条件; ! dminus(2)+dplus(2)=0; !第二个目标约束条件 @for(variable:@gin(x));! 结束 最终结果为 x1=2,x2=4,dplus(1)=100,最优利润为
-
一张图告诉你网站运营中的关键指标 | I. 网站运营内容指标数据分析
-
人工智能技术在摄影设备中的应用及分析实例
-
Apache Spark 的基本概念和在大数据分析中的应用
-
网络在单细胞转录组数据分析中的应用