单细胞测序—拟时序分析综合
拟时序分析（Pseudotime Analysis）在单细胞测序（Single-cell RNA-seq）中是一个重要的分析步骤，主要用于研究细胞在发育过程或其他生物学过程中所经历的状态变化。与传统的时间序列不同，拟时序分析不依赖于实际的时间信息，而是通过单细胞转录组数据来推测出细胞状态的动态变化轨迹。以下是进行拟时序分析的几个主要原因：
- 揭示细胞分化路径：在细胞分化过程中，不同的细胞会逐渐从一个祖细胞状态转变为不同的分化状态。拟时序分析可以帮助研究者识别出这些分化路径，并且理解在这些路径中哪些基因起了关键作用。
- 理解细胞状态的动态变化：拟时序分析允许研究者推断出细胞状态的连续变化，而不仅仅是静态的分类。这对于理解细胞如何响应环境变化或内部信号非常有用。
- 识别关键调控基因：通过分析细胞在拟时序中的基因表达变化，研究者可以识别出哪些基因在特定时间点或分化阶段起了关键作用，从而深入理解基因调控网络。
- 预测未记录的中间状态：在某些情况下，实际的生物实验可能无法捕捉到所有的细胞状态。拟时序分析可以通过现有的单细胞数据，预测出未被记录的中间状态，从而提供更完整的生物学过程图景。
- 比较不同条件下的细胞发育：通过对比不同条件下的拟时序轨迹，研究者可以理解不同实验处理（如药物作用、基因敲除等）如何影响细胞的发育路径和时间顺序。
通过拟时序分析，可以获得比传统方法更细致的关于细胞命运决定过程的见解，这对于理解复杂的生物过程、疾病机理，以及开发新的治疗策略具有重要意义
1 数据导入
- 准备单细胞RNA测序数据分析环境。
- 加载并过滤数据，排除不需要的细胞类型。
- 为后续的拟时序分析创建必要的文件结构和工作目录。
rm(list = ls())
library(cowplot)
library(dplyr)
library(data.table)
library(tidyverse)
library(paletteer)
library(vcd)
library(limma)
library(gplots)
library(Seurat)
library(clustree)
library(ggplot2)
library(reshape2)
library(dplyr)
library(monocle)
library(RColorBrewer)
library(ggpubr)
library(ggsignif)
library(patchwork)
library(tidydr)
library(ggforce)
library(ggrastr)
library(viridis)
library(gridExtra)
library(ggnewscale)
library(org.Mm.eg.db)
library(ComplexHeatmap)
library(ClusterGVis)
load(file = 'first_sce3.Rdata')
table(sce$celltype)
sce <- sce[,!(sce$celltype %in% c( "Carcer","CAF5","CAF6"))]
dir.create("monocle2")
setwd("monocle2/")2 构建CellDataSet对象、筛选、降维、排序
代码语言：r复制# 1.构建cds对象 ---------------------------------------------------------------
####################### 定义函数
run_monocle2 <- function(scRNAsub) {
scRNAsub <- sce
# 将counts矩阵转换为sparseMatrix
data <- as(as.matrix(scRNAsub@assays$RNA@counts), 'sparseMatrix')
# 将meta数据转换为AnnotatedDataFrame
pd <- new('AnnotatedDataFrame', data = scRNAsub@meta.data)
# 创建基因名表并转换为AnnotatedDataFrame
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
# 创建CellDataSet对象
mycds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
expressionFamily = negbinomial.size())
# 估计size factors和dispersions
mycds <- estimateSizeFactors(mycds)
mycds <- estimateDispersions(mycds, cores = 16, relative_expr = TRUE)
disp_table <- dispersionTable(mycds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
mycds <- setOrderingFilter(mycds, unsup_clustering_genes$gene_id)
pData(mycds)$orig.ident
> CONTROL GEM
> 1809 2974
# diff_test_res <- differentialGeneTest(mycds,
# fullModelFormulaStr = " celltype ~ orig.ident",
# reducedModelFormulaStr = " ~ orig.ident",
# relative_expr=TRUE, cores=16)
# ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
# ordering_genes
# mycds <- setOrderingFilter(mycds, ordering_genes)
#
# 寻找差异表达基因
diff.wilcox <- FindAllMarkers(scRNAsub)
all.markers <- diff.wilcox %>% dplyr::select(gene, everything()) %>% subset(p_val < 0.05)
diff.genes <- subset(all.markers, p_val_adj < 0.01)$gene
#设置代表性基因并绘图
mycds <- setOrderingFilter(mycds, diff.genes)
#降维
mycds <- reduceDimension(mycds , max_components = 4,
residualModelFormulaStr = "~orig.ident", #去除样本影响
reduction_method = 'DDRTree')
#排序
mycds <- orderCells(mycds )
save(mycds,file = "mycds.monocle2.Rdata")
# 返回处理过的CellDataSet对象和差异基因绘图
list(mycds = mycds)
}
# 使用示例
result <- run_monocle2(sce)
mycds <- result$mycds 代码解释
代码语言：r复制data <- ...
pd <- ...
fData <- ...
fd <- ...在使用 Monocle2 进行单细胞 RNA 测序数据分析时，数据的格式需要符合特定的要求，以便能够利用 Monocle2 的功能。这几行代码将数据从 Seurat 对象的结构转换为适合 Monocle2 的结构。
- 在 Seurat 对象中，基因表达数据通常以稀疏矩阵 (
sparseMatrix) 的形式存储，但也可以是其他形式。为了确保数据在 Monocle2 中的高效处理，这一步先将表达矩阵转换为常规的矩阵形式 (as.matrix)，然后再将其转换为稀疏矩阵形式 (sparseMatrix)。 - Seurat 对象中的元数据（meta data）通常以
data.frame的形式存储，包含每个细胞的附加信息，如细胞类型、样本来源等。在 Monocle2 中，这些元数据需要以AnnotatedDataFrame的形式提供。
AnnotatedDataFrame 是 Bioconductor 数据包中的一种特殊数据结构，能够更好地管理数据和元数据的关联。Monocle2 使用 AnnotatedDataFrame来存储细胞的注释信息，以便在后续分析中访问这些元数据。
通过调用 new() 函数，代码创建了一个 AnnotatedDataFrame 类的实例 pd，并将 scRNAsub@meta.data中的数据作为该对象的内容。
- fData是一个包含基因名称的
data.frame，它存储了每个基因的基本信息。在 Monocle2 中，这些信息必须以 AnnotatedDataFrame 的形式提供，
# 估计size factors和dispersions
mycds <- estimateSizeFactors(mycds)
mycds <- estimateDispersions(mycds, cores = 16, relative_expr = TRUE)- estimateSizeFactors(mycds)用于估计每个细胞的大小因子 (size factors)，以校正不同细胞间的测序深度差异。大小因子用于标准化基因表达数据，使得不同细胞的数据可以进行直接比较。
- estimateDispersions(mycds, cores = 16, relative_expr = TRUE)这个函数用于估计每个基因的分散性 (dispersion)，即基因表达水平的变异程度。分散性是判断基因在细胞群体中差异表达的重要指标。cores = 16 指定使用16个计算核心来加速计算。relative_expr = TRUE 指定使用相对表达量进行计算，以减少样本间的技术噪声。
#数据过滤
disp_table <- dispersionTable(mycds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
mycds <- setOrderingFilter(mycds, unsup_clustering_genes$gene_id)- dispersionTable(mycds)生成一个分散性表格，包含了所有基因的平均表达量和分散性信息。这个表格用于筛选在细胞群体中具有高变异性的基因，这些基因通常是细胞间异质性的重要标志。
- subset(disp_table, mean_expression >= 0.1)从分散性表格中筛选出平均表达量大于或等于 0.1 的基因。筛选这些基因是为了避免在分析中使用表达量过低的基因，因为这些基因可能贡献的信息较少。
- setOrderingFilter(mycds, unsup_clustering_genes$gene_id):筛选基因，这些基因将用于后续的无监督聚类和拟时序分析。
diff.wilcox <- FindAllMarkers(scRNAsub)
all.markers <- diff.wilcox %>% dplyr::select(gene, everything()) %>% subset(p_val < 0.05)
diff.genes <- subset(all.markers, p_val_adj < 0.01)$gene- FindAllMarkers（）用于查找在不同细胞群体间显著差异表达的基因。这个函数可以基于不同的统计测试（如 Wilcoxon 排序检验、ROC 分析等）来进行差异基因分析。
默认情况下，FindAllMarkers使用 Wilcoxon 排序检验来计算每个基因在每个细胞群体中与其他群体相比的差异表达程度。它会返回一个数据框，其中每行对应一个基因，包含基因名称、细胞群体、p 值等信息。
- 使用
dplyr包中的select函数，从diff.wilcox数据框中选择特定的列。gene指定了需要提取的列，这里它确保基因名称列被包括在内。everything()则会选择数据框中的所有其他列。这里筛选出所有 p 值小于 0.05 的基因，意味着这些基因在至少一个细胞群体中表现出显著的差异表达。 - 从 all.markers 中筛选出 p 值调整后（p_val_adj）小于 0.01 的基因。p 值调整是为了校正多重检验问题，从而减少假阳性的可能性。p 值调整后的基因可以更可靠地反映真实的差异表达。提取筛选后的基因名称列。最终结果是一个向量，包含了所有显著性更高（p 值调整后小于 0.01）的差异表达基因的名称
#设置代表性基因并绘图
mycds <- setOrderingFilter(mycds, diff.genes)
#降维
mycds <- reduceDimension(mycds , max_components = 4,
residualModelFormulaStr = "~orig.ident", #去除样本影响
reduction_method = 'DDRTree')
#排序
mycds <- orderCells(mycds )
save(mycds,file = "mycds.monocle2.Rdata")
# 返回处理过的CellDataSet对象和差异基因绘图
list(mycds = mycds) 这段代码是在 Monocle2 中进行降维、细胞排序和保存结果的过程。它的目的是利用筛选出的差异表达基因，按照拟时序的方式对细胞进行排序，从而推断出细胞状态的动态变化。
- setOrderingFilter()：用于在 CellDataSet 对象中设置用于排序的基因集合。在 Monocle2 中，排序基因决定了细胞如何在拟时序空间中排列。通常，这些基因是差异表达显著的基因，因为它们可以很好地捕捉到细胞状态的变化。
- reduceDimension()：这是 Monocle2 中的一个核心函数，用于对细胞进行降维处理。在单细胞 RNA 测序数据中，通常存在数千个基因，而降维的目的是将高维数据压缩到一个低维空间，以便可视化和分析。max_components = 4指定降维后保留的最大成分数量（或维度）。设置为 4，表示降维后的数据将保留 4 个主要成分，这些成分捕捉了数据中的主要变异。
- residualModelFormulaStr = "~orig.ident"：用于指定在降维时要校正的因素。orig.ident 通常代表样本或实验批次信息。通过在降维时去除样本效应，确保分析结果主要反映细胞间的生物学差异，而不是技术噪声。
- reduction_method = 'DDRTree'：
DDRTree（Discriminative Dimensionality Reduction with Trees）是 Monocle2 中的一种特殊降维方法，它不仅降低了数据的维度，还试图捕捉细胞状态之间的分支和拟时序轨迹。这种方法适用于拟时序分析，因为它可以揭示细胞状态如何随时间或分化进程变化。 - orderCells()：这个函数是 Monocle2 中拟时序分析的关键步骤之一。它基于前一步的降维结果，对细胞进行排序，确定它们在拟时序轨迹中的位置。排序的结果通常表示细胞状态的动态变化过程，比如细胞分化或其他生物学过程。
- 最后，将处理后的 CellDataSet 对象作为列表中的一个元素返回。这样做的目的是方便将该对象与其他可能的分析结果（如绘图对象）一起返回，便于后续操作或进一步分析。
3 定义根结点、重新排序
代码语言：r复制GM_state <- function(cds,x){
if (length(unique(pData(cds)$State)) > 1){
T0_counts <- table(pData(mycds)$State, pData(mycds)$RNA_snn_res.0.8)[,x]
return(as.numeric(names(T0_counts)[which
(T0_counts == max(T0_counts))]))
} else {
return (1)
}
}
# 使用示例
mycds <- orderCells(mycds, root_state = GM_state(mycds,"7"))代码解释
这段代码的目的是根据细胞状态信息确定拟时序分析中的“根状态”（root_state），并使用这个根状态来重新对细胞进行排序。这在拟时序分析中是关键的一步，因为根状态通常代表最早的或未分化的细胞状态，而其他细胞状态将在拟时序轨迹中相对于这个根状态进行排序。
- GM_state 函数:用于确定细胞的“根状态”（即最初状态或未分化状态）。在拟时序分析中，选择一个合适的根状态有助于正确地理解细胞分化或其他生物学过程的顺序。cds: 传入的 CellDataSet 对象（即 `mycds）; x: 指定某个簇的编号，比如 "7"，表示要考虑的细胞群体。
- 函数内部逻辑:
- if (length(unique(pData(cds)$State)) > 1): 这段代码检查 cds中的细胞状态数量。如果状态超过一个（即存在多个状态），则继续处理；否则，默认返回状态 1。
- T0_counts <- table(...),x:使用 table 函数生成一个二维表格，表格的行表示不同状态，列表示不同的簇（基于 RNA_snn_res.0.8 聚类结果）。这一步生成了一个关于某个簇（
x指定的簇，比如"7"）的细胞数目统计表。 - as.numeric(names(T0_counts)...),这一步找出在指定簇 x 中，哪个状态的细胞数目最多，并返回对应的状态编号。这个状态被认为是这个簇的“根状态”。
- return(...): 返回根状态的编号。
- 默认返回值,如果 cds只有一个状态，函数会返回 1作为默认根状态。
- orderCells(...):通过指定
root_state参数，可以确定细胞排序的起点，即轨迹的起始位置。调用 GM_state 函数，确定簇 "7"的根状态，然后将其作为 root_state 传递给 orderCells()函数。这意味着希望使用簇 "7"中细胞数目最多的状态作为拟时序轨迹的起点。
4 拟时序可视化（4张图）
代码语言：r复制######################作图
process_single_cell_data <- function(mycds) {
data_df <- t(reducedDimS(mycds)) %>%
as.data.frame() %>%
select_('Component 1' = 1, 'Component 2' = 2) %>%
rownames_to_column("Cells") %>%
mutate(pData(mycds)$State,
pData(mycds)$Pseudotime,
pData(mycds)$orig.ident,
pData(mycds)$celltype)
colnames(data_df) <- c("cells","Component_1","Component_2","State",
"Pseudotime","orig.ident","celltype")
reduced_dim_coords <- reducedDimK(mycds)
ca_space_df <- Matrix::t(reduced_dim_coords) %>%
as.data.frame() %>%
select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
mutate(sample_name = rownames(.), sample_state = rownames(.))
p1=ggplot()
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime),size =0.5)
scale_color_viridis()
theme_bw()
#theme_dr(arrow = grid::arrow(length = unit(0, "inches"))) #坐标轴主题修改
theme(
axis.line = element_line(color = "black", size = 0.5),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.title = element_text(size=15),
)
ggsave(p1,filename = "seurat_Pseudotime.pdf",width = 5,height = 3)
p3=ggplot()
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2, color = celltype),size =0.5)
scale_color_npg()
theme_bw()
#theme_dr(arrow = grid::arrow(length = unit(0, "inches"))) #坐标轴主题修改
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.line = element_line(color = "black", size = 0.5),
axis.title = element_text(size=15),
) facet_wrap(~celltype, nrow = 1)
w=length(unique(data_df$celltype))*3
ggsave(p3,filename = "seurat_celltype2.pdf",width = w,height = 3)
colnames(data_df)
p4=ggplot()
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2, color = celltype),size =0.5)
scale_color_npg()
theme_bw()
#theme_dr(arrow = grid::arrow(length = unit(0, "inches"))) #坐标轴主题修改
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.line = element_line(color = "black", size = 0.5),
axis.title = element_text(size=15),
) facet_wrap(~orig.ident, nrow = 1)
w=length(unique(data_df$orig.ident))*3
ggsave(p4,filename = "seurat_celltype3.pdf",width = w,height = 3)
# 返回两个数据框
list(data_df = data_df, ca_space_df = ca_space_df)
}
# 使用示例
result <- process_single_cell_data(mycds)
data_df <- result$data_df
ca_space_df <- result$ca_space_df代码解释
list(data_df = data_df, ca_space_df = ca_space_df)：
- 返回值：返回两个数据框，data_df包含降维后的细胞信息，ca_space_df 包含主成分图的降维坐标。
- 使用列表的好处是你可以一次性返回多个不同类型或用途的对象。在这个函数中，生成了两个不同的数据框（
data_df和ca_space_df），通过将它们放在列表中返回，可以方便地将它们一起传递给后续的处理步骤。
process_single_cell_data函数实现了对单细胞数据的伪时间分析和降维可视化，并生成了多种散点图，分别展示了伪时间、细胞类型和样本来源在低维空间中的分布。
函数的主要步骤包括：
- 提取降维后的数据并转换为数据框。
- 根据伪时间、细胞类型和样本类型生成散点图。
- 使用 ggplot2 保存生成的图像为 PDF 文件。
- 返回处理后的数据框，方便后续分析
p1 seurat_Pseudotime.pdf
伪时间散点图
p2 seurat_celltype.pdf
细胞类型散点图
p3 seurat_celltype2.pdf
分面显示不同细胞类型的散点图
p4 seurat_celltype3.pdf
分面显示不同样本来源的散点图
5 时序相关差异基因和富集分析
代码语言：r复制###################时序相关差异基因
clusterAndVisualize <- function(mycds, num_clusters = 4, gene_num = 200, enrich_db = "org.Mm.eg.db",
organism = "mmu", pvalue_cutoff = 0.05, topn = 5, seed = 123456,
output_pdf = "gh.pdf", pdf_height = 5, pdf_width = 12) {
# 差异基因测试
diff_test <- differentialGeneTest(mycds, cores = 16, fullModelFormulaStr = "~sm.ns(Pseudotime)")
diff_genes <- diff_test %>% arrange(qval) %>% head(gene_num) %>% dplyr::select(gene_short_name)
# 聚类并可视化
df <- plot_pseudotime_heatmap2(mycds[diff_genes[,1],],
num_clusters = num_clusters,
cores = 1,
use_gene_short_name = TRUE,
show_rownames = TRUE)
# 富集分析
enrich <- enrichCluster(object = df,
OrgDb = enrich_db,
type = "BP",
organism = organism,
pvalueCutoff = pvalue_cutoff,
topn = topn,
seed = seed)
# 采样标记基因
markGenes <- sample(unique(df$wide.res$gene), 25, replace = FALSE)
# 创建PDF输出
pdf(file = output_pdf, height = pdf_height, width = pdf_width, onefile = FALSE)
# 可视化聚类结果和富集分析结果
visCluster(object = df,
plot.type = "both",
column_names_rot = 45,
show_row_dend = FALSE,
markGenes = sample(df$wide.res$gene, gene_num, replace = FALSE),
markGenes.side = "right",
annoTerm.data = enrich,
go.col = rep(jjAnno::useMyCol("calm", n = num_clusters), each = 5),
line.side = "left")
# 关闭PDF设备
dev.off()
}
#使用示例
clusterAndVisualize(mycds, num_clusters = 4, gene_num = 400, enrich_db = "org.Mm.eg.db",organism = "mmu")代码解释
clusterAndVisualize 函数用于对单细胞测序数据中的拟时序分析结果进行聚类和可视化，并执行富集分析
函数参数
mycds：拟时序分析后的 CellDataSet 对象，包含了细胞的表达矩阵、元数据等信息。num_clusters：聚类的数量，用于将差异基因分成 num_clusters 个簇。gene_num：用于聚类和可视化的基因数量。enrich_db：用于富集分析的数据库，默认为"org.Mm.eg.db"，小鼠的基因数据库。organism：用于富集分析的物种，默认为"mmu"，即小鼠。pvalue_cutoff：富集分析的 P 值阈值。topn：在富集分析中显示的前topn个富集项。seed：用于随机过程中的种子，以保证结果的可重复性。output_pdf：输出 PDF 文件的名称。pdf_height和pdf_width：输出 PDF 的高度和宽度。
diff_test <- ...
diff_genes <- ...这两行代码主要涉及差异基因分析，用于识别那些在拟时序过程中显著变化的基因。
- differentialGeneTest对
mycds对象中的基因进行差异表达分析，目的是找出哪些基因在细胞的拟时序轨迹中发生显著变化。fullModelFormulaStr = "~sm.ns(Pseudotime)": 这是一个线性模型的公式，指定了如何解释基因表达的变化。返回一个包含测试结果的数据框，其中每一行代表一个基因，并包含该基因的统计检验结果，如 p 值（p-value）、校正后的 p 值（q-value）等。 - 对上一步生成的差异基因测试结果进行处理，选择前 gene_num 个最显著的差异表达基因。arrange(qval)：根据 qval（校正后的 p 值）对结果进行升序排序。
qval越小，基因的表达变化越显著，因此这种排序方式有助于找到那些在拟时序中变化最显著的基因。select(gene_short_name)：从结果中只选择gene_short_name列。最终结果 diff_genes 是一个包含 gene_num个基因符号的数据框，表示这些基因在拟时序分析中表达显著变化。这些基因随后可以用于进一步分析，如聚类分析或功能富集分析。
df <- plot_pseudotime_heatmap2(...)plot_pseudotime_heatmap2 函数生成一个热图(并不直接生成热图)，展示指定基因在拟时序轨迹中的表达情况，并根据这些基因的表达模式进行聚类。这种可视化有助于识别基因表达在细胞状态变化过程中的动态变化。
- mycds[diff_genes,1,] 提取的是这些基因的第一列（基因名称），即只选择这些基因的表达数据进行热图绘制。num_clusters = num_clusters:定义了要将基因分成的聚类数量。热图中的基因会根据它们的表达模式被分配到不同的聚类中，通常通过聚类分析（如 k-means 或 hierarchical clustering）实现。
df这个变量保存了 plot_pseudotime_heatmap2函数的输出。根据命名规则和上下文，这个输出可能是一个数据框，包含了基因表达和聚类的信息，也可能是热图的相关数据结构，用于进一步分析或绘图。
代码语言：r复制enrich <- enrichCluster(...)参数解释:
object = df: 输入的数据框，通常包含基因表达数据。OrgDb = enrich_db: 指定用于富集分析的基因注释数据库（如org.Mm.eg.db用于小鼠）。type = "BP": 指定富集分析的类型，这里是生物过程（BP）。organism = organism: 物种信息，例如小鼠（mmu）。pvalueCutoff = pvalue_cutoff: 指定显著性水平的阈值，筛选显著的富集结果。topn = topn: 返回前topn个显著的富集结果。seed = seed: 设置随机种子，以确保结果的可重复性。
在 enrichCluster 函数中设置 seed 是为了确保分析结果的可重复性。具体来说，以下是为什么在这样的函数中设置 seed 的原因：
- 在许多生物信息学分析中，特别是涉及富集分析或聚类分析的步骤，可能会包含随机性。例如，在富集分析中，可能会进行一些随机抽样、分配或算法初始化。这些步骤的输出可能会受到随机性的影响。
- 设置
seed可以使得这些随机过程在每次运行时产生相同的结果，从而确保分析的稳定性和可重复性。
markGenes <- sample(...)sample函数: 从 df$wide.res$gene 中不重复地随机抽取 25 个基因。unique(df$wide.res$gene): 获取所有独特的基因名。replace = FALSE`: 表示采样时不放回。
pdf(...)
visCluster(...)
dev.off()参数解释:
object = df: 输入的数据对象，通常是聚类或热图的数据。plot.type = "both": 指定绘制类型，如同时绘制热图和柱状图。column_names_rot = 45: 设置列名旋转角度。show_row_dend = FALSE: 指定是否显示行聚类树。markGenes: 指定要标记的基因。markGenes.side = "right": 标记基因显示在图形的右侧。annoTerm.data = enrich: 富集分析结果数据，用于在热图上进行注释。go.col: 设置不同聚类的颜色。line.side = "left": 将线条显示在图形的左侧。
gh.pdf
图片解读
这张图展示了基因表达模式在不同细胞簇中的变化，并将这些模式与特定的生物学过程（通过富集分析得到的结果）关联起来。以下是对这张图的详细解读：
1. 热图概述
- 列：每一列代表一个细胞，按照伪时间（pseudotime）的顺序排列。伪时间是一种衡量细胞在生物学过程（例如分化）中的进展程度的指标。伪时间的颜色梯度显示在右侧（从蓝色到红色）。
- 行：每一行代表一个基因，颜色表示基因的表达水平（蓝色表示低表达，红色表示高表达）。基因根据它们的表达模式被分成不同的簇。
2. 基因簇（C1, C2, C3, C4）
- C1, C2, C3, C4：这些是不同的基因簇（从C1到C4），簇内的基因在细胞中的表达模式相似。
- 基因数量：每个簇中的基因数量标注在簇旁边（例如，C1中有67个基因）。
- 表达趋势：每个簇旁边的小折线图显示了该簇基因在伪时间上的平均表达趋势。例如，C1中的基因表达逐渐增加，而C2中的基因表达逐渐减少。
3. 基因名称
- 热图右侧列出了每个簇中代表性的基因名称。这些基因是该簇中特征最明显的基因。
4. 富集的生物学过程
- 富集的GO术语：在每个基因簇（C1, C2, C3, C4）旁边的方框中列出了与该簇相关的Gene Ontology（GO）术语或生物学过程。这些过程是该簇基因可能参与的生物学过程。
- 例如
- C1 簇中富集了“细胞-基质粘附”和“胶质细胞生成”等过程。
- C2 簇与“葡萄糖分解过程”和“NADH再生”等过程相关。
- “n”旁边的数字表示在该过程中起作用的基因数量。
5. 颜色刻度
- Z-score刻度（左侧）：表示每个基因在细胞中的标准化表达水平（Z-score）。蓝色表示低于平均表达，红色表示高于平均表达。
- 伪时间刻度（右侧）：表示细胞在伪时间上的进展，从0（早期阶段，蓝色）到30（晚期阶段，红色）。
解读
- 这张图展示了不同基因在细胞进展过程中（伪时间）如何改变其表达模式。通过对基因进行聚类，我们可以发现哪些基因在相似的时间点共同上调或下调。此外，通过富集分析，我们可以将这些基因与特定的生物学过程关联起来，从而推测这些基因簇可能在细胞发育、分化或其他生物学功能中发挥的作用。
6 基因表达随时间变化的图
代码语言：r复制plotGeneOverPseudotime <- function(data, gene, output_file1, output_file2, colors_orig_ident) {
counts = Biobase::exprs(mycds)
phe <- pData(mycds)
dim(counts)
#z=as.data.frame(counts[gene_to_cluster,])
z=as.data.frame(counts)
ac=phe[,c("orig.ident","celltype","Pseudotime")]#这里是对热图进行分组的信息，从colnames(phe)中选
head(ac)
A=merge(ac,t(z),by = "row.names")
rownames(A)=A$Row.names
A=A[,-1]
# 绘制第一个图
p1 <- ggplot(data = A, aes(x = Pseudotime, y = .data[[gene]]))
geom_point(aes(color = celltype), alpha = 0.7, size = 0.5)
scale_color_npg() # 替换为实际的celltype值和颜色
new_scale_color()
geom_smooth(method = "loess", se = FALSE, aes(color = orig.ident))
scale_color_manual(values = colors_orig_ident, aesthetics = "colour", name = "orig.ident")
theme_bw()
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.line = element_line(color = "black", size = 0.5),
axis.title = element_text(size = 15)
)
ggsave(p1, filename = output_file1, width = 5, height = 3)
# 绘制第二个图，带facet_wrap
p2 <- ggplot(data = A, aes(x = Pseudotime, y = .data[[gene]]))
geom_point(aes(color = celltype), alpha = 0.7, size = 0.5)
scale_color_npg() # 替换为实际的celltype值和颜色
new_scale_color()
geom_smooth(method = "loess", se = FALSE, aes(color = orig.ident))
scale_color_manual(values = colors_orig_ident, aesthetics = "colour", name = "orig.ident")
facet_wrap(~orig.ident)
theme_bw()
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.line = element_line(color = "black", size = 0.5),
axis.title = element_text(size = 15)
)
# 计算宽度
w <- length(unique(data$orig.ident)) * 3
ggsave(p2, filename = output_file2, width = w, height = 3)
}
# 示例用法
plotGeneOverPseudotime(mycds, "S100a8", "gene_Pseudotime.pdf", "gene_Pseudotime_orig.ident.pdf",
colors_orig_ident = c("CONTROL" = "skyblue", "GEM" = "pink"))定义了一个名为 plotGeneOverPseudotime 的函数，用于绘制特定基因在伪时间（Pseudotime）上的表达趋势，并将其结果保存为PDF文件。这个函数包括了两种可视化方式：一个简单的点图和一个带有 facet_wrap 分组的图。以下是对这段代码的详细解释：
函数参数：
data：输入的数据集，一般是一个包含细胞信息和基因表达数据的对象。gene：指定要绘制的基因名称。output_file1和output_file2：两个PDF输出文件的路径，用于保存生成的图。colors_orig_ident：一个命名向量，用于指定不同orig.ident组别的颜色。
代码解释
代码语言：r复制counts = Biobase::exprs(mycds)
phe <- pData(mycds)counts：提取包含基因表达数据的矩阵，行表示基因，列表示细胞。phe：提取细胞的表型数据（metadata），包括 orig.ident、celltype和 Pseudotime 等信息。
z=as.data.frame(counts)
ac=phe[,c("orig.ident","celltype","Pseudotime")]
A=merge(ac,t(z),by = "row.names")
rownames(A)=A$Row.names
A=A[,-1]ac：提取与伪时间分析相关的表型信息，包括orig.ident（实验组别）、celltype（细胞类型）和Pseudotime（伪时间）。A：将表型信息与转置后的表达数据合并，生成一个综合数据框A。
gene_Pseudotime.pdf
gene_Pseudotime_orig.ident.pdf
7 差异分析分支节点
参考链接
https://cole-trapnell-lab.github.io/monocle-release/docs/#analyzing-branches-in-single-cell-trajectories
https://www.jianshu.com/p/9995cd707002
代码语言：r复制#########################差异分析分支节点
# 定义函数
plot_cell_trajectory(mycds, color_by = "State",size=0.01)
ggsave(filename = "plot_cell_trajectory.pdf",width = 5,height = 5)
analyze_and_visualize_BEAM <- function(mycds, branch_point_BEAM = 1, branch_point_heatmap = 1,
num_clusters = 4, cores = 8, qval_threshold = 1e-4,
enrich_db = "org.Mm.eg.db", organism = "mmu", pvalueCutoff = 0.05, topn = 5,
seed = 123456, output_filename = "branch-enrich.pdf") {
# 运行BEAM
BEAM_res <- BEAM(mycds, branch_point = branch_point_BEAM, cores = cores, progenitor_method = 'duplicate')
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
# 筛选基因并绘制热图
selected_genes <- row.names(subset(BEAM_res, qval < qval_threshold))
df <- plot_genes_branched_heatmap2(mycds[selected_genes,],
branch_point = branch_point_heatmap,
num_clusters = num_clusters,
cores = cores,
use_gene_short_name = TRUE,
show_rownames = TRUE)
# 富集分析
enrich <- enrichCluster(object = df,
OrgDb = enrich_db,
type = "BP",
organism = organism,
pvalueCutoff = pvalueCutoff,
topn = topn,
seed = seed)
# 随机选择标记基因
markGenes = sample(unique(df$wide.res$gene), 25, replace = FALSE)
# 可视化
pdf(output_filename, height = 6, width = 12, onefile = FALSE)
visCluster(object = df,
plot.type = "both",
column_names_rot = 45,
show_row_dend = FALSE,
markGenes = markGenes,
markGenes.side = "right",
annoTerm.data = enrich,
go.col = rep(jjAnno::useMyCol("calm", n = 4), each = 5),
go.size = "pval",
line.side = "left")
dev.off()
# 返回结果对象（如果需要）
list(heatmap_df = df, enrich_results = enrich, marked_genes = markGenes)
}
# 示例用法
result <- analyze_and_visualize_BEAM(mycds,branch_point_BEAM = 1,branch_point_heatmap = 1,
num_clusters = 4,
enrich_db = "org.Mm.eg.db", organism = "mmu")补充：指定基因在分支上的表达差异（只显示branch,不显示pre_branch）
代码语言：r复制test_genes <- row.names(subset(fData(mycds), gene_short_name %in% c("Fam92a", "Tecr", "Cdon")))
p_test <- plot_genes_branched_pseudotime(mycds[test_genes,],
branch_point = 1,
color_by = "State",
ncol = 1)
ggsave(p_test, filename = "test_gene_branch.pdf", width = 10, height = 8) test_gene_branch.pdf
对比branch-enrich.pdf理解
代码解释
代码语言：r复制plot_cell_trajectory（...）这两行代码先绘制细胞轨迹图，将轨迹按 "State" 着色
plot_cell_trajectory.pdf
analyze_and_visualize_BEAM <- function(...)定义了一个名为 analyze_and_visualize_BEAM 的函数，用于在单细胞轨迹分析中进行分支点的差异基因分析和富集分析，并生成相应的可视化结果。
参数解释:
mycds: 输入的单细胞 RNA-seq 数据对象，包含基因表达数据和细胞元数据。branch_point_BEAM = 1: 指定用于 BEAM 分析的分支点编号，确定在哪个细胞分支点进行差异分析,默认值为1。branch_point_heatmap = 1：指定用于绘制基因表达热图的分支点编号，用于展示该分支点上基因的表达模式。num_clusters = 4: 指定在热图中进行基因聚类的数量。基因根据表达模式被分为指定数量的聚类，并在热图中进行可视化。cores = 8: 指定用于并行计算的 CPU 核心数量，以加速 BEAM 分析和聚类过程。qval_threshold = 1e-4: 指定 q 值的阈值，用于筛选显著性基因。只有 q 值小于该阈值的基因才会被用于热图绘制。enrich_db = "org.Mm.eg.db": 指定用于富集分析的基因注释数据库。此处的数据库适用于小鼠 (Mus musculus) 的基因。如果分析的是其他物种的数据，需要使用对应物种的数据库。organism = "mmu": 指定物种代号，"mmu"表示小鼠。如果分析其他物种的数据，需要替换为相应的物种代号。pvalueCutoff = 0.05: 设置富集分析中显著性的 p 值阈值。只有 p 值小于该阈值的条目才会被认为是显著富集的。topn = 5: 指定在富集分析结果中展示的每个基因集群的前n个最显著的功能类别数量。seed = 123456：设置随机数生成的种子值，以确保分析的可重复性，特别是在涉及随机过程的步骤中（例如标记基因的随机选择）。output_filename = "branch-enrich.pdf": 指定输出 PDF 文件的文件名，该文件将保存热图和富集分析结果的可视化内容。
BEAM_res <- BEAM(...)
BEAM_res <- BEAM_res[order(BEAM_res$qval),] 注：什么是 BEAM 分析？
BEAM（Branch Expression Analysis Modeling） 是 Monocle 软件中的一种分析方法，专门用于处理单细胞 RNA-seq 数据。它用于识别在不同细胞命运分支之间差异表达的基因。简单来说，BEAM 分析可以帮助我们理解在细胞发育过程中，哪些基因在某个分支点上发生了显著变化，从而导致细胞向不同的方向分化。
- 参数说明 ：
1 branch_point = branch_point_BEAM: 指定进行 BEAM 分析的分支点（branch point），这个参数决定了在哪个细胞命运分支点上进行基因差异表达分析。
2 progenitor_method = 'duplicate': 指定使用哪种方法来处理 progenitor 细胞，'duplicate' 是一种方法，用于复制 progenitor 细胞的信息。
- order(BEAM_res$qval): 对 BEAM 分析的结果按 q 值（假阳性发现率）进行排序。q 值越低，表示基因在分支点上的差异表达越显著。通过排序，可以优先分析显著性较高的基因。
BEAM_res返回的是一个数据框，其中包含了每个基因在指定分支点上的差异表达统计结果。
selected_genes <- ...
df <- plot_genes_branched_heatmap2(...)- 筛选 q 值小于
qval_threshold的基因，然后使用plot_genes_branched_heatmap2绘制这些基因的分支热图。热图展示了不同分支上基因的表达情况，并按指定的集群数量进行聚类。
enrich <- enrichCluster(...)
markGenes = sample(...)- 对聚类结果进行 GO (Gene Ontology) 生物过程 (BP) 的富集分析，找到在每个集群中富集的生物学过程。seed 参数用于保证随机选择标记基因的一致性。
- 从聚类结果中随机选择 25 个基因，作为标记基因，后续将在热图中标注这些基因
pdf(...)
visCluster(...)
dev.off()
list(heatmap_df = df, enrich_results = enrich, marked_genes = markGenes)- 可视化、返回包含热图数据框、富集分析结果和标记基因的列表
branch-enrich.pdf
图片解读
曲线图，具体反映了在不同分支（branch1 和 branch2）上，各个基因聚类（C1, C2, C3, C4）的表达变化模式。以下是对这些曲线图的详细解释：
基因聚类（C1, C2, C3, C4）
- 每个基因聚类（如C1、C2、C3、C4）代表在分支点分析中被识别为具有相似表达模式的一组基因。也就是说，在分支点（branch point）附近，这些基因的表达模式具有共同的趋势。
曲线图
- 横轴：表示从
branch1到branch2的分化轨迹。左端代表branch1，右端代表branch2。 - 纵轴：表示基因的表达量。曲线的高度反映了基因在不同分支中的相对表达水平。
曲线的含义
- C1（曲线向下，然后再向上）：
- 在
branch1中，这些基因的表达水平较高（曲线起始较高）。 - 在过渡区域（
Pre-branch），表达水平有所下降（曲线下降）。 - 到达
branch2时，表达水平再次上升（曲线回升）。 - Gene Size: 1042：表示该聚类包含了1042个基因。
- 在
- C2（曲线变化较小）：
- 这些基因在
branch1和branch2中的表达水平相对稳定，但在过渡区域略有波动。 - Gene Size: 171：表示该聚类包含了171个基因。
- 这些基因在
- C3（曲线在
branch1中波动较大，在branch2中趋于稳定）：- 在
branch1中，这些基因的表达水平有波动，但整体较高。 - 到达
branch2时，表达水平逐渐趋于稳定。 - Gene Size: 290：表示该聚类包含了290个基因。
- 在
- C4（曲线先下后上，变化较大）：
- 在
branch1中，基因表达水平先下降。 - 到达过渡区域后，表达水平再次升高，在
branch2中达到峰值。 - Gene Size: 711：表示该聚类包含了711个基因。
- 在
总结
左侧的这些曲线图直观展示了各个基因聚类在两个分支间的表达变化模式。这些变化趋势反映了不同基因在细胞分化过程中如何响应分支选择，帮助理解特定基因在不同分化路径中的角色和功能。
热图，展示了单细胞分化轨迹上的基因表达模式，重点放在两个分支 branch1 和 branch2 之间的差异。
热图部分（中间部分）
- X轴：表示分支
branch1和branch2，它们代表了不同的分化轨迹或细胞命运。热图的列代表单个细胞，按分支顺序排列。 - Y轴：表示被分析的基因，这些基因在不同分支上的表达水平有所不同。基因分为不同的聚类（C1、C2、C3、C4），即在不同分支间具有相似表达模式的基因。
- 颜色：表示基因的Z-score（标准分数），即相对于基因平均表达量的表达水平。红色表示基因在特定分支中的高表达，蓝色表示低表达。
右侧部分
- 功能富集分析：显示了每个基因聚类所富集的生物学过程（BP）。
颜色图例
- Z-score 颜色图例：右侧的颜色条表示Z-score的范围。红色代表高表达，蓝色代表低表达。
- branch：用颜色（红色、灰色和蓝色）标识了细胞在不同分支（
Cell fate 1，Pre-branch，Cell fate 2）中的归属。
特别说明
借用官网对类似图片的描述
“列是伪时间中的点，行是基因，伪时间的开始位于热图的中间。当您从热图的中间向右阅读时，您正在通过伪时间跟踪一个谱系。当你向左阅读时，另一个。这些基因是分层聚类的，因此您可以可视化具有相似谱系依赖性表达模式的基因模块。”
8 不同分支上指定基因绘图
代码语言：r复制plot_gene_branched_pseudotime <- function(mycds, genes, gene_name ,branch_point = 1, color_by = "celltype",
ncol = 1, cell_size = 2, width = 5, height = 3,
output_filename = "gene_branched_Pseudotime.pdf",
colors_celltype = NULL, colors_branch = NULL) {
# 绘制基因表达伪时间图
kk <- plot_genes_branched_pseudotime(mycds[genes,],
branch_point = branch_point,
color_by = color_by,
ncol = ncol,
cell_size = cell_size)
# 将ggplot对象转换为数据框
kk <- ggplot_build(kk)
B <- kk$plot$data
# 选择需要的列
B2 <- B[, c("Cell", "expression", "orig.ident", "f_id", "Pseudotime", "celltype", "Branch")]
# 转换数据框为宽格式
df <- dcast(B2, Cell orig.ident Pseudotime celltype Branch ~ f_id, value.var = "expression")
# 检查genes中是否有多个基因，并只选择第一个（如果有多个，需要修改以适应多个基因）
gene_name <- gene_name # 假设我们只处理一个基因，如果有多个需要调整
# 绘制ggplot
p1 <- ggplot(data = df, aes(x = Pseudotime, y = get(gene_name))) # 使用get()来处理字符串作为列名
geom_point(aes(color = celltype), alpha = 0.7, size = 0.5)
scale_color_npg() # 替换为实际的celltype值和颜色
new_scale_color()
geom_smooth(method = "loess", se = FALSE, aes(color = Branch))
scale_color_manual(values = c("skyblue","pink"), name = "Branch", aesthetics = "colour")
theme_bw() labs(x="Pseudotime",y =gene_name ,title = "")
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.1, "lines"),
axis.ticks = element_blank(),
axis.line = element_line(color = "black", size = 0.5),
axis.title = element_text(size = 15))
# 保存图片
ggsave(p1, filename = output_filename, width = width, height = height)
# 可选：返回ggplot对象以供进一步使用
return(p1)
}
# 示例用法
plot_gene_branched_pseudotime(mycds, genes = c("S100a10", "S100a8"),gene_name= "S100a10" ,branch_point = 3)plot_genes_branched_pseudotimes是内置函数，plot_gene_branched_pseudotime是重新定义的函数。
代码解释
kk <- ggplot_build(kk): 这行代码将 ggplot 对象 kk 转换为一个列表，其中包含 ggplot内部的数据和布局信息。
B <- kk$plot$data： 从 ggplot 对象中提取实际用于绘图的数据框，并将其存储在 B 中。
B2 <- B, c("Cell", "expression", "orig.ident", "f_id", "Pseudotime", "celltype", "Branch"): 从数据框 B 中选择特定的列，并将这些列存储到一个新的数据框 B2 中。
df <- dcast(B2, Cell orig.ident Pseudotime celltype Branch ~ f_id, value.var = "expression")“：将数据框 B2 转换为宽格式，并将其存储在 df 中。
注：
问：我只关注gene_name= "S100a10",为什么还要设置genes = c("S100a10", "S100a8")？
答：即使只关注 S100a10 这个基因。但因为 plot_genes_branched_pseudotime() 函数需要输入一组基因（即使只打算绘制其中一个) ，即plot_genes_branched_pseudotime(mycdsgenes,, ...)
gene_branched_Pseudotime.pdf
注：图中只有branch上的信息，而没有pre_branch上的信息。
补充：总体展示降维聚类分群
代码语言：r复制DimPlot(sce, reduction = "umap", group.by = "celltype",
split.by = 'orig.ident',
label = T,pt.size = 0.1,label.size = 3,
repel = T,label.box = T)
scale_colour_manual(values = pal_d3("category20")(20),
aesthetics = c("colour", "fill"))
ggsave('umap-by-celltype-ggsci.pdf',height = 4,width=8)umap-by-celltype-ggsci.pdf
