Monocle2 | 单细胞测序的拟时序分析（细胞轨迹分析）

2021-12-06 本文已影响0人夕颜00

背景介绍

伪时间是衡量单个细胞在细胞分化等过程中取得了多大进展的指标。在许多生物学过程中，细胞并不是完全同步的。在细胞分化等过程的单细胞表达研究中，捕获的细胞在分化方面可能分布广泛。也就是说，在同一时间捕获的细胞群中，有些细胞可能已经很长时间了，而有些细胞甚至还没有开始这个过程。当您想要了解在细胞从一种状态转换到另一种状态时所发生的调节更改的顺序时，这种异步性会产生主要问题。跟踪同时捕获的细胞间的表达可以产生对基因动力学一个大致的认识，该基因表达的明显变异性将非常高。Monocle根据每个cell在学习轨迹上的进展对其进行排序，从而缓解了由于异步而产生的问题。Monocle不是跟踪表达式随时间变化的函数，而是跟踪沿轨迹变化的函数，我们称之为伪时间。伪时间是一个抽象的分化单位：它只是一个cell到轨迹起点的距离，沿着最短路径测量。轨迹的总长度是由细胞从起始状态移动到结束状态所经历的总转录变化量来定义的。

还是老规矩，跟着官网学习Monocle2的使用方法。官网如下：

http://cole-trapnell-lab.github.io/monocle-release/docs/

一、Abstract

单细胞基因表达研究使人们能够在复杂的生物过程和高度异质性的细胞群中描述转录调控。这些研究有助于发现识别特定亚型细胞的基因，或标记生物过程中的中间状态，以及在两种不同的细胞命运之间过渡态的基因。在许多单细胞研究中，单个细胞是以不同步的方式执行基因表达过程。实际上，每个细胞都是正在研究的转录过程中的一个瞬间。Monocle包是分析单细胞测序的工具。Monocle引入了在伪时间(拟时间)内对单个细胞排序的策略，利用单个细胞的非同步进程，将它们置于与细胞分化等生物学过程相对应的轨迹上。****Monocle利用先进的机器学习技术(反向图嵌入)从单细胞数据中学习显式的主图来对细胞进行排序，这可以强大而准确地解决复杂的生物过程。Monocle也可以进行聚类(即使用t-SNE和密度峰值聚类)。Monocle也可以进行差异基因表达测试，使人们能够识别在不同状态下差异表达的基因，沿着生物过程以及不同的细胞命运时基因表达的变化。Monocle是专为单细胞RNA-Seq研究设计的，但也可以用于其他分析。

二、Introduction

Monocle 2包括新的和改进的算法，用于对细胞进行分类和计数，执行细胞亚群之间的差异表达分析，以及重建细胞轨迹。Monocle 2也被重新设计，可以很好地完成包含数万个或更多细胞的大型单细胞RNA-Seq数据分析。

Monocle perform three main types of analysis:

Clustering, classifying, and counting cells. Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them.
Constructing single-cell trajectories. In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions.
Differential expression analysis. Characterizing new cell types and states begins with comparing them to other, better understood cells. Monocle includes a sophisticated but easy to use system for differential expression.

首先，Monocle 2使用一种简单的、无偏的和高度可扩展的统计程序来选择具有轨迹进展特征的基因。然后，它采用了一类流形学习算法，旨在在高维单细胞RNA-seq数据中嵌入一个主图。以前的方法是通过启发式分析细胞之间的成对距离来推断分支结构，而Monocle 2可以使用这张图来直接识别发育的命运决定。我们已经通过广泛的基准测试证明，Monocle 2优于其他工具，如Wishbone，而不需要用户指定轨迹的结构。

三、Monocle2算法

Monocle 2的算法是再2017年发表在nature methods 上，文章链接如下

https://www.nature.com/articles/nmeth.4402

