歡迎關(guān)注 生信寶典 公眾號(hào)，閱讀系列文章

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

Rstudio基礎(chǔ)

R語(yǔ)言是比較常用的統(tǒng)計(jì)分析和繪圖語(yǔ)言，擁有強(qiáng)大的統(tǒng)計(jì)庫(kù)、繪圖庫(kù)和生信分析的Bioconductor庫(kù)，是學(xué)習(xí)生物信息分析的必備語(yǔ)言之一。

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

Rstudio版本

Rsdutio分為桌面版和服務(wù)器版，桌面版可以在單機(jī)使用，服務(wù)器版可以從瀏覽器訪問(wèn)供多人使用。

服務(wù)器版安裝好之后，訪問(wèn)地址為<服務(wù)器IP:8787>，用戶名和密碼為L(zhǎng)inux用戶的用戶名和密碼。

Rstudio安裝

R安裝

Linux下安裝

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

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

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

# --enable-R-shlib 需要設(shè)置，使得其他程序包括Rstudio可以使用R的動(dòng)態(tài)庫(kù) # --prefix指定軟件安裝目錄，需使用絕對(duì)路徑 ./configure --prefix=/home/ehbio/R/3.4.0 --enable-R-shlib# 也可以使用這個(gè)命令，共享系統(tǒng)的blas庫(kù)，提高運(yùn)輸速度 #./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下安裝服務(wù)器版

• 安裝參考 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

• 安裝完之后的檢測(cè)、啟動(dòng)和配置

  sudo rstudio-server verify-installation #查看是否安裝正確 sudo rstudio-server start ## 啟動(dòng) sudo rstudio-server status ## 查看狀態(tài) sudo rstudio-server stop ## 停止 ifconfig | grep 'inet addr' ## 查看服務(wù)端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日志目錄，方便查看錯(cuò)誤信息：/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/)雙擊安裝，需要使用管理員權(quán)限，其它無(wú)需要注意的。

Rstudio 使用

Windows下桌面版直接雙擊打開(kāi)即可使用，Linux服務(wù)器版訪問(wèn)地址為<服務(wù)器IP:8787>，用戶名和密碼為L(zhǎng)inux用戶的用戶名和密碼。

Rstudio 界面

Rstudio中新建或打開(kāi)文件

如果是桌面版，直接就可以訪問(wèn)”我的電腦”去打開(kāi)之前寫過(guò)的腳本。如果是服務(wù)器版，可直接訪問(wèn)服務(wù)器上寫過(guò)的腳本。Rstudio右下1/4部分可以切換目錄，點(diǎn)擊more，設(shè)置工作目錄?？梢陨蟼鞅镜氐哪_本到對(duì)應(yīng)目錄打開(kāi)。

Rstudio還有其它的功能，不過(guò)這些對(duì)我們初學(xué)使用已經(jīng)足夠了。后續(xù)會(huì)根據(jù)需要推出基于ggplot2作圖的R入門介紹。

熱圖繪制

熱圖是做分析時(shí)常用的展示方式，簡(jiǎn)單、直觀、清晰?？梢杂脕?lái)顯示基因在不同樣品中表達(dá)的高低、表觀修飾水平的高低等。任何一個(gè)數(shù)值矩陣都可以通過(guò)合適的方式用熱圖展示。

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

上一篇講述了Rstudio的使用作為R寫作和編譯環(huán)境的入門，后面的命令都可以拷貝到Rstudio中運(yùn)行，或?qū)懗梢粋€(gè)R腳本，使用Rscript heatmap.r運(yùn)行。我們還提供了Bash的封裝，在不修改R腳本的情況下，改變參數(shù)繪制出不同的圖形。

生成測(cè)試數(shù)據(jù)

繪圖首先需要數(shù)據(jù)。通過(guò)生成一堆的向量，轉(zhuǎn)換為矩陣，得到想要的數(shù)據(jù)。

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

注意：運(yùn)算符的優(yōu)先級(jí)。

> 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轉(zhuǎn)為矩陣 (matrix)，再轉(zhuǎn)為數(shù)據(jù)框 (data.frame)。

# ncol: 指定列數(shù) # byrow: 先按行填充數(shù)據(jù) # ?matrix 可查看函數(shù)的使用方法 # as.data.frame的as系列是轉(zhuǎn)換用的 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

雖然方法比較繁瑣，但一個(gè)數(shù)值矩陣已經(jīng)獲得了。

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

• 讀入字符串

