SCS【10】单细胞转录组之差异表达分析 (Monocle 3)

点击关注，桓峰基因

桓峰基因公众号推出单细胞系列教程，有需要生信分析的老师可以联系我们！首选看下转录分析教程整理如下：

Topic 6. 克隆进化之 Canopy

Topic 7. 克隆进化之 Cardelino

Topic 8. 克隆进化之 RobustClone

SCS【1】今天开启单细胞之旅，述说单细胞测序的前世今生

SCS【2】单细胞转录组之 cellranger

SCS【3】单细胞转录组数据 GEO 下载及读取

SCS【4】单细胞转录组数据可视化分析 (Seurat 4.0)

SCS【5】单细胞转录组数据可视化分析 (scater)

SCS【6】单细胞转录组之细胞类型自动注释 (SingleR)

SCS【7】单细胞转录组之轨迹分析 (Monocle 3) 聚类、分类和计数细胞

SCS【8】单细胞转录组之筛选标记基因 (Monocle 3)

SCS【9】单细胞转录组之构建细胞轨迹 (Monocle 3)

SCS【10】单细胞转录组之差异表达分析 (Monocle 3)

今天来说说单细胞转录组数据的细胞轨迹分析，学会这些分析结果，距离发文章就只差样本的选择了，有创新性的样本将成为文章的亮点，并不是分析内容了！

这期继续介绍 Monocle 3 软件包用于差异表达分析。

前言

单细胞转录组测序(scRNA-seq)实验使我们能够发现新的细胞类型，并帮助我们了解它们是如何在发育过程中产生的。Monocle 3包提供了一个分析单细胞基因表达实验的工具包。

Monocle 3可以执行三种主要类型的分析:

聚类、分类和计数细胞。单细胞RNA-Seq实验允许发现新的(可能是罕见的)细胞亚型。
构建单细胞轨迹。在发育、疾病和整个生命过程中，细胞从一种状态过渡到另一种状态。Monocle 3可以发现这些转变。
差异表达分析。对新细胞类型和状态的描述，首先要与其他更容易理解的细胞进行比较。Monocle 3包括一个复杂的，但易于使用的表达系统。

工作流程图如下：

差异基因表达分析是RNA-Seq实验中的一项常见任务。Monocle可以帮助你找到细胞组之间表达差异的基因，并评估这些变化的统计意义。Monocle 3包含一个强大的系统，可以找到在不同类型的细胞中不同的基因，在不同的发育时间点收集的基因，或者以不同的方式被干扰的基因。

在Monocle中有两种差异分析方法:

回归分析:使用fit_models()，可以评估每个基因是否依赖于诸如时间、治疗等变量。
图自相关分析:使用 graph_test()，可以找到沿着轨迹或在簇之间变化的基因。

Monocle还具有寻找差异表达基因的协同调节模块的专门功能。Monocle还允许交互式地查询轨迹的特定 cluster或 partition (例如分支点)，以寻找其中变化的基因。

数据读取

expression_matrix <- readRDS("packer_embryo_expression.rds")
cell_metadata <- readRDS("packer_embryo_colData.rds")
gene_annotation <- readRDS("packer_embryo_rowData.rds")cds <- new_cell_data_set(expression_matrix, cell_metadata = cell_metadata, gene_metadata = gene_annotation)

数据集处理

## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 100)## Step 2: Remove batch effects with cell alignment
cds <- align_cds(cds, alignment_group = "batch")## Step 3: Reduce the dimensions using 'UMAP', 'tSNE', 'PCA', 'LSI', 'Aligned'
cds <- reduce_dimension(cds, reduction_method = "UMAP")## Step 4: Cluster the cells
cds <- cluster_cells(cds)## Step 5: Learn a graph
cds <- learn_graph(cds)## Step 6: Order cells
cds <- order_cells(cds)

Regression analysis

如何使用 Monocle 来查找根据几种不同标准差异表达的基因。对 cell_data_set 对象中的所有基因执行差异表达分析可能需要几分钟到几个小时，具体时间取决于分析的复杂程度。为了使小图简单快速，我们将使用小组基因。不过，请放心，Monocle可以分析数千个基因，甚至在大型实验中也可以分析，这使得有助于发现正在研究的生物过程中动态调节的基因。

让我们从一小组在纤毛神经元中很重要的基因开始:

