R語言學習 - 熱圖美化

2021-02-21 生信寶典
熱圖美化

上一期的繪圖命令中，最後一行的操作抹去了之前設定的橫軸標記的旋轉，最後出來的圖比較難看。上次我們是這麼寫的

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_5a   5.958073   5.843652   3.225465   4.886184     3.411362b  19.630582  20.376791  20.744580  18.534027    20.638288c 100.351299  99.849900 102.197343  98.583629    99.540488d 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 IDa   5.958073   5.843652   3.225465   4.886184     3.411362  ab  19.630582  20.376791  20.744580  18.534027    20.638288  bc 100.351299  99.849900 102.197343  98.583629    99.540488  cd 600.000000 700.000000 800.000000 900.000000 10000.000000  d

data_m <- melt(data, id.vars=c("ID"))
head(data_m)

 ID variable      value1  a    Grp_1   5.9580732  b    Grp_1  19.6305823  c    Grp_1 100.3512994  d    Grp_1 600.0000005  a    Grp_2   5.8436526  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_5a   6.61047  20.946720 100.133106 600.000000     5.267921b  20.80792  99.865962 700.000000   3.737228    19.289715c 100.06930 800.000000   6.252753  21.464081    98.607518d 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_5a 2.927986 4.455933 6.660112 9.231221  2.647987b 4.446780 6.656296 9.453271 2.244043  4.342677c 6.659201 9.645658 2.858529 4.489548  6.638183d 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="")


data <- data[apply(data,1,var)!=0,]

data

 Grp_1 Grp_2 Grp_3 Grp_4   Grp_5a   6.6  20.9 100.1 600.0     5.2b  20.8  99.8 700.0   3.7    19.2c 100.0 800.0   6.2  21.4    98.6d 900.0   3.3  20.3 101.1 10000.0


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_5a -0.5456953 -0.4899405 -0.1811446  1.7679341 -0.5511538b -0.4940465 -0.2301542  1.7747592 -0.5511674 -0.4993911c -0.3139042  1.7740182 -0.5936858 -0.5483481 -0.3180801d -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_5a   6.6  20.9 100.0 100.0   5.2b  20.8  99.8 100.0   3.7  19.2c 100.0 100.0   6.2  21.4  98.6d 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_5a   6.6  20.9 100.1 600.0     5.2b  20.8  99.8 700.0   3.7    19.2c 100.0 800.0   6.2  21.4    98.6d 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


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   value1   a    Grp_1   8.3012   b    Grp_1   24.843   c    Grp_1  101.454   d    Grp_1 2794.355   a    Grp_2   24.846   b    Grp_2  101.457   c    Grp_2 2794.358   d    Grp_2   5.7189   a    Grp_3  101.4510  b    Grp_3 2794.3511  c    Grp_3   8.30112  d    Grp_3   24.8413  a    Grp_4 2794.3514  b    Grp_4   5.71815  c    Grp_4   24.8416  d    Grp_4  101.4517  a    Grp_5   5.71818  b    Grp_5   24.8419  c    Grp_5  101.4520  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))


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繪製到現在就差不多了，但總是這麼畫下去也會覺得有點累，有沒有辦法更簡化呢？ 且聽下回分解。

