第一，這一節(jié)講什么？

? 這一節(jié)講述一些關(guān)于R語言處理數(shù)據(jù)框的基本操作，會詳細(xì)演示上一節(jié)常用函數(shù)的實(shí)際使用方法

? 根據(jù)我平時(shí)工作的經(jīng)驗(yàn)，只將一些最常用的操作介紹給大家

第二，這一節(jié)如何講？

? 我會提供一些示例數(shù)據(jù)，然后大家根據(jù)我的講述，一一朝下操作

? 在操作過程中，我會強(qiáng)調(diào)自己認(rèn)為一些比較重要的點(diǎn)

第三，這一節(jié)的目的？

? 演示，當(dāng)獲得最基本的差異分析結(jié)果后，如何快速篩選出自己的目標(biāo)基因

數(shù)據(jù)格式說明

生信分析的大部分結(jié)果是表格形式，這個(gè)表格文件的格式一般為txt文檔。

這些txt文檔的列通常為tab制表符分割，如果使用Notepad++打開文檔的話，你可能會看到類似下面的結(jié)果

image

不同列之間的→即上面所說tab制表符，純文本形式表示為\t

制表符一般是不可見，同樣不可見的常用字符還有每行結(jié)束時(shí)的換行符\n, 以及空格

Notepad++是一個(gè)常用的文本編輯器，在視圖設(shè)置里可以將所有字符顯示出來，例如在上圖中，空格和制表符分別被顯示為·和→

稍微啰嗦一下，為什么生信的結(jié)果多以txt文檔和csv文檔(以,分割的表格)？

原因有兩個(gè)：

1. excel文檔有行數(shù)上限，對2003版最大行數(shù)是65536行，對2007以上版本，最大行數(shù)是1048676行；生信分析結(jié)果數(shù)據(jù)量較大，有時(shí)候會有限制

2. excel不是純文本格式（主要原因），生信分析基本都在Linux操作系統(tǒng)下完成，而excel文件在Linux系統(tǒng)下是無法直接處理的

   寫到這里，我實(shí)在忍不住想要吐槽一下微軟，這真是一個(gè)神奇的公司……生信人員必須掌握的技能表中，dos2unix一定是名列前茅的……感覺有點(diǎn)莫名其妙的同學(xué)請?zhí)^這段吐槽哈~~~

對于大多數(shù)非生信專業(yè)同學(xué)，可以鼠標(biāo)右鍵選擇excel方式打開；或者直接將文件后綴修改為xls，然后就可以直接用excel打開了。但是請大家務(wù)必注意，excel會自動(dòng)將時(shí)間類的字符進(jìn)行轉(zhuǎn)化，例如在excel中輸入Sep 2字符會自動(dòng)變?yōu)?code>2-Sep。

測試數(shù)據(jù)說明

數(shù)據(jù)存放在百度網(wǎng)盤test_data目錄下，鏈接: https://pan.baidu.com/s/1iNo6_ni9mDaNzadE__ZNMg, 提取碼: y3bm（上次準(zhǔn)備的lncRNA相關(guān)數(shù)據(jù)，有的同學(xué)反映對lncRNA不是很熟悉，所以我更改為gene的數(shù)據(jù)類型）

Case-vs-Normal.xls是人細(xì)胞的差異表格示例，Case和Normal都是三個(gè)生物學(xué)重復(fù)

表格里面提供了差異比較的倍數(shù)和P值，以及每一個(gè)基因的GO和KEGG注釋信息

這種表格應(yīng)該是非常常見的，一般在生信公司做完基本的RNAseq分析后，都可以拿到類似的數(shù)據(jù)

那么接下來，我便給大家演示一下，如何通過R快速篩選出自己感興趣的基因

細(xì)心的同學(xué)可能會發(fā)現(xiàn)，測試數(shù)據(jù)里面有一些比較奇怪的數(shù)據(jù)

image

? 有的單元格是空的，有的單元格內(nèi)容為NA，這些都表示此處無內(nèi)容

? 還有的單元格數(shù)據(jù)為Inf，這代表無窮大

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

數(shù)據(jù)框讀取的時(shí)候，大家可能會比較熟悉read.table函數(shù)。但是我比較推薦read.delim函數(shù)，讀取我上面說的生信格式文件會更加方便。

在下面的案例中，可以看到，read.table函數(shù)會報(bào)錯(cuò)。

另外一點(diǎn)，雖然read.delim函數(shù)默認(rèn)參數(shù)sep='\t',header = T（列和列之間的分隔符，第一行是否為表頭），但是我每次都會自己手動(dòng)寫出來。這樣會增加自己對于即將處理的數(shù)據(jù)的熟悉度，算是個(gè)人一個(gè)小的工作習(xí)慣。

