劉小澤寫于18.12.10

生信必備三大件：生物、統(tǒng)計(jì)、技術(shù)，我想要借助R來學(xué)習(xí)統(tǒng)計(jì)學(xué)知識(shí)，因?yàn)槠綍r(shí)使用R比較頻繁，理解起來應(yīng)該也會(huì)更快一些。先來一次R的知識(shí)迭代，更新下知識(shí)庫

先看一些比較基礎(chǔ)的R操作

這其中有一些是我之前沒有學(xué)習(xí)到的，第一次get感覺超級(jí)有用

• 我們一般都會(huì)查詢陌生的函數(shù)如何使用，會(huì)用到help()。默認(rèn)狀態(tài)下，help()函數(shù)只會(huì)在載入加載過的程序包中搜索，如果想在所有包中進(jìn)行搜索，又不想先加載它，可以用help("xx", try.all.packages=TRUE)；如果已經(jīng)知道函數(shù)在哪個(gè)包中，可以用help("xx", package="yy")

• 數(shù)據(jù)直接編輯[這個(gè)東西很有意思]：借助edit()函數(shù)，我們可以直接在R中進(jìn)行數(shù)據(jù)修改（這個(gè)功能我其實(shí)一直在尋找，但最近才發(fā)現(xiàn)，這樣就不用導(dǎo)出后修改然后再導(dǎo)入R了）。
  例如，我們要修改mtcars這個(gè)數(shù)據(jù)集，但是為了避免修改原有的數(shù)據(jù)集，可以新保存到一個(gè)變量如new <- edit(mtcars) ，然后修改完點(diǎn)quit就好；
  如果想在現(xiàn)有數(shù)據(jù)集基礎(chǔ)上修改，不想新建，那么可以用fix（mtcars） ；
  創(chuàng)建一個(gè)新的數(shù)據(jù)集，再向其中填充數(shù)據(jù)new <- edit(data.frame()) ，默認(rèn)創(chuàng)建兩列，然后拖動(dòng)右下角就出現(xiàn)更多的列；

  Windows是可以直接使用edit的，Mac的話需要先安裝XQuartz這個(gè)工具

• 如果想研究數(shù)據(jù)集中的某個(gè)變量，比如mtcars中的mpg，可以使用mtcars$mpg，但是如果還需要?jiǎng)e的變量，那還是多敲幾遍mtcars$；其實(shí)可以先將數(shù)據(jù)集激活，放到環(huán)境中，什么時(shí)候想提取，就直接輸變量名就好了，比如先attach(mtcars)，就是將mtcars當(dāng)做當(dāng)前數(shù)據(jù)集，然后想調(diào)用mpg變量直接輸入mpg就好；不想用了就detach(mtcars)

• 不同數(shù)據(jù)，不同目的，不同分析：

  屬性數(shù)據(jù)

  例如mtcars中的cyl數(shù)據(jù)（記錄氣缸數(shù)）是屬性變量。先激活mtcars數(shù)據(jù)集，然后使用table(cyl)得到cyl中的三個(gè)值：4，6，8和相應(yīng)的頻數(shù)；還可以barplot畫直方圖barplot(table(cyl)) ，直接barplot(cyl)也能畫，但是分組信息就體現(xiàn)不出來，失去了意義

  數(shù)值型數(shù)據(jù)

  這一類在統(tǒng)計(jì)分析中經(jīng)常用到，例如mtcars中的mpg變量就是數(shù)值型數(shù)據(jù)；另外不同函數(shù)反映的意義也是不同的。

  作圖可以畫直方圖hist(mpg)，箱線圖boxplot()；

  統(tǒng)計(jì)分析可以分析均值mean()；計(jì)算截掉5%的均值mean(mpg, trim=.05) ；按分組變量cyl來計(jì)算mpg的分組均值tapply(mpg, cyl, mean) ；計(jì)算cyl為6的mpg的均值mean(mpg[cyl == 6 ]) ；計(jì)算常用分位數(shù)（包括極小、極大、中位數(shù)、兩個(gè)四分位數(shù)）quantile(mpg)或者fivenum(mpg) 或者summary(mpg)；計(jì)算四分位數(shù)的極差IQR(mpg)；標(biāo)準(zhǔn)差sd(mpg)；中位絕對(duì)離差(median absolute deviation) mad(mpg)

  探索二元關(guān)系

  1 擬合線性回歸，例如

  > z <- lm(cyl~mpg)
  > z
  # 結(jié)果
  Call:
  lm(formula = cyl ~ mpg)
  
  Coefficients:
  (Intercept)          mpg  
      11.2607      -0.2525 
  # 線性回歸截距11.26，斜率-0.25
  

  2 相關(guān)系數(shù)：考察回歸擬合的好壞，例如相關(guān)系數(shù)R可以用cor(cyl, mpg) 表示，更常用的R² 就能計(jì)算出來為0.726，表示數(shù)據(jù)變化的72.6%可以用氣缸數(shù)（cyl）與每加侖的英里數(shù)（mpg）表示

  3 殘差分析：殘差是估計(jì)值與觀測(cè)值的偏差

  z <- lm(cyl ~ mpg) #將回歸分析的結(jié)果作為對(duì)象保存在lm.res中
  lm.resids <- resid(z) # 提取殘差向量
  plot(lm.resids) # 畫個(gè)散點(diǎn)圖
  hist(lm.resids) # 畫個(gè)直方圖，是否為鐘形？
  qqnorm(lm.resids) # QQ圖(Quantile-Quantile Plots)是否落在直線上？
  # 都是為了檢驗(yàn)數(shù)據(jù)集分布假設(shè)是否有效
  

