刘小泽写于2020.5.18
上次看了热图如何去掉聚类树的同时保留聚类的顺序 ，有朋友留言给出了其他的解决方法，非常感谢！文中会给出答案
另外，这次还是带着数据，去看看pheatmap其他的有用操作
上数据！
test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
# 20行（基因）10列（样本）
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")
# 和上次不同的是，我们这次设置一个极值（为了后面演示scale）
test[10,10]=100
> rownames(test)[10];colnames(test)[10]
[1] "Gene10"
[1] "Test10"
上数据
画图
第一个原始图
可以看到Gene10 Test10这里的100大大超过了一般水平，导致画的图受它的影响极大，变得看不出整体的变化趋势
pheatmap(test)
第一个原始图
第二个：增加z-score参数的图
设置一个参数scale = 'row' 就会对行进行z-score标准化
什么是z-score？ 看个公式就明白了(表达量-均值)/标准差，目的就是不让数据受极值的影响太大。另外z-score一般还会设置一个最大最小值（例如
-2 ~ 2、-1 ~ 1）
Z score = (expression G - mean expression G1..Gn)/SD G1..Gn where G is any gene in the pathway and G1..Gn represent the aggregate measure of all of the genes.
Cite：Cheadle C, Vawter MP, Freed WJ, Becker KG. Analysis of microarray data using Z score transformation. J Mol Diagn. 2003;5(2):73‐81. doi:10.1016/S1525-1578(10)60455-2
scale做了什么
pheatmap(test,scale = 'row')
增加z-score参数的图
第三个图：手动zscore
如果想要理解某个参数做了什么，就自己模拟这个参数吧
我们这里也来模拟一下这个参数做了什么【并且用热心朋友推荐的方法去掉聚类树】
cal_z_score <- function(x){
(x - mean(x)) / sd(x)
}
diy1 <- t(apply(test, 1, cal_z_score))
# 很巧妙，设置高度为0，就不显示了
pheatmap(diy1,treeheight_row = 0,treeheight_col = 0)
手动zscore
第四个图：聚类是怎么聚的呢？
使用参数的话，很简单：
pheatmap(diy1,clustering_method = 'complete',cluster_rows = T, cluster_cols = T,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean")
这个method是直接使用的hlust计算的，全部聚类方法包括：
"ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
距离使用dist计算，默认使用欧式距离，全部距离包括：
"euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski"
提取聚类结果：
p1 <- pheatmap(diy1, silent = TRUE) # 只输出结果，不画图
> names(p1)
[1] "tree_row" "tree_col" "kmeans" "gtable"
# 做个丑陋的聚类图看看
library(dendextend)
p1$tree_row %>%
as.dendrogram() %>%
plot(horiz = TRUE)
# 聚类基因名
> rownames(diy1)[p1$tree_row[["order"]]]
[1] "Gene1" "Gene3" "Gene2" "Gene8" "Gene4" "Gene6" "Gene9" "Gene5" "Gene7" "Gene10" "Gene13" "Gene17"
[13] "Gene14" "Gene11" "Gene16" "Gene20" "Gene12" "Gene19" "Gene15" "Gene18"
自己聚类是这样的：
my_hclust_gene <- hclust(dist(diy1), method = "complete")
as.dendrogram(my_hclust_gene) %>%
plot(horiz = TRUE)
自己聚类是这样的
第五个图：把树切分
pheatmap(diy1,cutree_rows = 2,cutree_cols = 2)
把树切分
第六张图：自己给切分后的基因上色
# cutree_rows = 2就做了下面这行代码的事情
gene_anno <- cutree(tree = as.dendrogram(my_hclust_gene), k = 2)
gene_anno <- data.frame(cluster = ifelse(gene_anno == 1, "cluster1", "cluster2"))
> table(gene_anno)
gene_anno
cluster 1 cluster 2
9 11
pheatmap(diy1, annotation_row = gene_anno)
自己给切分后的基因上色
修改颜色
需要设置一个list，如果有列的注释信息，也是添加到这个my_colour的list中
my_colour = list(
cluster = c(cluster1 = "#e89829", cluster2 = "#cc4ee0")
)
pheatmap(diy1,
annotation_colors = my_colour,
annotation_row = gene_anno,
treeheight_row = 0)
修改颜色
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论