author：小杜的生信筆記
DEseq2官方网址：http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html
DESeq2是最常用的差异分析的方法，在上一个教程我们分享了limma差异分析 | 简书，差异分析--limma包 - 知乎 等平台，我们后续也会分享常用的差异分析方法，并进行一个比较。
DESeq2包是为高�维计数数据的归一化、可视化和差分分析而设计的。 它利用经验贝叶斯技术估计对数折叠变化和离散的先验，并计算这些量的后验估计。
image.png
1 导入所需包
#if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
#BiocManager::install("DESeq2")
library(DESeq2)
library(dplyr)
2 导入Counts数据矩阵
## 导入Counts数据矩阵
countdata <- read.table("Countdata.txt", header = T,row.names = 1)
countdata[1:10,1:6]
> countdata[1:10,1:6]
Control_01 Control_02 Control_03 Treat_01 Treat_02 Treat_03
gene01 11 11 11 10 11 10
gene02 14 14 14 13 14 13
gene03 14 14 14 13 14 13
gene04 7 7 7 6 7 6
gene05 11 11 11 11 11 10
gene06 14 14 14 13 14 12
gene07 13 13 13 12 13 12
gene08 11 11 10 10 11 9
gene09 11 11 11 10 11 9
gene10 10 12 12 12 12 10
3 数据过滤
## 过滤在所有重复样本中小于1的基因
countdata = countdata[rowMeans(countdata) > 1,]
4 数据样本信息
样本注释信息自己在本地后制作后进行导入
coldata <- read.csv("countdata.csv", header = T, row.names = 1)
head(coldata)
#############
> head(coldata)
condition
Control_01 CK
Control_02 CK
Control_03 CK
Treat_01 Treat
Treat_02 Treat
Treat_03 Treat
检查数据Counts文件与coldata数据是否匹配
all(rownames(coldata) %in% colnames(countdata))
all(rownames(coldata) == colnames(countdata))
当返回TRUE时，表明两个数匹配。
> all(rownames(coldata) %in% colnames(countdata)) #
[1] TRUE
> all(rownames(coldata) == colnames(countdata))
[1] TRUE
5 差异分析
## 制作差异矩阵
dds <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ condition)
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,"all_diff.csv") # 导出所有的差异文件
DESeq输出的差异文件
> head(diff)
log2 fold change (MLE): condition Treat vs CK
Wald test p-value: condition Treat vs CK
DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
gene01 10.6594 -0.01071206 0.788745 -0.0135811 0.989164 0.999747
gene02 13.6626 0.00959682 0.713528 0.0134498 0.989269 0.999747
gene03 13.6626 0.00959682 0.713528 0.0134498 0.989269 0.999747
gene05 10.8224 0.03339267 0.783215 0.0426354 0.965992 0.999747
gene06 13.4731 -0.03001613 0.716954 -0.0418662 0.966605 0.999747
gene07 12.6616 0.00389414 0.735478 0.0052947 0.995775 0.999747
筛选差异基因，使用P值和LogFC进行筛选，这里我们就不多做解释，详情可以看我们前面的推文。
# 设置筛选标准
foldChange = 1
padj = 0.05
#
diffsig <- diff[(diff$pvalue < padj & abs(diff$log2FoldChange) > foldChange),]
dim(diffsig)
write.csv(diffsig, "All_diffsig.csv")
注：上调和下调的差异基因的筛选，我们在这里不在讲述，请查看我们前面的推文即可
6 差异基因热图
6.1 标准化Count值
我们后面的热图使用Count值进行分析，以及后续的分析也可以使用count值进行分析。我们的Count值是基因的原始表达量，因此需要进行标准化出来，我们这里使用vst函数进行标准化处理。
## 标准化(标准化Counts值)
vsd <- vst(dds, blind = FALSE)
normalizeExp <- assay(vsd)
绘制差异基因热图
### 差异表达热图
df <- read.csv("All_diffsig.csv", header = T)
head(df)
df02 <- as.character(df$X)
## 标准化(标准化Counts值)
vsd <- vst(dds, blind = FALSE)
normalizeExp <- assay(vsd)
##
## 差异基因的Count
diff_expr <- normalizeExp[df02,]
head(diff_expr)
## 差异热图
library(pheatmap)
annotation_col <- data.frame(Group = factor(c(rep("CK",3), rep("Treat",3))))
rownames(annotation_col) <- colnames(diff_expr)
pdf("差异基因热图.pdf", height = 8, width = 12)
p <- pheatmap(diff_expr,
annotation_col = annotation_col,
color = colorRampPalette(c("green","black","red"))(50),
cluster_cols = F,
show_rownames = F,
show_colnames = F,
scale = "row", ## none, row, column
fontsize = 12,
fontsize_row = 12,
fontsize_col = 6,
border = FALSE)
print(p)
dev.off()
image.png
获取本章节数据和代码：关注微信公众号：小杜的生信筆記（ID：Du_Bioinformatics），回复关键词：DESeq2差异分析
“小杜的生信筆記” 公众号、知乎、简书、B站平台，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！
01.差异分析|使用limma包 - 简书
02. 差异分析--limma包 - 知乎
03.R语言差异分析 | limma包 (公众号)