文章中的核心理论为：每个细胞都可以表示为高维空间中的一个点，在高维空间中，每个维对应着一个有序基因的表达水平。高维数据首先通过几种降维方法，如PCA(默认)、扩散映射等，投射到低维空间。Monocle 2然后在自动选择的一组数据质心上构造一棵生成树。然后，该算法将细胞移动到它们最近的树的顶点，更新顶点的位置以适应细胞，学习新的生成树，并迭代地继续这个过程，直到树和细胞的位置已经收敛。在这个过程中，Monocle 2保持了高维空间和低维空间之间的可逆映射，从而既学习了轨迹，又降低了数据的维数。****一旦Monocle 2学会了树，用户就会选择一个tip作为根。计算每个单元的伪时间作为其沿树到根的测地线距离，并根据主图自动分配其分枝。因为monocle2学习树结构，与其他方法相比，分支结构自动出现。当它更新细胞位置并细化树时，monocle2简化了轨迹的结构，修剪了小的分支，这样最终的轨迹只保留了描述细胞状态显著差异的分支。

图片

四、Installing Monocle2

Monocle需要在R语言环境下运行。官网给出的都是R-3.4比较老的版本，现在都是R-3.6或者R-4.0以上版本。我们可以直接在Bioconductor中安装Monocle 2。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("monocle")

查看系统中安装的这个包的版本文档

browseVignettes("monocle")

图片

当然也可以在GitHub上安装最新版的Monocle2。

install.packages("devtools")
devtools::install_github("cole-trapnell-lab/monocle-release@develop")

有时我们添加的特性需要你安装某些额外的软件包。当您尝试上面的命令时，可能会看到错误。您可以通过输入(例如)来在错误消息中安装包。

BiocManager::install(c("DDRTree", "pheatmap"))

最后检查一下是否安装成功

library(monocle)

四、创建CellDataSet

方法一: 将Seurat object中数据提取来创建

data <- as(as.matrix(scRNA.Osteoclastic@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = scRNA.Osteoclastic@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
monocle_cds <- newCellDataSet(data,
                              phenoData = pd,
                              featureData = fd,
                              lowerDetectionLimit = 0.5,
                              expressionFamily = negbinomial.size())

FPKM/TPM值通常是对数正态分布的，而UMIs或读计数使用负二项更好地建模。要处理计数数据，请将负二项分布指定为newCellDataSet的expressionFamily参数：

稀疏矩阵用negbinomial.size()，

FPKM值用tobit()，

logFPKM值用gaussianff()

图片

方法二: 直接读取表达矩阵来创建

library(data.table)
##读取数据
data <- fread("fpkm.txt",data.table = F,header = T)
pd <-  fread("metadata.txt",data.table = F,header = T)
fd <- fread("gene_annotations.txt",data.table = F,header = T)
##创建
pd <- new("AnnotatedDataFrame", data = pd)
fd <- new("AnnotatedDataFrame", data = fd)
HSMM <- newCellDataSet(as.matrix(data),
                       phenoData = pd, featureData = fd,
                       expressionFamily = tobit())
###如果数据量大，建议转化为稀疏矩阵
HSMM <- newCellDataSet(as(as.matrix(data), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       expressionFamily = tobit())

方法三: 将Seurat对象直接转化为CellDataSet对象

importCDS(scRNA.Osteoclastic)

如果我们想要从Seurat对象或SCESet中导入所有的插槽，我们可以设置参数'import_all'为TRUE。#(默认为FALSE或只保留最小数据集)

五、估计size factor和离散度

size facotr帮助我们标准化细胞之间的mRNA的差异。

离散度值可以帮助我们进行后续的差异分析。


monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)

六、过滤低质量的细胞

大多数单细胞工作流程至少会包含一些由死细胞或空孔组成的库。同样重要的是要删除doublets:由两个或多个细胞意外生成的库。这些细胞可以破坏下游步骤，如伪时间排序或聚类。要知道一个特定的基因有多少个表达，或者一个给定的细胞有多少个基因表达，通常是很方便的。Monocle提供了一个简单的函数来计算这些统计数据。

monocle_cds <- detectGenes(monocle_cds, min_expr = 0.1)
print(head(fData(monocle_cds)))

图片

七、细胞分类

Monocle官网教程提供了4个分类方法：

**Classifying cells by type **
Clustering cells without marker genes
Clustering cells using marker genes
Imputing cell type

1、Classifying cells by type

Monocle提供了一个简单的系统，根据您选择的标记基因的表达来标记细胞。您只需提供一组函数，Monocle就可以使用这些函数对每个单元格进行注释。例如，您可以为几种单元格类型中的每种类型提供一个函数。这些函数接受每个单元格的表达式数据作为输入，并返回TRUE以告诉Monocle某个单元格满足函数定义的条件。所以你可以有一个功能，对表达成肌细胞特异性基因的细胞来说是真的，另一个功能对成纤维细胞特异性基因来说是真的，等等。


MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5"))
ANPEP_id <- row.names(subset(fData(HSMM),
    gene_short_name == "ANPEP"))

cth <- newCellTypeHierarchy()
cth <- addCellType(cth, "Myoblast", classify_func =
    function(x) { x[MYF5_id,] >= 1 })
cth <- addCellType(cth, "Fibroblast", classify_func = function(x)
{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })
HSMM <- classifyCells(HSMM, cth, 0.1)
table(pData(HSMM)$CellType)
###画图
pie <- ggplot(pData(HSMM),
aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1)
pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())

2、Clustering cells without marker genes

Monocle提供了一种算法，可以用来计算“未知”单元格的类型。该算法在函数clusterCells中实现，根据全局表达式轮廓将单元格分组在一起。这样的话，如果你的细胞表达了许多特定于成肌细胞的基因，但恰好缺少MYF5，我们仍然可以识别它为成肌细胞。clusterCells可以在无监督的方式下使用，也可以在半监督的“”模式下使用，这允许使用一些专家知识来辅助算法。我们先来看看无监督模式。

第一步是决定哪些基因用于细胞聚类。我们可以使用所有的基因，但是我们会包含很多基因，这些基因的表达水平不足以提供有意义的信号。包括它们只会给系统增加噪音。我们可以根据平均表达水平筛选基因，我们还可以选择在细胞间异常可变的基因。这些基因往往对细胞状态有很高的信息量。

HSMM=monocle_cds
disp_table <- dispersionTable(HSMM)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
plot_ordering_genes(HSMM)

图片

setOrderingFilter函数标记了在随后对clusterCells的调用中用于聚类的基因，尽管我们可以根据需要提供其他基因列表。plot_ordering_genes函数显示了基因表达的变异性（离散性）如何依赖于细胞间的平均表达。红线表示基于这种关系的单片体对色散的预期。我们标记用于聚类的基因显示为黑点，而其他基因显示为灰点。

plot_pc_variance_explained(HSMM, return_all = F)

图片

Monocle使用t-SNE来给细胞聚类可视化

HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 6,
                reduction_method = 'tSNE', verbose = T)
