【生信学习第一天】DEseq2 差异表达基因计算

一、介绍

• 分析来自 RNA-seq 的计数数据的一项基本任务是检测差异表达的基因。计数数据以表格的形式呈现，其中报告了每个样本已分配给每个基因的序列片段的数量。其他检测类型也有类似的数据，包括比较 ChIP-Seq、HiC、shRNA 筛选和质谱分析。一个重要的分析问题是与条件内的变异性相比，条件之间的系统变化的量化和统计推断。

• DESeq2是DEseq的升级版，但是DEseq2只适用于有生物学重复的试验，而DEseq既可以做有生物学重复也可以做无重复（或部分重复的）试验。

• DESeq2 包提供了使用负二项式广义线性模型测试差异表达的方法；离散度和对数倍数变化的估计包含数据驱动的先验分布。此小插图解释了包的使用并演示了典型的工作流程。

• 文中代码主要参考了 小明的数据分析笔记本 在B站发布的视频，对其进行了添加注释和修改。自学过程中观看了小明老师的很多教程和视频，推荐大家关注。

• 获取文中代码和示例文件及结果 回复 ：DESeq2

二、安装 DEseq2和dplyr

#安装DEseq2包需要借助BioManager包来安装，DEseq2安装包比较大，安装时间比较慢。
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")
BiocManager::install("dplyr")

三、计算差异表达

• 大家起初使用R包计算差异表达基因的时候，都会不太注意软件的运行原理，这里推荐这篇博客科学网—[转载]DESeq2原理不简单 - 胡耿的博文 (sciencenet.cn)，原理讲述的比较通俗易懂，因此这里我直接放上可以运行的代码和示例数据，通过回复关键词获取代码，供大家学习和交流。

##计算和过滤 reads count
# 获取工作目录路径
getwd()
#设置工作目录路径
setwd("C:/Users/11753/Desktop/")
# 读入原始的reads count文件，注意：一定是未经过标准化的row reads count
mycounts<-read.csv("merge_js153_guy11_readsconut_v2.csv",row.names = 1)
# 大概看一下文件内容和格式
head(mycounts)
dim(mycounts)
# 过滤掉只含有数目很少reads count的gene。
mycounts_1<-mycounts[rowSums(mycounts) != 0,]
dim(mycounts_1)
# 过滤数据存储到.csv文件中
mymeta<-read.csv("mymeta.csv",stringsAsFactors = T)
mymeta
colnames(mycounts_1) == mymeta$id


## 安装和计算差异表达基因
#载入R包
library(DESeq2)
# 获取reads count 矩阵
dds <- DESeqDataSetFromMatrix(countData=mycounts_1,
                              colData=mymeta,
                              design=~dex)
# 计算差异表达
dds <- DESeq(dds)
res <- results(dds)
# 查看结果
head(res)
class(res)
# 将结果文件转换成dataframe格式，用于后续的R包计算
res_1<-data.frame(res)
class(res_1)
head(res_1)

# 载入dptyr包，添加差异表达 up, down , not change
library(dplyr)
res_1 %>%
  mutate(group = case_when(
    log2FoldChange >= 2 & padj <= 0.05 ~ "UP",
    log2FoldChange <= -2 & padj <= 0.05 ~ "DOWN",
    TRUE ~ "NOT_CHANGE"
  )) -> res_2

table(res_2$group)
# 将结果写入.csv文件中
write.csv(res_2,file="diff_expr_result.csv",
          quote = F)
## 获取reads count标准化后的矩阵，对于差异表达基因筛选有参考作用
# 获取规范化计数并将其写入文件
nc = counts(dds,normalized=TRUE)

# 将其转换为数据帧以具有适当的列名。
dt = data.frame("id"=rownames(nc),nc)

# Save the normalize data matrix.
write.table(dt, file="norm-matrix-js153_guy11_deseq22.txt", sep="\t", row.name=FALSE, col.names=TRUE,quote=FALSE)

# 画一个样品的相关性 PCA plot
colData(dds)
rld <- rlog(dds)
plotPCA(rld,intgroup=c("id","sizeFactor"))

上一篇推文，介绍了chrome同步助手用于访问谷歌搜索和谷歌学术，但是不是很稳定，今天找到了新的一款Edge插件，IGG谷歌助手，仅需要一个qq邮箱，验证后就可使用，获取方式见前面。不过我还是推荐使用镜像。

获取新的Edge谷歌助手插件 回复：IGG

参考文献

• DESeq2 — 上海交大超算平台用户手册 文档 (sjtu.edu.cn)

• R语言DESeq2包做转录组RNAseq差异表达分析的一个简单小例子_哔哩哔哩_bilibili

• 科学网—[转载]DESeq2原理不简单 - 胡耿的博文 (sciencenet.cn)

• Gene-level differential expression analysis with DESeq2 | Introduction to DGE - ARCHIVED (hbctraining.github.io)

---

本人在读研和读博期间，收集了大量的生信书籍和科研写作书籍，30多本精选的书籍，约有500 Mb文件，可以按照需要是否获取。

**关注公众号 晓亮自学生信笔记，回复：生信资料，即可免费获取。**资料均为科研的交流学习使用，禁止商用，若有侵权请联系我删除。