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