HSMM <- clusterCells(HSMM, num_clusters = 2)
plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

图片

指定10个cluster

HSMM <- clusterCells(HSMM, num_clusters = 10)
plot_cell_clusters(HSMM)

图片

3、Clustering cells using marker genes
首先，我们将选择一组不同的基因用于细胞聚类。在我们挑选高表达和高度可变的基因之前。现在，我们将挑选与标记共变异的基因。从某种意义上说，我们将建立一个大的基因列表来作为标记，这样即使一个细胞没有MYF5，它也可以被其他基因识别为成肌细胞。

marker_diff <- markerDiffTable(HSMM[expressed_genes,],
            cth,
            residualModelFormulaStr = "~Media + num_genes_expressed",
            cores = 1)

函数markerDiffTable接受CellDataSet和CellTypeHierarchy，并根据提供的函数将所有单元格分类为类型。然后在识别不同基因类型之间差异表达的基因之前，它去除了所有“未知”和“不明确”的功能。同样，您可以提供要从该测试中排除的影响的剩余模型。然后函数返回一个测试结果的数据帧，您可以使用它来选择要用于聚类的基因。通常情况下，最好挑选最适合每种细胞类型的前10或20个基因。这确保了聚类基因不受一种细胞类型标记的支配。如果可能的话，通常需要为每种类型创建一个平衡的标记面板。Monocle提供了一个方便的功能，可以根据基因在每种类型中的表达受限程度来对基因进行排序。

candidate_clustering_genes <-
    row.names(subset(marker_diff, qval < 0.01))
marker_spec <-
  calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)
head(selectTopMarkers(marker_spec, 3))

semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)
HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)
plot_ordering_genes(HSMM)

