劉小澤寫于2020.5.18
上次看了熱圖如何去掉聚類樹的同時保留聚類的順序 ，有朋友留言給出了其他的解決方法，非常感謝！文中會給出答案
另外，這次還是帶著數(shù)據(jù)，去看看pheatmap其他的有用操作

上數(shù)據(jù)！

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 = "")

# 和上次不同的是，我們這次設(shè)置一個極值（為了后面演示scale）
test[10,10]=100
> rownames(test)[10];colnames(test)[10]
[1] "Gene10"
[1] "Test10"

上數(shù)據(jù)

畫圖

第一個原始圖

可以看到Gene10 Test10這里的100大大超過了一般水平，導(dǎo)致畫的圖受它的影響極大，變得看不出整體的變化趨勢

pheatmap(test)

第一個原始圖

第二個：增加z-score參數(shù)的圖

設(shè)置一個參數(shù)scale = 'row' 就會對行進(jìn)行z-score標(biāo)準(zhǔn)化

什么是z-score？ 看個公式就明白了(表達(dá)量-均值)/標(biāo)準(zhǔn)差，目的就是不讓數(shù)據(jù)受極值的影響太大。另外z-score一般還會設(shè)置一個最大最小值（例如-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參數(shù)的圖

第三個圖：手動zscore

如果想要理解某個參數(shù)做了什么，就自己模擬這個參數(shù)吧

我們這里也來模擬一下這個參數(shù)做了什么【并且用熱心朋友推薦的方法去掉聚類樹】

cal_z_score <- function(x){
    (x - mean(x)) / sd(x)
}
diy1 <- t(apply(test, 1, cal_z_score))
# 很巧妙，設(shè)置高度為0，就不顯示了
pheatmap(diy1,treeheight_row = 0,treeheight_col = 0)

手動zscore

第四個圖：聚類是怎么聚的呢？

使用參數(shù)的話，很簡單：

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計算，默認(rèn)使用歐式距離，全部距離包括：

"euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski"

提取聚類結(jié)果：

p1 <- pheatmap(diy1, silent = TRUE) # 只輸出結(jié)果，不畫圖
> 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)

自己給切分后的基因上色

修改顏色

需要設(shè)置一個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)

修改顏色

歡迎關(guān)注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!