monocle2轨迹分析

创建CellDataSet对象

CellDataSet对象需要三个输入文件，分别是：

• exprs：基因表达的数字矩阵，其中行是基因，列是细胞
• phenoData：行是细胞，列是细胞的属性（诸如细胞类型、培养条件、捕获日期）
• featureData：行是基因，列是基因的属性

创建CellDataSet对象的代码如下所示：

library(tidyverse)
library(monocle)

library(HSMMSingleCell)
data(HSMM_expr_matrix)
data(HSMM_gene_annotation)
data(HSMM_sample_sheet)
pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)

# First create a CellDataSet from the relative expression levels
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),   
                       phenoData = pd, 
                       featureData = fd,
                       lowerDetectionLimit=0.1,
                       expressionFamily=tobit(Lower=0.1))
# Next, use it to estimate RNA counts
rpc_matrix <- relative2abs(HSMM, method = "num_genes")

head(rpc_matrix)

# Now, make a new CellDataSet using the RNA counts
HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
                       phenoData = pd, 
                       featureData = fd,
                       lowerDetectionLimit=0.5,
                       expressionFamily=negbinomial.size())

数据分布的选择

Monocle可以很好地处理相对表达数据和基于计数的测量（例如UMI）。一般来说，它与转录计数数据（尤其是UMI数据）配合得最好。无论您的数据类型是什么，都必须为其指定适当的分布。FPKM/TPM值通常为对数正态分布，而UMI或读取计数最好使用负二项分布进行建模。

当我们使用的是Census mRNA计数值，我们改变了lowerDetectionLimit的值，以反映新的表达量表。重要的是，我们还将expressionFamily设置为negbinomal()，这告诉Monocle在某些下游统计测试中使用负二项分布。如果不改变这两个选项，可能会在以后产生问题，所以在使用Census 统计时，一定不要忘记它们。

Seurat对象转换为CellDataSet对象

phenotype_data <- SELECTCELL@meta.data |> dplyr::rename(Sample = sampleinf)
exprs_data  <- SELECTCELL@assays$RNA@counts
feature_data <- data.frame(gene_short_name=rownames(exprs_data),row.names = rownames(exprs_data))

# UMIs or read counts are better modeled with the negative binomial.
pd <- new("AnnotatedDataFrame", data = phenotype_data)
fd <- new("AnnotatedDataFrame", data = feature_data)
CDS <- newCellDataSet(as.matrix(exprs_data), 
                       phenoData = pd, 
                       featureData = fd,
                       expressionFamily=negbinomial.size())
pData(CDS)
fData(CDS)
exprs(CDS)

Estimate size factors and dispersions

• size factors：尺寸因子有助于我们对细胞间回收的 mRNA 差异进行标准化
• dispersions：帮助我们稍后执行差异表达式分析

HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)

Filtering low-quality cells

Constructing Single Cell Trajectories

Trajectory step 1: choose genes that define a cell's progress

通过差异表达选择基因

首先，我们必须决定我们将使用哪些基因来定义细胞通过肌肉生成的过程。我们最终想要一组基因，在我们正在研究的过程中，随着肌肉生成的过程而增加（或减少）表达。

获得一组随着肌肉生成的过程表达量增加或减少的基因，一种很容易想到的方法是比较在肌肉生成的过程开始时收集的细胞和在结束时收集的细胞，找到差异表达的基因。或者说是从生长培养基到分化培养基的转换而差异表达的所有基因。

HSMM <- detectGenes(HSMM, min_expr = 0.1)
expressed_genes <- row.names(subset(fData(HSMM), num_cells_expressed >= 10))
diff_test_res <- differentialGeneTest(HSMM[expressed_genes,], fullModelFormulaStr = "~Media")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))

一旦选择了需要排序的基因，就可以将其存储在monocle的对象中

HSMM <- setOrderingFilter(HSMM, ordering_genes)
plot_ordering_genes(HSMM)

通过无监督学习选择基因

monocle提供了一个称为dpfeature的程序选择基因，进行细胞排序

选择高分散度(dispersion)的基因

选择已知的marker基因

Trajectory step 2: reduce data dimensionality

