R语言DESeq2基因差异表达分析

经过表达定量后，我们已经得到了基因的表达量矩阵，差异表达分析通常是RNA-seq分析的第一步。

差异基因表达分析通常都是在R中，常用的有DESeq2，edgeR，limma等几种，这次主要介绍用DESeq2来进行差异表达分析。

需要准备的数据：基因表达定量矩阵（counts）及分组文件

安装

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

使用

#导入数据
count <- read.csv("RB_count.csv",header=T,row.names = 1)
myDesign <- data.frame(sampleID=colnames(count),condition=c("BL","BL","BL","WL","WL","WL"))#可自行准备导入或构建数据框
colnames(count)==myDesign$sampleID

#数据过滤
countData <- count[apply(count, 1, sum) > 0 , ]

#DEseq2运行
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
  countData = round(countData),
  colData = myDesign,
  design = ~ condition)

dds$condition <- relevel(dds$condition,ref="WL")  #设置对照组
dds <- DESeq(dds)
res <- results(dds)
res <- cbind(sample1="BL",sample2="WL",as.data.frame(res))

res$padj[is.na(res$padj)] <- 1 #将padj为NA的全部赋值为1
res005 <- res[res$padj<0.05,]#提取padj<0.05

#也可注释上调或下调
res <- data.frame(res)
library(dplyr)
res %>%
  mutate(group=case_when(log2FoldChange>=1&padj<=0.05~"UP",
                         log2FoldChange<= -1&padj<=0.05~"Down",
                         TRUE~"Not_Change")) -> res2
table(res2$group)

#写出文件
write.csv(res,"BL_WL_DEGs_all.csv")
write.csv(res005,"BL_WL_DEGs_padj005.csv")