了解下R的對(duì)象

對(duì)象屬性

R運(yùn)行都是靠對(duì)象，所有的對(duì)象都有兩個(gè)內(nèi)在屬性：類型和長度

類型包括四種：數(shù)值型、字符型、復(fù)數(shù)型、邏輯型（T/F/NA），用函數(shù)mode()查看。另外不管什么類型的數(shù)據(jù)，缺失值都是用NA（Not Available）表示，不是數(shù)值用NaN(Not a Number)表示；

• 數(shù)值型：數(shù)值太大用指數(shù)形式表示，如N <- 2e20， Inf表示正無窮

• 字符型：輸入時(shí)要加上雙引號(hào)，如果要在其中繼續(xù)引用雙引號(hào)的話，可以用\進(jìn)行轉(zhuǎn)義；或者直接用單引號(hào)

  > x <- "Double quotes \" delimitate R's strings." # 雙引號(hào)
  > x <- 'Double quotes " delimitate R\'s strings.' # 單引號(hào)
  >x
  [1] "Double quotes \" delimitate R's strings." # 在R中顯示的狀態(tài)[注意這里在R中顯示的還帶著轉(zhuǎn)義符，其實(shí)顯示出來是沒有的]
  > cat(x)
  Double quotes " delimitate R's strings. # 實(shí)際上的狀態(tài)
  

長度是對(duì)象中元素的數(shù)目，用length()查看

對(duì)象類別

只有數(shù)據(jù)框和列表支持多種對(duì)象并存；因子只有兩種類型

向量是一個(gè)變量的取值；因子是一個(gè)分類變量；數(shù)組是一個(gè)k維的數(shù)據(jù)表；矩陣是數(shù)組的特例【數(shù)組或者矩陣的所有元素都是一種】；數(shù)據(jù)框是一個(gè)或幾個(gè)向量/因子構(gòu)成，并且它們必須等長，可以是不同的類型；列表可以包含任何類型的對(duì)象（也可以列表包含列表）

瀏覽對(duì)象

利用ls()可以查看當(dāng)前在內(nèi)存中的對(duì)象，但是這個(gè)函數(shù)只列出了對(duì)象名，并且是所有的；

想要查看名稱中含有某個(gè)指定字符（如x）的對(duì)象，可以指定pattern：ls(pat = "x");

想要看以某個(gè)字母開頭的對(duì)象，可以利用ls(pat = "^x") ;

如果想看所有對(duì)象的詳細(xì)信息呢？ls.str()

刪除對(duì)象

• rm(x,y)刪除對(duì)象x和y
• rm(list=ls()) 刪除所有對(duì)象
• rm(list=ls(pat="^x")) 刪除所有x開頭的對(duì)象

了解下R的向量

建立向量

