轉(zhuǎn)錄組測(cè)序的最直接目的赏酥，就是設(shè)法尋找組間顯著表達(dá)變化的基因赖淤，解釋基因表達(dá)水平的變化對(duì)生物功能的影響医增。例如慎皱，腫瘤患者與正常人群相比哪些編碼蛋白或非編碼RNA分子發(fā)生了失調(diào)，這些失調(diào)分子是否是引發(fā)或加速腫瘤進(jìn)程的潛在因素叶骨？逆境脅迫下的植物體中哪些基因表達(dá)顯著激活茫多，這些激活的基因是否有利于植物應(yīng)對(duì)高溫、干旱忽刽、鹽脅迫等不利環(huán)境天揖。
那么夺欲，如何基于轉(zhuǎn)錄組測(cè)序獲得的定量表達(dá)值，識(shí)別差異表達(dá)變化的基因或其它非編碼RNA分子呢今膊？
查閱相關(guān)文獻(xiàn)尋找思路的習(xí)慣是一定要養(yǎng)成的些阅，文獻(xiàn)的材料方法或者正文部分都會(huì)有相關(guān)的描述，使用了哪個(gè)工具（例如DESeq2）作差異表達(dá)分析斑唬，閾值是什么（例如倍數(shù)變化≥2市埋，p<0.05等）。本篇就以R包DESeq2的差異表達(dá)基因分析方法為例做個(gè)簡(jiǎn)單演示恕刘，這是目前使用頻率最高的鑒定差異分子的R包之一缤谎。
1 關(guān)于DESeq2包的安裝
對(duì)于DESeq2的安裝也很簡(jiǎn)單，一般情況下褐着，直接通過(guò)Bioconductor安裝DESeq2就可以了坷澡。
考慮到19年下半年的時(shí)候，DESeq2做過(guò)重大更新和優(yōu)化献起，核心算法沒(méi)變洋访，但是運(yùn)行速率相對(duì)于先前的版本提升了幾十倍！小編不清楚最新版本的DESeq2是否已經(jīng)添加至Bioconductor中谴餐，因此如果想安裝最新版本提高計(jì)算效率姻政，而又不確定Bioconductor的DESeq2是否是新版本的情況下，建議直接鏈接至作者的github庫(kù)中直接安裝岂嗓。
##以下兩種安裝 DESeq2 的方法二選一
#（1）常規(guī)方法 bioconductor 安裝
#install.packages('BiocManager') #需要首先安裝 BiocManager汁展，如果尚未安裝請(qǐng)先執(zhí)行該步
BiocManager::install('DESeq2')
#（2）在 GitHub 中獲取最新版本的方式
#install.packages('devtools') #需要首先安裝 devtools，如果尚未安裝請(qǐng)先執(zhí)行該步
devtools::install_github('mikelove/DESeq2@ae7c6bd')
#其它備注項(xiàng)：若中間提示有其它依賴(lài) R 包的舊版包沖突的話厌殉，先刪除舊包再安裝新的
remove.packages('xxx')
BiocManager::install('xxx')
2 DESeq2分析差異表達(dá)基因的一般過(guò)程
DESeq2是一種基于負(fù)二項(xiàng)式分布的方法食绿，使用局部回歸推算均值和方差，通過(guò)離散度和倍數(shù)變化的收縮率估計(jì)以提高穩(wěn)定性公罕。定量分析關(guān)注的更多是差異表達(dá)的“強(qiáng)度”器紧，而非“存在”。
以下是DESeq2分析差異表達(dá)基因的一般過(guò)程楼眷。
2.1 準(zhǔn)備數(shù)據(jù)和文件讀取
首先準(zhǔn)備基因表達(dá)值矩陣铲汪。
本文的所有測(cè)試數(shù)據(jù)和R代碼，可在文末獲取
-
“control_treat.count.txt”罐柳，是6個(gè)測(cè)序樣本的基因表達(dá)值矩陣掌腰，包括3個(gè)處理組（treat）和3個(gè)對(duì)照組（control）。第1列是基因名稱(chēng)张吉，注意不能有重復(fù)值齿梁。
示例文件樣式，行是基因，列是樣本勺择，內(nèi)容為基因表達(dá)值
然后將準(zhǔn)備好的基因表達(dá)值矩陣讀到R中创南，并預(yù)定義樣本分組信息。
#讀取基因表達(dá)矩陣
dat <- read.delim('control_treat.count.txt', row.names = 1, sep = '\t', check.names = FALSE)
#指定分組因子順序
#注意要保證表達(dá)矩陣中的樣本順序和這里的分組順序是一一對(duì)應(yīng)的
coldata <- data.frame(condition = factor(rep(c('control', 'treat'), each = 3), levels = c('control', 'treat')))
2.2 計(jì)算差異倍數(shù)列表
整體過(guò)程非常簡(jiǎn)單酵幕，因?yàn)樽髡咭呀?jīng)將多步過(guò)程整合到一個(gè)函數(shù)中扰藕，因此使用起來(lái)非常的方便。大致上芳撒，DESeq2兩步就完成分析了邓深。
- 第1步是構(gòu)建DESeqDataSet對(duì)象，包括對(duì)表達(dá)值的標(biāo)準(zhǔn)化以及存儲(chǔ)輸入值和中間結(jié)果等笔刹。
- 第2步芥备，函數(shù)DESeq()是一個(gè)包含因子大小估計(jì)、離散度估計(jì)舌菜、負(fù)二項(xiàng)模型擬合萌壳、Wald統(tǒng)計(jì)等多步在內(nèi)的過(guò)程，結(jié)果將返回至DESeqDataSet對(duì)象日月。最后獲得各基因的差異倍數(shù)變化和顯著性p值袱瓮，用于后續(xù)的差異基因篩選。
##DESeq2 默認(rèn)流程
library(DESeq2)
#第一步爱咬，構(gòu)建 DESeqDataSet 對(duì)象
dds <- DESeqDataSetFromMatrix(countData = dat, colData = coldata, design= ~condition)
#第二步尺借，計(jì)算差異倍數(shù)并獲得 p 值
#備注：parallel = TRUE 可以多線程運(yùn)行，在數(shù)據(jù)量較大時(shí)建議開(kāi)啟
dds1 <- DESeq(dds, fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)
#注意精拟，需將 treat 在前燎斩，control 在后，意為 treat 相較于 control 中哪些基因上調(diào)/下調(diào)
res <- results(dds1, contrast = c('condition', 'treat', 'control'))
res
差異分析結(jié)果保存在“res”中蜂绎，包含了基因id栅表、標(biāo)準(zhǔn)化后的基因表達(dá)值、log2轉(zhuǎn)化后的差異倍數(shù)（Fold Change）值师枣、顯著性p值以及校正后p值（默認(rèn)FDR校正）等主要信息怪瓶。
如果期望將該表格輸出到本地，轉(zhuǎn)化為數(shù)據(jù)框結(jié)構(gòu)后直接write.table()就可以了践美。
#輸出表格至本地
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
write.table(res1, 'control_treat.DESeq2.txt', col.names = NA, sep = '\t', quote = FALSE)
2.3 篩選差異表達(dá)基因
后續(xù)通過(guò)該表劳殖，即可自定義閾值篩選差異表達(dá)基因了。
- 例如這里期望根據(jù)|log2FC|≥1以及padj<0.01篩選拨脉，并通過(guò)“up”和“down”分別區(qū)分上、下調(diào)的基因宣增。
##篩選差異表達(dá)基因
#首先對(duì)表格排個(gè)序玫膀，按 padj 值升序排序，相同 padj 值下繼續(xù)按 log2FC 降序排序
res1 <- res1[order(res1$padj, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
#log2FC≥1 & padj<0.01 標(biāo)識(shí) up爹脾，代表顯著上調(diào)的基因
#log2FC≤-1 & padj<0.01 標(biāo)識(shí) down帖旨，代表顯著下調(diào)的基因
#其余標(biāo)識(shí) none箕昭，代表非差異的基因
res1[which(res1$log2FoldChange >= 1 & res1$padj < 0.01),'sig'] <- 'up'
res1[which(res1$log2FoldChange <= -1 & res1$padj < 0.01),'sig'] <- 'down'
res1[which(abs(res1$log2FoldChange) <= 1 | res1$padj >= 0.01),'sig'] <- 'none'
#輸出選擇的差異基因總表
res1_select <- subset(res1, sig %in% c('up', 'down'))
write.table(res1_select, file = 'control_treat.DESeq2.select.txt', sep = '\t', col.names = NA, quote = FALSE)
#根據(jù) up 和 down 分開(kāi)輸出
res1_up <- subset(res1, sig == 'up')
res1_down <- subset(res1, sig == 'down')
write.table(res1_up, file = 'control_treat.DESeq2.up.txt', sep = '\t', col.names = NA, quote = FALSE)
write.table(res1_down, file = 'control_treat.DESeq2.down.txt', sep = '\t', col.names = NA, quote = FALSE)
2.4 差異表達(dá)火山圖示例
已經(jīng)識(shí)別了顯著差異表達(dá)的基因，最后期望通過(guò)火山圖將它們表示出來(lái)解阅，火山圖是文獻(xiàn)中常見(jiàn)統(tǒng)計(jì)圖表之一落竹。
- ggplot2為此提供了優(yōu)秀的作圖方案，參考以下示例货抄。
##ggplot2 差異火山圖
library(ggplot2)
#默認(rèn)情況下述召，橫軸展示 log2FoldChange，縱軸展示 -log10 轉(zhuǎn)化后的 padj
p <- ggplot(data = res1, aes(x = log2FoldChange, y = -log10(padj), color = sig)) +
geom_point(size = 1) + #繪制散點(diǎn)圖
scale_color_manual(values = c('red', 'gray', 'green'), limits = c('up', 'none', 'down')) + #自定義點(diǎn)的顏色
labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', title = 'control vs treat', color = '') + #坐標(biāo)軸標(biāo)題
theme(plot.title = element_text(hjust = 0.5, size = 14), panel.grid = element_blank(), #背景色蟹地、網(wǎng)格線积暖、圖例等主題修改
panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') + #添加閾值線
geom_hline(yintercept = 2, lty = 3, color = 'black') +
xlim(-12, 12) + ylim(0, 35) #定義刻度邊界
p