ciliated_genes <- c("che-1", "hlh-17", "nhr-6", "dmd-6", "ceh-36", "ham-1")
cds_subset <- cds[rowData(cds)$gene_short_name %in% ciliated_genes, ]

Monocle 中的差异分析工具非常灵活。Monocle的工作原理是为每个基因拟合一个回归模型。可以指定这个模型来考虑实验中的各种因素(时间、治疗等)。例如，在胚胎数据中，细胞被收集在不同的时间点。我们可以测试上述基因的表达是否会随着时间的推移而发生变化，首先对每个基因拟合一个广义线性模型:

y是有一个随机变量对应于基因I的表达值, xt是收集每个细胞的时间(以分钟为单位)，βt 捕捉时间对表达的影响，和β0 是截距项。我们可以通过将这个模型拟合到每一个基因上，来识别随时间变化的基因，然后测试是否βt明显不同于零。为此，我们首先调用fit_models()函数:

gene_fits <- fit_models(cds_subset, model_formula_str = "~embryo.time")

Gene_fits是包含每个基因一行的tibble。模型列包含广义线性模型对象，每个对象的目的是使用上面的方程解释基因在细胞中的表达。参数model_formula_str应该是一个指定模型公式的字符串。在测试中使用的模型公式可以包括colData表中作为列存在的任何术语，包括Monocle在其他分析步骤中添加的那些列。例如，如果您使用cluster_cells，您可以使用_cluster或partition(分别)作为模型公式来测试cluster和partition之间的差异基因。还可以包含多个变量，例如~embryo。时间+批次，这对于减去不需要的效果非常有帮助。

现在让我们看看这些基因中哪些有时间依赖性的表达。首先，我们使用argument_table()函数从每个模型中提取一个系数表:

fit_coefs <- coefficient_table(gene_fits)

请注意，表格中每个基因模型的每一项都包含一行。我们通常不关心截距项β0 ，所以我们可以很容易地提取时间项:

emb_time_terms <- fit_coefs %>%filter(term == "embryo.time")

现在，让我们找出那些具有显著时间成分的基因。argument_table()测试在Wald检验下，每个系数是否显著不同于零。默认情况下，coefficient_table()使用Benjamini和Hochberg方法调整这些p值以进行多重假设检验。这些调整的值可以在q_value列中找到。我们可以对结果进行过滤，并控制误发现率如下:

emb_time_terms %>%filter(q_value < 0.05) %>%select(gene_short_name, term, q_value, estimate)

我们可以看到6个基因中的5个随着时间的变化而显著变化。

Monocle还提供了一些简单的方法来绘制按差异分析中使用的因素分组的一小组基因的表达。这有助于将上面的测试显示的差异形象化。其中一种绘制方法是“小提琴”。

plot_genes_violin(cds_subset, group_cells_by = "embryo.time.bin", ncol = 2) + theme(axis.text.x = element_text(angle = 45,hjust = 1))

Controlling for batch effects and other factors

gene_fits <- fit_models(cds_subset, model_formula_str = "~embryo.time + batch")
fit_coefs <- coefficient_table(gene_fits)
fit_coefs %>%filter(term != "(Intercept)") %>%select(gene_short_name, term, q_value, estimate)

Evaluating models of gene expression

这些模型在“解释”基因表达方面有多好? 我们可以使用evaluate_fits()函数来计算每个模型的适合度:

evaluate_fits(gene_fits)

我们是否应该在我们的基因表达模型中包含批次? Monocle提供了一个compare_models()函数，可以帮助您做出决定。比较模型采用两个模型，并返回它们之间的似然比检验结果。每当您向模型添加术语时，它都会提高拟合性。但我们应该总是尽可能使用最简单的模型来解释我们的数据。似然比检验帮助我们确定拟合性的改善是否大到足以证明我们引入的额外术语的复杂性。像这样运行compare_models():

time_batch_models <- fit_models(cds_subset, model_formula_str = "~embryo.time + batch",expression_family = "negbinomial")
time_models <- fit_models(cds_subset, model_formula_str = "~embryo.time", expression_family = "negbinomial")
compare_models(time_batch_models, time_models) %>%select(gene_short_name, q_value)

