🔹在做单细胞转录组（scRNA-seq）时，我们拿到的往往是一张 “快照” ——成千上万个细胞在某个时刻的基因表达。
问题是：
👉 这些细胞是怎么一步步变化的？
👉 谁是“新生”的干细胞？谁又正走向成熟？
这时就要用到一种特别的方法：拟时序分析（Pseudotime Analysis）。
🌱 为什么要做拟时序？
还原 细胞发育轨迹：像拍一部细胞成长的“延时摄影”。
揭示 状态转变路径：比如免疫细胞如何被激活。
发现 关键调控基因：哪些转录因子在分支点发挥作用。
一句话总结：它能帮我们把静态的快照，拼成一条动态的“故事线”。
🛠️ 拟时序是怎么做的？
降维 📉
把几千个基因的高维表达数据，压缩成可视化的低维空间（PCA、UMAP）。连线成图 🔗
看哪些细胞更像邻居，画出“相似性网络”。找出轨迹 ➡️
在图上拟合一条或多条路径，就像给细胞排了队。标上时间轴 ⏳
设定“起点细胞”，依次给每个细胞一个拟时序值。挖掘关键基因 🔍
看哪些基因在轨迹上逐渐上升、下降，或者在哪些分支点起到“方向盘”作用。🔧 常见的工具
Monocle：最经典的拟时序工具（从 Monocle2 到 Monocle3）。
Slingshot：擅长发现分支轨迹。
PAGA (Scanpy)：更偏向于整体图结构。
SCORPIUS：适合直线型轨迹。
✅Monocle3介绍：Monocle 3 是分析这类数据的常用工具，它的三大功能是：
❖聚类识别细胞亚群；
❖构建发育轨迹（拟时序分析）；
❖做差异表达分析，帮助理解新发现的细胞类型/状态。
本期分享利用monocle3绘制单细胞拟时序分析。
library(monocle3)library(Seurat)library(SeuratObject)library(SeuratWrappers)library(Matrix)library(data.table)
前面，我们通过Seurat得到cell-type的聚类分簇信息，得到rds文件，那我们用之前的rds文件继续往下开展拟时序分析即可。
需要准备3个文件。
expression_matrix：Seurat对象里面的行是gene，列是cell的表达矩阵
cell_metadata:Seurat对象的meta.data信息
gene_mata：里面的gene信息
SRR <- readRDS("SRR.rds")# 原始 count 矩阵（推荐用于 Monocle3）get_counts_all <-function(obj, assay = "RNA", layers =NULL) {if (is.null(layers)) layers <- Layers(obj[[assay]])mats <- lapply(layers, function(ly) GetAssayData(obj, assay = assay, layer = ly))# 统一基因顺序（取交集，避免维度不齐）feats <- Reduce(intersect, lapply(mats, rownames))mats <- lapply(mats, function(m) m[feats, , drop=FALSE])do.call(cbind, mats)}expression_matrix <- get_counts_all(SRR, assay = "RNA",layers = c("counts.SRR11616456", "counts.SRR11616459"))cell_metadata=data.frame(SRR@meta.data)gene_annotation=data.frame(expression_matrix[,1])gene_annotation[,1]=row.names(gene_annotation)colnames(gene_annotation)=c("gene_short_name")##构建Monocle3cds对象cds<-new_cell_data_set(expression_matrix,cell_metadata=cell_metadata,gene_metadata=gene_annotation)#预处理，相当于Seurat流程normalize过程#counts数据选择norm_method=c("log")#data数据选择norm_method=c("none")cds<-preprocess_cds(cds,num_dim=50,norm_method=c("none"))#去除批次效应,多个样本时可以通过此方式去除批次，#根据自己的分类会去批次效应cds<-align_cds(cds,alignment_group="orig.ident")##降维，默认是"Umap"方式，scores是指5个核心cds<-reduce_dimension(cds,cores=5)##聚类分群，分辨率调小cds<-cluster_cells(cds,resolution=0.0000001)
##拟时序cds<-learn_graph(cds)# 可视化轨迹图plot_cells(cds,color_cells_by ="seurat_clusters", # 使用您之前注释的细胞类型label_groups_by_cluster =FALSE,label_leaves =FALSE,label_branch_points =TRUE, # 标记分支点graph_label_size =3)
坐标轴 (UMAP 1 & UMAP 2)
横纵坐标是 UMAP 的两个降维分量，主要用于把高维的基因表达数据映射到二维平面，便于可视化。
每个点代表一 单个细胞。
颜色（color_cells_by = "seurat_clusters"）
每个点的颜色对应 Seurat 聚类分群结果（之前在 Seurat 里注释的 cluster ID），也就是不同的细胞群。
比如图上看到不同颜色的“团块”，就是不同的细胞类型或亚群。
轨迹线 (learn_graph)
黑色的线是
learn_graph()算出来的 发育轨迹骨架（principal graph）。它表示单细胞在拟时序上的可能发育路径。
分支点 (label_branch_points = TRUE)
图上标注的数字（比如 1, 12, 14）就是 Monocle3 推断出的 branch points（分支点）。
分支点代表在细胞发育过程中可能出现的 分化选择，比如一类祖细胞往两个不同方向分化。
plot_cells(cds,color_cells_by ="partition", # 使用您之前注释的细胞类型label_groups_by_cluster =FALSE,label_leaves =FALSE,label_branch_points =TRUE, # 标记分支点graph_label_size =3)
color_cells_by = "partition"：每种颜色是一块 partition，表示在图学习（
learn_graph）时被判定为彼此不连通的细胞集合。Monocle3 会在每个 partition 内单独学习一条（或多条）轨迹；不同 partition 之间不推断发育连续性。黑色折线：
learn_graph()学到的 principal graph（轨迹骨架）。黑圈数字（如 1）：分支点（branch point），提示在该 partition 内可能发生分化分歧。
你图里可见：红色那一大片是一个大的 partition，内部有轨迹和分支点；右侧绿色、右上黄色、左下青色等是彼此独立的 partition——它们各自有一条独立的小轨迹/主干，但不会与红色那片连成一条全局路径。
这意味着什么？
拟时序是“分区内”定义的：当你用
order_cells()设定起点后，只会在该 partition 内得到有效的 pseudotime；其它 partition 会是 NA，或需要分别设起点。不能跨 partition 讲谱系：若生物学上你认为某两块应该连通，但被分成了不同 partition，Monocle3 默认不会连起来。
#指定细胞起点绘图----myselect<-function(cds,select.classify,my_select){cell_ids<-which(colData(cds)[,select.classify]==my_select)closest_vertex<-cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertexclosest_vertex<-as.matrix(closest_vertex[colnames(cds),])root_pr_nodes<-igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]root_pr_nodes}cds<-order_cells(cds,root_pr_nodes=myselect(cds,select.classify='seurat_clusters',my_select="1"))##使用Seurat的UMAP信息，这样可以与Seurat对象的细胞分布保持一致cds.embed<-cds@int_colData$reducedDims$UMAPint.embed<-Embeddings(SRR,reduction="umap")int.embed<-int.embed[rownames(cds.embed),]cds@int_colData$reducedDims$UMAP<-int.embedlibrary(ggsci)library(ggplot2)##不同细胞类型拟时序数值##拟时序值越高表示细胞分化程度越高，这里仅为演示，并非真实分化情况# 生成 16 种颜色colors16 <- c(pal_jco("default")(10), pal_npg("nrc")(6))plot_cells(cds,color_cells_by="pseudotime",show_trajectory_graph=T)+plot_cells(cds,color_cells_by="seurat_clusters",label_cell_groups=FALSE,label_leaves=FALSE,label_branch_points=FALSE,graph_label_size=1)+scale_color_manual(values = colors16)
用你指定的 Seurat cluster=“1” 作为起点计算 pseudotime；
把 Monocle3 的 UMAP 替换为 Seurat 的 UMAP，保证两边坐标一致；
画两幅图：左边按 pseudotime 着色，右边按 seurat_clusters 着色。你现在把两次
plot_cells()用+叠在一起了
Diff_genes <- graph_test(cds, neighbor_graph = "principal_graph", cores = 6)在“主轨迹图（principal_graph）”上做 graph_test（基于 Moran’s I 的空间自相关检验），找出沿轨迹表达有系统变化的基因。
cores=6并行加速。结果数据框里常见列：
morans_I（效应强度，越大越随轨迹变化）、p_value、q_value。
Diff_genes_sig <- Diff_genes %>%top_n(n =10, morans_I) %>% # 取 Moran’s I 前10的基因pull(gene_short_name) %>%as.character()
选取 Moran’s I 最高的 10 个基因的简称作为候选基因
范围：一般在 -1 到 1 之间。
接近 1：基因在相邻细胞上表达高度相似 → 说明它沿轨迹呈现连续变化趋势，很可能和细胞状态转变相关。
接近 0：基因在轨迹上分布比较随机，没有明显的拟时序趋势。
接近 -1：基因在相邻细胞上表现相反模式（少见，多数生物学上意义不大）。
plot_genes_in_pseudotime(cds[Diff_genes_sig, ],color_cells_by = "seurat_clusters",min_expr = 0.5, ncol = 2, cell_size = 1.5) +scale_color_manual(values = colors16)
按 pseudotime 画出这些基因的表达曲线（平滑趋势），并按
seurat_clusters给细胞着色（例如细胞类型）。min_expr=0.5：只显示表达量≥0.5 的数据，减少噪声。ncol=2：每页两列分面；cell_size=1.5：点大小。scale_color_manual(...)：自定义颜色
#-----选择特定基因展示---genes<-c("POSTN","FGFR3")genes_cds<-cds[rowData(cds)$gene_short_name%in%genes,]plot_genes_in_pseudotime(genes_cds,color_cells_by="seurat_clusters",min_expr=0.5,cell_size=1.5)+scale_color_manual(values = colors16)
数据仅供学习使用，无实际意义。好啦，本期分享就到这结束啦，感谢大家关注支持，代码可复制也可打赏后获得R文件及示例数据，祝大家科研顺利~
