歡迎關注 生信寶典 公眾號，閱讀系列文章

http://mp.weixin.qq.com/s/lKrhvYrwn93esC6MA3bHWw

Rstudio基礎

R語言是比較常用的統計分析和繪圖語言，擁有強大的統計庫、繪圖庫和生信分析的Bioconductor庫，是學習生物信息分析的必備語言之一。

Rstudio是編輯、運行R語言的最為理想的工具之一，支持純R腳本、Rmarkdown (腳本文檔混排)、Bookdown (腳本文檔混排成書)、Shiny (交互式網絡應用)等。

Rstudio版本

Rsdutio分為桌面版和服務器版，桌面版可以在單機使用，服務器版可以從瀏覽器訪問供多人使用。

服務器版安裝好之后，訪問地址為<服務器IP:8787>，用戶名和密碼為Linux用戶的用戶名和密碼。

Rstudio安裝

R安裝

Linux下安裝

Rstudio安裝前需要安裝R，如果使用的是新版的操作系統。直接可以用sudo apt-get installl r-base 或者yum install r-base來安裝。

若系統版本老，或沒有根用戶權限，則需要下載編譯源碼安裝，最新地址為https://cran.r-project.org/src/base/R-3/R-3.4.0.tar.gz。

具體編譯方式為 (Linux下軟件安裝見 http://blog.genesino.com/2016/06/bash1):

# --enable-R-shlib 需要設置，使得其他程序包括Rstudio可以使用R的動態庫 # --prefix指定軟件安裝目錄，需使用絕對路徑 ./configure --prefix=/home/ehbio/R/3.4.0 --enable-R-shlib# 也可以使用這個命令，共享系統的blas庫，提高運輸速度 #./configure --prefix=/home/ehbio/R/3.4.0 --enable-R-shlib --with-blas --with-lapack make make install

Windows下安裝

下載https://cran.r-project.org/bin/windows/雙擊就可以了。

Rstudio安裝

Linux下安裝服務器版

• 安裝參考 https://www.rstudio.com/products/rstudio/download-server/

  wget https://download2.rstudio.org/rstudio-server-rhel-1.0.136-x86_64.rpm sudo yum install --nogpgcheck rstudio-server-rhel-1.0.136-x86_64.rpm

• 安裝完之后的檢測、啟動和配置

  sudo rstudio-server verify-installation #查看是否安裝正確 sudo rstudio-server start ## 啟動 sudo rstudio-server status ## 查看狀態 sudo rstudio-server stop ## 停止 ifconfig | grep 'inet addr' ## 查看服務端ip地址 sudo rstudio-server start ## 修改配置文件后重啟 sudo rstudio-server active-sessions ## 列出活躍的sessions: sudo rstudio-server suspend-session <pid> ## 暫停session sudo rstudio-server suspend-all ##暫停所有session

• Rstudio日志目錄，方便查看錯誤信息：/var/log/rstudio-server/

• 配置文件：

  • /etc/rstudio/rserver.conf
  www-port=8787 (default) www-address=0.0.0.0 (default) rsession-ld-library-path=/opt/local/lib:/opt/local/someapp/lib rsession-which-r=/usr/local/bin/R
  • /etc/rstudio/rsession.conf

  • Timeout

    [user] session-timeout-minutes=30 [@powerusers] session-timeout-minutes=0

Windows下安裝桌面版

下載之后 (https://www.rstudio.com/products/rstudio/download2/)雙擊安裝，需要使用管理員權限，其它無需要注意的。

Rstudio 使用

Windows下桌面版直接雙擊打開即可使用，Linux服務器版訪問地址為<服務器IP:8787>，用戶名和密碼為Linux用戶的用戶名和密碼。

Rstudio 界面

Rstudio中新建或打開文件

如果是桌面版，直接就可以訪問”我的電腦”去打開之前寫過的腳本。如果是服務器版，可直接訪問服務器上寫過的腳本。Rstudio右下1/4部分可以切換目錄，點擊more，設置工作目錄。可以上傳本地的腳本到對應目錄打開。

Rstudio還有其它的功能，不過這些對我們初學使用已經足夠了。后續會根據需要推出基于ggplot2作圖的R入門介紹。

熱圖繪制

熱圖是做分析時常用的展示方式，簡單、直觀、清晰。可以用來顯示基因在不同樣品中表達的高低、表觀修飾水平的高低等。任何一個數值矩陣都可以通過合適的方式用熱圖展示。

本篇使用R的ggplot2包實現從原始數據讀入到熱圖輸出的過程，并在教程結束后提供一份封裝好的命令行繪圖工具，只需要提供矩陣，即可一鍵繪圖。

上一篇講述了Rstudio的使用作為R寫作和編譯環境的入門，后面的命令都可以拷貝到Rstudio中運行，或寫成一個R腳本，使用Rscript heatmap.r運行。我們還提供了Bash的封裝，在不修改R腳本的情況下，改變參數繪制出不同的圖形。

生成測試數據

繪圖首先需要數據。通過生成一堆的向量，轉換為矩陣，得到想要的數據。

data <- c(1:6,6:1,6:1,1:6, (6:1)/10,(1:6)/10,(1:6)/10,(6:1)/10,1:6,6:1,6:1,1:6, 6:1,1:6,1:6,6:1) [1] 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 6.0 5.0 4.0 3.0 2.0 1.0 1.0 [20] 2.0 3.0 4.0 5.0 6.0 0.6 0.5 0.4 0.3 0.2 0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 [39] 0.3 0.4 0.5 0.6 0.6 0.5 0.4 0.3 0.2 0.1 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 [58] 3.0 2.0 1.0 6.0 5.0 4.0 3.0 2.0 1.0 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 [77] 2.0 1.0 1.0 2.0 3.0 4.0 5.0 6.0 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 2.0 [96] 1.0

注意：運算符的優先級。

> 1:3+4 [1] 5 6 7 > (1:3)+4 [1] 5 6 7 > 1:(3+4) [1] 1 2 3 4 5 6 7

Vector轉為矩陣 (matrix)，再轉為數據框 (data.frame)。

# ncol: 指定列數 # byrow: 先按行填充數據 # ?matrix 可查看函數的使用方法 # as.data.frame的as系列是轉換用的 data <- as.data.frame(matrix(data, ncol=12, byrow=T)) V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 1 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 2 6.0 5.0 4.0 3.0 2.0 1.0 1.0 2.0 3.0 4.0 5.0 6.0 3 0.6 0.5 0.4 0.3 0.2 0.1 0.1 0.2 0.3 0.4 0.5 0.6 4 0.1 0.2 0.3 0.4 0.5 0.6 0.6 0.5 0.4 0.3 0.2 0.1 5 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 6 6.0 5.0 4.0 3.0 2.0 1.0 1.0 2.0 3.0 4.0 5.0 6.0 7 6.0 5.0 4.0 3.0 2.0 1.0 1.0 2.0 3.0 4.0 5.0 6.0 8 1.0 2.0 3.0 4.0 5.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 # 增加列的名字 colnames(data) <- c("Zygote","2_cell","4_cell","8_cell","Morula","ICM","ESC","4 week PGC","7 week PGC","10 week PGC","17 week PGC", "OOcyte")# 增加行的名字 # 注意paste和paste0的使用 rownames(data) <- paste("Gene", 1:8, sep="_")# 只顯示前6行和前4列 head(data)[,1:4] Zygote 2_cell 4_cell 8_cell Gene_1 1.0 2.0 3.0 4.0 Gene_2 6.0 5.0 4.0 3.0 Gene_3 0.6 0.5 0.4 0.3 Gene_4 0.1 0.2 0.3 0.4 Gene_5 1.0 2.0 3.0 4.0 Gene_6 6.0 5.0 4.0 3.0

雖然方法比較繁瑣，但一個數值矩陣已經獲得了。

還有另外2種獲取數值矩陣的方式。

• 讀入字符串

# 使用字符串的好處是不需要額外提供文件 # 簡單測試時可使用，寫起來不繁瑣，又方便重復 # 尤其適用于在線提問時作為測試案例 > txt <- "ID;Zygote;2_cell;4_cell;8_cell + Gene_1;1;2;3;4 + Gene_2;6;5;4;5 + Gene_3;0.6;0.5;0.4;0.4"# 習慣設置quote為空，避免部分基因名字或注釋中存在引號，導致讀入文件錯誤。 # 具體錯誤可查看 http://blog.genesino.com/collections/R_tips/ 中的記錄 > data2 <- read.table(text=txt,sep=";", header=T, row.names=1, quote="") > head(data2)Zygote X2_cell X4_cell X8_cell Gene_1 1.0 2.0 3.0 4.0 Gene_2 6.0 5.0 4.0 5.0 Gene_3 0.6 0.5 0.4 0.4

可以看到列名字中以數字開頭的列都加了X。一般要盡量避免行或列名字以數字開頭，會給后續分析帶去一些困難；另外名字中出現的非字母、數字、下劃線、點的字符都會被轉為點，也需要注意，盡量只用字母、下劃線和數字。

# 讀入時，增加一個參數`check.names=F`也可以解決問題。 # 這次數字前沒有再加 X 了 > data2 <- read.table(text=txt,sep=";", header=T, row.names=1, quote="", check.names = F) > head(data2)Zygote 2_cell 4_cell 8_cell Gene_1 1.0 2.0 3.0 4.0 Gene_2 6.0 5.0 4.0 5.0 Gene_3 0.6 0.5 0.4 0.4

• 讀入文件

與上一步類似，只是改為文件名，不再贅述。

> data2 <- read.table("filename",sep=";", header=T, row.names=1, quote="")

轉換數據格式

數據讀入后，還需要一步格式轉換。在使用ggplot2作圖時，有一種長表格模式是最為常用的，尤其是數據不規則時，更應該使用 （這點，我們在講解箱線圖時再說）。

# 如果包沒有安裝，運行下面一句，安裝包 #install.packages(c("reshape2","ggplot2"))library(reshape2) library(ggplot2)# 轉換前，先增加一列ID列，保存行名字 data$ID <- rownames(data)# melt：把正常矩陣轉換為長表格模式的函數。工作原理是把全部的非id列的數值列轉為1列，命名為value；所有字符列轉為variable列。# id.vars 列用于指定哪些列為id列；這些列不會被merge，會保留為完整一列。 data_m <- melt(data, id.vars=c("ID")) head(data_m) ID variable value 1 Gene_1 Zygote 1.0 2 Gene_2 Zygote 6.0 3 Gene_3 Zygote 0.6 4 Gene_4 Zygote 0.1 5 Gene_5 Zygote 1.0 6 Gene_6 Zygote 6.0 7 Gene_7 Zygote 6.0 8 Gene_8 Zygote 1.0 9 Gene_1 2_cell 2.0 10 Gene_2 2_cell 5.0 11 Gene_3 2_cell 0.5 12 Gene_4 2_cell 0.2 13 Gene_5 2_cell 2.0 14 Gene_6 2_cell 5.0 15 Gene_7 2_cell 5.0 16 Gene_8 2_cell 2.0

分解繪圖

數據轉換后就可以畫圖了，分解命令如下：

# data_m: 是前面費了九牛二虎之力得到的數據表 # aes: aesthetic的縮寫，一般指定整體的X軸、Y軸、顏色、形狀、大小等。 # 在最開始讀入數據時，一般只指定x和y，其它后續指定 p <- ggplot(data_m, aes(x=variable,y=ID)) # 熱圖就是一堆方塊根據其值賦予不同的顏色，所以這里使用fill=value, 用數值做填充色。 p <- p + geom_tile(aes(fill=value)) # ggplot2為圖層繪制，一層層添加，存儲在p中，在輸出p的內容時才會出圖。 p## 如果你沒有使用Rstudio或其它R圖形版工具，而是在遠程登錄的服務器上運行的交互式R，需要輸入下面的語句，獲得輸出圖形 （圖形存儲于R的工作目錄下的Rplots.pdf文件中）。## 如何指定輸出，后面會講到。 #dev.off()

熱圖出來了，但有點不對勁，橫軸重疊一起了。一個辦法是調整圖像的寬度，另一個是旋轉橫軸標記。

# theme: 是處理圖美觀的一個函數，可以調整橫縱軸label的選擇、圖例的位置等。 # 這里選擇X軸標簽45度。 # hjust和vjust調整標簽的相對位置，具體見 <https://stackoverflow.com/questions/7263849/what-do-hjust-and-vjust-do-when-making-a-plot-using-ggplot>。 # 簡單說，hjust是水平的對齊方式，0為左，1為右，0.5居中，0-1之間可以取任意值。vjust是垂直對齊方式，0底對齊，1為頂對齊，0.5居中，0-1之間可以取任意值。p <- p + theme(axis.text.x=element_text(angle=45,hjust=1, vjust=1)) p

設置想要的顏色。

# 連續的數字，指定最小數值代表的顏色和最大數值賦予的顏色 # 注意fill和color的區別，fill是填充，color只針對邊緣 p <- p + scale_fill_gradient(low = "white", high = "red") p

調整legend的位置。

# postion可以接受的值有 top, bottom, left, right, 和一個坐標 c(0.05,0.8) (左上角，坐標是相對于圖的左下角計算的） p <- p + theme(legend.position="top")

調整背景和背景格線以及X軸、Y軸的標題。

p <- p + xlab("samples") + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank()) p

合并以上命令，就得到了下面這個看似復雜的繪圖命令。

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

圖形存儲

圖形出來了，就得考慮存儲了，

# 可以跟輸出文件不同的后綴，以獲得不同的輸出格式 # colormode支持srgb （屏幕）和cmyk （打印，部分雜志需要，看上去有點褪色的感覺）格式 ggsave(p, filename="heatmap.pdf", width=10,height=15, units=c("cm"),colormodel="srgb")

至此，完成了簡單的heatmap的繪圖。但實際繪制時，經常會碰到由于數值變化很大，導致顏色過于集中，使得圖的可讀性下降很多。因此需要對數據進行一些處理，具體的下次再說。

R基礎概念和困惑點

插播一點R的基礎概念和疑難解釋。

R的交互式運行

在bash命令行下輸入大寫字母R即可啟動交互式界面

ct@ehbio:~$ RR version 3.3.1 (2016-06-21) -- "Bug in Your Hair" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-redhat-linux-gnu (64-bit)用'help.start()'通過HTML瀏覽器來看幫助文件。 用'q()'退出R.> 1:3 [1] 1 2 3 > a <- 1:3 > a [1] 1 2 3 # 使用source在交互式界面運行寫好的R腳本 > source('scrtpt.r') > quit() Save workspace image? [y/n/c]: n# ctrl+d也可退出R

R基本語法

獲取幫助文檔，查看命令或函數的使用方法、事例或適用范圍

>>> ?command >>> ??command #深度搜索或模糊搜索此命令>>> example(command) #得到命令的例子

R中的變量

> # 數字變量 > a <- 10 > a [1] 10 > > # 字符串變量 > a <- "abc" > a [1] "abc" > > # 邏輯變量 > a <- TRUE > > a [1] TRUE > > b <- T > > b [1] TRUE > > d <- FALSE > > d [1] FALSE > # 向量 > > a <- vector(mode="logical", length=5) > a [1] FALSE FALSE FALSE FALSE FALSE > > a <- c(1,2,3,4) # 判斷一個變量是不是vector > is.vector(a) [1] TRUE > > # 矩陣 > > a <- matrix(1:20,nrow=5,ncol=4,byrow=T) > a[,1] [,2] [,3] [,4] [1,] 1 2 3 4 [2,] 5 6 7 8 [3,] 9 10 11 12 [4,] 13 14 15 16 [5,] 17 18 19 20 > > is.matrix(a) [1] TRUE > > dim(a) #查看或設置數組的維度向量 [1] 5 4 > > # 錯誤的用法 > dim(a) <- c(4,4) Error in dim(a) <- c(4, 4) : dims [product 16]與對象長度[20]不匹配 > > # 正確的用法 > a <- 1:20 > dim(a) <- c(5,4) #轉換向量為矩陣 > a[,1] [,2] [,3] [,4] [1,] 1 6 11 16 [2,] 2 7 12 17 [3,] 3 8 13 18 [4,] 4 9 14 19 [5,] 5 10 15 20 > > print(paste("矩陣a的行數", nrow(a))) [1] "矩陣a的行數 5" > print(paste("矩陣a的列數", ncol(a))) [1] "矩陣a的列數 4" > > #查看或設置行列名 > rownames(a) NULL > rownames(a) <- c('a','b','c','d','e') > a[,1] [,2] [,3] [,4] a 1 6 11 16 b 2 7 12 17 c 3 8 13 18 d 4 9 14 19 e 5 10 15 20# R中獲取一系列的字母 > letters[1:4] [1] "a" "b" "c" "d" > colnames(a) <- letters[1:4] > aa b c d a 1 6 11 16 b 2 7 12 17 c 3 8 13 18 d 4 9 14 19 e 5 10 15 20 > # is系列和as系列函數用來判斷變量的屬性和轉換變量的屬性 # 矩陣轉換為data.frame > is.character(a) [1] FALSE > is.numeric(a) [1] TRUE > is.matrix(a) [1] TRUE > is.data.frame(a) [1] FALSE > is.data.frame(as.data.frame(a)) [1] TRUE

R中矩陣運算

# 數據產生 # rnorm(n, mean = 0, sd = 1) 正態分布的隨機數 # runif(n, min = 0, max = 1) 平均分布的隨機數 # rep(1,5) 把1重復5次 # scale(1:5) 標準化數據 > a <- c(rnorm(5), rnorm(5,1), runif(5), runif(5,-1,1), 1:5, rep(0,5), c(2,10,11,13,4), scale(1:5)[1:5]) > a[1] -0.41253556 0.12192929 -0.47635888 -0.97171653 1.09162243 1.87789657[7] -0.11717937 2.92953522 1.33836620 -0.03269026 0.87540920 0.13005744 [13] 0.11900686 0.76663940 0.28407356 -0.91251181 0.17997973 0.50452258 [19] 0.25961316 -0.58052230 1.00000000 2.00000000 3.00000000 4.00000000 [25] 5.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 [31] 2.00000000 10.00000000 11.00000000 13.00000000 4.00000000 -1.26491106 [37] -0.63245553 0.00000000 0.63245553 1.26491106 > a <- matrix(a, ncol=5, byrow=T) > a[,1] [,2] [,3] [,4] [,5] [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 [3,] 0.8754092 0.1300574 0.1190069 0.7666394 0.28407356 [4,] -0.9125118 0.1799797 0.5045226 0.2596132 -0.58052230 [5,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.00000000 [7,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 [8,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106# 求行的加和 > rowSums(a) [1] -0.6470593 5.9959284 2.1751865 -0.5489186 15.0000000 0.0000000 40.0000000 [8] 0.0000000## 注意檢查括號的配對 > a <- a[rowSums(abs(a)!=0,] 錯誤: 意外的']' in "a <- a[rowSums(abs(a)!=0,]"# 去除全部為0的行 > a <- a[rowSums(abs(a))!=0,]# 另外一種方式去除全部為0的行 > #a[rowSums(a==0)<ncol(a),] > a[,1] [,2] [,3] [,4] [,5] [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 [3,] 0.8754092 0.1300574 0.1190069 0.7666394 0.28407356 [4,] -0.9125118 0.1799797 0.5045226 0.2596132 -0.58052230 [5,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 [6,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 [7,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106# 矩陣運算，R默認針對整個數據進行常見運算# 所有值都乘以2> a * 2[,1] [,2] [,3] [,4] [,5] [1,] -0.8250711 0.2438586 -0.9527178 -1.9434331 2.18324487 [2,] 3.7557931 -0.2343587 5.8590704 2.6767324 -0.06538051 [3,] 1.7508184 0.2601149 0.2380137 1.5332788 0.56814712 [4,] -1.8250236 0.3599595 1.0090452 0.5192263 -1.16104460 [5,] 2.0000000 4.0000000 6.0000000 8.0000000 10.00000000 [6,] 4.0000000 20.0000000 22.0000000 26.0000000 8.00000000 [7,] -2.5298221 -1.2649111 0.0000000 1.2649111 2.52982213# 所有值取絕對值，再取對數 （取對數前一般加一個數避免對0或負值取對數） > log2(abs(a)+1)[,1] [,2] [,3] [,4] [,5] [1,] 0.4982872 0.1659818 0.5620435 0.9794522 1.0646224 [2,] 1.5250147 0.1598608 1.9743587 1.2255009 0.0464076 [3,] 0.9072054 0.1763961 0.1622189 0.8210076 0.3607278 [4,] 0.9354687 0.2387621 0.5893058 0.3329807 0.6604014 [5,] 1.0000000 1.5849625 2.0000000 2.3219281 2.5849625 [6,] 1.5849625 3.4594316 3.5849625 3.8073549 2.3219281 [7,] 1.1794544 0.7070437 0.0000000 0.7070437 1.1794544# 取出最大值、最小值、行數、列數 > max(a) [1] 13 > min(a) [1] -1.264911 > nrow(a) [1] 7 > ncol(a) [1] 5# 增加一列或一行 # cbind: column bind > cbind(a, 1:7)[,1] [,2] [,3] [,4] [,5] [,6] [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 1 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 2 [3,] 0.8754092 0.1300574 0.1190069 0.7666394 0.28407356 3 [4,] -0.9125118 0.1799797 0.5045226 0.2596132 -0.58052230 4 [5,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 5 [6,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 6 [7,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106 7 > cbind(a, seven=1:7)seven [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 1 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 2 [3,] 0.8754092 0.1300574 0.1190069 0.7666394 0.28407356 3 [4,] -0.9125118 0.1799797 0.5045226 0.2596132 -0.58052230 4 [5,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 5 [6,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 6 [7,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106 7# rbind: row bind > rbind(a,1:5)[,1] [,2] [,3] [,4] [,5] [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 [3,] 0.8754092 0.1300574 0.1190069 0.7666394 0.28407356 [4,] -0.9125118 0.1799797 0.5045226 0.2596132 -0.58052230 [5,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 [6,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 [7,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106 [8,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000# 計算每一行的mad (中值絕對偏差，一般認為比方差的魯棒性更強，更少受異常值的影響，更能反映數據間的差異) > apply(a,1,mad) [1] 0.7923976 2.0327283 0.2447279 0.4811672 1.4826000 4.4478000 0.9376786# 計算每一行的var (方差) # apply表示對數據（第一個參數）的每一行 (第二個參數賦值為1) 或每一列 (2)操作 # 最后返回一個列表 > apply(a,1,var) [1] 0.6160264 1.6811161 0.1298913 0.3659391 2.5000000 22.5000000 1.0000000# 計算每一列的平均值 > apply(a,2,mean) [1] 0.4519068 1.6689045 2.4395294 2.7179083 1.5753421# 取出中值絕對偏差大于0.5的行 > b = a[apply(a,1,mad)>0.5,] > b[,1] [,2] [,3] [,4] [,5] [1,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 [3,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 [4,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 [5,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106# 矩陣按照mad的大小降序排列 > c = b[order(apply(b,1,mad), decreasing=T),] > c[,1] [,2] [,3] [,4] [,5] [1,] 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 [2,] 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 [3,] 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 [4,] -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106 [5,] -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243> rownames(c) <- paste('Gene', letters[1:5], sep="_") > colnames(c) <- toupper(letters[1:5]) > cA B C D E Gene_a 2.0000000 10.0000000 11.0000000 13.0000000 4.00000000 Gene_b 1.8778966 -0.1171794 2.9295352 1.3383662 -0.03269026 Gene_c 1.0000000 2.0000000 3.0000000 4.0000000 5.00000000 Gene_d -1.2649111 -0.6324555 0.0000000 0.6324555 1.26491106 Gene_e -0.4125356 0.1219293 -0.4763589 -0.9717165 1.09162243# 矩陣轉置 > expr = t(c) > exprGene_a Gene_b Gene_c Gene_d Gene_e A 2 1.87789657 1 -1.2649111 -0.4125356 B 10 -0.11717937 2 -0.6324555 0.1219293 C 11 2.92953522 3 0.0000000 -0.4763589 D 13 1.33836620 4 0.6324555 -0.9717165 E 4 -0.03269026 5 1.2649111 1.0916224# 矩陣值的替換 > expr2 = expr > expr2[expr2<0] = 0 > expr2Gene_a Gene_b Gene_c Gene_d Gene_e A 2 1.877897 1 0.0000000 0.0000000 B 10 0.000000 2 0.0000000 0.1219293 C 11 2.929535 3 0.0000000 0.0000000 D 13 1.338366 4 0.6324555 0.0000000 E 4 0.000000 5 1.2649111 1.0916224# 矩陣中只針對某一列替換 # expr2是個矩陣不是數據框，不能使用列名字索引 > expr2[expr2$Gene_b<1, "Gene_b"] <- 1 Error in expr2$Gene_b : $ operator is invalid for atomic vectors # str是一個最為常用、好用的查看變量信息的工具，尤其是對特別復雜的變量， # 可以看清其層級結構，便于提取數據 > str(expr2)num [1:5, 1:5] 2 10 11 13 4 ...- attr(*, "dimnames")=List of 2..$ : chr [1:5] "A" "B" "C" "D" .....$ : chr [1:5] "Gene_a" "Gene_b" "Gene_c" "Gene_d" ...# 轉換為數據庫，再進行相應的操作 > expr2 <- as.data.frame(expr2) > str(expr2) 'data.frame': 5 obs. of 5 variables:$ Gene_a: num 2 10 11 13 4$ Gene_b: num 1.88 1 2.93 1.34 1$ Gene_c: num 1 2 3 4 5$ Gene_d: num 0 0 0 0.632 1.265$ Gene_e: num 0 0.122 0 0 1.092 > expr2[expr2$Gene_b<1, "Gene_b"] <- 1 > expr2 > expr2Gene_a Gene_b Gene_c Gene_d Gene_e A 2 1.877897 1 0.0000000 0.0000000 B 10 1.000000 2 0.0000000 0.1219293 C 11 2.929535 3 0.0000000 0.0000000 D 13 1.338366 4 0.6324555 0.0000000 E 4 1.000000 5 1.2649111 1.0916224

R中矩陣篩選合并

# 讀入樣品信息 > sampleInfo = "Samp;Group;Genotype + A;Control;WT + B;Control;WT + D;Treatment;Mutant + C;Treatment;Mutant + E;Treatment;WT + F;Treatment;WT" > phenoData = read.table(text=sampleInfo,sep=";", header=T, row.names=1, quote="") > phenoDataGroup Genotype A Control WT B Control WT D Treatment Mutant C Treatment Mutant E Treatment WT F Treatment WT# 把樣品信息按照基因表達矩陣中的樣品信息排序，并只保留有基因表達信息的樣品 # match() returns a vector of the positions of (first) matches ofits first argument in its second. > phenoData[match(rownames(expr), rownames(phenoData)),]Group Genotype A Control WT B Control WT C Treatment Mutant D Treatment Mutant E Treatment WT# ‘%in%’ is a more intuitive interface as a binary operator, whichreturns a logical vector indicating if there is a match or not forits left operand.# 注意順序，%in%比match更好理解一些 > phenoData = phenoData[rownames(phenoData) %in% rownames(expr),] > phenoDataGroup Genotype A Control WT B Control WT C Treatment Mutant D Treatment Mutant E Treatment WT# 合并矩陣 # by=0 表示按照行的名字排序 # by=columnname 表示按照共有的某一列排序 # 合并后多出了新的一列Row.names > merge_data = merge(expr, phenoData, by=0, all.x=T) > merge_dataRow.names Gene_a Gene_b Gene_c Gene_d Gene_e Group Genotype 1 A 2 1.87789657 1 -1.2649111 -0.4125356 Control WT 2 B 10 -0.11717937 2 -0.6324555 0.1219293 Control WT 3 C 11 2.92953522 3 0.0000000 -0.4763589 Treatment Mutant 4 D 13 1.33836620 4 0.6324555 -0.9717165 Treatment Mutant 5 E 4 -0.03269026 5 1.2649111 1.0916224 Treatment WT> rownames(merge_data) <- merge_data$Row.names > merge_data Row.names Gene_a Gene_b Gene_c Gene_d Gene_e Group Genotype A A 2 1.87789657 1 -1.2649111 -0.4125356 Control WT B B 10 -0.11717937 2 -0.6324555 0.1219293 Control WT C C 11 2.92953522 3 0.0000000 -0.4763589 Treatment Mutant D D 13 1.33836620 4 0.6324555 -0.9717165 Treatment Mutant E E 4 -0.03269026 5 1.2649111 1.0916224 Treatment WT# 去除一列；-1表示去除第一列 > merge_data = merge_data[,-1] > merge_dataGene_a Gene_b Gene_c Gene_d Gene_e Group Genotype A 2 1.87789657 1 -1.2649111 -0.4125356 Control WT B 10 -0.11717937 2 -0.6324555 0.1219293 Control WT C 11 2.92953522 3 0.0000000 -0.4763589 Treatment Mutant D 13 1.33836620 4 0.6324555 -0.9717165 Treatment Mutant E 4 -0.03269026 5 1.2649111 1.0916224 Treatment WT# 提取出所有的數值列 > merge_data[sapply(merge_data, is.numeric)]Gene_a Gene_b Gene_c Gene_d Gene_e A 2 1.87789657 1 -1.2649111 -0.4125356 B 10 -0.11717937 2 -0.6324555 0.1219293 C 11 2.92953522 3 0.0000000 -0.4763589 D 13 1.33836620 4 0.6324555 -0.9717165 E 4 -0.03269026 5 1.2649111 1.0916224

str的應用

str: Compactly display the internal *str*ucture of an R object, a
diagnostic function and an alternative to ‘summary (and to some
extent, ‘dput’). Ideally, only one line for each ‘basic’
structure is displayed. It is especially well suited to compactly
display the (abbreviated) contents of (possibly nested) lists.
The idea is to give reasonable output for any R object. It
calls ‘args’ for (non-primitive) function objects.

str用來告訴結果的構成方式，對于不少Bioconductor的包，或者復雜的R函數的輸出，都是一堆列表的嵌套，str(complex_result)會輸出每個列表的名字，方便提取對應的信息。

# str的一個應用例子 > str(list(a = "A", L = as.list(1:100)), list.len = 9) List of 2$ a: chr "A"$ L:List of 100..$ : int 1..$ : int 2..$ : int 3..$ : int 4..$ : int 5..$ : int 6..$ : int 7..$ : int 8..$ : int 9.. [list output truncated]# 利用str查看pca的結果，具體的PCA應用查看http://mp.weixin.qq.com/s/sRElBMkyR9rGa4TQp9KjNQ> pca_result <- prcomp(expr) > pca_result Standard deviations: [1] 4.769900e+00 1.790861e+00 1.072560e+00 1.578391e-01 2.752128e-16Rotation:PC1 PC2 PC3 PC4 PC5 Gene_a 0.99422750 -0.02965529 0.078809521 0.01444655 0.06490461 Gene_b 0.04824368 -0.44384942 -0.885305329 0.03127940 0.12619948 Gene_c 0.08258192 0.81118590 -0.451360828 0.05440417 -0.35842886 Gene_d -0.01936958 0.30237826 -0.079325524 -0.66399283 0.67897952 Gene_e -0.04460135 0.22948437 -0.002097256 0.74496081 0.62480128 > str(pca_result) List of 5$ sdev : num [1:5] 4.77 1.79 1.07 1.58e-01 2.75e-16$ rotation: num [1:5, 1:5] 0.9942 0.0482 0.0826 -0.0194 -0.0446 .....- attr(*, "dimnames")=List of 2.. ..$ : chr [1:5] "Gene_a" "Gene_b" "Gene_c" "Gene_d" ..... ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...$ center : Named num [1:5] 8 1.229 3 0.379 0.243..- attr(*, "names")= chr [1:5] "Gene_a" "Gene_b" "Gene_c" "Gene_d" ...$ scale : logi FALSE$ x : num [1:5, 1:5] -6.08 1.86 3.08 5.06 -3.93 .....- attr(*, "dimnames")=List of 2.. ..$ : chr [1:5] "A" "B" "C" "D" ..... ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...- attr(*, "class")= chr "prcomp"# 取出每個主成分解釋的差異 > pca_result$sdev [1] 4.769900e+00 1.790861e+00 1.072560e+00 1.578391e-01 2.752128e-16

R的包管理

# 什么時候需要安裝包 > library('unExistedPackage') Error in library("unExistedPackage") : 不存在叫‘unExistedPackage’這個名字的程輯包# 安裝包 > install.packages("package_name") # 指定安裝來源 > install.packages("package_name", repo="http://cran.us.r-project.org")# 安裝Bioconductor的包 > source('https://bioconductor.org/biocLite.R') > biocLite('BiocInstaller') > biocLite(c("RUVSeq","pcaMethods"))# 安裝Github的R包 > install.packages("devtools") > devtools::install_github("JustinaZ/pcaReduce")# 手動安裝, 首先下載包的源文件（壓縮版就可），然后在終端運行下面的命令。 ct@ehbio:~$ R CMD INSTALL package.tar.gz# 移除包 >remove.packages("package_name")# 查看所有安裝的包 >library()# 查看特定安裝包的版本 > installed.packages()[c("DESeq2"), c("Package", "Version")]Package Version "DESeq2" "1.14.1" > # 查看默認安裝包的位置 >.libPaths()# 調用安裝的包 >library(package_name)#devtools::install_github("hms-dbmi/scde", build_vignettes = FALSE) #install.packages(c("mvoutlier","ROCR")) #biocLite(c("RUVSeq","pcaMethods","SC3","TSCAN","monocle","MultiAssayExperiment","SummarizedExperiment")) #devtools::install_github("satijalab/seurat")

熱圖美化

上一期的繪圖命令中，最后一行的操作抹去了之前設定的橫軸標記的旋轉，最后出來的圖比較難看。

上次我們是這么寫的

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 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 # 標準化數據，并轉換為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繪制到現在就差不多了，但總是這么畫下去也會覺得有點累，有沒有辦法更簡化呢？ 且聽下回分解。

熱圖繪制 - pheatmap

繪制熱圖除了使用ggplot2，還可以有其它的包或函數，比如pheatmap::pheatmap (pheatmap包中的pheatmap函數)、gplots::heatmap.2等。

相比于ggplot2作heatmap, pheatmap會更為簡單一些，一個函數設置不同的參數，可以完成行列聚類、行列注釋、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="") 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 pheatmap::pheatmap(data, filename="pheatmap_1.pdf")