plot_pc_variance_explained(HSMM, return_all = F)

HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 3,
  norm_method = 'log',
  reduction_method = 'tSNE',
  residualModelFormulaStr = "~Media + num_genes_expressed",
  verbose = T)
HSMM <- clusterCells(HSMM, num_clusters = 2)
plot_cell_clusters(HSMM, 1, 2, color = "CellType")

4、Imputing cell type
注意，我们已经减少了成肌细胞群中“污染”成纤维细胞的数量，反之亦然。但是“未知”细胞呢？如果为clusterCells提供CellTypeHierarcy，Monocle将使用它对整个簇进行分类，而不仅仅是单个细胞。基本上，clusterCells的工作原理与以前完全相同，只是在构建集群之后，它会统计每个集群中每种细胞类型的频率。当一个簇由超过一定百分比（在本例中为10%）的特定类型组成时，簇中的所有单元格都被设置为该类型。如果一个集群由多个单元类型组成，则整个集群都被标记为“不明确”。如果没有高于阈值的单元类型，则集群被标记为“未知”。因此，Monocle可以帮助您在缺少标记数据的情况下估算每个细胞的类型。

HSMM <- clusterCells(HSMM,
              num_clusters = 2,
              frequency_thresh = 0.1,
              cell_type_hierarchy = cth)
plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

pie <- ggplot(pData(HSMM), aes(x = factor(1), fill =
    factor(CellType))) +
geom_bar(width = 1)
pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())

八、构建轨迹

在发育过程中，为了对刺激做出反应，在整个生命周期中，细胞从一种功能“状态”过渡到另一种功能“状态”。处于不同状态的细胞表达不同的基因，产生蛋白质和代谢物的动态重复，从而完成它们的工作。当细胞在不同的状态间移动时，会经历一个转录重组的过程，一些基因被沉默，而另一些则被激活。这些瞬态通常很难描述，因为在更稳定的端点状态之间净化细胞可能很困难或不可能。单细胞RNA-Seq可以让你看到这些状态而不需要纯化。然而，要做到这一点，我们必须确定每个细胞的可能状态范围。
Monocle介绍了利用RNA-Seq进行单细胞轨迹分析的策略。Monocle不是通过实验将细胞净化成离散状态，而是使用一种算法来学习每个细胞作为动态生物过程的一部分必须经历的基因表达变化序列。一旦掌握了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞放置在轨迹中的适当位置。然后，你可以使用Monocle的差异分析工具来寻找在轨迹过程中受调控的基因，如“寻找随假时间变化的基因”一节所述。如果这个过程有多个结果，Monocle将重建一个“分支”轨迹。这些分支与细胞的“决定”相对应，Monocle提供了强有力的工具来识别受它们影响并参与制造它们的基因。在“分析单细胞轨迹中的分支”部分中，可以看到如何分析分支。Monocle依靠一种称为反向图嵌入的机器学习技术来构造单细胞轨迹.

The ordering workflow

Step 1: choosing genes that define progress
Step 2: reducing the dimensionality of the data
Step 3: ordering the cells in pseudotime

Step 1: 选择定义过程的基因

推断单细胞轨迹是一个机器学习问题。第一步是选择Monocle将用作机器学习方法输入的基因。这叫做特征选择，它对轨迹的形状有很大的影响。在单细胞RNA-Seq中，低水平表达的基因通常非常嘈杂，但有些基因包含有关细胞状态的重要信息。单核细胞通过检测这些基因在细胞群体中的表达模式来对细胞进行排序。Monocle寻找以“有趣的”（即不只是嘈杂的）方式变化的基因，并用它们来构造数据。Monocle为您提供了多种工具来选择基因，这些基因将产生一个健壮、准确和具有生物学意义的轨迹。你可以使用这些工具来执行一个完全“无监督”的分析，在这个分析中，Monocle不知道你认为哪个基因是重要的。或者，你可以利用marker gene知识来定义生物学进展，从而形成Monocle的轨迹。我们认为这种模式是“半监督”的分析。

Monocle官网教程提供了4个选择方法：

选择发育差异表达基因
选择clusters差异表达基因
选择离散程度高的基因
自定义发育marker基因

