[R]https://www.r-project.org/是一套完備的專業(yè)數(shù)據(jù)分析和圖形展示的統(tǒng)計(jì)環(huán)境和編程語言。R是免費(fèi)的，并且支持所有常用的系統(tǒng)?；赗的[Bioconductor]http://www.bioconductor.org提供有大量的生物科學(xué)領(lǐng)域的數(shù)據(jù)分析軟件包。
本教程目的在于讓沒有接觸過R的人快速入門。 在閱讀的過程中如果遇到不理解的地方，可以不用太糾結(jié)，繼續(xù)往后讀，可能讀了后邊的一些內(nèi)容，就會(huì)理解前邊的一些內(nèi)容。 然而，對(duì)于提供的代碼，最還還是要都輸入到R控制臺(tái)去看看效果。
安裝R和R包http://www.itdecent.cn/p/41690c94954f
1、R之初體驗(yàn)
在R控制臺(tái)逐行輸入如下代碼，不要想太多為什么，先看看是什么效果：
1
1:10
1 + 2
"2"
1 + "2"
1:10 + 1
c(1,3,13) + 1
c(1,3,13) * 5
pi
sin(pi)
log2(2)
log10(10)
log(2)
log(2,2)
log2(2,2)
round(pi,2)
(i+1)*(i+2)
(1i+3)*(2i+2)
sqrt(2)
2^(1/2)
(-1)^(1/2)
(0i - 1)^(1/2)
as.complex(-1)
x <- 1
x
y <- sin(x)
y
(z <- asin(y)) #加括號(hào)的時(shí)候會(huì)同時(shí)在屏幕輸出對(duì)象的值
z
2.1 R的一般語法結(jié)構(gòu)
object_name <- function(arguments)
R里的一切都可以看作是對(duì)象(object)，將一個(gè)或多個(gè)對(duì)象作為參數(shù)(arguments)，經(jīng)過一系列運(yùn)算輸出特定的值的對(duì)象，就是函數(shù)(function)。 這里的 <nowiki><-</nowiki> 是賦值符，就是把后邊的結(jié)果取個(gè)名字叫object_name。 這里的object_name本身不是object，就像你的名字和你的關(guān)系一樣。 這里的function可以是運(yùn)算符的形式，也可以省略（也就是直接在控制臺(tái)輸入對(duì)象或者其名稱時(shí)，就會(huì)返回相應(yīng)的結(jié)果，比如對(duì)象的值或者摘要）。
另外，除了用賦值運(yùn)算符以外，還可以用assign函數(shù)，如下：
assign("object_name", function(arguments))
需要區(qū)別的是這里的"object_name"是用引號(hào)引起來的，這和不引的時(shí)候是有區(qū)別的，這就牽涉到R的數(shù)據(jù)類型（data type）的問題。
2.2 R數(shù)據(jù)類型
R的數(shù)據(jù)類型（data type）有數(shù)值(Numeric)、字符(Character) 、復(fù)數(shù)（Complex）邏輯（Logical） 等。
使用mode查看
數(shù)值(Numeric)
x <- c(1, 3, 5)
x
is.numeric(x)
(y <- as.character(x))
is.numeric(y)
字符(Character)
x <- c("1", "2", "3") # 創(chuàng)建一個(gè)字符型的向量
x
is.character(x) #檢查是否是字符型
as.numeric(x) #轉(zhuǎn)化成為數(shù)值型
as.numeric("A") #強(qiáng)制把一些字符轉(zhuǎn)化成為數(shù)值會(huì)產(chǎn)生NA值
is.numeric(as.numeric("A"))
as.character(as.numeric("A"))
is.numeric(as.character(as.numeric("A")))
復(fù)數(shù)（Complex）
sqrt(-1) #數(shù)值型的-1是沒法開根號(hào)的
mode(-1)
as.complex(-1) #轉(zhuǎn)成復(fù)數(shù)形式
sqrt(as.complex(-1)) #可以開根號(hào)了
sqrt(0i+1) #加上 xi （x為數(shù)值），就是復(fù)數(shù)的形式
sqrt(3i-4)
is.complex(3i-4)
邏輯（Logical）
TRUE # TRUE 和 FALSE就是邏輯值
5 < 10
5 > 10
1:10 < 5 #返回邏輯值
sum(1:10 < 5) #可以求和，TRUE就是1，F(xiàn)ALSE就是0
as.character(1:10 < 5)
as.numeric(1:10 < 5) #邏輯值轉(zhuǎn)為數(shù)值
mode(c(TRUE, FALSE, NA))
mode(c(TRUE, FALSE, "NA")) #引號(hào)引起來的時(shí)候就是一個(gè)字符型
mode(c(TRUE, FALSE, 1))
mode(c(TRUE, FALSE, 1))
mode(c(1, 2, "1"))
mode(c(1, 2, "1"))
mode(c(1, 2, 0i))
0i
mode(c(1i, 2i, "1"))
一個(gè)向量只能是一種類型，在使用c()函數(shù)創(chuàng)建向量的時(shí)候，如果將不同類型的值放在一起，會(huì)有一個(gè)轉(zhuǎn)化，這種優(yōu)先級(jí)是這樣的：字符型 > 復(fù)數(shù) > 數(shù)值 > 邏輯 。一個(gè)向量只有一種類型，在列表讀入的時(shí)候，也會(huì)遇到這種同一列里可能有幾種類型的情況，這時(shí)候也會(huì)有轉(zhuǎn)化，比如一列數(shù)字里如果具有含有標(biāo)點(diǎn)的字符，這時(shí)候就會(huì)被當(dāng)作字符型的讀入。 而R在默認(rèn)的情情況下會(huì)將字符轉(zhuǎn)為因子型的，這樣如果之間使用as.numeric轉(zhuǎn)換成數(shù)值型的，就會(huì)發(fā)生錯(cuò)誤，就需要先轉(zhuǎn)成字符型的才行。
這時(shí)候就可以使用如下的命令來關(guān)掉R的這功能，避免出現(xiàn)沒有預(yù)料的情況。
options(stringsAsFactors=FALSE)
2.3 R對(duì)象類型
R常用的對(duì)象類型有向量(vector)、因子(factor)、數(shù)據(jù)框(data frame)、矩陣(array)、列表(list)、函數(shù)(function)等。
向量(vector)
按順序排列的數(shù)值、字符、復(fù)數(shù)或者邏輯值的集合，就是向量。單獨(dú)的一個(gè)這些類型的值出現(xiàn)的時(shí)候，可以看作是長(zhǎng)度為一的向量。
c(1,4,7)
x <- c("cat", "dog", "mouse", "human")
x[1] #取出第一個(gè)元素，要注意的是R和Python等不一樣的地方是，index是從1開始的而不是從0開始的。
x[-2] #去掉第二個(gè)元素，這里-號(hào)的意義和Python里也是不一樣的。
x[-c(1,3)]
x[10] #R里超過了向量的長(zhǎng)度取值也不會(huì)報(bào)錯(cuò)，而是返回NA值
x[3] <- "mice"
x # 可以看到，是很容易對(duì)向量里的值進(jìn)行替換的
x[c(1,2,1,3,2,2,4,4)] #再體會(huì)一下R向量取元素的方法
因子(factor)
含有分組信息的向量，比較特殊。
x <- c("馬", "羊", "牛", "雞", "鴨")
y <- factor(c("家畜", "家畜", "家畜", "家禽", "家禽"))
mode(y)
as.numeric(y)
as.character(y)
(z <- factor(c("1", "2", "4", "5", "7")))
as.numeric(z)
as.numeric(as.character(z)) # 這就是上文中講到的問題，需要注意讀入文本框的過程中可能把數(shù)字當(dāng)字符
數(shù)據(jù)框(data frame)
二維的結(jié)構(gòu)，每列是一個(gè)向量?？梢钥醋魇且环N特殊的列表，即每個(gè)元素的長(zhǎng)度一致。
options(stringsAsFactors=FALSE)
(x <- data.frame(A=1:10, B=letters[1:10])) #創(chuàng)建一個(gè)數(shù)據(jù)框
x[1, ] #數(shù)據(jù)框中的值的提取
x[ , 1]
x[1, 1]
x[1:3, 1:2]
x[1:3, 1:2] <- data.frame(LETTERS[2:4], 3:5) #可用同樣的行列的數(shù)據(jù)框來替換數(shù)據(jù)框里同樣行列的數(shù)據(jù)
x
x[ , "A"] #也可以直接用列名
x[ ,3] #體會(huì)這和向量里是不一樣的，不能超出列表的長(zhǎng)度來取值
x[1 ,3]
x[1 , ][3]
as.numeric(x[1 , ])[3] #as.numeric 將只有一行的data frame轉(zhuǎn)成了vector
as.numeric(x[1:2 , ]) #而對(duì)于兩行的data frame卻不能轉(zhuǎn)了
矩陣(matrices)
含有同一種數(shù)據(jù)類型的二維結(jié)構(gòu)，其本質(zhì)是vector，但是具有行和列的屬性。
x <- matrix(1:20, ncol=4, nrow=5) #創(chuàng)建一個(gè)具有4列5行的矩陣
x[1， 3] #矩陣取出元素的方法
x[4:2, c(2,4)] #體會(huì)一下R里的靈活取元素的方法
x[4:2, c(2,4)] <- "OK" #和向量或者數(shù)據(jù)框一樣可以替換
x # 留意現(xiàn)在數(shù)據(jù)類型的變化
x[4:2, c(2,4)] <- (1:6)^2 #同樣長(zhǎng)度的向量的替換，長(zhǎng)度不等時(shí)重復(fù)使用
x #變成字符型之后就變不回來了。
x[10] # 這里就能看出矩陣其實(shí)是向量
數(shù)組(array)
矩陣可以看作二維的數(shù)組，因此推廣一下理解更高維數(shù)組
# Create an array.
a <- array(c('green','yellow'),dim = c(3,3,2))
print(a)
列表(list)
列表可以看作更為一般的向量，只是其元素可以是不同的數(shù)據(jù)類型，甚至是不同的對(duì)象類型，比如list的元素可以是list。
x=list(A=1:3, B=letters[1:10])
x
函數(shù)(function)
將一系列代碼放在一起，就成為函數(shù)，函數(shù)可以預(yù)設(shè)參數(shù)
f <- function(){
print("A")
print("Z")
1+2
}
f
f() #記住函數(shù)的調(diào)用方式
f <- function(x="A", y="B", a=1, b=2){ #可以先設(shè)置默認(rèn)值
print(x)
print(y)
a+b
}
f()
f("a", "b", 3, 5) #參數(shù)可以按順序?qū)?
f(b=10, x="OK") #也可以不安順序但是指明誰是誰
R的函數(shù)分為向量化的和非向量化的，向量化的函數(shù)對(duì)參數(shù)中的向量逐個(gè)元素地進(jìn)行運(yùn)算，逐一返回運(yùn)算結(jié)果組成新的向量。運(yùn)算符也是一種函數(shù)，只是用法和一般的函數(shù)有別。在使用函數(shù)時(shí)需要區(qū)分向量化的和非向量化的區(qū)別。 一般來說，常用的一些函數(shù)用習(xí)慣了就知道怎么用了，確實(shí)不需要大腦里經(jīng)常去想哪個(gè)是向量化的哪個(gè)不是。
我們可以看到，我們上邊區(qū)分了數(shù)據(jù)類型和對(duì)象類型，其實(shí)這兩者很難去區(qū)分，數(shù)據(jù)也是對(duì)象。
3.1常規(guī)取子集
前邊已經(jīng)說到了對(duì)象取子集，這里展開說一下。
R中對(duì)象取子集的基本語法
my_object[row] #一維對(duì)象取子集，如向量、因子、列表
my_object[row, col] #二維對(duì)象取子集，如矩陣、數(shù)據(jù)框
my_object[row, col, dim] #三維對(duì)象取子集，如三維數(shù)組
#以上語法可推廣到多維數(shù)組
=== 三種可能方法 ===
取子集語法中的row、col、dim、等等可能有如下三種值： \
1. 使用位置向量，向量前邊加"-"表示不包括 \
2. 使用等長(zhǎng)的邏輯向量 \
3. 使用元素名稱向量，當(dāng)元素?zé)o名稱時(shí)就不好用了 \
X <- 1:10 #創(chuàng)建一個(gè)向量
names(X) <- LETTERS[1:10] #為每個(gè)元素命名, LETTERS為內(nèi)置的26個(gè)字母按順序組成的向量
X
X[1:4] # 使用位置取第 1-4 元素
X[1, 2, 4] #
X[c(1, 2, 4)] # 體會(huì)一下與上一條命令的區(qū)別
X[-c(1:4)] # 不包括 第1-4個(gè)元素。
X[-1:4] # 體會(huì)一下與上一條命令的區(qū)別
X[c(1,3,5,8,6,5,4,2,50)]
X[-c(1,3,5,8,6,5,4,2,50)]
X >5 # 邏輯向量
X[X > 5] #使用邏輯向量取子集
X[c("B", "F", "J")] #使用元素名稱向量取子集
X <- data.frame(N=1:6, S=c("a", "b", "a", "c", "d", "b")) #創(chuàng)建一個(gè)數(shù)據(jù)庫(kù)
X
X[c(1, 3, 6), ] #使用位置向量
X[X$N < 5 & X$N > 1, ] #使用邏輯向量
X[X$S %in% c("a", "c"), ] #使用邏輯向量
X[ , c("S", "N")] #通過元素名，這里是列名
row.names(X) <- LETTERS[1:6] #更改行名
X[c("D", "E", "A"), ] #使用行名
Y <- list(A=1:3, B=letters[1:5], C=sample(letters, 3)) #創(chuàng)建一個(gè)列表, sample用于隨機(jī)抽取
Y[1] #返回的是一個(gè)列表
Y[c("A","B")]
Y[1:2] #
$符號(hào)的使用
對(duì)于列表（包括數(shù)據(jù)框）可以使用$加元素的名稱來取出相應(yīng)的元素。
Y <- list(A=1:3, B=letters[1:5], C=sample(letters, 3)) #創(chuàng)建一個(gè)列表, sample用于隨機(jī)抽取
Y$A
Y$A[2]
Y$B[2:3]
[[...]]的使用
[[...]]可用于取出單一的元素，可遞歸
X <- 1:10 #創(chuàng)建一個(gè)向量
names(X) <- LETTERS[1:10] #為每個(gè)元素命名, LETTERS為內(nèi)置的26個(gè)字母按順序組成的向量
X
X[[1]]
X[[1:2]] #遞歸失敗
Y <- list(A=1:3, B=letters[1:5], C=sample(letters, 3)) #創(chuàng)建一個(gè)列表, sample用于隨機(jī)抽取
Y[1] #返回的是一個(gè)列表
Y[[1]] #返回的是一個(gè)向量(元素)
Y[c("A","B")] #取子集
Y[[c("A","B")]] #遞歸失敗，"A"元素里沒有"B"元素
Y[1:2] #
Y[[1:2]] #遞歸成功，Y第一個(gè)元素里的第二個(gè)元素被取出
Z <- list(a = list(b = 9, c = "hello", d=list(e="3rd", f=list(g="4th"))), h = 1:5)
Z
Z[[c("a", "b")]]
Z[[c("a", "b", "c")]]
Z[[c("a", "d", "f")]]
Z[[c("a", "d", "f", "g")]]
Z$a$d$f$g
3.2 給元素的對(duì)象賦值
在取子集時(shí)賦值，即能改變對(duì)象中相應(yīng)元素的值
X <- 1:10 #創(chuàng)建一個(gè)向量
names(X) <- LETTERS[1:10] #為每個(gè)元素命名, LETTERS為內(nèi)置的26個(gè)字母按順序組成的向量
X
X[5] <- 1000
X
X["B"] <- "b"
X
Y <- list(A=1:3, B=letters[1:5], C=sample(letters, 3)) #創(chuàng)建一個(gè)列表, sample用于隨機(jī)抽取
Y
Y[1:2] <- 1
Y
Y$B <- "abc"
Y
4.1 R的運(yùn)算
向量運(yùn)算符
*算術(shù)運(yùn)算：加(+)、減(-)、乘()、除以(/)、指數(shù)(^)、整除(%/%)、取模(%%)
**比較運(yùn)算：相等()、不等(!=)、大于(>)、小于(<)、大于等于(>=)、小于等于(<=)
**邏輯運(yùn)算：或(|)、且(&)、非(!)
算術(shù)運(yùn)算針對(duì)數(shù)值型或者復(fù)數(shù)型的對(duì)象進(jìn)行加減乘數(shù)的基本算術(shù)，比較運(yùn)算、邏輯運(yùn)算結(jié)果為邏輯型。
10+2*3+4/5-2*6
10 + 2*3 + 4/5 -2*6 + 2^3
10.3/2
11.3 %/% 2
11.3 %% 2
1 == 2
1 <= 2
1 <= 2 & 1 == 2
1 <= 2 | 1 == 2
!(1 <= 2 & 1 == 2)
向量運(yùn)算
如上提到的運(yùn)算符都可以用于向量和向量之間，兩個(gè)向量之前相同位置上的數(shù)值之間進(jìn)行運(yùn)算，運(yùn)算結(jié)果按順序返回生成新的向量。當(dāng)兩個(gè)向量長(zhǎng)度不等時(shí)，較短的向量將重新從第一個(gè)元素開始重復(fù)使用，如此直到較長(zhǎng)的向量所有元素都完成。這樣當(dāng)一個(gè)向量與長(zhǎng)度為1的向量進(jìn)行運(yùn)算時(shí)，相當(dāng)于每個(gè)元素都與這個(gè)值做同樣的運(yùn)算。
1:5 + 10:14
1:5 + 10:19
1:5 + 10:23 #當(dāng)長(zhǎng)度之間不是整數(shù)倍的時(shí)候會(huì)提示，但是依然有運(yùn)算結(jié)果
更多運(yùn)算符
&& (且，非向量化的，用于長(zhǎng)度為1的向量)；
|| (或，非向量化的，用于長(zhǎng)度為1的向量)；
%in% (前一個(gè)向量的元素是否出現(xiàn)在后一個(gè)向量里,對(duì)前一個(gè)向量中每個(gè)元素依次判斷返回向量結(jié)果)
<- (向左賦值)
-> (向右賦值)
TRUE && c(TRUE, FALSE) # && 和 || 一般用于 if() 中
TRUE && c(FALSE, TRUE)
FALSE || c(TRUE, FALSE)
FALSE || c(FALSE, TRUE)
c(1,3,5) %in% c(2,3,5,6)
c(2,3,5,6) %in% c(1,3,5)
y <- 1 -> x
y
x
5.1 常用函數(shù)
獲取幫助
help(help) # 查看幫助信息，能查看已載入的包的，需要注意的是對(duì)一些R的關(guān)鍵字需要用引號(hào)引起來，比如"if" "[[" 等
help.search("help") # 對(duì)所有的包進(jìn)行搜索，包括未載入的, 留意這里的參數(shù)必須是字符型的
? help # help()的快捷辦法，
? "if"
?? help
學(xué)會(huì)使用help函數(shù)是至關(guān)重要的，每次遇到不理解的地方，就查看一下help，積少成多，慢慢就成高手了。這里只提供一些函數(shù)的最基本的用法，更為廣泛的用法都可以通過help來查看詳細(xì)的參數(shù)說明。 大多數(shù)情況也也沒有必要記住那些參數(shù)，因?yàn)楹芏鄥?shù)一般情況下用不到。
工作目錄
getwd() # 查看當(dāng)前的工作目錄
dir() # 查看當(dāng)前工作目錄下文件和文件夾
setwd('/BJPROJ/RNA/rna_test/XXXXXX') # 設(shè)置工作目錄為 /BJPROJ/RNA/rna_test/XXXXXX
setwd("D:/XXXX")
setwd("D:\\XXXX") # 注意windows下的文件夾路徑的兩種寫法。
數(shù)學(xué)運(yùn)算
(x <- 1:10)
(y <- sample(1:100,10))
(z <- sample(1:100,10)) #先創(chuàng)建三個(gè)向量，用于展示
sum(x) # 數(shù)值型向量的所有元素的總和
prod(x) # 數(shù)值型向量的所有元素的乘積
mean(x) # 數(shù)值型向量的所有元素的平均值
sd(x) # 標(biāo)準(zhǔn)差
median(x) # 中位數(shù)
max(x) # 最大值
min(x) # 最小值
var(x) # 方差 variance of x and the covariance or correlation of x and y if these are vectors
var(x, y) # 協(xié)方差
cov(x, y) # 協(xié)方差
cor(x, y) # 皮爾遜相關(guān)性系數(shù)
cor(x, y, method = "spearman") # spearman相關(guān)性系數(shù)
max(x, y, z)
pmax(x, y, z) #體會(huì)這種向量化的運(yùn)算與max函數(shù)的區(qū)別
pmin(x, y, z)
which.max(z) #返回向量中最大值的index
which.min(y)
sum(1:10) #求和
cumsum(1:10) #依次累加，看看效果
prod(1:10) #求積
cumprod(1:10) #累乘
字符處理
(x <- "ATTATGCATCGGC")
(y <- c("CGATTACGGC", "ATGAATGGC", "ATGTCGC", "ATNAGGGC"))
nchar(y) # 字符串中的字符數(shù)，這也是向量化的函數(shù)
length(y) # 向量y的長(zhǎng)度，即元素個(gè)數(shù)
substr(y, 2,3) #向量化操作，對(duì)每個(gè)元素都進(jìn)行同樣操作，返回對(duì)應(yīng)順序結(jié)果
substring(x,4) #從第4個(gè)字符開始取至結(jié)尾
strsplit(x, "T") # 以T為間隔講字符串分割，返回結(jié)果為列表的形式
strsplit(y, "T") # 以T為間隔講字符串分割，返回結(jié)果為列表的形式
grep("CG", y) # 返回含特定字符的元素的序號(hào)(index)
grep("CG", y, value = TRUE) #返回含特定字符的元素
grepl("CG", y) # 向量化運(yùn)算，返回邏輯型向量
grepl("^CG", y) # ^表示起始
grepl("C.G", y) # . 表示任意值
grepl("A.G", y) # . 表示任意值
? regex # 查看R正則表達(dá)式的用法
gsub("CG","_CG_", y)
(z <- letters[1:5])
paste(x, z, sep="_") #向量化的函數(shù)，將兩個(gè)向量的元素粘在一起，以特定字符相隔
paste(z, collapse="_") #設(shè)置collapse參數(shù)，將向量中所有的元素粘在一起，以特定字符相隔
sprintf
這個(gè)函數(shù)用于將一些變量格式化并加入字符串中。
sprintf("%s is %f feet tall\n", "Sven", 7.1)
sprintf("%s is %f,: %d", "A", 99, 50)
sprintf("%f", pi) #以下展示不同形式的pi,不用都記住，知道有這么回事就好
sprintf("%.3f", pi)
sprintf("%1.0f", pi)
sprintf("%5.1f", pi)
sprintf("%05.1f", pi)
sprintf("%+f", pi)
sprintf("% f", pi)
sprintf("%-10f", pi)
sprintf("%e", pi)
sprintf("%E", pi)
sprintf("%g", pi)
sprintf("%g", 1e6 * pi)
sprintf("%.9g", 1e6 * pi)
sprintf("%G", 1e-6 * pi)
(A <- sprintf("%03d", seq(0,100,5)))
print(A)
6.1 讀寫數(shù)據(jù)與基本作圖
這里我們使用rpkm.xls和id2symbol.tsv作為例子。下載后放到工作目錄或者把工作目錄設(shè)成存放此文件的文件夾。
rpkm <- read.delim("rpkm.xls", header=TRUE, sep='\t') # 這里的rpkm.xls 文件實(shí)際是純文本的，并不是EXCEL格式
head(rpkm) #查看前幾行
Lines <- readLines("rpkm.xls") # 逐行讀入
head(Lines) #注意與rpkm的區(qū)別
RPKM <- rpkm
names(RPKM) <- letters[1:10] #更改表頭，也就是列標(biāo)題
head(RPKM)
rm(RPKM) #移除RPKM
write.table(rpkm, 'rpkm.tsv', row.names=FALSE, col.names=TRUE, sep='\t', quote=F) #保存數(shù)據(jù)框到表格，sep='\t' 這個(gè)參數(shù)在聚群上使用時(shí)需要注意，集群上默認(rèn)的是空格，因此，需要使用制表符的時(shí)候，一定要加上。
plot(x=c(1,2,5) ,y= c(4,8,6)) #散點(diǎn)圖，x和y為長(zhǎng)度相等的向量，同一位置數(shù)值為一個(gè)點(diǎn)的橫軸坐標(biāo)。
plot(rpkm$T1,rpkm$T2) #兩個(gè)參數(shù)，不指明x ， y的時(shí)候，前邊的為x，后邊為y
plot(y=log(rpkm$T1+0.1), x=log(rpkm$T2+0.1), type='h',col='#FF0000') #type指明圖的類型
hist(log2(rpkm$T1+0.1)) #直方圖
plot(density(log2(rpkm$T1+0.1))) #密度曲線圖，先用density函數(shù)，再用plot函數(shù)
str(density(log2(rpkm$T1+0.1))) #看看dengsity函數(shù)返回的對(duì)象的元素
pie(c(A=1,b=2,C=3), # 餅圖
labels = c("ok","byeye","p"), # labels指定扇形的標(biāo)記
col=c("red","green","blue")) # 指定顏色，顏色的選擇可參考：http://colorbrewer2.org/
png('pie_test.png', type="cairo-png") #圖像的保存為png格式，這里默認(rèn)輸出到當(dāng)前工作目錄
pie(c(A=1,b=2,C=3), c("ok","byeye","p"), col=c("red","green","blue"))
dev.off()
pdf('./pie_test.pdf') #圖像保存為pdf格式，文件名可以用絕對(duì)路徑，改變輸出到特定的目錄
pie(c(A=1,b=2,C=3), c("ok","byeye","p"), col=c("red","green","blue"))
dev.off()
apply序列
R里盡量避免for循環(huán)，運(yùn)行速度比較慢。需要對(duì)一系列值做同樣的操作的時(shí)候，可以使用apply序列的函數(shù)。
rpkm <- read.delim("rpkm.xls", header=TRUE, sep='\t')
rpkm$mean_Ts <- apply(rpkm[ , 2:4], 1, mean) #apply 對(duì)數(shù)據(jù)框逐行運(yùn)算，mean可以換成別的方程，可以自己定義
head(rpkm) #看一下計(jì)算的結(jié)果
apply(rpkm[ , -1], 2, mean) #對(duì)數(shù)據(jù)框的每列進(jìn)行運(yùn)算，第二個(gè)參數(shù)填2
lapply(1:10, function(x) x^2)
lapply(as.list(1:10), function(x) x^2) #lapply用于向量或列表，對(duì)列表每一個(gè)函數(shù)做同樣的操作，依次返回結(jié)果，結(jié)果為列表
sapply(1:10, function(x) x^2)
sapply(as.list(1:10), function(x) x^2) #sapply用于向量或列表，對(duì)列表每一個(gè)函數(shù)做同樣的操作，依次返回結(jié)果，結(jié)果為向量
rpkm$mean_W <- mapply(mean, rpkm$W1, rpkm$W2, rpkm$W3) #mapply是sapply的多變量版本，返回的值是向量
head(rpkm) #注意上邊的這種用法是錯(cuò)誤的
rpkm$mean_W <- mapply(function(x,y,z) mean(c(x,y,z)), x=rpkm$W1, y=rpkm$W2, z=rpkm$W3)
head(rpkm) #注意mean函數(shù)的參數(shù)形式
未分類
2:9 # : 用于產(chǎn)生間隔為1的等差數(shù)列, 從前一個(gè)數(shù)開始，最大不超過后一個(gè)數(shù)
2.1:9.2
seq(2, 13, 2) # 產(chǎn)生等差數(shù)列，從第一個(gè)參數(shù)開始，以第三個(gè)參數(shù)為差，不大于第二個(gè)數(shù)
seq(2, 13, length=6) # 使用length參數(shù)，產(chǎn)生指定長(zhǎng)度的等差數(shù)列，差通過計(jì)算得到
rep(c(1,3,4), 3) # 重復(fù)向量指定的次數(shù)
rep(c(1,3,4), each=3) # 每個(gè)元素依次進(jìn)行重復(fù)
unique(c(1,2,3,1,2,2,4,5)) # 返回唯一的值
combn(1:5, 2) # 求組合,
combn(1:5, 3)
# rpkm <- read.delim("rpkm.xls", header=TRUE, sep='\t')
ID2Symbol <- read.delim("id2symbol.tsv")
RPKM <- merge(ID2Symbol, rpkm, by="geneID", all.y=TRUE) # merge 函數(shù)根據(jù)數(shù)據(jù)框的某一列將兩個(gè)數(shù)據(jù)框聯(lián)合在一起
head(RPKM) # 用來查看前邊的幾行
tail(RPKM) # 用來查看后邊的幾行
dim(RPKM) # 查看行數(shù)和列數(shù)
names(ID2Symbol) <- c("ID", "Symbol")
RPKM <- merge(ID2Symbol, rpkm, by.x="ID", by.y="geneID") #
head(RPKM)
dim(RPKM)
str(rpkm) #查看R對(duì)象的結(jié)構(gòu)(structure), 可以看到每個(gè)元素的類型
x <- 1:10
stopifnot(x <= 10) #stopifnot函數(shù)用于判斷一個(gè)邏輯向量是否都是TRUE
stopifnot(x <= 5)
all(x <= 5) #判斷邏輯向量是否都是TRUE
any(x <=5) #判斷是否有一個(gè)為TRUE
which(x<=8 & x >=5) #返回TRUE的值的index
7.1 R中的流程控制
常用的就是流程控制為for, if, else, while等，R里的流程控制語法中注意括號(hào)的用法。
for(i in 1:10){ # for對(duì)所有的元素依次執(zhí)行
if(i < 5){ # if做出判斷
print(sprintf("%s < 5", i))
} else if(i < 8) { # else if 部分可以沒有，也可以有多個(gè)
print(sprintf("%s >= 5 & %s < 8", i, i))
} else { #
print(sprintf("%s >= 8", i))
}
}
for(i in 1:10){
if(i < 5){ #
next # 跳出當(dāng)前循環(huán)小節(jié)
print(sprintf("%s < 5", i))
} else if(i < 8) { #
print(sprintf("%s >= 5 & %s < 8", i, i))
} else { #
print(sprintf("%s >= 8", i))
break # 終止整個(gè)for循環(huán)
}
}
i=1
while(i <= 20){ # 滿足條件的時(shí)候執(zhí)行，需預(yù)先定義i，且循環(huán)中i應(yīng)變化，否則即是死循環(huán)
if(i < 5){ #
print(sprintf("%s < 5", i))
i = i + 2
} else if(i < 8) { #
print(sprintf("%s >= 5 & %s < 8", i, i))
i = i + 1
} else { #
print(sprintf("%s >= 8", i))
i = i + 3
}
}
8.1 R包安裝與運(yùn)用
R之所以強(qiáng)大，是因?yàn)橛写罅康腞包，通過安裝新的R包獲得新的功能。R包根據(jù)不同的來源有不同的安裝方法:
R包安裝
CRAN中的R包
install.packages(c("ggplot2", "plyr", "reshape", "gplots", "pheatmap")) #可以看到這也是個(gè)向量化的函數(shù)， 可以同時(shí)安裝多個(gè)包
install.packages("devtools") #安裝github中包需要用到devtools這個(gè)包
#BioConductor中的R包
source("https://bioconductor.org/biocLite.R") #不行的話試著改成http://bioconductor.org/biocLite.R
biocLite("DESeq2") #安裝DESeq2， 其它包仿此
#github中的R包，一些CRAN中R包的最新開發(fā)中版本可以在github中找到
library(devtools) #載入devtools
install_github("hadley/ggplot2") #安裝github中的開發(fā)中版本
#其它來源的版本
#一般都會(huì)說明如何安裝，也可以下載下來安裝，這種情況比較少，具體問題具體再說。
install.packages("路徑到文件或者網(wǎng)絡(luò)連接到文件")
R包的運(yùn)用
library(gglot2) # 載入R包，載入之后就能使用其中函數(shù)和數(shù)據(jù)了
browseVignettes("ggplot2") # 查看手冊(cè)
search() # 查看已經(jīng)加載的R包以及對(duì)象
detach("package:pheatmap") # 卸載已加載的R包，在需要換一個(gè)版本的包的時(shí)候可以用到
接下來會(huì)對(duì)一些包當(dāng)中的常用函數(shù)做介紹，這些包的功能比這里介紹的要強(qiáng)大地多，通過考核之后可以自己去仔細(xì)研究。
reshape包
library(reshape) #如果沒有安裝這個(gè)包，自己安裝一下
# rpkm <- read.delim("rpkm.xls", header=TRUE, sep='\t')
RPKM <- melt(rpkm, id=1) #轉(zhuǎn)換數(shù)據(jù)的格式，自己體會(huì)一下效果， id設(shè)置不變的列，這里是第一列不變
head(RPKM)
head( melt(rpkm, id=c(1,3,5))) #1,3,5列不變，其它的melt
names(RPKM) <- c("geneID", "sample", "rpkm") # 注意作為列名的rpkm和對(duì)象名rpkm的區(qū)別
Rpkm <- cast(RPKM, geneID~sample, value="rpkm") # 體會(huì)一下效果
head(Rpkm)
plyr包
library(plyr)
stat <- ddply(RPKM, .(sample), summarize, Mean=mean(rpkm), SD=sd(rpkm), SEM=sd(rpkm)/sqrt(length(rpkm)))
#使用ddply來根據(jù)一列或多列來求另一列按這兩列的值來分類的各種值，可自定義函數(shù)。
# 這里是根據(jù)樣本一列，求一個(gè)樣本中所有基因的表達(dá)均值等等的
head(stat)
# 大家可以自己想一下，怎么分基因求幾個(gè)樣本中的均值和標(biāo)準(zhǔn)差等等的
Stat <- ddply(RPKM, .(sample), Mean=mean(rpkm), SD=sd(rpkm), SEM=sd(rpkm)/sqrt(length(rpkm)))
head(Stat) # 沒有summarize是不可以的
pheatmap包
library(pheatmap)
test_data <- rpkm[apply(rpkm[ , 2:10], 1, sd) !=0, ] #有點(diǎn)復(fù)雜，自己慢慢領(lǐng)會(huì)，領(lǐng)會(huì)不了也沒關(guān)系，這里重點(diǎn)在pheatmap
pheatmap(test_data[1:100, 2:10], scale="row", cluster_cols = FALSE) #做熱圖
pdf("heatmap.pdf", height=7, width=7)
pheatmap(test_data[1:100, 2:10], scale="row", cluster_cols = FALSE) #做熱圖
dev.off()
pdf("heatmap2.pdf", height=10, width=7, onefile=FALSE) # 行數(shù)比較多，增加高度
pheatmap(test_data[1:100, 2:10], scale="row", cluster_cols = FALSE, labels_row = test_data$geneID, fontsize_row = 5) #labels_row可以手動(dòng)設(shè)置行的標(biāo)記 ， fontsize_row 設(shè)置字號(hào)，這里變小
dev.off()
ggplot2畫圖
ggplot2是很強(qiáng)大的作圖R包，學(xué)了ggplot2之后似乎都可以把基本的畫圖函數(shù)忘掉了。
以下是一些例子，可以試試結(jié)合各種元素畫出不同風(fēng)格和外觀的圖。
library(ggplot2)
library(reshape)
library(RColorBrewer)
rpkm_mel <- melt(rpkm, id=1) # ? melt
rpkm_mel$group <- substr(rpkm_mel$variable, 1, 1) #選取樣本名的首字母作為組名（示例用，隨便選的）
rpkm_mel$group <- factor(rpkm_mel$group, levels=c("W","T")) #改變group的順序
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value)))
g #使用g來喊出圖來
g <-
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable)) # 增加colour
g
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable, fill=group)) #增加fill
g #使用g來喊出圖來
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable, fill=group)) + #增加fill
scale_fill_manual(values=c("red","green")) # 手動(dòng)設(shè)置填充顏色
g #使用g來喊出圖來
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable, fill=group)) + #增加fill
scale_fill_manual(values=c("red","green")) + # 手動(dòng)設(shè)置填充顏色
scale_x_discrete(labels=c(1:9)) #設(shè)置橫坐標(biāo)
g #使用g來畫出圖來
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_boxplot(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable, fill=group))+
xlab("samples")+
ylab("log RPKM") +
ggtitle("boxplot of gene expression") +
scale_x_discrete(labels=c(1:9)) +
scale_y_continuous(limits=c(-10,10), breaks=seq(-10,10,by=2.5))+
scale_fill_manual(values=c("red","green")) +
theme_bw() +
theme(legend.position='top', # ? theme查看更多的theme參數(shù)
legend.direction="horizontal",
axis.text.x=element_text(size=20)
)
ggsave("boxplot.pdf", plot=g)
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_histogram(data=rpkm_mel, aes(x=log(value), fill=variable), position='dodge',binwidth=1) + #控制bin的寬度
facet_grid(.~group) +
xlab("log RPKM")+
ggtitle("histogram of gene expression") +
scale_fill_manual(values=brewer.pal(9, 'Set1'), name='sample') + #用 ? RColorBrewer 查看這個(gè)顏色的信息
theme_bw() +
theme(legend.position='top',
legend.direction="horizontal",
axis.text.x=element_text(size=20)
)
ggsave("histogram.pdf", plot=g)
g <- #將后邊的ggplot返回的對(duì)象“命名”為g
ggplot() +
geom_point(data=rpkm_mel, aes(x=variable, y=log(value), colour=variable, fill=group), position='jitter')+ v#jitter 會(huì)在坐標(biāo)上隨機(jī)加上一些偏離
xlab("sample")+
ylab("log RPKM") +
ggtitle("scatterplot of gene expression") +
scale_x_discrete(labels=c(1:9)) +
scale_y_continuous(limits=c(-10,10), breaks=seq(-10,10,by=2.5))+
scale_fill_manual(values=c("red","green")) +
theme_bw() +
theme(legend.position='top',
legend.direction="horizontal",
axis.text.x=element_text(size=20)
)
ggsave("scatterplot.pdf", plot=g, useDingbats=FALSE)
散點(diǎn)圖中包括圓圈(shape=1,10,13,16,19,20等)時(shí)，在一些可能閱讀器中可能會(huì)顯示為q，這時(shí)設(shè)置：useDingbats=FALSE
R腳本
除了如上在R控制臺(tái)運(yùn)行R的代碼以外，還能把代碼放到一個(gè)文件里，作為腳本的形式運(yùn)行,比如這個(gè)test.R腳本。
print("Hi RNA")
cat("Hi Hong")
Rscript test1.R
另外，還能將一些變量作為參數(shù)的形式傳入，只需要用commandArgs函數(shù)即可。(trailingOnly = FALSE)
args <- commandArgs(TRUE)
person1 <- args[1]
person2 <- args[2]
#Say Hellow to Person 1
cat(sprintf("Hi %s \n", person1))
#Say Hellow to Person 2
cat(sprintf("Hi %s \n", person2))
Rscript test2.R Hong RNA
