DESeq2 差异表达分析

Question

学习目标

• 了解如何准备用于 pseudobulk 差异表达分析的单细胞 RNA-seq 原始计数数据
• 利用 DESeq2 工具对特定细胞类型聚类进行 pseudobulk 差异表达分析
• 创建函数以遍历不同细胞类型的 pseudobulk 差异表达分析

本课程基于 2019 Bioconductor tutorial on scRNA-seq pseudobulk DE analysis 。特别是，许多数据整理步骤均来自这个教程。

DESeq2 差异表达分析

在鉴定了 scRNA-seq 簇的细胞类型之后，我们通常希望在特定细胞类型内的条件之间执行差异表达分析。虽然 Seurat 中存在执行此分析的函数，但这些分析的 p 值通常会被夸大，因为每个细胞都被视为样本。我们知道，样本中的单个细胞并不是彼此独立的，因为它们是从相同的动物/样本中分离出来的，来自相同的环境。

如果我们把细胞当作样本，那么我们真正研究的不是群体间的变异，而是个体之间的变异。因此，我们只能在个人的层面上作出结论，而不能在群体的层面上作出结论。通常，我们想要研究的是哪些基因对群体水平（而不是个体水平）的某条件下很重要，所以我们需要从不同的生物/样本（而不是从不同的细胞）中获取样本。为此，当前的最佳做法是使用 pseudobulk 方法，该方法涉及以下步骤：

1. 将子集替换为感兴趣的细胞类型以执行 DE 分析。
2. 提取 QC 过滤后的原始计数用于 DE 分析
3. 将计数和元数据聚合到样本级别
4. 进行 DE 分析（每个条件至少需要两个生物重复才能执行分析，但建议进行更多重复)。

我们将使用与其余工作流相同的数据集，现在已将其多路分解为单个样本，以便使用复制来进行差异表达分析。我们将把它作为 SingleCellExperient 对象导入。

注意：要从我们在单细胞分析工作流程结束时创建的 Seurat 对象中提取子集并提取细胞，我们可以使用类似于以下代码：

# Bring in Seurat object
seurat <- readRDS("path/to/seurat.rds")

# Extract raw counts and metadata to create SingleCellExperiment object
counts <- seurat@assays$RNA@counts 

metadata <- seurat@meta.data

# Set up metadata as desired for aggregation and DE analysis
metadata$cluster_id <- factor(seurat@active.ident)

# Create single cell experiment object
sce <- SingleCellExperiment(assays = list(counts = counts), 
                           colData = metadata)

# Identify groups for aggregation of counts
groups <- colData(sce)[, c("cluster_id", "sample_id")]

探索数据集

在这次分析中，我们将使用 Kang et al, 2017 的单细胞 RNA-seq 数据集，我们曾在单细胞 RNA-seq 分析工作流程的其余部分使用过。然而，对于差异表达分析，我们使用的是具有 8 个对照样本和 8 个干扰素刺激样本的非混合计数数据。这与 scRNA-seq 分析的其余部分形成对比，该分析使用了取自 8 名狼疮患者的Pooled Peripheral Blood Mononuclear Cells (PBMCs)，分为 single pooled control and a single pooled interferon-stimulated

注意：如果可能，您应该始终从 scRNA-seq 工作流开始使用非混合样本。

我们从 ExperimentHub R 包中获取了拆分成 8 个单独样本的原始计数数据集，如下所述：http://biocworkshops2019.bioconductor.org.s3-website-us-east-1.amazonaws.com/page/muscWorkshop__vignette/​。

Metadata

除了原始数据之外，我们还需要收集有关数据的信息；这称为 metadata。经常有一如果我们拿到数据就直接开始上手分析，但对这些数据的来源样本一无所知，那就没有太大的意义了。

下面提供了我们的数据集的一些相关 metadata：

• 这些文库是用 10X Genomics 第 2 版化学方法制备的
• 样品在 Illumina NextSeq 500 上测序
• 将 8 例狼疮患者的 PBMC 样本各分成两等份，然后进行 demultiplexed
  • 用 100U/mL recombinant IFN-β刺激一份 PBMC 6h。
  • 另一份不做处理。
  • 6 小时后，将每个条件下的 8 个样本混合在两个最终池（刺激细胞和对照细胞）中。
  • 对照和刺激混合样本分别鉴定了 12,138 和 12,167 个细胞（去除二倍体后)。
  • 使用工具 Demuxlet 对样本进行 demultiplexed
