差異分析|DESeq2完成配對樣本的差異分析

本文為群中小夥伴進行的一次差異分析探索的記錄。

前段時間拿到一個RNA-seq測序數據（病人的癌和癌旁樣本，共5對）及公司做的差異分析結果（1200+差異基因），公司告知用的是配對樣本的DESeq分析。

考慮到平時limma和DESeq2包進行差異分析時沒有特別註明是否配對，這配對和非配對有啥區別呢？

於是分別嘗試使用limma和DESeq2包的非配對分析，發現得到的差異基因和公司的差距很大。我查了好多關於RNA-seq配對分析的資料，發現幾乎沒有這方面的帖子。

詢問公司DESeq配對分析的代碼，公司說保密不能給，此外公司還告知現在的配對樣本的分析都改用了DESeq2。好吧，那就只能自己動手，以下為探索過程的一個記錄。

# 加載包
library(openxlsx)
library(DESeq2)
library(limma)
library(edgeR)

1.讀入原始數據及分組信息
rowdata <- read.xlsx("Count_data.xlsx",sheet = 1,rowNames = T)
gset <- rowdata[rowMeans(rowdata)>0,] # 剔除表達量低的基因

group_list <- c(rep("case",5), rep("control",5)) 
group_list <- factor(group_list,levels = c("control","case"))


## 2.1表達矩陣
data = gset

## 2.2分組矩陣
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)

## 2.3差比表達矩陣構建，並過濾
DGElist <- DGEList( counts = data, group = group_list)
keep <- rowSums(cpm(DGElist) > 0.5 ) >=2 
table(keep) 
# FALSE 2443 TRUE 15909

DGElist_QC <- DGElist[keep, ,keep.lib.sizes=FALSE] 
dim(DGElist) 
# 18352 10 過濾前

dim(DGElist_QC)
# 15909 10過濾後

## 2.4歸一化基因表達分布

DGElist_norm <- calcNormFactors(DGElist_QC, method = "TMM")
DGElist_norm$samples$norm.factors # 查看每個樣品歸一化後的係數
DGElist_QC$samples$norm.factors #未歸一化之前的樣本係數，都是1

## 2.5 limma包進行voom函數
v <- voom(DGElist_norm, design, plot = TRUE, normalize = "quantile")
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(contrasts = c('case-control'), levels = design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## 2.6提取差異矩陣
nrDEG_limma_voom = topTable(fit2, coef = "case-control", n = Inf)
nrDEG_limma_voom = na.omit(nrDEG_limma_voom)
nrDEG_limma_voom = nrDEG_limma_voom[order(nrDEG_limma_voom$logFC),] 

## 2.7 定義差異基因
nrDEG <- nrDEG_limma_voom
nrDEG$Group = "notsignificant"

logFC_cutoff <- 0.6 # 定義差異基因的標準，自定義

nrDEG$Group [which( (nrDEG$P.Value < 0.05) & (nrDEG$logFC > logFC_cutoff) )] = "upregulated"
nrDEG$Group [which( (nrDEG$P.Value < 0.05) & (nrDEG$logFC < -logFC_cutoff) )] = "downregulated"

table(nrDEG$Group)
可以看到只有67個下調的33個上調的，火山圖不好看，而且根本沒法繼續做GO和KEGG分析。
OK，嘗試使用DESeq2包的非配對差異分析。


## 3.1表達矩陣
data = apply(gset, 2, as.integer) ## DESeq2分析需要是整數
row.names(data) <- row.names(gset)

## 3.2分組矩陣
condition = group_list
coldata <- data.frame(row.names = colnames(data), condition)

dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = coldata,
                              design = ~condition)

dds$condition<- relevel(dds$condition, ref = "control") # 指定哪一組作為對照組

## 3.3差異表達矩陣

dds <- DESeq(dds)  
nrDEG_DESeq2 <- as.data.frame(results(dds))
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2$log2FoldChange),] 

## 3.4定義差異基因
nrDEG <- nrDEG_DESeq2
nrDEG$Group = "notsignificant"
logFC_cutoff <- 0.6

nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange < -logFC_cutoff) )] = "downregulated"

table(nrDEG$Group)
 
可以看到常規的DESeq2分析比limma voom分析多了一些差異基因，但是和公司給的1200+的差異基因還是差遠了。

發現差異之後開始了檢索和求助之旅，查了很多帖子，也求助了一些大神，似乎很少人注意過DESeq2包做配對的差異分析。
討論群中小夥伴貼了limma進行配對分析的方式，於是我去查閱了DESeq2的說明書，可以看到說明書就這麼一句話：
剩下的事情就簡單了，依此修改後，DESeq2包成功做出了配對差異分析，復現了公司的結果。好了，下面就是使用DESeq2包完成配對差異分析的代碼了，自取！
## 4.1表達矩陣

data = apply(gset, 2, as.integer) ## DESeq2分析需要是整數
row.names(data) <- row.names(gset)

## 4.2分組矩陣，配對分析與常規分析最大的區別就在分組矩陣

condition = group_list
# 配對分析要加上這段代碼，知道誰和誰是一對，比如1,1是一對，5,5是一對
subject <- factor(c(1,2,3,4,5,1,2,3,4,5))  

coldata <- data.frame(row.names = colnames(data), condition)

# 注意在design中加上配對信息
dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = coldata,
                              design = ~subject +condition) 

dds$condition<- relevel(dds$condition, ref = "control") 

## 4.3差異表達矩陣，還是和常規分析一樣

dds <- DESeq(dds)
nrDEG_DESeq2 <- as.data.frame(results(dds))
rld <- rlog(dds)
# 這裡我還提取了標準化後的表達矩陣，可以用於後續的熱圖繪製等等
normal_gset <- assay(rld) 
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2$log2FoldChange),] 

## 4.4定義差異基因

nrDEG <- nrDEG_DESeq2
nrDEG$Group = "notsignificant"

logFC_cutoff <- 0.6
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange < -logFC_cutoff) )] = "downregulated"

table(nrDEG$Group)
此時終於成功得到了1200+的差異基因，通過對比公司的分析結果和我做的結果，幾乎完成了復現。有一點區別在於我做了一些低表達基因的過濾，導致log2Foldchang和pvalue略微有些變化，但這區別可以忽略不計。
總結來說，由於算法的不同，不同差異分析的R包得到的差異基因數量不完全一致。重要的是，針對配對的樣本，如果不進行配對分析而用常規的差異分析，這樣的結果可能會大不相同。因此，在分析數據的時候，一定要明白實驗設計。

最後，我還發現有意思的一個情況。在進行clusterProfilerR包的GSEA分析的時候，我用非配對分析得到的log2FoldChang出的GSEA結果，和配對分析得到的log2FoldChang出的GSEA結果幾乎是一致的，儘管兩種分析方法得到的log2FoldChang有很大差別。
PS：公眾號後臺回復「入群」，備註  姓名+單位/學校+研究方向  即邀請您進交流群！
◆ ◆ ◆  ◆ ◆
精心整理（含圖版）|R語言生信分析，可視化，你要的全拿走，建議收藏！