前三种都是无监督分析方法，细胞发育轨迹生成完全不受人工干预；最后一种是半监督分析方法，可以使用先验知识辅助分析。

##使用clusters差异表达基因
deg.cluster <- FindAllMarkers(scRNA.Osteoclastic)
diff.genes <- subset(deg.cluster,p_val_adj<0.05)$gene
HSMM <- setOrderingFilter(HSMM, diff.genes)
plot_ordering_genes(HSMM)
##使用seurat选择的高变基因
var.seurat <- VariableFeatures(scRNA.Osteoclastic)
HSMM <- setOrderingFilter(HSMM, var.genes)
plot_ordering_genes(HSMM)
##使用monocle选择的高变基因
disp_table <- dispersionTable(HSMM)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
HSMM <- setOrderingFilter(HSMM, disp.genes)
plot_ordering_genes(HSMM)

图片

理想情况下，我们希望尽可能少地使用正在研究的系统生物学的先验知识。我们希望从数据中发现重要的排序基因，而不是依赖于文献和教科书，因为这可能会在排序中引入偏见。我们将从一种更简单的方法开始，但是我们通常推荐一种更复杂的方法，称为“dpFeature”。


diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
              fullModelFormulaStr = "~Media")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
plot_ordering_genes(HSMM_myo)

Step 2: 降维

下一步，我们将把空间缩小到一个二维空间，当Monocle对细胞进行排序时，我们将能够很容易地可视化和解释这些空间。

HSMM <- reduceDimension(HSMM, max_components = 2,
    method = 'DDRTree')

Step 3: 按照轨迹排序细胞

HSMM <- orderCells(HSMM)
plot_cell_trajectory(HSMM, color_by = "seurat_clusters")

图片

plot_cell_trajectory(HSMM, color_by = "State")

图片

plot_cell_trajectory(HSMM, color_by = "Pseudotime")

图片

如果你的树上有一堆状态，那就很难弄清楚每个状态都落在树上的什么地方了。有时，可以方便地“分面”轨迹图，以便更容易地查看每个状态的位置.

plot_cell_trajectory(HSMM, color_by = "State") +
  facet_wrap(~State, nrow = 1)

图片

如果你没有时间序列，你可能需要利用你对系统的生物学知识，根据特定标记基因的表达位置来设置根。例如，在这个实验中，一个高度增殖的祖细胞群正在产生两种类型的有丝分裂后细胞。因此，根部应该有表达高水平增殖标记的细胞。我们可以使用抖动图来判断哪种状态对应于快速扩散：

blast_genes <- row.names(subset(fData(HSMM),
                                gene_short_name %in% c("CCNB2", "MYOD1", "MYOG")))
plot_genes_jitter(HSMM[blast_genes,],
                  grouping = "State",
                  min_expr = 0.1)

图片

HSMM_expressed_genes <-  row.names(subset(fData(HSMM),
                                          num_cells_expressed >= 10))
HSMM_filtered <- HSMM[HSMM_expressed_genes,]
my_genes <- row.names(subset(fData(HSMM_filtered),
                             gene_short_name %in% c("YWHAB", "GAPDH", "TNNC1")))
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")

图片

plot_genes_in_pseudotime(cds_subset, color_by =  "State")

图片