• 经过聚类和标记鉴定，鉴定出以下细胞类型
  • B cells
  • CD4 T cells
  • CD8 T cells
  • NK cells
  • FCGR3A+ Monocytes
  • CD14+ Monocytes
  • Dendritic cells
  • Megakaryocytes

注意：在单细胞工作流程中，我们还确定了一些其他细胞类型，但是我们将继续使用此数据集和在分析中确定的细胞类型。

设置 R 环境

差异表达分析的做准备，我们需要设置项目和目录结构，加载必要的库，并引入原始计数的单细胞 RNA-seq 基因表达数据。打开 RStudio 并创建一个名为“DE_Analysis_scrnaseq”的新 R 项目。然后，创建以下目录：

DE_analysis_scrnaseq/
├── data
├── results
└── figures

下载数据

将 RData 对象下载到数据文件夹中：scRNA-seq raw counts（ https://www.dropbox.com/s/dgdooocb1a03dvk/scRNA-seq_input_data_for_DE.rds?dl=1）

新建脚本

打开一个新的 Rscript 文件，并以一些注释开始，以指示该文件将包含的内容：

# May 2020

# Single-cell RNA-seq analysis - Pseudobulk DE analysis with DESeq2

保存脚本为： DE_analysis_scrnaseq.R

加载库

引入特定细胞类型的原始计数数据后，我们将使用来自各种程序包的工具将数据整理为所需的格式，然后将单细胞的原始计数聚合到样本级别。然后，我们将使用 DESeq2 对感兴趣的条件进行差异表达分析。要了解有关 DESeq2 方法和分析步骤解构的更多信息，可参考（ https://hbctraining.github.io/DGE_workshop_salmon/schedule/）。

加载我们将用于分析的库

# Load libraries
library(scater)
library(Seurat)
library(tidyverse)
library(cowplot)
library(Matrix.utils)
library(edgeR)
library(dplyr)
library(magrittr)
library(Matrix)
library(purrr)
library(reshape2)
library(S4Vectors)
library(tibble)
library(SingleCellExperiment)
library(pheatmap)
library(apeglm)
library(png)
library(DESeq2)
library(RColorBrewer)

加载 RData(RDS) 对象

我们正在使用的数据集已作为 RData 对象保存到 RDS 文件。我们可以使用 readRDS() 函数读取它。

# Read in the dataset
sce <- readRDS("data/scRNA-seq_input_data_for_DE.rds")

RData 对象是一个单细胞实验对象，它是使用 SingleCellExperient 包生成的专用列表类型。这些对象具有以下结构：

我们可以使用 SingleCellExperient 包中的函数来提取不同的组件。首先我们可以查看一下实验数据的计数和元数据。

# Explore the raw counts for the dataset

## Check the assays present
assays(sce)

## Explore the raw counts for the dataset
dim(counts(sce))

counts(sce)[1:6, 1:6]

sc_DE_countdata.png

我们看到原始计数数据是一个一个接一个的基因稀疏矩阵，有超过 35,000 行（基因）和近 30,000 列（细胞)。

注意：不要对这个数据集运行 head() ，因为它仍然会显示数千列，所以我们只查看了前六行和前六列。

接下来，我们可以了解一下每个细胞的元数据。

## Explore the cellular metadata for the dataset
dim(colData(sce))

head(colData(sce))

对于每个细胞，我们都有关于相关条件(ctrl 或 stim)、样本 ID 和细胞类型的信息。我们将使用此信息来执行感兴趣的任何特定细胞类型的条件之间的差异表达分析。

获取样本中细胞间聚合的必要指标

首先，我们需要确定数据集中存在的群集数量和群集名称。

# Named vector of cluster names
kids <- purrr::set_names(levels(sce$cluster_id))
kids

# Total number of clusters
nk <- length(kids)
nk

# Named vector of sample names
sids <- purrr::set_names(levels(sce$sample_id))

# Total number of samples 
ns <- length(sids)
ns

要执行样本水平差异表达分析，我们需要生成样本水平的元数据。为此，我们将以匹配样本 ID 的因子级别的顺序，对单个细胞元数据中的样本进行重新排序，然后只从与该样本对应的第一个细胞中提取样本信息。