# 使用字符串的好處是不需要額外提供文件 # 簡(jiǎn)單測(cè)試時(shí)可使用，寫起來(lái)不繁瑣，又方便重復(fù) # 尤其適用于在線提問(wèn)時(shí)作為測(cè)試案例 > 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"# 習(xí)慣設(shè)置quote為空，避免部分基因名字或注釋中存在引號(hào)，導(dǎo)致讀入文件錯(cuò)誤。 # 具體錯(cuò)誤可查看 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

可以看到列名字中以數(shù)字開(kāi)頭的列都加了X。一般要盡量避免行或列名字以數(shù)字開(kāi)頭，會(huì)給后續(xù)分析帶去一些困難；另外名字中出現(xiàn)的非字母、數(shù)字、下劃線、點(diǎn)的字符都會(huì)被轉(zhuǎn)為點(diǎn)，也需要注意，盡量只用字母、下劃線和數(shù)字。

# 讀入時(shí)，增加一個(gè)參數(shù)`check.names=F`也可以解決問(wèn)題。 # 這次數(shù)字前沒(méi)有再加 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="")

轉(zhuǎn)換數(shù)據(jù)格式

數(shù)據(jù)讀入后，還需要一步格式轉(zhuǎn)換。在使用ggplot2作圖時(shí)，有一種長(zhǎng)表格模式是最為常用的，尤其是數(shù)據(jù)不規(guī)則時(shí)，更應(yīng)該使用 （這點(diǎn)，我們?cè)谥v解箱線圖時(shí)再說(shuō)）。

# 如果包沒(méi)有安裝，運(yùn)行下面一句，安裝包 #install.packages(c("reshape2","ggplot2"))library(reshape2) library(ggplot2)# 轉(zhuǎn)換前，先增加一列ID列，保存行名字 data$ID <- rownames(data)# melt：把正常矩陣轉(zhuǎn)換為長(zhǎng)表格模式的函數(shù)。工作原理是把全部的非id列的數(shù)值列轉(zhuǎn)為1列，命名為value；所有字符列轉(zhuǎn)為variable列。# id.vars 列用于指定哪些列為id列；這些列不會(huì)被merge，會(huì)保留為完整一列。 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

分解繪圖

數(shù)據(jù)轉(zhuǎn)換后就可以畫圖了，分解命令如下：

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

熱圖出來(lái)了，但有點(diǎn)不對(duì)勁，橫軸重疊一起了。一個(gè)辦法是調(diào)整圖像的寬度，另一個(gè)是旋轉(zhuǎn)橫軸標(biāo)記。

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

設(shè)置想要的顏色。

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

調(diào)整legend的位置。

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

調(diào)整背景和背景格線以及X軸、Y軸的標(biāo)題。

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

合并以上命令，就得到了下面這個(gè)看似復(fù)雜的繪圖命令。

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

圖形存儲(chǔ)

圖形出來(lái)了，就得考慮存儲(chǔ)了，

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

至此，完成了簡(jiǎn)單的heatmap的繪圖。但實(shí)際繪制時(shí)，經(jīng)常會(huì)碰到由于數(shù)值變化很大，導(dǎo)致顏色過(guò)于集中，使得圖的可讀性下降很多。因此需要對(duì)數(shù)據(jù)進(jìn)行一些處理，具體的下次再說(shuō)。

R基礎(chǔ)概念和困惑點(diǎn)

插播一點(diǎn)R的基礎(chǔ)概念和疑難解釋。

R的交互式運(yùn)行

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

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()'通過(guò)HTML瀏覽器來(lái)看幫助文件。 用'q()'退出R.> 1:3 [1] 1 2 3 > a <- 1:3 > a [1] 1 2 3 # 使用source在交互式界面運(yùn)行寫好的R腳本 > source('scrtpt.r') > quit() Save workspace image? [y/n/c]: n# ctrl+d也可退出R

R基本語(yǔ)法

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

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

R中的變量

> # 數(shù)字變量 > 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) # 判斷一個(gè)變量是不是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) #查看或設(shè)置數(shù)組的維度向量 [1] 5 4 > > # 錯(cuò)誤的用法 > dim(a) <- c(4,4) Error in dim(a) <- c(4, 4) : dims [product 16]與對(duì)象長(zhǎng)度[20]不匹配 > > # 正確的用法 > a <- 1:20 > dim(a) <- c(5,4) #轉(zhuǎn)換向量為矩陣 > 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的行數(shù)", nrow(a))) [1] "矩陣a的行數(shù) 5" > print(paste("矩陣a的列數(shù)", ncol(a))) [1] "矩陣a的列數(shù) 4" > > #查看或設(shè)置行列名 > 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系列函數(shù)用來(lái)判斷變量的屬性和轉(zhuǎn)換變量的屬性 # 矩陣轉(zhuǎn)換為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中矩陣運(yùn)算