雖然有點丑，但一步就出來了。

在heatmap美化篇提到的數據前期處理方式，都可以用于pheatmap的畫圖。此外Z-score計算在pheatmap中只要一個參數就可以實現。

pheatmap::pheatmap(data, scale="row", filename="pheatmap_1.pdf")

有時可能不需要行或列的聚類，原始展示就可以了。

pheatmap::pheatmap(data, scale="row", cluster_rows=FALSE, cluster_cols=FALSE, filename="pheatmap_1.pdf")

給矩陣 (data)中行和列不同的分組注釋。假如有兩個文件，第一個文件為行注釋，其第一列與矩陣中的第一列內容相同 (順序沒有關系)，其它列為第一列的不同的標記，如下面示例中(假設行為基因，列為樣品)的2,3列對應基因的不同類型 (TF or enzyme)和不同分組。第二個文件為列注釋，其第一列與矩陣中第一行內容相同，其它列則為樣品的注釋。

row_anno = data.frame(type=c("TF","Enzyme","Enzyme","TF"), class=c("clu1","clu1","clu2","clu2"), row.names=rownames(data)) row_anno type class a TF clu1 b Enzyme clu1 c Enzyme clu2 d TF clu2 col_anno = data.frame(grp=c("A","A","A","B","B"), size=1:5, row.names=colnames(data)) col_anno grp size Grp_1 A 1 Grp_2 A 2 Grp_3 A 3 Grp_4 B 4 Grp_5 B 5 pheatmap::pheatmap(data, scale="row", cluster_rows=FALSE, annotation_col=col_anno, annotation_row=row_anno, filename="pheatmap_1.pdf")