# Generate sample level metadata

## Determine the number of cells per sample
table(sce$sample_id)

## Turn named vector into a numeric vector of number of cells per sample
n_cells <- as.numeric(table(sce$sample_id))

## Determine how to reoder the samples (rows) of the metadata to match the order of sample names in sids vector
m <- match(sids, sce$sample_id)

## Create the sample level metadata by combining the reordered metadata with the number of cells corresponding to each sample.
ei <- data.frame(colData(sce)[m, ], 
                  n_cells, row.names = NULL) %>% 
                select(-"cluster_id")
ei

在将细胞聚合到样本级之前，如果尚未执行此步骤，我们希望确保移除质量较差的细胞。一般来说，我们建议对质量控制指标进行更严格，实际的探索，并对过滤阈值进行更细致的选择，如此前所述 scRNA-seq—质量控制 ；然而，为了更快地进行差异表达分析，我们将只使用 BioConductor 教程中执行的 Scater 程序包中的函数来删除计数异常值和低计数基因。

# Perform QC if not already performed
dim(sce)

# Calculate quality control (QC) metrics
sce <- calculateQCMetrics(sce)

# Get cells w/ few/many detected genes
sce$is_outlier <- isOutlier(
        metric = sce$total_features_by_counts,
        nmads = 2, type = "both", log = TRUE)

# Remove outlier cells
sce <- sce[, !sce$is_outlier]
dim(sce)

## Remove lowly expressed genes which have less than 10 cells with any counts
sce <- sce[rowSums(counts(sce) > 1) >= 10, ]

dim(sce)

现在，我们准备将计数聚合到样本级别。从本质上讲，我们取的是每种细胞类型中每个样本的计数总和。

Count aggregation to sample level

# Aggregate the counts per sample_id and cluster_id

# Subset metadata to only include the cluster and sample IDs to aggregate across
groups <- colData(sce)[, c("cluster_id", "sample_id")]

# Aggregate across cluster-sample groups
pb <- aggregate.Matrix(t(counts(sce)), 
                       groupings = groups, fun = "sum") 

class(pb)

dim(pb)

pb[1:6, 1:6]

这个聚合的输出是一个稀疏矩阵，当我们快速查看时，我们可以看到它是一个基于细胞类型的基因-样本矩阵。例如，在 B 细胞中，样本 ctrl101 的 NOC2L 基因有 12 个相关计数。

要在每个细胞类型的基础上执行 DE 分析，我们需要通过几种方式来处理我们的数据。我们需要做以下几个步骤：

1. 按细胞类型拆分数据
2. 变换矩阵，使基因成为行名，样本成为列名

我们将按细胞类型划分数据；但是，并非所有样本都包含每种细胞类型的细胞。要确定每种细胞类型存在哪些样本，我们可以运行以下命令：

# Not every cluster is present in all samples; create a vector that represents how to split samples
splitf <- sapply(stringr::str_split(rownames(pb), 
                                    pattern = "_",  
                                    n = 2), `[`, 1)

现在，我们可以将矩阵转换为一个列表，该列表被分成每个群集的计数矩阵，然后对每个数据框进行转换，这样行就是基因，列就是样本。

# Turn into a list and split the list into components for each cluster and transform, so rows are genes and columns are samples and make rownames as the sample IDs
pb <- split.data.frame(pb, 
                       factor(splitf)) %>%
        lapply(function(u) 
                set_colnames(t(u), 
                             stringr::str_extract(rownames(u), "(?<=_)[:alnum:]+")))

class(pb)

# Explore the different components of list
str(pb)

还可以检查每个群集的每个样本的计数：

# Print out the table of cells in each cluster-sample group
options(width = 100)
table(sce$cluster_id, sce$sample_id)

用 DESeq2 进行基因的差异表达分析

我们将使用 DESeq2 进行 DE 分析，下面的流程图中用绿色显示了使用 DESeq2 的分析步骤。DESeq2 首先将计数数据归一化，以消除样本之间文库大小和 RNA 组成的差异。然后，我们将使用归一化计数在基因和样本水平上为 QC 绘制一些曲线图。最后一步是使用 DESeq2 包中的适当函数来执行差异表达式分析。在接下来的课程中，我们将深入讨论这些步骤中的每一个步骤，但有关 DESeq2 的更多细节和有用建议可以在我们的材料中找到，这些材料详细介绍了 bulk RNA-seq 数据和 DESeq2 vignette 的工作流程。 https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html ​