其實(shí)還有一個(gè)參數(shù)是check.names=F，后面會講

> df <- read.table('Case-vs-Normal.xls',sep='\t',header = T)
Warning messages:
1: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string
2: In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  number of items read is not a multiple of the number of columns
  line 10 did not have 10 elements
           
> df <- read.delim('Case-vs-Normal.xls',sep='\t',header = T)

為什么此處read.table使用會出錯(cuò)，我暫時(shí)不做過多解釋。如果有人比較糾結(jié)這個(gè)事情的話，可以嘗試一下下面的代碼。

df <- read.table('Case-vs-Normal.xls',sep='\t',header = T, quote = '')

點(diǎn)擊數(shù)據(jù)框預(yù)覽，發(fā)現(xiàn)怎么有的列名變了？比如原來是FPKM_Case-1，現(xiàn)在變?yōu)?code>FPKM_Case.1。大家注意一下，R在讀取數(shù)據(jù)框是，會自動(dòng)將中橫杠-轉(zhuǎn)換為點(diǎn).，添加check.names=F，可以避免這種情況。

image

所以，在數(shù)據(jù)框讀取時(shí)，我自己一般習(xí)慣的讀取方式為

df <- read.delim('Case-vs-Normal.xls',sep='\t',header = T, check.names=F)

你可以再預(yù)覽一下，看看是不是正常了？

還有啰嗦的一個(gè)點(diǎn)是，大家仔細(xì)看看一下數(shù)據(jù)框預(yù)覽的結(jié)果就會發(fā)現(xiàn)，不同列中空值結(jié)果的表示是不一樣的

在foldChange和pval兩列中，無論原始數(shù)據(jù)是NA還是空的，此處都是NA

GO和GO_description兩列中，空值沒有任何元素

但是在KEGG列第一行中，空值表示為NA

大家記住

1.純數(shù)字列，無論原始數(shù)據(jù)是什么形式，讀取的結(jié)果，空值都是NA；字符串列，空值是和其原始數(shù)據(jù)形式統(tǒng)一

2.GO和GO_description兩列中的空值，其實(shí)本質(zhì)上空字符串，代碼表示為""

3.""和NA在進(jìn)行空值剔除時(shí)會有影響，這個(gè)后面碰到了再講

image

數(shù)據(jù)框基本操作

快速瀏覽

讀取完表格之后，我們通常會查看一下數(shù)據(jù)框的基本信息

1.通過summary函數(shù)查看整體情況

summary(df)

image

仔細(xì)看一下結(jié)果，可以發(fā)現(xiàn)summary返回的結(jié)果分為兩種

一種是針對純數(shù)字的列，例如foldChange列，返回值分別是該列最小值、第1四分位數(shù)、中位值、均值、第3四分位數(shù)、最大值等信息，同時(shí)還統(tǒng)計(jì)了空值NA的個(gè)數(shù)

在這個(gè)案例中，由于foldChange列中存在Inf，所以最大值和均值均為無窮大

另一種是針對含有字符串的列，例如GO列，返回值為各個(gè)字符串的頻數(shù)統(tǒng)計(jì)，并且是由大到小依次展示

在這個(gè)案例中，出現(xiàn)次數(shù)最多的空值，有4727次；其次是GO:0016021，有377次

summary只顯示6行結(jié)果，多余的內(nèi)容都用other來代替了

2.如果想要快速瀏覽數(shù)據(jù)，可以通過head(df)查看數(shù)據(jù)庫前6行，或者head(df, 10)查看前10行

由于結(jié)果過多，就不展示圖片了

3.查看數(shù)據(jù)框維度

> dim(df) # 查看行列數(shù)
[1] 20030    13
> nrow(df) # 查看行數(shù)
[1] 20030
> ncol(df) # 查看列數(shù)
[1] 13

我一般就使用dim（），dim(df)[1]是行數(shù)，dim(df)[2]是列數(shù)，不同記那么多函數(shù)名稱了

代碼后面的#表示注釋內(nèi)容，

注釋內(nèi)容是不會被運(yùn)行的

4.行列名查看

rownames(df) # 查看行名
colnames(df) # 查看列名

剛剛我們讀取表格的時(shí)候，只用header參數(shù)定義的列名，

所以行名默認(rèn)是從1開始的數(shù)字

切片、提取子集

在R里面，數(shù)據(jù)框操作的基本原則為：（我試著寫著順口溜哈）

逗號(,)間隔行和列，行的操作必有逗(,)

