edgeR分析RNAseq数据

生信喵

背景
       分析原始count数据，除了DEseq2包还有一个R包熟为人知，就是edgeR包。今天带来的是教程和实战，使用的数据还是前一篇推文中的。可以关注公众号回复"20221009"获取。
      
教程地址：
文档：

1. http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/designmatrices.html

复制代码

英文：

1. http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

复制代码

中文：

1. http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html

复制代码

1. rm(list=ls())
   
2. library(limma)
   
3. library(Glimma)
   
4. library(edgeR)
   
5. setwd('/home/xhs/jyxy/14/')

复制代码

一、读取数据
1.1 方法一：读取 featureCounts 表达矩阵

1. countdata <- read.csv("CountMatrix.csv",header = T,row.names = 1)
   
2. head(countdata)
   
3. #过滤表达矩阵
   
4. nrow(countdata)
   
5. countdata <- countdata[rowSums(countdata)>10,]
   
6. nrow(countdata) #21281
   
7. 
   
8. #读取样品信息
   
9. coldata <- read.csv("sample.csv",header = T)
   
10. head(coldata)
    
11. 
    
12. #生成DGElist对象
    
13. dge <- DGEList(counts = countdata, genes = row.names(countdata),group = as.factor(coldata$dex))
    
14. dge
    
15. 
    
16. #添加样品信息
    
17. dge$samples$cell <- as.factor(coldata$cell)
    
18. 
    
19. #备份一份
    
20. dge2 <- dge

复制代码

1.2 方法二：读取 htseq 多个文件表达信息

1. setwd("./htseq/")
   
2. files <- dir(pattern = "tsv")
   
3. files
   
4. 
   
5. #查看文件内容
   
6. read.delim(files[1], nrow=5,header = F)
   
7. #生成DGElist对象
   
8. x <- readDGE(files = files,columns = c(1,2))
   
9. x
   
10. #查看属性
    
11. class(x)
    
12. x$samples
    
13. x$counts
    
14. dim(x)
    
15. samplenames <- substring(colnames(x),1,nchar(colnames(x)))
    
16. samplenames
    
17. 
    
18. colnames(x) <- samplenames
    
19. 
    
20. x$samples$group <- as.factor(coldata$dex)
    
21. x$samples$cell <- as.factor(coldata$cell)
    
22. x$samples

复制代码

二、DGElist 对象探索

1. barplot(colSums(dge$counts),las=2,col = rep(rainbow(4),each=2))
   
2. opar <- par(no.readonly = TRUE)
   
3. par(mar=c(7.1,4.1,4.1,2.1))
   
4. barplot(colSums(dge$counts),las=2,col = rep(rainbow(4),each=2))
   
5. par(opar)

复制代码

      
       不同样品 reads count

三、对数据计算标准化因子

1. # Normalization method: "TMM","TMMwsp","RLE","upperquartile","none"
   
2. dge.tmm <- calcNormFactors(dge,method = "TMM")
   
3. dge.tmm$samples
   
4. plotMDS(dge.tmm)
   
5. 
   
6. dge.tmmwsp <- calcNormFactors(dge,method = "TMMwsp")
   
7. dge.tmmwsp$samples
   
8. plotMDS(dge.tmmwsp)
   
9. 
   
10. dge.upperquartile <- calcNormFactors(dge,method = "upperquartile")
    
11. dge.upperquartile$samples
    
12. plotMDS(dge.upperquartile)
    
13. 
    
14. dge.rle <- calcNormFactors(dge,method = "RLE") # DESeq2, cuffdiff
    
15. dge.rle$samples
    
16. plotMDS(dge.rle)
    
17. 
    
18. dge.none <- calcNormFactors(dge,method = "none")
    
19. dge.none$samples
    
20. plotMDS(dge.none)
    
21. 
    
22. dge <- calcNormFactors(dge,method = "TMM")
    
23. dge
    
24. dge$samples
    
25. 
    
