edgeR差异表达分析

1、安装包

library(edgeR)
setwd("E:/4转录组数据/")

#1. 导入原始counts数据和表型数据
dat <- read.table("counts_matrix_HCC31_symble.txt",header=T,stringsAsFactors=FALSE)
class(dat)
dat[1:3,1:3]
group <- read.csv("counts_matrix_condition_HCC31.csv",header=T,stringsAsFactors=FALSE)
head(group)
y<-DGEList(counts=dat[,2:3],group=group$condition,genes=dat[,1])
y #查看y
dim(y)

#2. 过滤表达量偏低的基因
#基因在至少1个样本中得count per million（cpm）要大于1
keep <- rowSums(cpm(y)>1) >= 1
y <- y[keep,]
dim(y)
y

#3. 进行标准化因子计算 默认使用TMM方法
y <- calcNormFactors(y)
y
#这里主要是通过图形的方式来展示实验组与对照组样本是否能明显的分开
#以及同一组内样本是否能聚的比较近 即重复样本是否具有一致性
#但是需要至少三个样本才可。
#plotMDS(y)

##对于没有重复的情况，可以自定义。
##方法一：
#4-2. 自定义离散值
#如果是人类数据, 且实验做的很好(无过多的其他因素影响), 设置为0.4
#如果是遗传上相似的模式物种（这里为小鼠）, 设置为0.1
y_bcv <- y
#因为本次的数据使用的是人的数据，所以使用0.4
#如果基因很少或者很多，还可以选择修改bcv。
##这个值越小，差异基因越多，所以根据需求啦。
bcv <- 0.05
et <- exactTest(y_bcv, dispersion = bcv ^ 2)
gene1 <- decideTestsDGE(et, p.value = 0.05, lfc = 0)
head(gene1)
summary(gene1)


detags <- rownames(y)[as.logical(gene1)]
plotSmear(et,gene1.tags=detags)
abline(h=c(-1,1),col="blue")

##输出结果
# 改一下gene1的名称
colnames(gene1) <- "Signifi"
# 组合将所需要的数据组成一个新的data.frame
results <- cbind(y$genes,y$counts,et$table,gene1)
head(results)
# 将新生成的results数据框写成一个excel数据表
write.csv(x = results,file = "HCC31_symble_edgeRsig.csv",row.names = F)

##对于没有重复的情况，可以自定义。
##方法2
##假设你已经知道了一些基因是不会发生变化，例如管家基因，那么我们就可以通过它们来预测dispersion.
#4-3. 根据已知一些不会发生改变的基因推测dispersion
nosig <- names(res@.Data[res@.Data ==0,])
gene_name <- sample(nosig,500,replace = FALSE)

y1_gene <- rownames(y1)
housekeeping <- which(y1_gene %in% gene_name)

y1 <- y

y1$samples$group <- 1

y0 <- estimateDisp(y1[housekeeping,],trend = "none",tagwise =FALSE)

y$common.dispersion <- y0$common.dispersion

design <- model.matrix(~group)

fit <- glmFit(y,design)

lrt <- glmLRT(fit)

genes <- decideTestDGE(lrt,p.value=0.05,lfc = 0)


##这个运行的时候会报错，具体什么原因还不知道。
##时间紧张，就先这样吧。
#4-1. 估计离散度，这是有重复的情况下。
y <- estimateCommonDisp(y,verbose=TRUE)
y <- estimateTagwiseDisp(y)
plotBCV(y)

#5. 通过检验来鉴别差异表达基因
et <- exactTest(y)
top <- topTags(et)
top

#6. 定义差异表达基因与基本统计
summary(de<-decideTestsDGE(et));# 默认选取FDR = 0.05为阈值

#7. 图形展示检验结果
detags <- rownames(y)[as.logical(de)]
plotSmear(et,de.tags=detags)
abline(h=c(-1,1),col="blue")