# 數(shù)據(jù)產(chǎn)生 # rnorm(n, mean = 0, sd = 1) 正態(tài)分布的隨機(jī)數(shù) # runif(n, min = 0, max = 1) 平均分布的隨機(jī)數(shù) # rep(1,5) 把1重復(fù)5次 # scale(1:5) 標(biāo)準(zhǔn)化數(shù)據(jù) > 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## 注意檢查括號(hào)的配對(duì) > a <- a[rowSums(abs(a)!=0,] 錯(cuò)誤: 意外的']' 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# 矩陣運(yùn)算，R默認(rèn)針對(duì)整個(gè)數(shù)據(jù)進(jìn)行常見(jiàn)運(yùn)算# 所有值都乘以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# 所有值取絕對(duì)值，再取對(duì)數(shù) （取對(duì)數(shù)前一般加一個(gè)數(shù)避免對(duì)0或負(fù)值取對(duì)數(shù)） > 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# 取出最大值、最小值、行數(shù)、列數(shù) > 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# 計(jì)算每一行的mad (中值絕對(duì)偏差，一般認(rèn)為比方差的魯棒性更強(qiáng)，更少受異常值的影響，更能反映數(shù)據(jù)間的差異) > apply(a,1,mad) [1] 0.7923976 2.0327283 0.2447279 0.4811672 1.4826000 4.4478000 0.9376786# 計(jì)算每一行的var (方差) # apply表示對(duì)數(shù)據(jù)（第一個(gè)參數(shù)）的每一行 (第二個(gè)參數(shù)賦值為1) 或每一列 (2)操作 # 最后返回一個(gè)列表 > apply(a,1,var) [1] 0.6160264 1.6811161 0.1298913 0.3659391 2.5000000 22.5000000 1.0000000# 計(jì)算每一列的平均值 > apply(a,2,mean) [1] 0.4519068 1.6689045 2.4395294 2.7179083 1.5753421# 取出中值絕對(duì)偏差大于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# 矩陣轉(zhuǎn)置 > 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# 矩陣中只針對(duì)某一列替換 # expr2是個(gè)矩陣不是數(shù)據(jù)框，不能使用列名字索引 > expr2[expr2$Gene_b<1, "Gene_b"] <- 1 Error in expr2$Gene_b : $ operator is invalid for atomic vectors # str是一個(gè)最為常用、好用的查看變量信息的工具，尤其是對(duì)特別復(fù)雜的變量， # 可以看清其層級(jí)結(jié)構(gòu)，便于提取數(shù)據(jù) > 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" ...# 轉(zhuǎn)換為數(shù)據(jù)庫(kù)，再進(jìn)行相應(yīng)的操作 > 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# 把樣品信息按照基因表達(dá)矩陣中的樣品信息排序，并只保留有基因表達(dá)信息的樣品 # 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# 提取出所有的數(shù)值列 > 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的應(yīng)用

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 ‘a(chǎn)rgs’ for (non-primitive) function objects.

str用來(lái)告訴結(jié)果的構(gòu)成方式，對(duì)于不少Bioconductor的包，或者復(fù)雜的R函數(shù)的輸出，都是一堆列表的嵌套，str(complex_result)會(huì)輸出每個(gè)列表的名字，方便提取對(duì)應(yīng)的信息。

# str的一個(gè)應(yīng)用例子 > 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的結(jié)果，具體的PCA應(yīng)用查看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"# 取出每個(gè)主成分解釋的差異 > pca_result$sdev [1] 4.769900e+00 1.790861e+00 1.072560e+00 1.578391e-01 2.752128e-16

R的包管理

# 什么時(shí)候需要安裝包 > library('unExistedPackage') Error in library("unExistedPackage") : 不存在叫‘unExistedPackage’這個(gè)名字的程輯包# 安裝包 > install.packages("package_name") # 指定安裝來(lái)源 > 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")# 手動(dòng)安裝, 首先下載包的源文件（壓縮版就可），然后在終端運(yùn)行下面的命令。 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" > # 查看默認(rèn)安裝包的位置 >.libPaths()# 調(diào)用安裝的包 >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")

熱圖美化

上一期的繪圖命令中，最后一行的操作抹去了之前設(shè)定的橫軸標(biāo)記的旋轉(zhuǎn)，最后出來(lái)的圖比較難看。

上次我們是這么寫的

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

為了使橫軸旋轉(zhuǎn)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))

