Differential expression analysis -- DESeq2

前言

DESeq2 是一款基于负二项分布的R包,对于 count data 进行处理,用于 RNA-seq 基因差异分析。

安装

RGUI中输入:

1
2
3
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2", version = "3.8")

矩阵准备

1. 载入表达矩阵

表达矩阵的数据结构:
Count_data

gene 列为基因 ensemble_id,后几个列名为样本名称,中间为 count(次数)

2. 设置分组信息

分组信息的结构:
Coldata

id 列为样本名称,后几列的列名为因素

核心步骤

1.载入数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
options(stringsAsFactors = FALSE)
# 载入 count 文件(由 htseq-count 计数,未标准化),sep 指定分隔符,col.names 指定列名
# 一个因素:种类 bec 和 lyec,三个生物学重复
ten <- read.table("E:/RNA-SEQ/4.0deseq/E10_count.txt", sep="\t", col.names = c("gene_id","bec1"))
ele <- read.table("E:/RNA-SEQ/4.0deseq/E11_count.txt", sep="\t", col.names = c("gene_id","bec2"))
fif <- read.table("E:/RNA-SEQ/4.0deseq/E15_count.txt", sep="\t", col.names = c("gene_id","bec3"))
six <- read.table("E:/RNA-SEQ/4.0deseq/E16_count.txt", sep="\t", col.names = c("gene_id","lyec1"))
eig <- read.table("E:/RNA-SEQ/4.0deseq/E18_count.txt", sep="\t", col.names = c("gene_id","lyec2"))
thi <- read.table("E:/RNA-SEQ/4.0deseq/E13_count.txt", sep="\t", col.names = c("gene_id","lyec3"))
# 合并各个表达矩阵,一次性合并报错了(估计是行数不统一问题)
all <- merge(ten, ele, by="gene_id")
all <- merge(all, fif, by="gene_id")
all <- merge(all, six, by="gene_id")
all <- merge(all, eig, by="gene_id")
all <- merge(all, thi, by="gene_id")
# 把 ensemble_id 作为行号
rownames(all) <- all[,1]
# 把 ensemble_id 列删除,因为已经由行号来表示
all = all[,-1]
# 删除整行 count 皆为零的行
all = all[rowSums(all) >0,]
# 删除 count 文件最后五行的结果统计
all = all[-1,]
all = all[-1,]
all = all[-1,]
all = all[-1,]
all = all[-1,]

设置分组信息

1
2
3
4
# rep 指定生物学重复数
condition <- factor(c(rep("bec",3),rep("lyec",3)), levels = c("bec","lyec"))
# 创建 colData
colData <- data.frame(row.names=colnames(all), condition)

构建 dds(DESeqDataSet)

1
2
3
4
5
# 载入包
library(DESeq2)
# 构建 dds
dds <- DESeqDataSetFromMatrix(all, colData, design= ~ condition)
dds <- DESeq(dds)

结果

1
2
3
4
5
6
7
8
9
res= results(dds)  # 或者 res = results(dds, contrast=c("condition", "bec", "lyec"))
# 根据 p-value 排序
res = res[order(res$pvalue),]
# 显示一共多少个 genes 上调和下调(FDR0.1)
summary(res)
# 保存为文件
write.csv(res,file="All_results.csv")
# padj小于0.05的 genes
table(res$padj<0.05)

提取差异基因

1
2
3
4
# padj<0.05 且 表达倍数取以2为对数后大于1或者小于-1的差异表达基因
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
# 保存差异表达基因的结果
write.csv(diff_gene_deseq2,file= "DEGs_BEC_vs_LYEC.csv")

gene symbol 注释

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 载入包,没安装的百度包名,进入 Bioconductor 安装
library('biomaRt')
library("curl")
# 重新载入差异表达基因,head=T:将首行作为列名
DEG <- read.table("E:/RNA-SEQ/4.0deseq/DEGs_BEC_vs_LYEC.csv", sep=",", head = T)
# 查看有哪些数据库可用
ensembl=useMart("ensembl")
# 查看需要用到的物种的数据集
listDatasets(ensembl)
# ggallus_gene_ensembl:鸡的数据集,数据库选择 ensemble
mart <- useDataset("ggallus_gene_ensembl", useMart("ensembl"))
# 将 id 赋值给 my_ensembl_gene_id
my_ensembl_gene_id<-DEG[,1]
# 注释结果
gg_symbols<- getBM(attributes=c('ensembl_gene_id','external_gene_name',"description"),filters = 'ensembl_gene_id', values = my_ensembl_gene_id, mart = mart)

合并差异基因和注释结果

被合并的两个表要有共同的列和列名

1
2
3
4
# 按 ensembl_gene_id 合并,如果表不大可用 Excel 插入共同的列和列名
DEG_ann<-merge(DEG,gg_symbols,by="ensembl_gene_id")
# 将合并结果保存为文件
write.csv(DEG_ann,file= "DEGs_BEC_vs_LYEC_ANN.csv")

下游分析

GO 和 KEGG 富集分析(也可以有蛋白网络互作分析)

参考

  1. RNA-seq(7): DEseq2筛选差异表达基因并注释(bioMart)
  2. 简单使用DESeq2/EdgeR做差异分析
  3. Bioconductor:DESeq2
  4. RNA_seq分析流程软件思维导图
-------------本文结束 感谢您的阅读-------------
暖一下
ZJohnson wechat
扫一扫,领红包!
0%