相關焦點

  • 我的R語言學習方法
    我為什麼要學習R語言？利用R語言做數據分析的工作；我怎麼學習R語言？在快速學習R語言基礎後，採用邊學邊做，不斷強化的策略學習和應用R語言；我學習R語言的什麼，我主要學習了R語言做數據整理，數據分析、數據建模和數據可視化這四方面的內容，並且對幫助我有效完成數據工作非常有幫助。第二點，我談一下自己在學習R語言的具體方法和做法。
  • R語言學習 - 熱圖美化 (數值標準化和調整坐標軸順序)
    溫故知新熱圖美化上一期的繪圖命令中，最後一行的操作抹去了之前設定的橫軸標記的旋轉，最後出來的圖比較難看。
  • 【R函數學習】R語言時間序列函數整理
    d=scan("sha.csv")sha=ts(d,start=1964,freq=1)plot.ts(sha) #繪製時序圖acf(sha,22) #繪製自相關圖，滯後期數22pacf(sha,22) #繪製偏自相關圖，滯後期數22corr=acf(sha,22) #保存相關係數cov=acf(sha,22,type = "covariance
  • 【C語言】什麼是 lvalue,什麼是 rvalue?
    隨著學習的深入，你會經常看到 lvalue 這個詞。一般出現在諸如各種、書籍中，更頻繁遇到的是在你的錯誤提示中：一般你會看到這個是因為你的代碼這麼寫：……        int i;        5 = i;……一些朋友想當然的就覺得 lvalue 指的就是賦值運算符左邊的那個值，而 rvalue 當然就是右邊那個值啦。
  • PPT中柱形圖與餅圖的美化方法
    但凡涉及數據處理，柱形圖與餅圖幾乎是使用最多的數據展示工具，其導入及使用的方法相對簡單，很多人相對欠缺的是對其進行適當的美化，使之簡約大方、突出重點
  • 教您用R語言畫直方圖!| 繽紛燦爛R語言 | 醫學方R語言高階課程
    橘黃色表示死亡，藍色表示生存，當然，如前所述，這張圖還缺少一個圖例。：接受一個數值向量或者一個矩陣；當為一個向量時，則該向量用來繪製左半邊的條形圖或者直方圖；當為一個矩陣時，該矩陣的第一列數據用來繪製左半邊的圖，第二列用來繪製右半邊的圖，而其他列（假如存在其他列）則被自動忽略了。
  • 如何學習好R語言?【全套R語音書籍+視頻下載】
    包括了閱讀經典的教材、代碼、論文、學習公開課。 - 通過牛人來學習。 包括同行的聚會、討論、大牛的博客、微博、twitter、RSS。 - 通過練習來學習。 包括代碼練習題、參加kaggle比賽、解決實際工作中的難題。 - 通過分享來學習。 包括自己寫筆記、寫博客、寫書、翻譯書，和同伴分享交流、培訓新人。# 全套R語音書籍下載。關注公眾號後，回覆：R語言
  • R語言時間序列函數大全(收藏!)
    acf(sha,22) #繪製自相關圖，滯後期數22pacf(sha,22) #繪製偏自相關圖，滯後期數22corr=acf(sha,22) #保存相關係數cov=acf(sha,22,type = 「covariance」) #保存協方差2、同時繪製兩組數據的時序圖d=read.csv(「double.csv
  • 如何入門R語言
    學歷歷程一開始，接觸R語言我也懵逼，畢竟除了小學學過的VB和大二為了完成線代作業學了一點點的Matlab，我就一點編程基礎都沒有了。那時候百度了挺多資料，發現都比較零散不成系統。好在最後買了《R語言實戰》，在這本書的幫助下才入了門。而後又陸續看了《R語言編程藝術》，《統計學習導論》等書，學習總算是越來越有感覺。
  • 跟著Nature microbiology學畫圖~R語言畫韋恩圖然後拼圖
    R語言裡比較常用的畫韋恩圖的包是VennDiagram，但是今天的內容涉及到拼圖，用VennDiagram畫圖後如何拼圖我暫時還不知道。image.png這裡遇到一個問題是文字標籤超出邊界了，如何調節文字標籤的位置暫時還不知道如何通過代碼實現，只能出圖後手動編輯了。
  • 跟著Nature microbiology學畫圖~R語言ggplot2畫散點圖
    期刊是 Nature microbiologyimage.png今天重複的圖片是Figure2中的氣泡圖read.csv("example_data/bubble_plot_example.csv",             header=T)df讀入的數據是寬格式，ggplot2作圖需要用長格式數據，對寬格式數據進行轉化df1<-reshape2::melt(df)df1最基本的散點圖
  • Excel美化二維折線圖!
    美化後的Excel二維折線圖：
  • 最全的 R 語言學習路線在這裡,讓你少走彎路!
    現在對R感興趣的人越來越多，很多人都想快速的掌握R語言，然而，由於目前大部分高校都沒有開設R語言課程，這就導致很多人不知道如何著手學習
  • C語言學習基礎必會
    C語言學習目標：將C語言簡單概述一下我將把我學到的知識和自己的理解以自己的話說出來。可能有很多問題邏輯也不是那麼的完整。
  • C語言中的奇技淫巧
    前言學習C語言的過程中，總會遇到很多令人眼前一亮的代碼，尤其是你寫了幾十行的代碼，別人只用了簡單幾行的遞歸就實現的功能
  • R語言讀取excel數據時的一點小建議
    比如在R語言中，讀取數據後，head或者tail一下，初步看下數據是否正確讀入。事情的原由是，群裡有位剛學R語言的小夥伴問道，他想畫相關性圖，但是不知為何數據中一直有NA，也就是我們常說的缺失值。如下圖：隨後，我讓他把數據發給我看下，我讀取數據後發現如下：
  • C語言 | 求圓周長 面積 表面積 體積
    「要成為絕世高手，並非一朝一夕，除非是天生武學奇才，但是這種人…萬中無一」這道理放在C語言學習上也一併受用
  • 不一定正確的多分組差異分析結果熱圖展現
    <- paste0(colnames(group_list)[i],colnames(tT_tmp))  deg_re <- cbind(deg_re,tT_tmp)  }deg_re <- deg_re[,-1]save(deg_re,file='deg.Rdata')Step04-pheatmap 文章是用logFC做熱圖
  • 《R語言實戰》菜雞筆記(八):回歸
    Income + Frost, data=data_1)qqPlot(fit_4, labels=row.names(data_1), id.method="identify", simulate=TRUE, main="Q-Q Plot") data_1["Nevada",]fitted(fit_4)["Nevada" ] residuals(fit_4)["Nevada"] rstudent
  • R語言讀取xlsx文件
    關於R語言讀取Excel文件，比較麻煩，我從來都反對直接讀取xlsx文件，因為爬蟲數據時，一般保存的格式都是csv文件，或者直接保存到資料庫裡面