單個(gè)數(shù)字處理列，知道名稱也能做

冒號:對應(yīng)連續(xù)值，間斷數(shù)據(jù)要c括(c())

列的提取樣式多，數(shù)據(jù)格式不一樣

1.提取某個(gè)指定元素

> df[5,6] # 提取第5行第6列元素
[1] 0
> df[8,9] # 提取第8行第9列元素
[1] 0.9832643

2.提取整行數(shù)據(jù)

#### 提取單行數(shù)據(jù)
df[2,] # 提取第2行

##### 提取多行數(shù)據(jù)
df[1:4,] # 提取1-4行
df[c(2,4,5),] # 提取2、4、5行

3.提取整列數(shù)據(jù)

對于數(shù)據(jù)框而言，提取列的方式有多種

#### 按照位置信息提取列
df[,1] # 提取第1列
df[1] # 提取第1列

#### 按照列名提取列
df['gene_id'] #提取gene_id列
df$gene_id #提取gene_id列

#### 提取多列
df[2:6] # 提取2-6列
df[c(5,8,10)] # 提取5、8、10列

但是上面不同方式提取出來的結(jié)果，其數(shù)據(jù)類型有所區(qū)別，可以通過class函數(shù)來查看相應(yīng)數(shù)據(jù)類型

> class(df[1]) 
[1] "data.frame"
> class(df[,1]) 
[1] "factor"
> class(df['gene_id']) 
[1] "data.frame"
> class(df$gene_id) 
[1] "factor"

對于普通用戶，僅進(jìn)行簡單的數(shù)據(jù)處理或繪圖，上面幾種方法可以隨便使用，

如果以后有高級分析需求，例如做差異分析等復(fù)雜處理，需要注意一下數(shù)據(jù)結(jié)構(gòu)的問題

這個(gè)我們以后碰到了再說

4.提取子集

df[1:4,2:6] #提取第1-4行,2-6列
df[c(1,5,10),c(2,9,10)] # 提取第1、5、10行,第2、9、10列
df[1:4,c('gene_id','pval')] # 提取'gene_id'和’pval'兩列的第1-4行

數(shù)據(jù)處理

數(shù)據(jù)框基本操作還有很多，就不過多講解了，

下面直接進(jìn)入實(shí)際操作——數(shù)據(jù)篩選，

將其他常見的數(shù)據(jù)框操作融合進(jìn)實(shí)際案例處理，

這樣能夠讓大家加深理解

去除NA

前面提到過，foldChange和pval兩列中有NA值

這說明差異比較的時(shí)候，Normal的表達(dá)值為0(分母為0，沒有意義)

所以我們首先要將這些基因剔除了，可以通過na.omit()函數(shù)處理

> dim(df)
[1] 20030    13
> dim(na.omit(df)) # 剔除含有NA的整行內(nèi)容
[1] 17161    13

注意，na.omit（）只剔除含有NA的行

上面提到那些GO和GO_description等列，雖然有空值，但是na.omit()無法處理這些數(shù)據(jù)

df_delNA <- na.omit(df)
head(df_delNA) # 可以看到，GO和GO_description等列還是有空值

如果執(zhí)行df <- na.omit(df)的話，會直接替換df數(shù)據(jù)本身

通常，我會將處理過后的數(shù)據(jù)保存進(jìn)新的變量，如df_delNA，

這樣的好處是不會改變原始數(shù)據(jù)，如果那里有問題了，便于重新處理數(shù)據(jù)

條件篩選

直接篩選

剔除完空行之后，我們希望獲得差異倍數(shù)較大（foldChange絕對值大于等于2），且極顯著(pval小于等于0.01)的基因

首先，進(jìn)行顯著性篩選

> head(df_delNA$pval <= 0.01)
[1] FALSE FALSE FALSE FALSE FALSE FALSE
> dim(df_delNA[df_delNA$pval <= 0.01,])
[1] 145  13

df_delNA$pval <= 0.01會將該列的數(shù)值同0.01一一比較，符合條件的話返回TRUE，否則返回FALSE

這些TRUE或FALSE組成的列表對應(yīng)不同行，可以將其作為行的篩選條件

最后，我們還是將篩選出來的結(jié)果保存為新的變量

df_delNA_0.01 <- df_delNA[df_delNA$pval <= 0.01,]

然后，對差異倍數(shù)進(jìn)行篩選

篩選條件foldChange絕對值大于等于2，可以拆分為foldChange大于等于2或foldChange小于等于-2

用代碼表示就是

df_delNA_0.01$foldChange >=2 | df_delNA_0.01$foldChange <=-2

在R里面，多個(gè)條件通過&(并且)和|(或)來表示

