基因芯片数据分析（八）：DESeq2差异分析实战案例

DoubleHelix

发布于 2019-12-13 10:38:21

3.5K0

发布于 2019-12-13 10:38:21

文章被收录于专栏：生物信息云生物信息云

我们在前2篇文章分别介绍了edgeR和DESe

基因芯片数据分析（六）：DESeq2包的基本原理我们接下来通过一个案例介绍利用edgeR进行差异分析。

基因芯片数据分析（七）：edgeR差异分析实战案例本文接着介绍DESeq2包进行差异分析。

包的安装和加载

# 包的安装和加载
BiocManager::install("DESeq2")
library("DESeq2")

读入数据

这里我们用的数据是一个原始的counts数据的Excel文件，和上一讲中用的数据一样（想运行案例，文末获取文件）

# 读入原始的counts数据
counts <- read.table("gene_counts.xls", sep = "\t", header = T, row.names = 1)

导入的数据我们通过head()函数查看前6行。

行名A1,A2,A3,B1,B2,B3为样本名，列名是基因名。

创建分组

设置实验组别，在基因芯片数据分析（六）：DESeq2包的基本原理这篇文章中我们介绍基本原理时，有一步需要选择参考样本，在实际分析中，我们可以自己选择参考样本，一般都是对照组作为参考样本，在DESeq2中，对照组是是，control样本，实验组是case。原始数据中有6个样本，对照和实验组各三个。所以我们可以创建一个分组向量。

# 创建分组
colData <- data.frame(row.names = c("A1","A2","A3","B1","B2","B3"),
                      condition =
                        factor(c("control","control","control","case","case","case"),
                               levels = c("control","case")))

构建DESeqDataSet对象

构建DESeqDataSet对象用于分析，colData指定我们的样本信息和分组列表，design = ~ condition表示分组信息安装colData的condition列分组。

#构建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(countData = counts, colData = colData,
                              design = ~ condition)

差异分析

# 函数分析差异
dds <- DESeq(dds)
# 计算标准化因子
sizeFactors(dds)
#提取差异表达结果
res <- results(dds)

class(res)
res <- as.data.frame(res)

结果处理与保存

# 添加一列
res <- cbind(rownames(res), res)
# 重命名列名
colnames(res) <- c("gene_id", "baseMean", "log2FoldChange", "lfcSE", "stat",
                   "pval", "padj")
#保存文件到本地
write.table(res, "case-vs-control-all-DESeq2.gene.xls",
            sep = "\t", col.names = TRUE, row.names = FALSE, quote = FALSE,
            na = "")

差异基因筛选

这里和前文基因芯片数据分析（七）：edgeR差异分析实战案例差不多，不多解释！

# 获取差异基因
resSig <- res[which(res$pval < 0.05 & abs(res$log2FoldChange) > 1),]
resSig[which(resSig$log2FoldChange > 0), "up_down"] <- "Up"
resSig[which(resSig$log2FoldChange < 0), "up_down"] <- "Down"
write.table(resSig, "case-vs-control-diff-pval-0.05-FC-2-DESeq2.gene.xls",
            sep = "\t", col.names = TRUE, row.names = FALSE, quote = FALSE, na = "")

这样我们就可以获取差异表达的基因了，可以用于后面的作图，比如火山图！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2019-12-06，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

本文分享自 MedBioInfoCloud 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

数据分析

登录后参与评论

0 条评论

热度