簡單的轉錄組差異基因表達分析 -- DESeq2

經(jīng)典的轉錄組差異分析通常會使用到三個工具limma/voom, edgeR和DESeq2揣云。今天我們就通過一個小規(guī)模的轉錄組測序數(shù)據(jù)來演示DESeq2的簡單流程际起。

? 文中使用的數(shù)據(jù)來自Standford 大學的一個擬南芥的small RNA-seq數(shù)據(jù)（https://bios221.stanford.edu/data/mobData.RData）

該數(shù)據(jù)包含6個樣本：SL236奸柬，SL260，SL237湿蛔，SL238棋嘲，SL239店溢，SL240叁熔，并分成了三組，分別是：

? MM="triple mutatnt shoot grafted onto triple mutant root"

? WM="wild-type shoot grafted onto triple mutant root"

? WW="wild-type shoot grafted onto wild-type root"

簡而言之床牧，WW組可以認為是實驗的對照組荣回，而MM和WM則是兩個處理組。

對于DESeq2的分析流程而言戈咳，我們需要輸入的數(shù)據(jù)包括：

表達矩陣（countData）
分組信息（colData）
擬合信息（design）：指明如何根據(jù)樣本的分組進行建模

? 下面就以mobData 中的數(shù)據(jù)為例簡單介紹DESeq2 的分析流程

載入數(shù)據(jù)及生成DESeqDataSet

? 由于mobData 中的行名沒有提供基因的ID心软，我們也不是為了探究生物學問題，就以mobData 的行數(shù)作為其ID

library(DESeq2)
library(dplyr)
load("data/mobData.RData")
head(mobData)
##      SL236 SL260 SL237 SL238 SL239 SL240
## [1,]    21    52     4     4    86    68
## [2,]    18    21     1     5     1     1
## [3,]     1     2     2     3     0     0
## [4,]    68    87   270   184   396   368
## [5,]    68    87   270   183   396   368
## [6,]     1     0     6    10     6    12
dim(mobData)
## 3000    6

row.names(mobData) <- as.character(c(1:dim(mobData)[1]))
mobGroups <- c("MM", "MM", "WM", "WM", "WW", "WW")
colData <- data.frame(row.names=colnames(mobData),
                      group_list=mobGroups)
dds <- DESeqDataSetFromMatrix(countData = mobData,
                              colData = colData,
                              design = ~group_list,
                              tidy = F)
dds
## class: DESeqDataSet 
## dim: 3000 6 
## metadata(1): version
## assays(1): counts
## rownames(3000): 1 2 ... 2999 3000
## rowData names(0):
## colnames(6): SL236 SL260 ... SL239 SL240
## colData names(1): group_list

DESeqDataSet 是DESeq2流程中儲存read counts和中間統(tǒng)計分析數(shù)據(jù)的對象著蛙，之后的分析都建立在該對象之上進行删铃。

DESeqDataSet 可以通過以下四種方法產(chǎn)生：

From transcript abundance files and tximport —— DESeqDataSetFromTximport

From a count matrix —— DESeqDataSetFromMatrix

From htseq-count files —— DESeqDataSetFromHTseqCount

From a SummarizedExperiment object —— DESeqDataSet

預處理

? 在進行差異分析之前，需要對樣本數(shù)據(jù)的表達矩陣進行預處理踏堡，包括：

去除低表達量基因
探索樣本分組信息 -- 有助于挖掘潛在的差異樣本

table(rowSums(counts(dds)==0))
## 
##    0    1    2    3    4    5    6 
## 1833  442  324  170   92   84   55
keep <- rowSums(counts(dds)==0) < 6
dds.keep <- dds[keep, ]
dim(dds.keep)
## [1] 2945    6

library(factoextra)
library(FactoMineR)
count.pca <- counts(dds.keep) %>% as.matrix %>% t %>% PCA(.,graph = F)
fviz_pca_ind(count.pca,
             col.ind = mobGroups, 
             addEllipses = F, 
             mean.point = F,
             legend.title = "Groups")

colnames(dds.keep) <- c(paste(colnames(dds.keep),mobGroups, sep =".")) 
counts(dds.keep) %>% as.matrix %>% t %>% dist %>% hclust %>% plot
colnames(dds.keep) <- colnames(dds)

? 通過PCA結果來看各組樣本分組情況還是不錯的猎唁，但hclust的聚類結果反映的分組就略微有點混雜了，可能要聚類計算的距離函數(shù)選用不當有關暂吉。

差異分析

使用DESeq()函數(shù)進行差異分析時胖秒，該函數(shù)干了以下三件事：

計算文庫大小，預測用于標準化的size factor
預測離散度（dispersions）
使用負二項分布的廣義線性回歸模型（GLM）進行擬合慕的，以及進行Wald test

The Wald test (also called the Wald Chi-Squared Test) is a way to find out if explanatory variables in a model are significant.

ref : https://www.statisticshowto.com/wald-test/

degobj <- DESeq(dds.keep);degobj
## class: DESeqDataSet 
## dim: 2945 6 
## metadata(1): version
## assays(4): counts mu H cooks
## rownames(2945): 1 2 ... 2999 3000
## rowData names(26): baseMean baseVar ... deviance maxCooks
## colnames(6): SL236 SL260 ... SL239 SL240
## colData names(2): group_list sizeFactor
countdds <- counts(degobj, normalized = T)
deg.res1 <- results(degobj, contrast=c("group_list","MM","WW"), pAdjustMethod = 'fdr')

? counts() 可以提取DESeq object中的表達矩陣，而results() 可以提取差異分析的結果挤渔，其中包括了：

樣本間的均值, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values.

