如何进行差异表达基因分析-基于DESeq2包

崔耳又又

差异表达分析是做RNA-seq分析必经的一步，这里给大家介绍一下。

1 DEseq2安装

在命令行中输入以下内容：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")

2 准备文件

在做差异表达基因分析时需要两个输入文件，一个是表达矩阵文件，一个是分组信息，将两个文件整理为“.txt”格式。

表达矩阵文件（第一列为基因ID），行为样本信息，列为count值（一般不用FPKM值）：

gene_id	NIP_1	NIP_2	NIP_3	MUT_1	MUT_2	MUT_3
Os01g0100100	643	589	591	644	778	650
Os01g0100200	9	7	7	13	5	13
Os01g0100400	187	171	130	189	182	184
Os01g0100500	1131	1084	1033	1226	1165	1160
Os01g0100600	280	332	299	310	415	406
Os01g0100700	2135	1883	1531	1913	1593	2304
Os01g0100800	400	389	288	305	325	405
Os01g0100900	1384	1724	1356	1629	1678	1480
Os01g0101200	81	69	71	64	86	83
Os01g0101300	57	32	46	29	48	62
Os01g0101600	1445	1190	1221	1437	1814	1345

分组文件，第一列为样本名，第二列为分组名：

sample	group
NIP_1	NIP
NIP_2	NIP
NIP_3	NIP
MUT_1	MUT
MUT_2	MUT
MUT_3	MUT

3 使用R语言进行分析，话不多说，直接上脚本

#---导入R包---
library(DESeq2)
library(dplyr)

#---导入文件---
countdata <- read.table("C:/Users/bf/Desktop/新建文件夹/shootOEvsNIP.txt", header = T,row.names = 1) #表达矩阵文件,注意更换为自己的文件
countdata = countdata[rowMeans(countdata) > 1,]
coldata  <- read.table("C:/Users/bf/Desktop/新建文件夹/sampleOE_NIP.txt", header = T, row.names = 1) #样本分组文件,注意更换为自己的文件
all(rownames(coldata) %in% colnames(countdata))  
all(rownames(coldata) == colnames(countdata))

#---制作差异矩阵---
dds <-  DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ group) 
dim(dds)

#---过滤---
dds <- dds[rowSums(counts(dds)) > 1,]  
nrow(dds)

#---差异比较---
dep <- DESeq(dds)
res <- results(dep)
diff = res
diff <- na.omit(diff)  ## 去除NA
dim(diff)

#---导出文件---
write.csv(diff,"C:/Users/bf/Desktop/shootOEvsNIP.xls") 

4 输出文件

输出文件格式如下：

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
Os01g0100100	1164.78	0.43	0.11	3.90	0.00	0.00
Os01g0100200	26.25	1.68	0.51	3.31	0.00	0.00
Os01g0100400	256.76	-0.09	0.24	-0.39	0.69	0.76
Os01g0100500	1463.07	0.02	0.12	0.17	0.87	0.90
Os01g0100600	380.22	-0.34	0.16	-2.17	0.03	0.05
Os01g0100700	2196.86	0.96	0.11	8.81	0.00	0.00
Os01g0100800	473.50	-0.11	0.17	-0.64	0.52	0.60
Os01g0100900	1786.58	-0.36	0.11	-3.31	0.00	0.00
Os01g0101200	686.24	-10.36	0.67	-15.54	0.00	0.00
Os01g0101300	68.55	0.02	0.32	0.07	0.94	0.96
Os01g0101600	1886.67	-4.21	0.13	-33.60	0.00	0.00
Os01g0101800	89.00	-0.27	0.32	-0.85	0.40	0.48
Os01g0102000	219.78	0.62	0.25	2.45	0.01	0.03

· gene_id: 基因或转录本的标识符。

· baseMean: 平均表达量。

· log2FoldChange: 对数变换后的 fold change，表示在不同条件之间的表达倍数变化。

· lfcSE: log2 fold change 的标准误差。

· tat: 统计检验值。

· pvalue: 未经校正的 p 值。

· padj: 经过多重假设检验校正后的调整 p 值（通常使用 FDR 校正），用于控制假阳性发现率。

微信扫一扫分享文章