volcano plot

針對(duì)測(cè)序數(shù)據(jù)和芯片數(shù)據(jù)，目前常用差異分析的R包有edgeR、limma、DESeq2，做一簡單比較，方便平時(shí)分析。內(nèi)容多為搬運(yùn)，主要方便下次尋找。

1. 三種分析方法的比較

三種分析比較

image.png

三種packages的比較

1.limma包做差異分析要求數(shù)據(jù)滿足正態(tài)分布或近似正態(tài)分布，如基因芯片、TPM格式的高通量測(cè)序數(shù)據(jù)。
2.通常認(rèn)為Count數(shù)據(jù)不符合正態(tài)分布而服從泊松分布。對(duì)于count數(shù)據(jù)來說，用limma包做差異分析，誤差較大
3.DESeq2、和 EdgeR都是基于count，然后兩個(gè)都是NB（negative binomial）但是在估計(jì)dispersion parameter的方法上面不一樣。
4.limma，edgeR，DESeq2三大包基本是做轉(zhuǎn)錄組差異分析的金標(biāo)準(zhǔn)，大多數(shù)轉(zhuǎn)錄組的文章都是用這三個(gè)R包進(jìn)行差異分析。
5.edgeR差異分析速度快，得到的基因數(shù)目比較多，假陽性高（實(shí)際不差異，結(jié)果差異）。DESeq2差異分析速度慢，得到的基因數(shù)目比較少，假陰性高（實(shí)際差異，結(jié)果不差異）。
6.需要注意的是制作分組信息的因子向量是，因子水平的前后順序，在R的很多模型中，默認(rèn)將因子向量的第一個(gè)水平看作對(duì)照組。

2.實(shí)戰(zhàn)

2.1數(shù)據(jù)準(zhǔn)備

rm(list = ls())
library("DESeq2")
library("limma")
library("edgeR")
expr = read.csv("mRNA_exprSet.csv",sep = ',',header=T)  
head(expr)

讀取的基因矩陣文件，行為基因名，列為樣本名

2.2 表達(dá)數(shù)據(jù)整理

# 對(duì)重復(fù)基因名取平均表達(dá)量，然后將基因名作為行名
expr = avereps(expr[,-1],ID = expr$X) # 自定義

# 去除低表達(dá)的基因
expr = expr[rowMeans(expr)>1,] # 自定義

# 表達(dá)矩陣分組(癌癥組織和癌旁組織）
library(stringr)
tumor <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) < 10]
normal <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) >= 10]

tumor_sample <- expr[,tumor]
normal_sample <- expr[,normal]

exprSet_by_group <- cbind(tumor_sample,normal_sample)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))

save(exprSet_by_group, group_list, file = 'exprSet_by_group_list.Rdata')

2.3 edgeR包進(jìn)行差異分析

# 表達(dá)矩陣
data = exprSet_by_group

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

# 差異表達(dá)矩陣
DGElist <- DGEList( counts = data, group = group_list)
## Counts per Million or Reads per Kilobase per Million
keep_gene <- rowSums( cpm(DGElist) > 1 ) >= 2 ## 自定義
table(keep_gene)
DGElist <- DGElist[ keep_gene, , keep.lib.sizes = FALSE ]

DGElist <- calcNormFactors( DGElist )
DGElist <- estimateGLMCommonDisp(DGElist, design)
DGElist <- estimateGLMTrendedDisp(DGElist, design)
DGElist <- estimateGLMTagwiseDisp(DGElist, design)

fit <- glmFit(DGElist, design)
results <- glmLRT(fit, contrast = c(-1, 1)) 
nrDEG_edgeR <- topTags(results, n = nrow(DGElist))
nrDEG_edgeR <- as.data.frame(nrDEG_edgeR)
head(nrDEG_edgeR)

# 提取基因差異顯著的差異矩陣
padj = 0.01 # 自定義
foldChange= 2 # 自定義
nrDEG_edgeR_signif  = nrDEG_edgeR[(nrDEG_edgeR$FDR < padj & 
                                     (nrDEG_edgeR$logFC>foldChange | nrDEG_edgeR$logFC<(-foldChange))),]
nrDEG_edgeR_signif = nrDEG_edgeR_signif[order(nrDEG_edgeR_signif$logFC),]
save(nrDEG_edgeR_signif,file = 'nrDEG_edgeR_signif.Rdata')

2.4 DESeq2包做差異表達(dá)

data = exprSet_by_group

# 分組矩陣
condition = factor(group_list)
coldata <- data.frame(row.names = colnames(data), condition)
dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = coldata,
                              design = ~condition)
dds$condition<- relevel(dds$condition, ref = "normal") # 指定哪一組作為對(duì)照組

# 差異表達(dá)矩陣
dds <- DESeq(dds)  
allDEG2 <- as.data.frame(results(dds))

# 提取基因差異顯著的差異矩陣
padj = 0.01 # 自定義
foldChange= 2 # 自定義
nrDEG_DESeq2_signif = allDEG2[(allDEG2$padj < padj & 
                                 (allDEG2$log2FoldChange>foldChange | allDEG2$log2FoldChange<(-foldChange))),]
nrDEG_DESeq2_signif = nrDEG_DESeq2_signif[order(nrDEG_DESeq2_signif$log2FoldChange),]
save(nrDEG_DESeq2_signif, file = 'nrDEG_DESeq2_signif.Rdata')

2.5 limma包分析過程

# 表達(dá)矩陣
data = exprSet_by_group

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

# 差異表達(dá)矩陣
DGElist <- DGEList( counts = data, group = group_list )
keep_gene <- rowSums( cpm(DGElist) > 1 ) >= 2 # 自定義
table(keep_gene)
DGElist <- DGElist[ keep_gene, , keep.lib.sizes = FALSE ]

DGElist <- calcNormFactors( DGElist )
v <- voom(DGElist, design, plot = TRUE, normalize = "quantile")
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(contrasts = c('tumor-normal'), levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

nrDEG_limma_voom = topTable(fit2, coef = 'tumor-normal', n = Inf)
nrDEG_limma_voom = na.omit(nrDEG_limma_voom)
head(nrDEG_limma_voom)

# 提取基因差異顯著的差異矩陣
padj = 0.01 # 自定義
foldChange= 2 # 自定義
nrDEG_limma_voom_signif = nrDEG_limma_voom[(nrDEG_limma_voom$adj.P.Val < padj & 
                                              (nrDEG_limma_voom$logFC>foldChange | nrDEG_limma_voom$logFC<(-foldChange))),]
nrDEG_limma_voom_signif = nrDEG_limma_voom_signif[order(nrDEG_limma_voom_signif$logFC),]
save(nrDEG_limma_voom_signif, file = 'nrDEG_limma_voom_signif.RDATA')

最后，參考文檔也很重要?？梢苑嗊@三個(gè)包的說明文檔

參考：
http://www.itdecent.cn/p/cf2ec58e5361
https://blog.csdn.net/weixin_43700050/article/details/98085127