自定義下顏色吧。

# <bias> values larger than 1 will give more color for high end. # Values between 0-1 will give more color for low end. pheatmap::pheatmap(data, scale="row", cluster_rows=FALSE, annotation_col=col_anno, annotation_row=row_anno, color=colorRampPalette(c('green','yellow','red'), bias=1)(50), filename="pheatmap_1.pdf")

heatmap.2的使用就不介紹了，跟pheatmap有些類似，而且也有不少教程。

不改腳本的熱圖繪制

繪圖時通常會碰到兩個頭疼的問題：
1. 需要畫很多的圖，唯一的不同就是輸出文件，其它都不需要修改。如果用R腳本，需要反復替換文件名，繁瑣又容易出錯。

每次繪圖都需要不斷的調整參數，時間久了不用，就忘記參數放哪了；或者調整次數過多，有了很多版本，最后不知道用哪個了。

為了簡化繪圖、維持腳本的一致，我用bash對R做了一個封裝，然后就可以通過修改命令好參數繪制不同的圖了。

先看一看怎么使用

首先把測試數據存儲到文件中方便調用。數據矩陣存儲在heatmap_data.xls文件中；行注釋存儲在heatmap_row_anno.xls文件中；列注釋存儲在heatmap_col_anno.xls文件中。

# tab鍵分割，每列不加引號 write.table(data, file="heatmap_data.xls", sep="\t", row.names=T, col.names=T,quote=F) # 如果看著第一行少了ID列不爽，可以填補下 system("sed -i '1 s/^/ID\t/' heatmap_data.xls")write.table(row_anno, file="heatmap_row_anno.xls", sep="\t", row.names=T, col.names=T,quote=F) write.table(col_anno, file="heatmap_col_anno.xls", sep="\t", row.names=T, col.names=T,quote=F)