26. #检查样本中的异常值，使用plotMDS()作图。
    
27. plotMDS(dge)
    
28. plotMDS(dge,labels = dge$samples$cell,col=rep(rainbow(4),each=2))

复制代码

      
       MDS 图

四、差异表达计算

1. #1 生成实验矩阵
   
2. model.matrix( ~ group,data = dge$samples)
   
3. design.matrix <- model.matrix( ~ cell + group ,data = dge$samples)
   
4. rownames(design.matrix) <- rownames(dge$samples)
   
5. design.matrix
   
6. 
   
7. #2 估计离散系数（dispersion）
   
8. dge <- estimateDisp(dge,design = design.matrix)
   
9. dge$common.dispersion
   
10. dge$tagwise.dispersion
    
11. 
    
12. #可以用BCV plot查看离散系数；
    
13. plotBCV(dge, cex = 0.8)

复制代码

      
       BCV plot 离散系数

1. # plot var and mean
   
2. plotMeanVar(dge, show.raw=TRUE, show.tagwise=TRUE, show.binned=TRUE)

复制代码

      
       var and mean 图

1. #3 exactTest检验
   
2. dge_res <- exactTest(dge)
   
3. dge_results <- as.data.frame(topTags(dge_res,n=nrow(dge$counts),sort.by = "logFC"))
   
4. View(dge_results)
   
5. 
   
6. #3  最大似然法检测
   
7. fit <- glmFit(y = dge, design = design.matrix)
   
8. lrt <- glmLRT(fit, coef=2)
   
9. lrt_results <- as.data.frame(topTags(lrt,n=nrow(dge$counts),sort.by = "logFC"))
   
10. View(lrt_results)
    
11. 
    
12. #4 筛选差异表达基因
    
13. write.csv(x = dge_results,file = "dge_results.csv")
    
14. 
    
15. 
    
16. #筛选出差异表达基因，logFC >=1或者<=-1,并且q值小于0.05的基因
    
17. res <- na.omit(dge_results) #21281
    
18. res0.05 <- res[(abs(res$logFC)>=1 & res$FDR <=0.05), ]
    
19. nrow(res0.05)#1531
    
20. 
    
21. #筛选出差异表达基因，logFC >=2,并且q值小于0.01的基因
    
22. res0.01 <- res[(abs(res$logFC)>=1 & res$FDR <=0.01), ]
    
23. nrow(res0.01)#1286
    
24. 
    
25. write.csv(res0.05,file = "res0.05.csv")
    
26. #找出差异表达最显著的基因，按log2差异倍数排序
    
27. head(res0.05[order(res0.05$logFC,decreasing = TRUE), ])
    
28. head(res0.05[order(res0.05$FDR,decreasing = TRUE), ])

复制代码

五、结果可视化
       X-axis (A) 表示两个样品的⼏何平均数；
       Y-axis (M) 表示两个样品的 fold-change；

      

1. #绘制MA plot
   
2. select.sign.gene = decideTestsDGE(dge_res, p.value = 0.01)
   
3. select.sign.gene_id = rownames(dge_res)[as.logical(select.sign.gene)]
   
4. plotSmear(dge_res, de.tags = select.sign.gene_id, cex = 0.5,ylim=c(-4,4))
   
5. abline(h = c(-2, 2), col = "blue")

复制代码

      
       MA 图

1. #绘制差异基因散点图（类似火山图）；
   
2. plotMD(dge_res,ylim=c(-10,10))
   
3. abline(h=c(-1, 1), col="green")

复制代码

      
       差异基因散点图

1. #热图
   
2. cpm <- cpm(dge$counts, log=FALSE)
   
3. geneid <- rownames(res0.05[order(res0.05$logFC,decreasing = TRUE), ])
   
4. cpm <- cpm[as.factor(geneid[1:50]),]
   
5. library(pheatmap)
   
6. pheatmap(cpm)

复制代码

      
       差异基因表达热图