DESeq2差异表达分析

分享是一种态度

在前文scRNA-seq marker identification(二)，我们我们提到了差异分析，下面我们来详细了解下

学习目标

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

本课程基于2019 Bioconductor tutorial on scRNA-seq pseudobulk DE analysis (http://biocworkshops2019.bioconductor.org.s3-website-us-east-1.amazonaws.com/page/muscWorkshop__vignette/)。特别是，许多数据整理步骤均来自这个教程。

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 <- [e-mail]

# Set up metadata as desired for aggregation and DE analysis
metadata$cluster_id <- factor([e-mail])

# 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)

未完......

---

注：以上内容来自哈佛大学生物信息中心(HBC)_的教学团队的生物信息学培训课程。原文链接：https://hbctraining.github.io/scRNA-seq/schedule/  点击 “阅读原文” 可直达

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注