两个模型中的第一个被称为完整模型。这个模型本质上是一种预测给定细胞中每个基因表达值的方法，既知道它是什么时候收集的，又知道它来自哪批细胞。第二个模型称为简化模型，它做同样的事情，但它只知道每个细胞被收集的时间。因为完整的模型有更多关于每个细胞的信息，它能更好地预测基因在每个细胞中的表达。Monocle必须为每个基因回答的问题是，完整模型的预测比简化模型的预测好多少。通过了解每个细胞的批次获得的改进越大，似然比检验的结果就越显著。

我们可以看到，所有基因的似然比检验都是显著的，这表明数据中存在大量的批处理效应。因此，我们有理由将批处理项添加到模型中。

Choosing a distribution for modeling gene expression

Monocle使用广义线性模型来捕捉基因的表达如何依赖于实验中的每个变量。这些模型要求您指定一个描述基因表达值的分布。大多数使用这种方法分析基因表达数据的研究使用负二项分布，这通常适用于测序读取或UMI计数数据。负二项式是许多RNA-seq分析包的核心，如DESeq2。

Monocle的fit_models()支持负二项分布和下表中列出的其他几个分布。默认值是“拟泊松”，它与负二项式非常相似。与负二项式相比，Quasipoisson的精度略低，但拟合速度快得多，这使得它非常适合具有数千个细胞的数据集。

expression_family有几个允许的值:

quasipoisson Quasi-poisson   ++  ++  Default for fit_models(). Recommended for most users.
negbinomial Negative binomial   +++ +   Recommended for users with small datasets (fewer than 1,000 cells).
poisson Poisson -   +++ Not recommended. For debugging and testing only.
binomial    Binomial    ++  ++  Recommended for single-cell ATAC-seq

Finding genes that change as a function of pseudotime

这种分析的核心目标是识别随着细胞沿轨迹发展而发生变化的基因。了解基因启动和关闭的顺序可以为新的发育模式提供信息。例如，Sharon和Chawla等人最近分析了拟时间依赖基因，得出了胰岛如何在胰腺中形成的全新模型。

让我们看看胚胎数据:

plot_cells(cds, color_cells_by = "cell.type", label_groups_by_cluster = FALSE, label_leaves = FALSE,label_branch_points = FALSE)

我们如何找到在不同路径上不同表达的基因?
我们如何找到那些被限制在轨迹起点的?
还是被排除在外?

再次，我们转向graph_test()，这一次传递给它neighbor_graph=“principal_graph”，这告诉它测试轨迹上相似位置的单元格是否有相关表达式:

ciliated_cds_pr_test_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
pr_deg_ids <- row.names(subset(ciliated_cds_pr_test_res, q_value < 0.05))

下面是graph_test()给出的一些非常重要的有趣基因:

plot_cells(cds, genes = c("hlh-4", "gcy-8", "dac-1", "oig-8"), show_trajectory_graph = FALSE,label_cell_groups = FALSE, label_leaves = FALSE)

和前面一样，我们可以将轨迹可变的基因收集到模块中:

gene_module_df <- find_gene_modules(cds[pr_deg_ids, ], resolution = c(10^seq(-6,-1)))

在这里，我们绘制了Packer和Zhu等人注释的每组细胞类型的模块总分:

cell_group_df <- tibble::tibble(cell = row.names(colData(cds)), cell_group = colData(cds)$cell.type)
agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
pheatmap::pheatmap(agg_mat, scale = "column", clustering_method = "ward.D2")

我们还可以将gene_module_df传递给plot_cells()。

plot_cells(cds, genes = gene_module_df %>%filter(module %in% c(27, 10, 7, 30)), label_cell_groups = FALSE, show_trajectory_graph = FALSE)

Monocle提供了另一个绘图功能，有时可以更清楚地看到基因沿单一路径的动态。您可以使用choose_cells()选择一个路径，也可以通过按cluster、细胞类型或其他限制于该路径的注释对细胞数据集进行子集设置。让我们选择一个这样的路径，AFD细胞:

AFD_genes <- c("gcy-8", "dac-1", "oig-8")
AFD_lineage_cds <- cds[rowData(cds)$gene_short_name %in% AFD_genes, colData(cds)$cell.type %in%c("AFD")]

函数plot_genes_in_pseudotime()取一小组基因，并显示它们作为拟时间函数的动态，你可以看到dac -1在其他两个基因之前被激活。

plot_genes_in_pseudotime(AFD_lineage_cds, color_cells_by = "embryo.time.bin", min_expr = 0.5)