最后的圖應(yīng)該是下邊樣子的。

上次的測(cè)試數(shù)據(jù)，數(shù)值的分布比較均一，相差不是太大，但是Gene_4和Gene_5由于整體的值低于其它的基因，從顏色上看，不仔細(xì)看，看不出差別。

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

輸出的結(jié)果是這個(gè)樣子的

圖中只有右上角可以看到紅色，其他地方就沒(méi)了顏色的差異。這通常不是我們想要的。為了更好的可視化效果，需要對(duì)數(shù)據(jù)做些預(yù)處理，主要有 對(duì)數(shù)轉(zhuǎn)換，Z-score轉(zhuǎn)換，抹去異常值，非線性顏色等方式。

對(duì)數(shù)轉(zhuǎn)換

假設(shè)下面的數(shù)據(jù)是基因表達(dá)數(shù)據(jù)，4個(gè)基因 (a, b, c, d)和5個(gè)樣品 (Grp_1, Grp_2, Grp_3, Grp_4)，矩陣中的值代表基因表達(dá)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")

對(duì)數(shù)轉(zhuǎn)換后的數(shù)據(jù)，看起來(lái)就清晰的多了。而且對(duì)數(shù)轉(zhuǎn)換后，數(shù)據(jù)還保留著之前的變化趨勢(shì)，不只是基因在不同樣品之間的表達(dá)可比 (同一行的不同列)，不同基因在同一樣品的值也可比 (同一列的不同行) (不同基因之間比較表達(dá)值存在理論上的問(wèn)題，即便是按照長(zhǎng)度標(biāo)準(zhǔn)化之后的FPKM也不代表基因之間是完全可比的)。

Z-score轉(zhuǎn)換

Z-score又稱為標(biāo)準(zhǔn)分?jǐn)?shù)，是一組數(shù)中的每個(gè)數(shù)減去這一組數(shù)的平均值再除以這一組數(shù)的標(biāo)準(zhǔn)差，代表的是原始分?jǐn)?shù)距離原始平均值的距離，以標(biāo)準(zhǔn)差為單位?？梢詫?duì)不同分布的各原始分?jǐn)?shù)進(jìn)行比較，用來(lái)反映數(shù)據(jù)的相對(duì)變化趨勢(shì)，而非絕對(duì)變化量。

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 # 標(biāo)準(zhǔn)化數(shù)據(jù)，并轉(zhuǎn)換為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轉(zhuǎn)換后，顏色分布也相對(duì)均一了，每個(gè)基因在不同樣品之間的表達(dá)的高低一目了然。但是不同基因之間就完全不可比了。

抹去異常值

粗暴一點(diǎn)，假設(shè)檢測(cè)飽和度為100，大于100的值都視為100對(duì)待。

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

雖然損失了一部分信息，但整體模式還是出來(lái)了。只是在選擇異常值標(biāo)準(zhǔn)時(shí)需要根據(jù)實(shí)際確認(rèn)。

非線性顏色

正常來(lái)講，顏色的賦予在最小值到最大值之間是均勻分布的。非線性顏色則是對(duì)數(shù)據(jù)比較小但密集的地方賦予更多顏色，數(shù)據(jù)大但分布散的地方賦予更少顏色，這樣既能加大區(qū)分度，又最小的影響原始數(shù)值。通?？梢愿鶕?jù)數(shù)據(jù)模式，手動(dòng)設(shè)置顏色區(qū)間。為了方便自動(dòng)化處理，我一般選擇用四分位數(shù)的方式設(shè)置顏色區(qū)間。

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")) # 獲取數(shù)據(jù)的最大、最小、第一四分位數(shù)、中位數(shù)、第三四分位數(shù) 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 # 在最小值和第一四分位數(shù)之間劃出6個(gè)區(qū)間，第一四分位數(shù)和中位數(shù)之間劃出6個(gè)區(qū)間，中位數(shù)和第三四分位數(shù)之間劃出5個(gè)區(qū)間，最后的數(shù)劃出5個(gè)區(qū)間 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 # 安照設(shè)定的區(qū)間分割數(shù)據(jù) # 原始數(shù)據(jù)替換為了其所在的區(qū)間的數(shù)值 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 # 雖然看上去還是數(shù)值，但已經(jīng)不是數(shù)字類型了 # 而是不同的因子了，這樣就可以對(duì)不同的因子賦予不同的顏色了 > 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

# 產(chǎn)生對(duì)應(yīng)數(shù)目的顏色 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逐個(gè)賦值 p <- p + scale_fill_manual(values=col) ggsave(p, filename="heatmap_nonlinear.pdf", width=8, height=12, units=c("cm"),colormodel="srgb")