? 使用results()函數(shù)時需要指明進行比較的樣本肮街，這里用 contrast=c("group_list","MM","WW") 提取MM組和WW組進行差異分析的結果。如果想要比較WM組和WW組判导，只要改變contrast=c("group_list","WM","WW") 即可嫉父。

# result default
results(object, contrast, name, lfcThreshold = 0,
    altHypothesis = c("greaterAbs", "lessAbs", "greater", "less"),
    listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE,
    alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun,
    format = c("DataFrame", "GRanges", "GRangesList"), test,
    addMLE = FALSE, tidy = FALSE, parallel = FALSE,
    BPPARAM = bpparam(), minmu = 0.5)

? 檢查結果中是否包含NA值

dim(deg.res1)
## [1] 2945    6
apply(deg.res1,2, function(x) sum(is.na(x)))
##       baseMean log2FoldChange          lfcSE           stat         pvalue 
##              0              0              0              0              0 
##           padj 
##           1142

? 這里padj中有1142個NA值是因為使用results()提取差異分析結果時沛硅，大于alpha值（這里是0.1）的矯正后p-value都會被當做是NA。因此绕辖，我們將這些padj值都設為1

deg.res1$padj[is.na(deg.res1$padj)] <- 1
table(is.na(deg.res1$padj))
## 
## FALSE 
##  2945
summary(deg.res1)
## 
## out of 2945 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 102, 3.5%
## LFC < 0 (down)     : 215, 7.3%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 4)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotMA(deg.res1)

使用MA-plot查看fold change與read counts的關系

排序后以log2FoldChange絕對值大于1摇肌， padj小于0.05為條件篩選顯著的差異表達基因

deg.ord1 <- deg.res1[order(deg.res1$padj, decreasing = F),]
deg.sig1 <- subset(deg.ord1,abs(deg.ord1$log2FoldChange)>1 & deg.ord1$padj<0.05)
deg.sig1
## log2 fold change (MLE): group_list MM vs WW 
## Wald test p-value: group_list MM vs WW 
## DataFrame with 217 rows and 6 columns
##              baseMean    log2FoldChange             lfcSE              stat
##             <numeric>         <numeric>         <numeric>         <numeric>
## 2500 159.774533763595 -4.11069945121082 0.432943287486228 -9.49477580557611
## 1212 112.711030407211  3.96456274705605 0.449008108800994  8.82960166942819
## 1917  303.58162963147  3.45479344900962 0.406885893259752  8.49081648255649
## 490   71.687159921068 -6.18926587678444 0.753506435554632 -8.21395224345857
## 1317 222.718984861955 -3.03292643435534  0.37190983516799 -8.15500464779421
## ...               ...               ...               ...               ...
## 1792 5.98829479907705 -4.93244157414254  1.77616534224116 -2.77701712607386
## 98   6699.69196582045  1.19245883058489 0.429894831042278   2.7738384936934
## 295  6699.58197403621  1.19259109051591 0.429948214187822  2.77380170718637
## 662  5.15806783175131 -5.09298667756621  1.83903942690847 -2.76937329512713
## 899  18.1185589826749 -3.43948745771197  1.24400277511213 -2.76485513257955
##                    pvalue                 padj
##                 <numeric>            <numeric>
## 2500 2.20685744092375e-21 3.97896396598552e-18
## 1212 1.05049159823894e-18 9.47018175812401e-16
## 1917 2.05190369014182e-17 1.23319411777524e-14
## 490  2.14024730609951e-16 9.64716473224354e-14
## 1317  3.4916644951535e-16 1.25909421695235e-13
## ...                   ...                  ...
## 1792  0.00548602882496584    0.046221074632773
## 98    0.00553991737467495   0.0462481504691228
## 295    0.0055405438166004   0.0462481504691228
## 662   0.00561642447454876   0.0466654992055826
## 899   0.00569480805304419   0.0468846526010898

? 至此，便篩選出了217個在MM組和WW組之間的顯著差異表達基因仪际。至于后續(xù)的可視化分析則是因課題而異了围小，等以后有空了再補坑吧！

? 有同學可能注意到树碱，雖然我們的樣本有多個組肯适，但在差異分析時進行的還是pairwise的分析，為什么我們不可以三個組一起分析呢成榜？學過ANOVA的同學應該都知道框舔，ANOVA就是可以應對這種多組差異分析的情況。但要注意的是赎婚，ANOVA只可以告訴我們對于某個基因在這三組中是否存在差異刘绣，想要找出是哪一組有其他組別有差異還是需要進行pairwise t-test之類的分析。所以挣输，在這里我們兩組兩組地進行分析正是出于這個考慮额港，并且有更方便我們解釋差異分析的結果，說明在A組基因的表達量相對于B組的是上調(diào)還是下調(diào)歧焦。另外移斩，本文的差異分析還是處于單因子水平（只有一個變量），至于多因子的差異分析以后研究透了再和大家進行分享绢馍。

一點想法：

? 最近想要整理這些轉錄組常規(guī)分析流程向瓷，一來是因為之前學習的時候很零碎，想借此機會系統(tǒng)性地整理這方面的知識舰涌。另一方面也是因為我感覺平時做濕實驗的時候都講究個protocol猖任，一些基本的實驗，像PCR這種我們只要根據(jù)具體實驗調(diào)整一些參數(shù)就都可以跟著protocol來做了瓷耙。所以朱躺，一些經(jīng)典的生信分析也可以有一個這樣的“protocol”，要用的時候調(diào)整參數(shù)就可以了搁痛。

參考文章：

https://web.stanford.edu/class/bios221/labs/rnaseq/lab_4_rnaseq.html

https://www.researchgate.net/post/How_to_interpret_results_of_DESeq2_with_more_than_two_experimental_groups

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq

https://bioconductor.org/packages/release/bioc/manuals/DESeq2/man/DESeq2.pdf

完长搀。