Skip to content

3.3 使用 DEseq2 提取 RNA seq 分析数据中的 DEG

DEseq2 是什么？

DESeq2 是一个 R 软件包，用于对 RNA 测序（RNA-seq）数据进行差异表达分析。RNA-seq 是一种量化转录组表达水平的技术，可用于研究不同条件和时间点下基因表达的变化。

image.png

DESeq2 的主要特点

  • 正规化：考虑到样本间读取总数的差异，对基因的表达量进行正规化。
  • 差异表达分析：统计评估两个或多个条件间基因表达的变化。
  • 方差估计：通过共享基因间的方差，提高表达量变化的估计。
  • 多重检验校正：为了同时评估大量基因，进行校正以减少假阳性风险。

DESeq2 差异表达分析

DESeq2 使用正规化因子（尺寸因子）考虑文库深度差异，对原始计数进行建模。接着，估计每个基因的方差，并对估计值进行缩小，以生成更准确的方差估计值，从而对计数进行建模。最后，DESeq2 拟合负二项分布模型，并使用 Wald 检验或似然比检验进行假设检验。

运行 DESeq2

本次 DESeq2 使用以下文件。请下载并放在桌面上：

  • counts_data.csv
  • sample_info.csv

需要加载以下库：

library(DESeq2)
library(tidyverse)

读取数据：

setwd("~/Desktop")

counts_data <- read.csv("counts_data.csv")
colData <- read.csv("sample_info.csv")

查看 counts_data 和 sample_info 的内容：

counts_data 是一个矩阵文件，行是基因名，列是样本名，每个单元格记录对应的读取计数。

> head(counts_data)
#                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
# ENSG00000000003        679        448        873        408       1138       1047
# ENSG00000000005          0          0          0          0          0          0
# ENSG00000000419        467        515        621        365        587        799
# ENSG00000000457        260        211        263        164        245        331
# ENSG00000000460         60         55         40         35         78         63
# ENSG00000000938          0          0          2          0          1          0
#                 SRR1039520 SRR1039521
# ENSG00000000003        770        572
# ENSG00000000005          0          0
# ENSG00000000419        417        508
# ENSG00000000457        233        229
# ENSG00000000460         76         60
# ENSG00000000938          0          0

colData 提供了每个样本的实验变量或条件等元数据。

> head(colData)
#            cellLine dexamethasone
# SRR1039508   N61311     untreated
# SRR1039509   N61311       treated
# SRR1039512  N052611     untreated
# SRR1039513  N052611       treated
# SRR1039516  N080611     untreated
# SRR1039517  N080611       treated

colData 示例输出：

  • cellLine: 表示样本的细胞株，用于区分不同细胞株的样本。
  • dexamethasone: 表示样本是否经过地塞米松处理。

基于这些数据，DESeq2 会按照 design = ~ dexamethasone 指定评估地塞米松处理对基因表达的影响。

使用 DESeqDataSetFromMatrix 函数创建 DESeq2 分析的数据集对象（DESeqDataSet）：

dds <- DESeqDataSetFromMatrix(
  countData = counts_data,
  colData = colData,
  design = ~dexamethasone
)

design 是解析的设计公式，指定要基于哪些变量构建统计模型。在此例中，design 指定为 dexamethasone（地塞米松处理的有无）以研究其对基因表达的影响。

在执行 DESeq2 前进行预过滤步骤。总读取数不少于 10 的基因会被保留：

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

执行 DESeq2 并保存结果：

dds <- DESeq(dds)

res <- results(dds)

write.csv(res, "deseq2_results.csv")

查看 DESeq2 结果

查看 res：

> res
# log2 fold change (MLE): dexamethasone untreated vs treated 
# Wald test p-value: dexamethasone untreated vs treated 
# DataFrame with 22369 rows and 6 columns
#                  baseMean log2FoldChange     lfcSE      stat    pvalue      padj
#                 <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
# ENSG00000000003  708.6022      0.3788474  0.173165  2.187781 0.0286856  0.138465
# ENSG00000000419  520.2979     -0.2037622  0.100756 -2.022329 0.0431424  0.183141
# ENSG00000000457  237.1630     -0.0340364  0.126493 -0.269078 0.7878699  0.929861
# ENSG00000000460   57.9326      0.1171825  0.301588  0.388551 0.6976080  0.894183
# ENSG00000000971 5817.3529     -0.4409590  0.258783 -1.703975 0.0883858  0.297174
# ...                   ...            ...       ...       ...       ...       ...
# ENSG00000273483   2.68957      -0.600436  1.084493 -0.553656 0.5798145        NA
# ENSG00000273485   1.28645      -0.194034  1.346631 -0.144089 0.8854304        NA
# ENSG00000273486  15.45254       0.113349  0.426061  0.266040 0.7902086  0.930772
# ENSG00000273487   8.16323      -1.017777  0.575826 -1.767508 0.0771432  0.271676
# ENSG00000273488   8.58448      -0.218082  0.570756 -0.382094 0.7023919  0.896539

列的含义：

  • baseMean: 基因的平均计数值。
  • log2FoldChange: 对数 2 倍变化值。
  • lfcSE: 对数 2 倍变化的标准误差。
  • stat: Wald 统计量。
  • pvalue: Wald 检验的 p 值。
  • padj: 校正后的 p 值。

查看 res 的摘要：

> summary(res)
# out of 22369 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up)       : 1502, 6.7%
# LFC < 0 (down)     : 1884, 8.4%
# outliers [1]       : 51, 0.23%
# low counts [2]     : 3903, 17%
# (mean count < 4)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results

summary(res) 的输出解读：

  • 22369 个非零总读取计数的基因。
  • 调整后的 p 值< 0.1 的基因信息。
  • LFC > 0（上调）: 1502 个基因（占 6.7%）。
  • LFC < 0（下调）: 1884 个基因（占 8.4%）。
  • outliers（离群值）: 51 个基因（占 0.23%）。
  • low counts（低计数）: 3903 个基因（占 17%）。

最后，使用 MAplot 可视化变化基因：

plotMA(res)

成功后会显示如下图表，蓝色点表示校正后 p 值小于指定阈值（默认为 0.1）的统计显著差异表达基因。

image.png