調(diào)整行的順序或列

如果想保持圖中每一行的順序與輸入的數(shù)據(jù)框一致，需要設(shè)置因子的水平。這也是ggplot2中調(diào)整圖例或橫縱軸字符順序的常用方式。

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繪制到現(xiàn)在就差不多了，但總是這么畫下去也會(huì)覺(jué)得有點(diǎn)累，有沒(méi)有辦法更簡(jiǎn)化呢？ 且聽(tīng)下回分解。

熱圖繪制 - pheatmap

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

相比于ggplot2作heatmap, pheatmap會(huì)更為簡(jiǎn)單一些，一個(gè)函數(shù)設(shè)置不同的參數(shù)，可以完成行列聚類、行列注釋、Z-score計(jì)算、顏色自定義等。那我們來(lái)看看效果怎樣。

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

雖然有點(diǎn)丑，但一步就出來(lái)了。

在heatmap美化篇提到的數(shù)據(jù)前期處理方式，都可以用于pheatmap的畫圖。此外Z-score計(jì)算在pheatmap中只要一個(gè)參數(shù)就可以實(shí)現(xiàn)。

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

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

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

給矩陣 (data)中行和列不同的分組注釋。假如有兩個(gè)文件，第一個(gè)文件為行注釋，其第一列與矩陣中的第一列內(nèi)容相同 (順序沒(méi)有關(guān)系)，其它列為第一列的不同的標(biāo)記，如下面示例中(假設(shè)行為基因，列為樣品)的2,3列對(duì)應(yīng)基因的不同類型 (TF or enzyme)和不同分組。第二個(gè)文件為列注釋，其第一列與矩陣中第一行內(nèi)容相同，其它列則為樣品的注釋。

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有些類似，而且也有不少教程。

不改腳本的熱圖繪制

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

每次繪圖都需要不斷的調(diào)整參數(shù)，時(shí)間久了不用，就忘記參數(shù)放哪了；或者調(diào)整次數(shù)過(guò)多，有了很多版本，最后不知道用哪個(gè)了。

為了簡(jiǎn)化繪圖、維持腳本的一致，我用bash對(duì)R做了一個(gè)封裝，然后就可以通過(guò)修改命令好參數(shù)繪制不同的圖了。

先看一看怎么使用

首先把測(cè)試數(shù)據(jù)存儲(chǔ)到文件中方便調(diào)用。數(shù)據(jù)矩陣存儲(chǔ)在heatmap_data.xls文件中；行注釋存儲(chǔ)在heatmap_row_anno.xls文件中；列注釋存儲(chǔ)在heatmap_col_anno.xls文件中。

# tab鍵分割，每列不加引號(hào) write.table(data, file="heatmap_data.xls", sep="\t", row.names=T, col.names=T,quote=F) # 如果看著第一行少了ID列不爽，可以填補(bǔ)下 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：指定是否計(jì)算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

一個(gè)回車就得到了下面的圖

字有點(diǎn)小，是因?yàn)閳D太大了，把圖的寬和高縮小下試試。

# -f: 指定輸入的矩陣文件 # -d：指定是否計(jì)算Z-score，<none> (否), <row> (按行算), <col> (按列算) # -P: 行注釋文件 # -Q: 列注釋文件 # -u: 設(shè)置寬度，單位是inch # -v: 設(shè)置高度，單位是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

橫軸的標(biāo)記水平放置

# -A: 0, X軸標(biāo)簽選擇0度 # -C: 自定義顏色，注意引號(hào)的使用，最外層引號(hào)與內(nèi)層引號(hào)不同，引號(hào)之間無(wú)交叉 # -T: 指定給定的顏色的類型；如果給的是vector (如下面的例子), 則-T需要指定為vector; 否則結(jié)果會(huì)很怪異，只有倆顏色。 # -t: 指定圖形的題目，注意引號(hào)的使用；參數(shù)中包含空格或特殊字符等都要用引號(hào)引起來(lái)作為一個(gè)整體。 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的參數(shù)還有一些，可以完成前面講述過(guò)的所有熱圖的繪制，具體如下：

***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的一個(gè)功能，s-plot可以繪制的圖的類型還有一些，列舉如下；在后面的教程中，會(huì)一一提起。

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腳本的，還需要?jiǎng)跓┐蠹覄?dòng)動(dòng)手，轉(zhuǎn)發(fā)此文章到朋友圈，并留言索取。

Reference

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

3

2

總結(jié)

以上是生活随笔為你收集整理的R绘制热图的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。