genes <- c("TNNT2", "TNNC1", "CDK1")
p1 <- plot_genes_jitter(HSMM[genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(HSMM[genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(HSMM[genes,], color_by = "State")
plotc <- p1|p2|p3

图片

九、差异分析

官方给出的差异分析有三大方法，我们重点关注第三个：****根据伪时间功能寻找差异基因

1、Basic Differential Analysis
2、Finding Genes that Distinguish Cell Type or State
3、Finding Genes that Change as a Function of Pseudotime

Monocle的主要工作是通过生物过程（如细胞分化）将细胞按顺序排列，而不知道要提前查看哪些基因。一旦这样做了，你就可以分析细胞，找到随着细胞进展而变化的基因。例如，你可以发现随着细胞“成熟”而显著上调的基因。让我们来看看一组对肌肉生成很重要的基因：

to_be_tested <- row.names(subset(fData(HSMM),
                                 gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
cds_subset <- HSMM[to_be_tested,]

diff_test_res <- differentialGeneTest(cds_subset,
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")

diff_test_res[,c("gene_short_name", "pval", "qval")]
plot_genes_in_pseudotime(cds_subset, color_by ="State")

这个sm.ns函数指出Monocle应该通过表达式值拟合自然样条曲线，以帮助它将表达式的变化描述为进程的函数。让我们加入基因注释，这样就很容易看出哪些基因是重要的。

图片

十、根据伪时间表达pattern聚类基因

** 研究时间序列基因表达研究时出现的一个常见问题是：“哪些基因遵循相似的动力学趋势”？Monocle可以通过对具有相似趋势的基因进行分组来帮助你回答这个问题，所以你可以分析这些基因组，看看它们有什么共同点。Monocle提供了一种方便的方法来可视化所有假时间依赖性基因。函数plot_pseudotime_heatmap采用CellDataSet对象（通常只包含重要基因的子集），并生成平滑的表达曲线，非常类似于plot_genes_in_pseudotime。然后，它对这些基因进行聚类，并使用pheatmap软件包进行绘图。这让你可以想象出基因模块在不同的时间内共同变化。**

marker_genes <- row.names(subset(fData(HSMM),
                                 gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
                                                        "ANPEP", "PDGFRA","MYOG",
                                                        "TPM1",  "TPM2",  "MYH2",
                                                        "MYH3",  "NCAM1", "TNNT1",
                                                        "TNNT2", "TNNC1", "CDK1",
                                                        "CDK2",  "CCNB1", "CCNB2",
                                                        "CCND1", "CCNA1", "ID1")))

diff_test_res <- differentialGeneTest(HSMM[marker_genes,],
                                      fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
plot_pseudotime_heatmap(HSMM[sig_gene_names,],
                        num_clusters = 6,
                        cores = 1,
                        show_rownames = T)

图片

十一、单细胞轨迹的“分支”分析

通常，单细胞轨迹包括分支。这些分支的产生是因为细胞执行不同的基因表达程序。在发育过程中，当细胞做出命运选择时，分支出现在轨迹中：一个发育谱系沿着一条路径前进，而另一个谱系产生第二条路径。Monocle包含分析这些分支事件的广泛功能。

以Steve Quake's的实验室进行的实验为例，Barbara Treutlein团队从发育中的小鼠肺中捕获了细胞。他们在发育早期捕获细胞，后来当肺中同时含有两种主要类型的上皮细胞（AT1和AT2）时，细胞就要决定变成AT1或AT2。Monocle可以将这个过程重构为一个分支轨迹，允许您对决策点进行非常详细的分析。下图显示了使用部分数据重建的Monocle轨迹。有一个分支，标记为“1”。当细胞从树的左上角经过树枝的早期发育阶段时，哪些基因会发生变化？哪些基因在分枝间有差异表达？为了回答这个问题，Monocle为您提供了一个特殊的统计测试：分支表达式分析建模，或BEAM。BEAM(Branched expression analysis modeling)是一种统计方法，用于寻找以依赖于分支的方式调控的基因。

plot_cell_trajectory(HSMM, color_by = "State")

图片

BEAM进行统计分析

BEAM_res <- BEAM(HSMM, branch_point = 1, cores = 8)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
plot_genes_branched_heatmap(HSMM[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1,
                            num_clusters = 4,
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)

该热图显示的是同一时间点两个谱系的变化。热图的列是伪时间的点，行是基因。从热图中间往右读，是伪时间的一个谱系；往左是另一个谱系。基因是被按照等级聚类的，所以你看到的基因表达模式和谱系的表达模式是非常相似的。

图片

We can plot a couple of these genes, such as Pdpn and Sftpb (both known markers of fate in this system), using the plot_genes_branched_pseudotime function, which works a lot like the plot_genes_in_pseudotime function, except it shows two kinetic trends, one for each lineage, instead of one.

genes <- row.names(subset(fData(HSMM),
                          gene_short_name %in% c( "MEF2C", "CCNB2", "TNNT1")))

plot_genes_branched_pseudotime(HSMM[genes,],
                               branch_point = 1,
                               color_by = "State",
                               ncol = 1)

图片

转载来自：
Monocle2 | 单细胞测序的拟时序分析（细胞轨迹分析）