然后用程序sp_pheatmap.sh繪圖。

# -f: 指定輸入的矩陣文件 # -d：指定是否計算Z-score，<none> (否), <row> (按行算), <col> (按列算) # -P: 行注釋文件 # -Q: 列注釋文件 ct@ehbio:~/$ sp_pheatmap.sh -f heatmap_data.xls -d row -P heatmap_row_anno.xls -Q heatmap_col_anno.xls

一個回車就得到了下面的圖

字有點小，是因為圖太大了，把圖的寬和高縮小下試試。

# -f: 指定輸入的矩陣文件 # -d：指定是否計算Z-score，<none> (否), <row> (按行算), <col> (按列算) # -P: 行注釋文件 # -Q: 列注釋文件 # -u: 設置寬度，單位是inch # -v: 設置高度，單位是inch ct@ehbio:~/$ sp_pheatmap.sh -f heatmap_data.xls -d row -P heatmap_row_anno.xls -Q heatmap_col_anno.xls -u 8 -v 12

橫軸的標記水平放置

# -A: 0, X軸標簽選擇0度 # -C: 自定義顏色，注意引號的使用，最外層引號與內層引號不同，引號之間無交叉 # -T: 指定給定的顏色的類型；如果給的是vector (如下面的例子), 則-T需要指定為vector; 否則結果會很怪異，只有倆顏色。 # -t: 指定圖形的題目，注意引號的使用；參數中包含空格或特殊字符等都要用引號引起來作為一個整體。 ct@ehbio:~/$ sp_pheatmap.sh -f heatmap_data.xls -d row -P heatmap_row_anno.xls -Q heatmap_col_anno.xls -u 8 -v 12 -A 0 -C 'c("white", "blue")' -T vector -t "Heatmap of gene expression profile"

