果子老師做過(guò)一個(gè)非常驚人的舉動(dòng)，用DESeq2處理1225例樣本的TCGA數(shù)據(jù)，在沒(méi)有使用DESeq多線程參數(shù)parallel的情況下，跑了將近40個(gè)小時(shí)。
那么問(wèn)題來(lái)了，在那么大的樣本量的情況下，應(yīng)該用DESeq2進(jìn)行數(shù)據(jù)處理嗎？我的結(jié)論是不應(yīng)該，DESeq2的適用場(chǎng)景是小樣本的差異表分析，降低假陽(yáng)性。當(dāng)你的樣本量足夠多的時(shí)候，我們其實(shí)有更好的選擇。
這里以果子老師的數(shù)據(jù)為例，來(lái)對(duì)比DESeq2的結(jié)果和我的分析結(jié)果進(jìn)行比較.
加載DESeq2結(jié)果
load(file="dds_very_long.Rdata")
library(DESeq2)
deseq2_result <- results(dds)
table(deseq2_result$padj < 0.01)
# FALSE TRUE
# 23997 25072
下面我分析時(shí)的數(shù)據(jù)預(yù)處理 部分，
options(stringsAsFactors = FALSE)
# 加載數(shù)據(jù)
load(file = "BRCA_RNASEQ_exprdf.Rdata")
# 提取表達(dá)量矩陣
expr_mt <- as.matrix(expr_df[,-1])
row.names(expr_mt) <- expr_df$gene_id
colnames(expr_mt) <- colnames(expr_df)[-1]
# 根據(jù)文庫(kù)大小標(biāo)準(zhǔn)化
expr_mt <- expr_mt / rep(colSums(expr_mt), each=nrow(expr_mt)) * 1e6
# 過(guò)濾地表達(dá)基因
expr_mt <- expr_mt[rowSums(expr_mt > 0) > (ncol(expr_mt) / 3), ]
# 統(tǒng)計(jì)癌癥和癌旁
TCGA_id <- colnames(expr_mt)
table(substring(TCGA_id,14,15))
### 我們發(fā)現(xiàn)了7個(gè)轉(zhuǎn)移的樣本，本次分析，我們關(guān)注的是癌癥和癌旁，先把轉(zhuǎn)移的樣本去掉
### 原發(fā)和轉(zhuǎn)移的對(duì)比作為家庭作業(yè)
TCGA_id <- TCGA_id[substring(TCGA_id,14,15)!="06"]
### 創(chuàng)建metadata
sample <- ifelse(substring(TCGA_id,14,15)=="01","cancer","normal")
sample <- factor(sample,levels = c("normal","cancer"),ordered = F)
metadata <- data.frame(TCGA_id,sample)
下一步，利用非參數(shù)檢驗(yàn)方法， wilcox.test，關(guān)于非參數(shù)檢驗(yàn)的緣起可以看「女士品茶」的第16章擺脫參數(shù)
威爾科克森注釋著計(jì)算t檢驗(yàn)和方法分析的公式，意識(shí)到這些不同尋常的極端數(shù)值會(huì)對(duì)結(jié)果產(chǎn)生極大的影響，導(dǎo)致“學(xué)生”的t檢驗(yàn)偏小。 ... 如果異常值體現(xiàn)了某種因素對(duì)系統(tǒng)數(shù)據(jù)的系統(tǒng)性污染，那么使用非參數(shù)方法只會(huì)讓事情變得更糟。
# wilcox.test差異分析 ---------------------------------------------------------
cancer_sample <- metadata[metadata$sample == "cancer", "TCGA_id"]
normal_sample <- metadata[metadata$sample == "normal", "TCGA_id"]
cancer_mt <- expr_mt[,colnames(expr_mt) %in% cancer_sample ]
normal_mt <- expr_mt[,colnames(expr_mt) %in% normal_sample ]
# 計(jì)算logFoldChanges
logFC <- log2(rowMeans(as.matrix(cancer_df)) / rowMeans(as.matrix(normal_df)))
library(future.apply)
plan(multiprocess)
p_values <- future_lapply(seq(nrow(cancer_df)), function(x){
res <- wilcox.test(x = cancer_mt[x,], y = normal_mt[x,] )
res$p.value
})
p <- unlist(p_values)
p.adj <- p.adjust(p, method = "fdr")
table(p.adj < 0.01)
# FALSE TRUE
# 10997 24030
我們得到了24,030個(gè)校正后p值小于0.01的基因，而DESeq2是25,072個(gè)。如果比較全部的基因的話，韋恩圖上可以發(fā)現(xiàn)，絕大部分基因都是相同的。
但是通常情況下，我們會(huì)更去關(guān)注一些變化比較大且p值顯著的基因，用這些基因去做下游的富集分析。所以，下一步就是看看后面富集分析結(jié)果兩者有什么區(qū)別。
我們用Y叔的clusterProfiler，去分析倍數(shù)變化4倍，矯正p值小于0.01的基因
# 提取基因
library(clusterProfiler)
library(org.Hs.eg.db)
org <- org.Hs.eg.db
diffgene1 <- row.names(expr_mt)[p.adj < 0.01 & abs(logFC) > 2]
diffgene1 <- substr(diffgene1, 1, 15)
diffgene2 <- row.names(deseq2_result)[deseq2_result$padj < 0.01 &
! is.na(deseq2_result$padj) &
abs(deseq2_result$log2FoldChange) > 2]
diffgene2 <- substr(diffgene2, 1, 15)
GO富集分析
library(clusterProfiler)
library(org.Hs.eg.db)
org <- org.Hs.eg.db
diffgene1 <- row.names(expr_mt)[p.adj < 0.01 & abs(logFC) > 2]
diffgene1 <- substr(diffgene1, 1, 15)
diffgene2 <- row.names(deseq2_result)[deseq2_result$padj < 0.01 &
! is.na(deseq2_result$padj) &
abs(deseq2_result$log2FoldChange) > 2]
diffgene2 <- substr(diffgene2, 1, 15)
ego1 <- enrichGO(diffgene1,
OrgDb = org,
keyType = "ENSEMBL",
ont = "BP"
)
ego2 <- enrichGO(diffgene2,
OrgDb = org,
keyType = "ENSEMBL",
ont = "BP"
)
merge_result <- merge_result(list(wilcox=ego1,DESeq2=ego2))
dotplot(merge_result,showCategory= 20 )
從點(diǎn)圖中，你可以認(rèn)為這兩個(gè)分析結(jié)果是一致。
綜上，當(dāng)你在樣本量足夠多（兩組都不少于10吧），其實(shí)沒(méi)有去用DESeq2這些復(fù)雜的工具，用基礎(chǔ)的統(tǒng)計(jì)學(xué)檢驗(yàn)方法就能得到很好的結(jié)果了。
在樣本量比較小的時(shí)候，用復(fù)雜的模型是無(wú)奈之舉，它有很多假設(shè)成分在，尤其是你還想從無(wú)重復(fù)的實(shí)驗(yàn)設(shè)計(jì)中算p值。當(dāng)你樣本量夠多的時(shí)候，用最簡(jiǎn)單的模型其實(shí)就會(huì)有很好的結(jié)果。
