你是不是也這樣——
拿到一個(gè)RNA-seq表達(dá)矩陣，想干幾件再正常不過(guò)的事：
——只要p值小于0.05的基因；
——按表達(dá)量從高到低排個(gè)序；
——把重復(fù)的基因行去掉；
——按分組求個(gè)平均表達(dá)量。
在Excel里，你得先點(diǎn)篩選按鈕，再點(diǎn)排序按鈕，再用數(shù)據(jù)透視表，一操作就是十分鐘，關(guān)鍵下回打開(kāi)很可能忘記當(dāng)時(shí)怎么篩的，老板一問(wèn)沒(méi)法復(fù)現(xiàn)。
今天師姐就給你一套 dplyr 五大金剛，讓你在R里面用幾行代碼，把這些洗數(shù)據(jù)的事干得明明白白、全程可追溯。
看完這篇，你會(huì)：
——用 filter 做任何你想要的篩選；
——用 arrange 排序排得明明白白；
——用 select 挑想要的列；
——用 distinct 一鍵去重；
——用 mutate + summarise 加新列、做分組統(tǒng)計(jì)。
01 先裝包，別拿base R硬剛
dplyr是tidyverse宇宙的核心成員。還沒(méi)有的同學(xué)先裝：
install.packages("dplyr")
library(dplyr)
師姐后面所有的演示，都基于一個(gè)你熟悉的迷你表達(dá)矩陣：
expr <- data.frame(
Gene = c("TP53", "BRCA1", "EGFR", "TP53", "MYC"),
Sample = c("Tumor1", "Tumor1", "Tumor2", "Tumor2", "Tumor1"),
Expression = c(120, 0, 200, 130, 450),
pvalue = c(0.01, 0.03, 0.001, 0.04, 0.06),
stringsAsFactors = FALSE
)
expr
假設(shè)這里面出了一些常見(jiàn)問(wèn)題：TP53重復(fù)了、有p值大于0.05的、有表達(dá)量為0的異常值。我們一個(gè)一個(gè)來(lái)處理。
02 filter —— 你想要什么行，它就給你什么行
場(chǎng)景一：只要顯著差異的基因（p < 0.05）
sig_genes <- filter(expr, pvalue < 0.05)
sig_genes
就這么簡(jiǎn)單，比你在Excel里點(diǎn)“篩選→數(shù)字篩選→小于”，再OK一下快好幾秒，而且完全可重復(fù)。
場(chǎng)景二：多條件組合
我要 “Tumor1” 樣本里，表達(dá)量大于100的基因：
filter(expr, Sample == "Tumor1", Expression > 100)
逗號(hào)就是“且”。如果想用“或”，就用 | 符號(hào)：
filter(expr, Gene == "TP53" | Gene == "MYC")
場(chǎng)景三：排除條件
去掉表達(dá)量為0的行：
filter(expr, Expression != 0)
學(xué)完這一招，你再也不用拿眼睛一行一行掃幾萬(wàn)行數(shù)據(jù)了。
03 arrange —— 想怎么排，就怎么排
場(chǎng)景：按表達(dá)量從高到低排序
arrange(expr, desc(Expression))
不加 desc 默認(rèn)是從小到大。想多級(jí)排序？輕松：
arrange(expr, Sample, desc(Expression))
先按樣本名字排序，同一個(gè)樣本內(nèi)部再按表達(dá)量降序。Excel里你得多步排序，R里一行搞定。
04 select —— 只挑你想要的列
場(chǎng)景：我只想保留基因名和表達(dá)量?jī)闪?，p值暫時(shí)不看了
select(expr, Gene, Expression)
還可以按條件挑列，比如只要數(shù)值型列：
select(expr, where(is.numeric))
或者剔除某列：
select(expr, -pvalue)
05 distinct —— 一鍵去重，比Excel刪除重復(fù)項(xiàng)快
場(chǎng)景：TP53重復(fù)了，每個(gè)基因只保留第一次出現(xiàn)的行
expr_unique <- distinct(expr, Gene, .keep_all = TRUE)
expr_unique
.keep_all = TRUE 意思是保留所有列。
06 mutate + summarise —— 加新列和分組統(tǒng)計(jì)，數(shù)據(jù)清洗的終極大招
先看 mutate，加一列l(wèi)og2轉(zhuǎn)換后的表達(dá)量：
expr <- mutate(expr, Log2Expr = log2(Expression + 1))
注意加1是防止有0取不了對(duì)數(shù)。
再看 summarise，求所有基因的平均表達(dá)量：
summarise(expr, MeanExpr = mean(Expression))
但這個(gè)不太實(shí)用，因?yàn)榘阉袠颖舅谢蚯罅藗€(gè)總平均。真正厲害的是分組統(tǒng)計(jì)。
分組統(tǒng)計(jì)：每個(gè)基因在不同樣本中的平均表達(dá)量
這次用管道符 %>% 把操作串起來(lái)讀，更符合思考順序：
expr %>%
group_by(Gene) %>%
summarise(AvgExpr = mean(Expression))
如果還想知道每組有多少個(gè)觀測(cè)值：
expr %>%
group_by(Gene) %>%
summarise(AvgExpr = mean(Expression), Count = n())
n() 是dplyr里統(tǒng)計(jì)行數(shù)的函數(shù)，非常常用。
今日總結(jié)：師姐送你一套清洗模板
以后拿到表達(dá)矩陣 df，標(biāo)準(zhǔn)清洗流程直接套：
library(dplyr)
clean_df <- expr %>%
filter(pvalue < 0.05) %>% # 只要顯著基因
distinct(Gene, .keep_all = TRUE) %>% # 去重
mutate(Log2Expr = log2(Expression + 1)) %>% # 加新列
arrange(desc(Expression)) # 排序
從今往后，80%的數(shù)據(jù)清洗交給這五大金剛，Excel只留給老板看最終結(jié)果。
下期預(yù)告：
數(shù)據(jù)洗干凈了，接下來(lái)當(dāng)然要畫(huà)圖。下期師姐帶你扔掉GraphPad，用R畫(huà)出第一張SCI級(jí)別的箱線圖。繪圖系列即將開(kāi)啟，千萬(wàn)別走丟！
這里是兩棲生物手冊(cè)，中科院生物醫(yī)學(xué)在讀博士，記錄生物醫(yī)學(xué)實(shí)驗(yàn)和生物信息學(xué)筆記，全平臺(tái)同名，正在努力成為貼心負(fù)責(zé)的RNA-seq數(shù)據(jù)分析者和實(shí)驗(yàn)技術(shù)分享者~