sp_pheatmap.sh的參數還有一些，可以完成前面講述過的所有熱圖的繪制，具體如下：

***CREATED BY Chen Tong (chentong_biology@163.com)***----Matrix file-------------- Name T0_1 T0_2 T0_3 T4_1 T4_2 TR19267|c0_g1|CYP703A2 1.431 0.77 1.309 1.247 0.485 TR19612|c1_g3|CYP707A1 0.72 0.161 0.301 2.457 2.794 TR60337|c4_g9|CYP707A1 0.056 0.09 0.038 7.643 15.379 TR19612|c0_g1|CYP707A3 2.011 0.689 1.29 0 0 TR35761|c0_g1|CYP707A4 1.946 1.575 1.892 1.019 0.999 TR58054|c0_g2|CYP707A4 12.338 10.016 9.387 0.782 0.563 TR14082|c7_g4|CYP707A4 10.505 8.709 7.212 4.395 6.103 TR60509|c0_g1|CYP707A7 3.527 3.348 2.128 3.257 2.338 TR26914|c0_g1|CYP710A1 1.899 1.54 0.998 0.255 0.427 ----Matrix file------------------Row annorarion file -------------- ------1. At least two columns-------------- ------2. The first column should be the same as the first column inmatrix (order does not matter)-------------- Name Clan Family TR19267|c0_g1|CYP703A2 CYP71 CYP703 TR19612|c1_g3|CYP707A1 CYP85 CYP707 TR60337|c4_g9|CYP707A1 CYP85 CYP707 TR19612|c0_g1|CYP707A3 CYP85 CYP707 TR35761|c0_g1|CYP707A4 CYP85 CYP707 TR58054|c0_g2|CYP707A4 CYP85 CYP707 TR14082|c7_g4|CYP707A4 CYP85 CYP707 TR60509|c0_g1|CYP707A7 CYP85 CYP707 TR26914|c0_g1|CYP710A1 CYP710 CYP710 ----Row annorarion file ------------------Column annorarion file -------------- ------1. At least two columns-------------- ------2. The first column should be the same as the first row in ---------matrix (order does not matter)-------------- Name Sample T0_1 T0 T0_2 T0 T0_3 T0 T4_1 T4 T4_2 T4 ----Column annorarion file --------------Usage:sp_pheatmap.sh optionsFunction:This script is used to do heatmap using package pheatmap.The parameters for logical variable are either TRUE or FALSE.OPTIONS:-f Data file (with header line, the first column is therowname, tab seperated. Colnames must be unique unless youknow what you are doing.)[NECESSARY]-t Title of picture[Default empty title]["Heatmap of gene expression profile"]-a Display xtics. [Default TRUE]-A Rotation angle for x-axis value (anti clockwise)[Default 90]-b Display ytics. [Default TRUE]-H Hieratical cluster for columns.Default FALSE, accept TRUE-R Hieratical cluster for rows.Default TRUE, accept FALSE-c Clustering method, Default "complete". Accept "ward.D", "ward.D2","single", "average" (=UPGMA), "mcquitty" (=WPGMA), "median" (=WPGMC) or "centroid" (=UPGMC)-C Color vector. Default pheatmap_default. Aceept a vector containing multiple colors such as <'c("white", "blue")'> will be transferred to <colorRampPalette(c("white", "blue"), bias=1)(30)>or an R function <colorRampPalette(rev(brewer.pal(n=7, name="RdYlBu")))(100)>generating a list of colors.-T Color type, a vetcor which will be transferred as described in <-C> [vector] ora raw vector [direct vector] or a function [function (default)].-B A positive number. Default 1. Values larger than 1 will give more colorfor high end. Values between 0-1 will give more color for low end. -D Clustering distance method for rows.Default 'correlation', accept 'euclidean', "manhattan", "maximum", "canberra", "binary", "minkowski". -I Clustering distance method for cols.Default 'correlation', accept 'euclidean', "manhattan", "maximum", "canberra", "binary", "minkowski". -L First get log-value, then do other analysis.Accept an R function log2 or log10. [Default FALSE]-d Scale the data or not for clustering and visualization.[Default 'none' means no scale, accept 'row', 'column' to scale by row or column.]-m The maximum value you want to keep, any number larger willlbe taken as this given maximum value.[Default Inf, Optional] -s The smallest value you want to keep, any number smaller willbe taken as this given minimum value.[Default -Inf, Optional] -k Aggregate the rows using kmeans clustering. This is advisable if number of rows is so big that R cannot handle their hierarchical clustering anymore, roughly more than 1000.Instead of showing all the rows separately one can cluster therows in advance and show only the cluster centers. The numberof clusters can be tuned here.[Default 'NA' which means nocluster, other positive interger is accepted for executingkmeans cluster, also the parameter represents the number ofexpected clusters.]-P A file to specify row-annotation with format described above.[Default NA]-Q A file to specify col-annotation with format described above.[Default NA]-u The width of output picture.[Default 20]-v The height of output picture.[Default 20] -E The type of output figures.[Default pdf, accepteps/ps, tex (pictex), png, jpeg, tiff, bmp, svg and wmf)]-r The resolution of output picture.[Default 300 ppi]-F Font size [Default 14]-p Preprocess data matrix to avoid 'STDERR 0 in cor(t(mat))'.Lowercase <p>.[Default TRUE]-e Execute script (Default) or just output the script.[Default TRUE]-i Install the required packages. Normmaly should be TRUE if this is your first time run s-plot.[Default FALSE]