Analyzing branches in single-cell trajectories

分析在轨迹分支点周围被调节的基因通常可以帮助我们深入了解控制细胞命运决定的遗传回路。Monocle可以萃取与系统中的命运决定相对应的分支点。这样做就像选择感兴趣的细胞(和分支点)一样简单choose_cells():

cds_subset <- choose_cells(cds)

然后在子集上调用graph_test()。这将识别出具有有趣的表达模式的基因，这些基因只在您所选择的轨迹区域内，从而为您提供更精确和相关的基因集。

subset_pr_test_res <- graph_test(cds_subset, neighbor_graph = "principal_graph",cores = 4)
pr_deg_ids <- row.names(subset(subset_pr_test_res, q_value < 0.05))

将这些基因分组到模块中可以揭示命运特定基因或那些在分支点之前或之后立即被激活的基因:

gene_module_df <- find_gene_modules(cds_subset[pr_deg_ids, ], resolution = 0.001)

我们将根据模块在轨迹上的相似度(使用hclust)来组织模块，这样更容易看到哪些模块先于其他模块启动:

agg_mat <- aggregate_gene_expression(cds_subset, gene_module_df)
module_dendro <- hclust(dist(agg_mat))
gene_module_df$module <- factor(gene_module_df$module, levels = row.names(agg_mat)[module_dendro$order])plot_cells(cds_subset, genes = gene_module_df, label_cell_groups = FALSE, show_trajectory_graph = FALSE)

Graph-autocorrelation analysis for comparing clusters

在L2蠕虫的数据中，我们发现了一些非常独特的神经元聚类:

# reload and reprocess the data as described in the 'Clustering and classifying
# your cells' section
expression_matrix <- readRDS("cao_l2_expression.rds")
cell_metadata <- readRDS("cao_l2_colData.rds")
gene_annotation <- readRDS("cao_l2_rowData.rds")# Make the CDS object
cds <- new_cell_data_set(expression_matrix, cell_metadata = cell_metadata, gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 100)
cds <- reduce_dimension(cds)
cds <- cluster_cells(cds, resolution = 1e-05)colData(cds)$assigned_cell_type <- as.character(partitions(cds))
colData(cds)$assigned_cell_type <- dplyr::recode(colData(cds)$assigned_cell_type,`1` = "Body wall muscle", `2` = "Germline", `3` = "Motor neurons", `4` = "Seam cells",`5` = "Sex myoblasts", `6` = "Socket cells", `7` = "Marginal_cell", `8` = "Coelomocyte",`9` = "Am/PH sheath cells", `10` = "Ciliated neurons", `11` = "Intestinal/rectal muscle",`12` = "Excretory gland", `13` = "Chemosensory neurons", `14` = "Interneurons",`15` = "Unclassified eurons", `16` = "Ciliated neurons", `17` = "Pharyngeal gland cells",`18` = "Unclassified neurons", `19` = "Chemosensory neurons", `20` = "Ciliated neurons",`21` = "Ciliated neurons", `22` = "Inner labial neuron", `23` = "Ciliated neurons",`24` = "Ciliated neurons", `25` = "Ciliated neurons", `26` = "Hypodermal cells",`27` = "Mesodermal cells", `28` = "Motor neurons", `29` = "Pharyngeal gland cells",`30` = "Ciliated neurons", `31` = "Excretory cells", `32` = "Amphid neuron",`33` = "Pharyngeal muscle")# Reload and re-process the data as described in the 'Clustering and
# classifying your cells' section In some clusters the majority of cells may
# have NA as the cell type so use this function with caution.
consensus <- function(x) {uniqx <- unique(na.omit(x))uniqx[which.max(tabulate(match(x, uniqx)))]
}l_partition <- unique(partitions(cds))
l_cell_type <- list()
for (i_partition in l_partition) {l_cell_type[as.character(i_partition)] <- consensus(colData(cds)[partitions(cds) ==i_partition, ][["cao_cell_type"]])
}
colData(cds)[["assigned_cell_type"]] <- as.character(l_cell_type[as.character(partitions(cds))])

只对神经元进行子集:

neurons_cds <- cds[, grepl("neurons", colData(cds)$assigned_cell_type, ignore.case = TRUE)]
plot_cells(neurons_cds, color_cells_by = "partition")