另外一種表示方式，可以直接用絕對值函數(shù)abs()來處理

abs(df_delNA_0.01$foldChange) >=2

這兩種的結(jié)果是一樣的

> dim(df_delNA_0.01[df_delNA_0.01$foldChange >= 2 | 
+                       df_delNA_0.01$foldChange <= -2, ])
[1] 58 13
> 
> dim(df_delNA_0.01[abs(df_delNA_0.01$foldChange) >= 2, ])
[1] 58 13

1.在R里面撰寫代碼的時(shí)候，如果一條命令太長了，可以將其分為多行來撰寫

雖然df_delNA_0.01$foldChange >= 2 | df_delNA_0.01$foldChange <= -2分為兩行撰寫了，

但是中間|會將其連接為一個(gè)整體

分行撰寫會增加代碼的可讀性

2.通過dim(df_delNA_0.01[abs(df_delNA_0.01$foldChange) >= 2, ])可以看出

在R里面，函數(shù)是可以直接疊加使用的

上面的代碼邏輯是

先abs()求該列的絕對值，然后判斷絕對值是否大于2

其次將判斷結(jié)果作為行篩選條件賦給數(shù)據(jù)框，

最后dim()計(jì)算最終篩選結(jié)果的行列數(shù)

先處理，再篩選

上面展示的是，我們根據(jù)表格中的原始數(shù)據(jù)，設(shè)置篩選條件，縮減表格

但是也有很多情況是，我們需要先對原始數(shù)據(jù)進(jìn)行一定的理，然后根據(jù)處理后的結(jié)果，再進(jìn)行篩選

比如我們希望拿到實(shí)驗(yàn)組中低豐度(FPKM小于1)，但是差異卻極顯著(pval小于0.01)的基因

那么首先，我們需要計(jì)算實(shí)驗(yàn)組的FPKM均值

實(shí)驗(yàn)組的數(shù)據(jù)在哪里呢？

> colnames(df_delNA)
 [1] "gene_id"          "FPKM_Case-1"      "FPKM_Case-2"     
 [4] "FPKM_Case-3"      "FPKM_Normal-1"    "FPKM_Normal-2"   
 [7] "FPKM_Normal-3"    "foldChange"       "pval"            
[10] "GO"               "GO_description"   "KEGG"            
[13] "KEGG_description" "caseMean"  

數(shù)一下，實(shí)驗(yàn)組(Case)位于2-4列，所以計(jì)算他們的均值就是

df_delNA$caseMean <- rowMeans(df_delNA[2:4])

1.df_delNA[2:4]直接提取數(shù)據(jù)框df_delNA第2-4列

但是我一般會寫成df_delNA[,2:4]

這樣，我會很清楚的知道，我在操作列

2.rowMeans函數(shù)可以計(jì)算數(shù)據(jù)框的行均值，對應(yīng)的colMeans可以計(jì)算列均值

類似的函數(shù)還有rowSums、colSums等

3.df_delNA$caseMean 會在數(shù)據(jù)框df_delNA中直接添加新列caseMean

> colnames(df_delNA)
 [1] "gene_id"          "FPKM_Case-1"      "FPKM_Case-2"     
 [4] "FPKM_Case-3"      "FPKM_Normal-1"    "FPKM_Normal-2"   
 [7] "FPKM_Normal-3"    "foldChange"       "pval"            
[10] "GO"               "GO_description"   "KEGG"            
[13] "KEGG_description" "caseMean"      

聰明的同學(xué)可能已經(jīng)發(fā)現(xiàn)，既然實(shí)驗(yàn)組的名稱里面有Case關(guān)鍵字，那么是否有快捷的方法？

上一節(jié)常用函數(shù)里面，我們提過一個(gè)grep函數(shù)，大家可以試一下

df_delNA$caseMean <- rowMeans(df_delNA[,grep('Case',colnames(df_delNA))])

grep可以進(jìn)行關(guān)鍵詞搜索，返回匹配的位置信息

如果忘了的同學(xué)可以回顧一下上一節(jié)的內(nèi)容

> grep('Case',colnames(df_delNA))
[1] 2 3 4

接下來篩選就和上面講的內(nèi)容差不多了

df_delNA_lowEx <- df_delNA[df_delNA$caseMean <=1, ]
df_delNA_lowEx_sig <- df_delNA_lowEx[df_delNA_lowEx$pval <= 0.01,]

最終，我們篩選到63個(gè)在實(shí)驗(yàn)組中表達(dá)豐度低，但卻具有顯著調(diào)控作用的基因

> dim(df_delNA_lowEx_sig)
[1] 63 14