sp_pheatmap.sh是我寫作的繪圖工具s-plot的一個功能，s-plot可以繪制的圖的類型還有一些，列舉如下；在后面的教程中，會一一提起。

Usage:s-plot optionsFunction:This software is designed to simply the process of plotting and help researchers focus more on data rather than technology.Currently, the following types of plot are supported.#### Bars s-plot barPlot s-plot horizontalBar s-plot multiBar s-plot colorBar#### Lines s-plot lines#### Dots s-plot pca s-plot scatterplot s-plot scatterplot3d s-plot scatterplot2 s-plot scatterplotColor s-plot scatterplotContour s-plot scatterplotLotsData s-plot scatterplotMatrix s-plot scatterplotDoubleVariable s-plot contourPlot s-plot density2d#### Distribution s-plot areaplot s-plot boxplot s-plot densityPlot s-plot densityHistPlot s-plot histogram#### Cluster s-plot hcluster_gg (latest) s-plot hcluster s-plot hclust (depleted)#### Heatmap s-plot heatmapS s-plot heatmapM s-plot heatmap.2 s-plot pheatmap s-plot pretteyHeatmap # obseleted s-plot prettyHeatmap#### Others s-plot volcano s-plot vennDiagram s-plot upsetView

為了推廣，也為了激起大家的熱情，如果想要sp_pheatmap.sh腳本的，還需要勞煩大家動動手，轉發此文章到朋友圈，并留言索取。

Reference

• http://blog.genesino.com/2017/06/R-Rstudio

3

2

總結

以上是生活随笔為你收集整理的R绘制热图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。