數(shù)值型向量

• 向量沒規(guī)律：c()

• 向量的規(guī)律比較簡單：seq()或“:”

  > 1:10
  # 1 2 3 4 5 6 7 8 910
  > 1:10-1
  # 0 1 2 3 4 5 6 7 8 9
  > 1:(10-1) # 注意有括號(hào)與上面沒有的區(qū)別
  # 1 2 3 4 5 6 7 8 9
  > z <- seq(1,5,by=0.5)
  # 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
  > z <- seq(1,10,length=11)
  # 1.0 1.9 2.8 3.7 4.6 5.5 6.4 7.3 8.2 9.1 10.0
  

• 向量的規(guī)律比較復(fù)雜：rep()

  > z <- rep(2:5,2)
  # 2 3 4 5 2 3 4 5
  > z <- rep(2:5,each=2)
  # 2 2 3 3 4 4 5 5
  

• 通過鍵盤輸入向量scan()

  > z <- scan()
  1: 1 2 3 4 5 6 
  7: 
  # Read 6 items
  > z
  # 1 2 3 4 5 6
  

字符型向量

注意引號(hào)在輸入時(shí)應(yīng)該寫作：\" ；paste（）可以連接多個(gè)參數(shù)成為字符串，其中如果有數(shù)值，那么數(shù)值會(huì)被強(qiáng)制轉(zhuǎn)為字符串；默認(rèn)空格分割各個(gè)字符，使用sep=自定義

邏輯型向量

TRUE、FALSE可以簡寫成T、F；如果轉(zhuǎn)換為數(shù)值，F(xiàn)ALSE為0，TRUE為1

因子型向量

利用factor（）建立

factor(x, levels = sort(unique(x), na.last = TRUE),
         labels = levels, exclude = NA, ordered = is.ordered(x))

levels 指定因子水平（缺省值是向量x中不同的值）；
labels 指定水平名稱；
ordered是邏輯型選項(xiàng)，表示因子水平是否有順序

• 字符型因子轉(zhuǎn)為數(shù)值型因子【反之亦然】

> a <- c("x","y","z","x")
> a
[1] "x" "y" "z" "x"
> levels(a) <- c(1,2,3)
> a
[1] 1 2 3 1
Levels: 1 2 3

• 產(chǎn)生規(guī)則的因子序列，利用gl(k,n) ，k是水平數(shù)，n是重復(fù)數(shù)；另外還有兩個(gè)選項(xiàng)，length指定總共數(shù)據(jù)個(gè)數(shù)，label指定每個(gè)水平名稱

  > gl(2,3)
  [1] 1 1 1 2 2 2
  Levels: 1 2
  > gl(2,3,length=10)
  [1] 1 1 1 2 2 2 1 1 1 2
  Levels: 1 2
  > gl(2,3,labels = c("x","y"))
  [1] x x x y y y
  Levels: x y
  

向量提取

這里簡單說下根據(jù)邏輯值提取

# 例如x向量是這樣的
> x <- c(20, 15, 11, 6)
> x>10 # 這個(gè)僅僅是判斷，返回的也是F/T
[1]  TRUE  TRUE  TRUE FALSE
> x[x>10] # 提取大于10的元素
[1] 20 15 11
> x[ x < 15 & x > 10] # 限定多個(gè)提取條件
[1] 11

> y <- x[!is.na(x)] # 找到x中的非缺失值
> z <- x[(!is.na(x))&(x>0)] # 非負(fù)非缺失值

了解下R的矩陣 Matrix

相比于數(shù)組，矩陣使用頻率更高，構(gòu)建矩陣使用matrix

> x <- matrix(1:3, 2, 3) # 默認(rèn)按列填充
> x
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    1    3
# 如果改成按行填充，將 byrow=TRUE就好

矩陣元素的替換

> x[,3] <- NA
> x
     [,1] [,2] [,3]
[1,]    1    3   NA
[2,]    2    1   NA
> x[is.na(x)] <- 1 # 缺失值用1代替
> x
     [,1] [,2] [,3]
[1,]    1    3    1
[2,]    2    1    1

矩陣統(tǒng)計(jì)運(yùn)算