神经元有许多子类型，所以可能不同的神经元簇对应不同的子类型。为了研究哪些基因在cluster中表达差异，我们可以使用上面讨论的回归分析工具。然而，Monocle提供了一种寻找UMAP或t-SNE空间中不同细胞组间差异基因的替代方法。函数graph_test()使用了来自空间自相关分析的统计数据，称为Moran’s I, Cao和Spielmann等人证明该统计数据可以有效地发现单细胞RNA-seq数据集中变化的基因。

你可以像这样运行graph_test():

pr_graph_test_res <- graph_test(neurons_cds, neighbor_graph = "knn", cores = 8)
pr_deg_ids <- row.names(subset(pr_graph_test_res, q_value < 0.05))

数据框pr_graph_test_res包含cell_data_set中每个基因的Moran’s I测试结果。如果您想按效应大小对基因进行排序，请按morans_Icolumn对该表进行排序，其范围从-1到+1。值为0表示没有影响，而+1表示完全正自相关，表明附近的细胞具有非常相似的基因表达值。远小于零的显著值通常很少。

阳性值表示基因在UMAP空间的一个焦点区域表达(例如，特定于一个或多个簇)。但是我们如何把基因和簇联系起来呢?

Finding modules of co-regulated genes

一旦有了一组基因，在簇中以某种有趣的方式发生变化，Monocle就提供了一种将它们分组为模块的方法。你可以调用find_gene_modules()，它实际上是在基因上运行UMAP(与细胞相反)，然后使用Louvain社区分析将它们分组到模块中:

gene_module_df <- find_gene_modules(neurons_cds[pr_deg_ids, ], resolution = 0.01)

数据帧gene_module_df包含每个基因的一行，并标识它所属的模块。要查看哪些模块在哪些集群或分区中表示，可以使用两种不同的可视化方法。第一种方法是制作一个简单的表，显示所有集群中每个模块中所有基因的聚合表达。Monocle为此提供了一个简单的实用函数aggregate_gene_expression:

cell_group_df <- tibble::tibble(cell = row.names(colData(neurons_cds)), cell_group = partitions(cds)[colnames(neurons_cds)])
agg_mat <- aggregate_gene_expression(neurons_cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
colnames(agg_mat) <- stringr::str_c("Partition ", colnames(agg_mat))pheatmap::pheatmap(agg_mat, cluster_rows = TRUE, cluster_cols = TRUE, scale = "column",clustering_method = "ward.D2", fontsize = 6)

有些模块高度特定于单元的某些分区，而其他模块则跨多个分区共享。注意，aggregate_gene_expression可以用于细胞和基因的任意分组。您不局限于查看来自find_gene_modules()、clusters()和partitions()的模块。

查看模块及其表达式的第二种方法是将gene_module_df直接传递给plot_cells()。如果有很多模块，就很难看到每个模块在哪里表示，所以我们只看它们的一个子集:

plot_cells(neurons_cds, genes = gene_module_df %>%filter(module %in% c(8, 28, 33, 37)), group_cells_by = "partition", color_cells_by = "partition",show_trajectory_graph = FALSE)

我们这期主要介绍差异表达分析，内容过多，一次完成有点太乱了，目前单细胞测序的费用也在降低，单细胞系列可算是目前的测序神器，有这方面需求的老师，联系桓峰基因，提供最高端的科研服务！

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

有想进生信交流群的老师可以扫最后一个二维码加微信，备注“单位+姓名+目的”，有些想发广告的就免打扰吧，还得费力气把你踢出去！

References:

Trapnell C. et. al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381-386 (2014).
Qiu, X. et. al. Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14, 979-982 (2017).
Cao, J. et. al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496-502 (2019).
Haghverdi, L. et. al. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat. Biotechnol. 36, 421-427 (2018).
McInnes, L., Healy, J. & Melville, J. UMAP: Uniform Manifold Approximation and Projection for dimension reduction. Preprint at https://arxiv.org/abs/1802.03426 (2018).
Traag, V.A., Waltman, L. & van Eck, N.J. From Louvain to Leiden: guaranteeing well-connected communities. Scientific Reportsvolume 9, Article number: 5233 (2019).
Levine, J. H. et. al. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell 162, 184-197 (2015).
Levine, J. H., et. al. Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. Cell 162, 184-197 (2015).