Xueyd's Blog

0%

edgeR差异分析

发表于 分类于
本文字数： 547 阅读时长 ≈ 1 分钟

总结edgeR差异分析。

1. 有重复的差异分析

1. 构建edgeR的DGEList对象，并过滤

1
2
3
4
5
6
7
8
library(edgeR)
exprSet <- read.csv("counts.csv", header=T)
group_list <- rep(c("Control", "FU"), each = 3)

# 保留非全0的基因
exprSet <- exprSet[rowSums(exprSet) > 0,]
# 构建edgeR的DGEList对象
dge <- DGEList(counts = exprSet,group = factor(group_list))

保留在至少在两个样本中cpm值大约1的基因

1
2
3
4
keep <- rowSums(cpm(dge)>1) >= 2
table(keep)
dge <- dge[keep, , keep.lib.sizes=FALSE]
dge$samples$lib.size <- colSums(dge$counts)

归一化基因表达分布

1
dge <- calcNormFactors(dge)

假设数据符合正态分布，构建线性模型，0代表x线性模型的截距为0。

1
2
3
design <- model.matrix(~0+factor(group_list))
rownames(design) <- colnames(dge)
colnames(design) <- levels(factor(group_list))

2. 差异表达分析

计算线性模型的参数，并拟合线性模型

1
2
3
4
dge <- estimateGLMCommonDisp(dge,design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)

进行差异分析

1
lrt <- glmLRT(fit, contrast=c(-1,1))   # 1,-1意味着后比前

3. 提取过滤差异分析结果

1
2
DEG <- topTags(lrt, n=nrow(dge))
DEG <- as.data.frame(DEG)

2. 无重复edgeR差异分析

1. 数据过滤

1
2
3
4
5
6
7
8
9
library(edgeR)

count <- read.csv(file = "counts.csv",row.names = 1,header = T)
########## treat vs NC ############
colnames(count) <- c("Control","treat")

# 过滤全为0的基因
count <- count[rowSums(count[,1:2])>0,]
count <- as.matrix(count)

2. 构建DGEList对象

1
2
3
4
5
6
7
8
9
group_list <- c("Control","treat")
dge <- DGEList(counts = count,
               group = factor(group_list))

# 过滤掉低表达基因
# 保留在至少在一个样本里有表达的基因(CPM > 1)
keep <- rowSums(cpm(dge)>1) >= 1
table(keep)
dge <- dge[keep, , keep.lib.sizes=FALSE]

3. 标准化

考虑到测序深度不同, 我们需要对其进行标准化, 避免文库大小不同导致的分析误差。
edgeR里默认采用TMM(trimmed mean of M-values) 对配对样本进行标准化。

1
dge <- calcNormFactors(dge)

4. 差异表达分析

根据经验给定一个值(BCV, square-root-dispersion). edgeR给的建议是, 如果你是人类数据, 且实验做的很好(无过多的其他因素影响), 设置为0.4, 如果是遗传上相似的模式物种, 设置为0.1, 如果是技术重复, 那么设置为0.01。
bcv值大，得到的差异基因数量就少，因此可以根据需求设置。

1
2
bcv <- 0.2
et <- exactTest(dge, dispersion=bcv^2)

用decideTestsDGE看下有多少基因上调, 多少基因下调. 设置p.value=0.05。

1
2
3
4
gene1 <-decideTestsDGE(et,p.value =0.05,lfc =0)
summary(gene1)
DEG <- topTags(et, n=nrow(dge))
DEG <- as.data.frame(DEG)

参考：

  1. https://blog.csdn.net/u012110870/article/details/102804557