Sample-level metadata

要执行 DE 分析，除了任何其他样本水平的元数据（例如批次、性别、年龄等）之外，我们还需要所有样本的元数据，包括群集 ID、样本 ID 和感兴趣的条件(Group_Id)。EI 数据框保存样本 ID 和条件信息，但是我们需要将该信息与群集 ID 结合起来。

首先，我们将为每个细胞类型群集创建一个样本名称组合向量。

# Get sample names for each of the cell type clusters

# prep. data.frame for plotting
get_sample_ids <- function(x){
        pb[[x]] %>%
                colnames()
}

de_samples <- map(1:length(kids), get_sample_ids) %>%
        unlist()

然后，我们可以获得与向量中的每个样本相对应的群集 ID。

# Get cluster IDs for each of the samples

samples_list <- map(1:length(kids), get_sample_ids)

get_cluster_ids <- function(x){
        rep(names(pb)[x], 
            each = length(samples_list[[x]]))
}

de_cluster_ids <- map(1:length(kids), get_cluster_ids) %>%
        unlist()

最后，让我们使用集群 ID 和相应的样本 ID 创建一个数据框。我们将把状态信息合并在一起。

# Create a data frame with the sample IDs, cluster IDs and condition

gg_df <- data.frame(cluster_id = de_cluster_ids,
                    sample_id = de_samples)

gg_df <- left_join(gg_df, ei[, c("sample_id", "group_id")]) 

metadata <- gg_df %>%
        dplyr::select(cluster_id, sample_id, group_id) 
        
metadata

对感兴趣的群集取子集

现在我们有了样本级别的元数据，我们可以使用 DESeq2 运行差异表达式分析。通常，我们希望对多个不同的群集执行分析，这样我们就可以将工作流设置为在任何群集上轻松运行。

为此，我们可以创建数据集中所有群集细胞类型 ID 的群集向量。然后，我们可以选择要对其执行 DE 分析的细胞类型。

查看一下集群细胞类型 ID：

# Generate vector of cluster IDs
clusters <- levels(metadata$cluster_id)
clusters

[1] "B cells"           "CD14+ Monocytes"   "CD4 T cells"       "CD8 T cells"      
[5] "Dendritic cells"   "FCGR3A+ Monocytes" "Megakaryocytes"    "NK cells"    

我们在数据集中看到了多种不同的免疫细胞类型。让我们对 B 细胞执行 DE 分析，它是我们向量中的第一个元素。从向量中提取 B 细胞：

clusters[1]

我们可以使用此输出对 B 细胞运行 DE 分析。首先，我们可以仅将元数据和计数设置为 B 细胞。

# Subset the metadata to only the B cells
cluster_metadata <- metadata[which(metadata$cluster_id == clusters[1]), ]
head(cluster_metadata)

# Assign the rownames of the metadata to be the sample IDs
rownames(cluster_metadata) <- cluster_metadata$sample_id
head(cluster_metadata)

# Subset the counts to only the B cells
counts <- pb[[clusters[1]]]

cluster_counts <- data.frame(counts[, which(colnames(counts) %in% rownames(cluster_metadata))])

# Check that all of the row names of the metadata are the same and in the same order as the column names of the counts in order to use as input to DESeq2
all(rownames(cluster_metadata) == colnames(cluster_counts))  

创建 DESeq2 对象

现在，我们可以创建 DESeq2 对象以准备运行 DE 分析。我们需要包括计数，元数据和设计公式以进行我们感兴趣的比较。在设计公式中，我们还应在元数据中包含我们想要回归其变化的任何其他列（例如批次，性别，年龄等）。我们只需要比较感兴趣的内容，它作为 group_id 存储在我们的元数据数据框中。

有关 DESeq2 工作流程和设计公式的更多信息，请参见 DESeq2 的材料（ https://hbctraining.github.io/DGE_workshop_salmon/schedule/）。

dds <- DESeqDataSetFromMatrix(cluster_counts, 
                              colData = cluster_metadata, 
                              design = ~ group_id)