接下来，我们将把空间缩小到二维空间，这样我们就可以在Monocle对细胞进行排序时轻松地可视化和解释。

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

Trajectory step 3: order cells along the trajectory

Monocle事先不知道该将树的哪个轨迹称为“开始”，因此我们经常不得不使用root_state参数调用orderCells来指定开始。如果数据包含时间序列，这里我们可以使用时间零点作为轨迹的开始。

如果没有时间序列，需要利用你对系统的生物学知识，根据某些标记基因的表达位置来确定根源。例如，一个高度增殖的祖细胞群体正在产生两种类型的有丝分裂后细胞。所以根部应该有表达高水平增殖标记物的细胞。我们可以使用jitter来找出哪个状态对应于快速增殖。

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

可以看到位于零时刻的细胞位于轨迹的一个顶端（红圈的地方）

接下来我们通过转态（State）着色

可以看到大多数零时刻的细胞位于State 1，因此我们将State1设置为轨迹的起点。

HSMM <- orderCells(HSMM, root_state = 1)

也可以通过函数识别包含从时间零点开始的大多数cell格的状态

GM_state <- function(cds){
  if (length(unique(pData(cds)$State)) > 1){
    T0_counts <- table(pData(cds)$State, pData(cds)$Hours)[,"0"]
    return(as.numeric(names(T0_counts)[which
          (T0_counts == max(T0_counts))]))
  } else {
    return (1)
  }
}
HSMM <- orderCells(HSMM, root_state = GM_state(HSMM))

定位细胞分化的起点后，我们就可以绘制细胞沿伪时间分化的轨迹

plot_cell_trajectory(HSMM, color_by = "Pseudotime")

Finding Genes that Change as a Function of Pseudotime

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

gene <- diff_test_res|> rownames_to_column("gene") |>
    arrange(qval)|> pull(gene)
HSMM_subset <- HSMM[gene[1:5],]
plot_genes_in_pseudotime(HSMM_subset, color_by = "Hours")

Clustering Genes by Pseudotemporal Expression Pattern

https://www.jianshu.com/p/a8c9ff0fe4a8

研究时间序列基因表达研究时出现的一个常见问题是：“哪些基因遵循类似的动力学趋势（kinetic trends）”？Monocle可以通过将具有相似趋势的基因分组来帮助你回答这个问题，这样你就可以分析这些基因组，看看它们有什么共同点。

函数plot_pseudotime_heatmap获取CellDataSet对象（通常只包含重要基因的子集），并生成与plot_genes_in_pseudotime类似的平滑表式曲线。然后，它对这些基因进行聚类，并使用pheatmap软件包对它们进行绘图。这可以让你可视化基因模块，它们在假时间内共同变化。

sig_gene_names <- row.names(subset(diff_test_res, qval < 0.05))
library(parallel)
png(file = "a.png")
plot_pseudotime_heatmap(HSMM[sig_gene_names[1:50],],
                        num_clusters = 4,
                        cores = 10,
                        show_rownames = T)
dev.off()

Analyzing Branches in Single-Cell Trajectories

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

Monocle为您提供了一种特殊的统计测试：branched expression analysis modeling, or BEAM.

BEAM将使用orderCells排序的CellDataSet和轨迹中分支点的名称作为输入。它返回每个基因的显著性得分表。得分显著的基因在其表达中具有分支依赖性。

BEAM_res <- BEAM(HSMM, branch_point = 1, cores = 1)

png(file = "a.png")
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)
dev.off()

这张热图同时显示了两个谱系的变化。它还要求您选择branch point。列是伪时间中的点，行是基因，伪时间的开始在热图的中间。当你从热图的中间向右阅读时，你通过伪时间跟随一个谱系。当你向左阅读时，另一个。这些基因是分层聚集的，因此您可以可视化具有类似谱系依赖表达模式的基因模块。

gene <- BEAM_res|> rownames_to_column("gene") |>
    arrange(qval)|> pull(gene)
plot_genes_branched_pseudotime(HSMM[gene[1:5],],
                       branch_point = 1,
                       color_by = "Pseudotime",
                       ncol = 1)