矩陣的排列是有方向性的，規(guī)定矩陣按列排列，一般不說明的時(shí)候，統(tǒng)計(jì)函數(shù)也是按列計(jì)算（但是可以用MARGIN來改變，等于1代表列，等于2代表行）

cov()與cor()分別計(jì)算矩陣的協(xié)方差矩陣和相關(guān)系數(shù)陣；

可以進(jìn)行標(biāo)準(zhǔn)化scale(x, center=T, scale=T) ；

按列求均值apply(x, MARGIN=2, FUN=mean)

了解下R的數(shù)據(jù)框 Data frame

雖然說數(shù)據(jù)框與矩陣很相似，也是二維表格，也是要求各個(gè)變量的觀測(cè)值長度相等，但是，在數(shù)據(jù)框中，行和列的意義是不一樣的， 其中列表示變量，行為觀測(cè)

提取

• 提取一列：數(shù)據(jù)框名$變量名
• 提取滿足條件子集：例如 subset(df, con1 == "treated" & con2 > 100)

添加

• 基本方法：df$val1 <- con df是數(shù)據(jù)框，val1是新變量，con是定義的val1元素
• 使用with：df$val1 <- with(df, con)
• transform一次增加多個(gè)變量：df <- transform(df, val1 = con, val2 = con2)

數(shù)據(jù)保存

一般用write.table()或者save()

• 保存為簡單的文本文件：
  write.table(df, file="/path/", row.names = F, quote = F)
  row.names =F 表示行名不寫入文件；quote = F表示變量名不放在雙引號(hào)中

• 保存為逗號(hào)分隔（csv）：用write.csv()

• 保存為R格式文件：
  save(df, file="/path/")

  或者保存全部save.image() ，它等價(jià)于save(list =ls(all=TRUE), file=".RData")

數(shù)據(jù)讀取

讀取數(shù)據(jù)集

R內(nèi)置的基本數(shù)據(jù)集有100多個(gè)（常為數(shù)據(jù)框和列表）。它們隨R的啟動(dòng)全部一次性自動(dòng)載入，通過命令data() 可以查看全部的數(shù)據(jù)集（也包含了通過library()加載的包中數(shù)據(jù)集）；使用data(package = "pkname") 可以列出包pkname中的所有數(shù)據(jù)集，但是可能還未被加載，確定要用的時(shí)候可以加載包

讀取Rdata

涉及多個(gè)數(shù)據(jù)集的分析時(shí)，最常使用load("/path/")

R的編程思維

• 控制結(jié)構(gòu)中，條件語句少不了，例如if (條件) 表達(dá)式1 else 表達(dá)式2 ，但是這樣有些冗余，可以利用ifelse(條件， 表達(dá)式1， 表達(dá)式2)

• 循環(huán)結(jié)構(gòu)中，知道終止條件用for（）；不知道運(yùn)行次數(shù)用while（）

• 一般將一組命令放在大括號(hào)中=》模塊化編程

  # 一般格式：函數(shù)名 = function(變量列表) 函數(shù)體
  # 例如一個(gè)畫圖函數(shù)
  > myfun <- function(x, y) {
       data <- read.table(x)
       plot(data$V1, data$V2, type="l")
       title(y)
  }
  

• 代碼要增加注釋，設(shè)置行前自動(dòng)縮進(jìn)

• 命令向量化（vectorization） ，就是將循環(huán)隱含在表達(dá)式中，例如條件語句可以用邏輯索引向量代替

  for (i in 1:length(x)){
       if (x[i] == b)
            x[i] <- 0
      else
          y[i] <- 1
  }
  # 這么長的函數(shù)，其實(shí)就是為了判斷x中元素是否為b，等于b就賦值0；不等就賦值1
  # 可以用下面代替
  x[x == b] <- 0
  x[x != b] <- 1
  

  向量化快的原因是：R中使用向量化會(huì)立即調(diào)用C語言，計(jì)算時(shí)比R快100倍

• 創(chuàng)建project來運(yùn)行項(xiàng)目，這樣一是可以方便讀取存儲(chǔ)；二是不同項(xiàng)目分類存放，方便查找更新

• 想要把腳本從頭到尾運(yùn)行一遍，用source()，記得添加絕對(duì)路徑

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

Welcome to our bioinfoplanet!