Single Cell——轨迹分析（Monocle3）
·
monocle3
Monocle 3包提供了分析单细胞基因表达实验的工具包。
Monocle3和Monocle2并没有本质上的区别，只是把降维图从DDRTree改成了UMAP。原因可能是包的作者认为UMAP比DDRTree降维更能反映高维空间的数据。
monocle3的三个主要功能：
- 分群、计数细胞
- 构建细胞轨迹
- 差异表达分析
细胞聚类及鉴定亚群
安装包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 一些依赖包
bioc_pkgs <- c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
'limma', 'S4Vectors', 'SingleCellExperiment',
'SummarizedExperiment')
for (bpkg in bioc_pkgs){
if (! require(bpkg,character.only=T) ) {
BiocManager::install(bpkg,ask = F,update = F)
require(bpkg,character.only=T)
}
}
# (可选)如果要在细胞聚类时设定分辨率参数(resolution)，就要安装一个python包
install.packages("reticulate")
reticulate::py_install("louvain")
# 现在安装3版本，还是要通过github
devtools::install_github('cole-trapnell-lab/monocle3')
# 重启Rstudio，检测是否成功
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
library(ggplot2)
library(dplyr)
加载数据
expression_matrix是一个行为基因，列为细胞样本的表达矩阵cell_metadata是一个数据框，行为细胞，列是细胞的属性（比如细胞类型、培养环境、培养时间等）gene_metadata是一个数据框，行为feature信息（比如基因），列是基因属性（比如GC含量）
必须要满足：
- expression_matrix的行与gene_metadata的行对应
- expression_matrix的列与cell_metadata的行对应
# 下载数据(如果已经有了，可以跳过)
expression_matrix <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("http://staff.washington.edu/hpliner/data/cao_l2_rowData.rds"))
# 加载、探索数据
expression_matrix <- readRDS(file = 'cao_l2_expression.rds')
cell_metadata <- readRDS(file = 'cao_l2_colData.rds')
gene_annotation <- readRDS(file = 'cao_l2_rowData.rds')
dim(expression_matrix)
expression_matrix[1:3,1:3]
cell_metadata[1:3,1:3]
head(gene_annotation)
将数据存储在cell_data_set对象中
# 创建CDS(Cell Data Set)对象
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
cds
数据预处理
# 对数据进行规范化和预处理
# 对数据进行规范化和预处理
cds <- preprocess_cds(cds, num_dim = 100) #大约2 mins
plot_pc_variance_explained(cds)
我们可以看到，使用超过100台PC只会捕获少量额外的变化，并且每增加一台PC都会使Monocle的下游步骤变慢。
降低维度并可视化细胞
可以使用t-SNE，这在单细胞rna测序中非常流行，或者使用UMAP。Monocle 3默认使用UMAP，它更适合于RNA-seq中的聚类和轨迹分析。调用reduce_dimension()将数据的维数降低到X， Y平面，以便我们可以轻松地绘制它，请：
# umap降维
cds <- reduce_dimension(cds, preprocess_method = "PCA",reduction_method = c("UMAP"))
plot_cells(cds)
图中的每个点都代表cds对象中的不同细胞，可以看到不同的细胞分成了不同的群，有的群有几千个细胞，有的群只有少数几个，我们这里直接使用测试数据中做好的细胞标记
# 数据中一共提供了30类细胞(细胞类型及数量见：)
table(colData(cds)$"cao_cell_type")
# 根据这个上色
plot_cells(cds, color_cells_by="cao_cell_type")
# (另外除了用cao_cell_type细胞类型，还能用下面的任意一个)
colnames(colData(cds))
可以看到图中许多细胞类型再UMAP结果中靠的很近
根据基因表达量进行映射：
#根据基因表达量进行映射
plot_cells(cds, genes=c("cpna-2", "egl-21", "ram-2", "inos-1"))
使用tSNE方法降维的话：
#进行tSNE降维
cds <- reduce_dimension(cds, reduction_method="tSNE")
# 然后对tSNE结果可视化
plot_cells(cds, reduction_method="tSNE", color_cells_by="cao_cell_type")
其实可以看到，tSNE结果的各个细胞亚群内部的关联性不如UMAP
检查、移除批次效应 => residual_model_formula_str参数
为了更加真实地反映差异基因的来源是生物学因素，而不是技术因素（比如细胞板批次、细胞培养差别、测序差别等），就首先要检查有没有所谓的”技术批次“对数据差异产生影响。
当进行降维分析时，应该经常检查批次效应，最好在colData(cds) 添加一列，其中注明细胞的批次信息（比如来自哪个细胞板），然后就能根据这个对细胞的批次进行可视化，结果最好是各个细胞亚群的不同批次混在一起，这样我们就能认为：降维所采用的”特异性“的特征并不来自批次
#对细胞的批次进行可视化
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)
如果真的担心不同批次会产生其他的影响，可以在降维之前对批次进行校正，保证后面的各个细胞群只来自一个批次，那么它们之间的分群差异就更可能是由于生物因素导致
# 假入要对批次进行校正
cds = preprocess_cds(cds, num_dim = 100, residual_model_formula_str = "~ plate")
# 校正后再降维
cds = reduce_dimension(cds)
# 再次检查批次效应
plot_cells(cds, color_cells_by="plate", label_cell_groups=FALSE)
细胞聚类
这是细胞分群过程中非常重要的一步。Monocle利用Louvain community detection 这个非监督聚类方法（它是phenoGraph算法的一部分），它和我们常见的Kmeans、hclust层次聚类等还不太一样，这个方法更倾向于找到一个网络中联系紧密的部分，而经常忽略节点的特性；而我们常见的聚类呢，它一般会忽视整体中各个部分的联系，而通过计算两个节点(目标)之间的距离（如欧式距离、曼哈顿距离、余弦相似度等） 找到相似的特性
按照cluster可视化
cds = cluster_cells(cds, resolution=c(10^seq(-6,-1)))
plot_cells(cds)
# 这个可视化函数如果不加任何参数，它默认对不同的cluster上色
# 用下面函数查看有哪些cluster：
cds@clusters$UMAP$clusters #或者
clusters(cds)
按照partition进行可视化
不过这样细分看的太混杂，于是 Alex Wolf et al 利用他们开发的 PAGA 算法将一些小的cluster聚合成一个大的partition
plot_cells(cds, color_cells_by="partition", group_cells_by="partition")
# 检查有多少partition
cds@clusters$UMAP$partitions
按照每个cluster对应的细胞类型
除了默认按照cluster编号进行标记不同的亚群，还可以按照每个cluster对应的细胞类型进行可视化：
#按照每个cluster对应的细胞类型进行可视化
plot_cells(cds, color_cells_by="cao_cell_type")
对label进行简化：
plot_cells(cds, color_cells_by="cao_cell_type", label_groups_by_cluster=FALSE)
找marker基因
确定了各种亚群以后，我们就想知道什么导致了它们的差异，这就是寻找marker基因的过程，利用top_markers()函数
#找marker基因
marker_test_res = top_markers(cds, group_cells_by="partition", reference_cells=1000, cores=8)
# 设置reference_cells是随机挑出来这些数量的细胞作为参照，然后让top_markers和参照集中的基因进行显著性检验；另外reference_cells还可以是来自colnames(cds)的细胞名
marker_test_res[1:4,1:4]
自行过滤：
top_specific_markers = marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>%
top_n(1, pseudo_R2)
# 基因id去重
top_specific_marker_ids = unique(top_specific_markers %>% pull(gene_id))
对每组的marker基因可视化
plot_genes_by_group(cds,
top_specific_marker_ids,
group_cells_by="partition",
ordering_type="maximal_on_diag",
max.size=3)
如果要多看几个marker，只需要修改一下top_n
top_specific_markers = marker_test_res %>%
filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>%
top_n(2, pseudo_R2)
top_specific_marker_ids = unique(top_specific_markers %>% pull(gene_id))
plot_genes_by_group(cds,
top_specific_marker_ids,
group_cells_by="partition",
ordering_type="cluster_row_col",
max.size=3)
对细胞类型进行注释
上面的测试数据中给出了细胞类型，所以可以进行聚类后的可视化，但实际上，我们通常只能看到不同的cluster数字和它们的分布（也就是类似于as.character(partitions(cds))的结果），因此利用marker基因去重新定义细胞类型是非常关键的一步。
# 先将partitions的分组由因子型转为字符型
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"="Unclassified neurons",
"4"="Seam cells",
"5"="Coelomocytes",
"6"="Pharyngeal epithelia",
"7"="Vulval precursors",
"8"="Non-seam hypodermis",
"9"="Intestinal/rectal muscle",
"10"="Touch receptor neurons",
"11"="Pharyngeal neurons",
"12"="Am/PH sheath cells",
"13"="NA",
"14"="Unclassified neurons",
"15"="flp-1(+) interneurons",
"16"="Canal associated neurons",
"17"="Pharyngeal gland",
"18"="Other interneurons",
"19"="Ciliated sensory neurons",
"20"="Ciliated sensory neurons",
"21"="Ciliated sensory neurons",
"22"="Ciliated sensory neurons",
"23"="Ciliated sensory neurons",
"24"="Ciliated sensory neurons",
"25"="Oxygen sensory neurons",
"26"="Ciliated sensory neurons",
"27"="Unclassified neurons",
"28"="Pharyngeal gland",
"29"="Ciliated sensory neurons",
"30"="Ciliated sensory neurons",
"31"="Ciliated sensory neurons",
"32"="Ciliated sensory neurons",
"33"="Pharyngeal muscle",
"34"="Failed QC")
plot_cells(cds, group_cells_by="partition", color_cells_by="assigned_cell_type")
# Seurat中用RenameIdents 进行细胞类型重定义
# 另外，想从中取子集、过滤的话
cds[,colData(cds)$assigned_cell_type != "Failed QC"]
基于Garnett的自动化注释
预处理
# step-1：首先根据上面得到的细胞类型assigned_cell_type，找top_marker
assigned_type_marker_test_res = top_markers(cds,
group_cells_by="assigned_cell_type",
reference_cells=1000,
cores=8)
# step-2：过滤（阈值自定义）
garnett_markers = assigned_type_marker_test_res %>%
filter(marker_test_q_value < 0.01 & specificity >= 0.5) %>%
group_by(cell_group) %>%
top_n(5, marker_score)
# step-3：去重复
garnett_markers = garnett_markers %>% group_by(gene_short_name) %>%
filter(n() == 1)
# step-4：生成marker文件
generate_garnett_marker_file(garnett_markers, file="./marker_file.txt")
- 方法一：自己使用Garnett
# 安装Garnett
#BiocManager::install(c("org.Mm.eg.db", "org.Hs.eg.db"))
#devtools::install_github("cole-trapnell-lab/garnett", ref="monocle3")
library(garnett)
# 安装物种注释库(这里的物种是C. elegans秀丽隐杆线虫)，目的是为了做基因ID转换
#BiocManager::install("org.Ce.eg.db")
library(org.Ce.eg.db)
#colData(cds)$garnett_cluster = cds@clusters$UMAP$clusters
# install.packages("igraph")
#devtools::install_version("igraph", version = "2.0.1")
library(igraph)
#class(cds)
##########?train_cell_classifier
# rlang::last_trace()
colData(cds)$garnett_cluster <- clusters(cds)
worm_classifier <- train_cell_classifier(cds = cds,
marker_file = "./marker_file.txt",
db=org.Ce.eg.db,
cds_gene_id_type = "ENSEMBL",
num_unknown = 50,
marker_file_gene_id_type = "SYMBOL",
cores = 8)
# 这样自己就做了一个分组数据集，官方也建议将高质量的训练集汇总到：https://github.com/cole-trapnell-lab/garnett/issues
#
cds = classify_cells(cds, worm_classifier,
db = org.Ce.eg.db::org.Ce.eg.db,
cluster_extend = TRUE,
cds_gene_id_type = "ENSEMBL")
plot_cells(cds,
group_cells_by="partition",
color_cells_by="cluster_ext_type")
- 方法二：使用别人提供的分组数据集
#用提供的分组数据集
ceWhole <- readRDS(url("https://cole-trapnell-lab.github.io/garnett/classifiers/ceWhole_20191017.RDS"))
cds <- classify_cells(cds, ceWhole,
db = org.Ce.eg.db,
cluster_extend = TRUE,
cds_gene_id_type = "ENSEMBL")
plot_cells(cds,
group_cells_by="partition",
color_cells_by="cluster_ext_type")
构建发育轨迹
在发育生物学中，细胞整个发育阶段中对刺激做出的不同相应导致了细胞功能状态的转化。不同状态的细胞表达不同的基因，从而产生不同的蛋白、代谢物参与不同的生命过程。伴随着细胞状态的转化，转录信息也在不断更换，有的基因起初沉默后来激活。但是这个事情很难去探索，因为要纯化出不同状态并且稳定的细胞基本是不现实的，而scRNA可以不用通过实验方法的纯化细胞，就能研究细胞生命的各个状态。
这也是Monocle最大的亮点—利用算法去学习每个细胞基因表达量的变化，从而对细胞状态进行推断，推断出整体的一个表达轨迹，那么就可以将每个细胞放在特定位置上；另外有可能一个过程有多个结果，那么在轨迹图上就会有多个分支，代表了细胞分化
载入数据，构建对象
expression_matrix <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_expression.rds"))
cell_metadata <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_colData.rds"))
gene_annotation <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_rowData.rds"))
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
预处理
cds <- preprocess_cds(cds, num_dim = 50)
cds <- align_cds(cds, alignment_group = "batch", residual_model_formula_str = "~ bg.300.loading + bg.400.loading + bg.500.1.loading + bg.500.2.loading + bg.r17.loading + bg.b01.loading + bg.b02.loading")
降维可视化
cds <- reduce_dimension(cds)#降维
plot_cells(cds, label_groups_by_cluster=FALSE, color_cells_by = "cell.type")#可视化
看不同基因在不同分支的表达量:
ciliated_genes = c("che-1",
"hlh-17",
"nhr-6",
"dmd-6",
"ceh-36",
"ham-1")
plot_cells(cds,
genes=ciliated_genes,
label_cell_groups=FALSE,#不在可视化中标记细胞群体
show_trajectory_graph=FALSE)#不显示轨迹图
细胞聚类
Monocle不认为一组数据中的所有细胞都来自同一个”祖先“，很多实验中，它们会有多个发育轨迹。Monocle会通过聚类来判断细胞是否应该归属同一个发育轨迹。之前介绍的cluster_cells()中，细胞可以按cluster细分，还可以按partition归为大类。
#聚类
cds <- cluster_cells(cds)#分组聚类
plot_cells(cds, color_cells_by = "partition")#可视化，并按分区着色
利用learn graph在每个partition中寻找主路径
cds <- learn_graph(cds)#学习细胞之间的图形结构
#可视化
plot_cells(cds,
color_cells_by = "cell.type",
label_groups_by_cluster=FALSE,#不显示群组
label_leaves=FALSE,#不显示叶节点
label_branch_points=FALSE)#不显示分支点
对细胞出现的先后进行排序
Monocle将基因表达量的变化定义为发育轨迹上生命过程的变化（体现为不同的细胞类型），而不是真正的时间变化（因此我们叫它”拟时序“，而不是”真时间分析“）
之前得到的轨迹总长度就是细胞从一个阶段的起始到终止所产生的转录水平变化总量；拟时序分析会沿着最短的路径，计算每个细胞和轨迹起始位点的距离
# 先将每个细胞根据”胚胎发育时间区间“进行上色，然后根据胚胎发育过程前后绘制节点
# 在cds@colData中提供了细胞对应的类型以及胚胎发育时间，这个才是真正的关键信息。有了这个信息才能作图。工具不是万能的，它不可能帮助我们去完成生物学中的推断，一定是我们自己先定义好，再交给它进行可视化而已
plot_cells(cds,
color_cells_by = "embryo.time.bin",
label_cell_groups=FALSE,
label_leaves=TRUE,
label_branch_points=TRUE,
graph_label_size=1.5)
- 图中的黑线就是整个架构（注意到整个图并非完全连接的，毕竟是按照partition聚类，每个partition中的细胞差异有点大）
- 浅灰色的圆圈表示”叶片“表示发育轨迹中的不同结局（可以认为拟时序分析的图是一个”根-茎-叶“结构），用参数label_leaves控制
- 黑色的圆圈表示分支节点，预示着其中的细胞会有不同发展方向，用参数label_branch_points控制
- 圈中数字大小表示出现时间的先后
从上面的图中我们就能知道出现时间比较早的细胞位置，我们需要对全部的细胞都排个先后顺序，因此要用到order_cells() ，不过这个函数需要一个”根节点“位置，一般是一个partition分配一个根节点。
简单高效挑选根节点
首先根据轨迹图中的节点将与它们最相近的细胞分成组，然后计算来自最早时间点的细胞所占组分，最后挑出包含早期细胞数量最多的节点，认为它是就是根节点
cds <- order_cells(cds)
plot_cells(cds,
color_cells_by = "pseudotime",
#在伪时间（pseudotime）轨迹上的位置进行着色
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5)
# 官方给出了一个函数，这里定义了一个time_bin，选择了最早的时间点区间。
get_earliest_principal_node <- function(cds, time_bin="130-170")
#定义函数
#time_bin：时间区间，默认值为 "130-170"
{
cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)
#选择时间区间为 time_bin 的细胞
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
#每个细胞在UMAP图中最接近的顶点的信息
closest_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,]))))]
#从UMAP主要图形中通过 计算次数最多的值 提取顶点
root_pr_nodes
#返回找到的根节点
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
#使用找到的根节点对细胞数据集进行排序
cds_sub <- choose_graph_segments(cds)
plot_cells(cds,
color_cells_by = "pseudotime",
#在伪时间（pseudotime）轨迹上的位置进行着色
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5)
利用3D的发育轨迹对上述内容做个概述
3D轨迹实际上就是降维时选前3个主成分 => max_components = 3，后续都和2D保持类似
cds_3d <- reduce_dimension(cds, max_components = 3)#降维到3
cds_3d <- cluster_cells(cds_3d)#聚类
cds_3d <- learn_graph(cds_3d)#学习图形结构
cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))#指定根节点排序
cds_3d_plot_obj <- plot_cells_3d(cds_3d, color_cells_by="partition")
cds_3d_plot_obj
Monocle差异分析
Monocle3提供了不同细胞类型之间寻找差异基因的方法，主要有两种：
- Regression analysis：利用fit_models()，用来评价基因表达是否会受到诸如时间、处理等的影响
- Graph-autocorrelation analysis：利用graph_test()，用来寻找一条轨迹或不同cluster中基因的差异
方法一：Regression analysis
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会对每个基因拟合一个回归模型，其中会加入实验中的一些因素（比如时间、处理等）。例如，在胚胎相关的单细胞数据中，细胞会在不同的时间点进行收集，这里就可以检测是否有基因在这段时间内的表达量发生了变化，利用fit_models
gene_fits = fit_models(cds_subset, model_formula_str = "~embryo.time")
# 其中model_formula_str就是要比较的分组对象，如果相获得不同的cluster或者partition的差异基因，就用model_formula_str = "~cluster"或者model_formula_str = "~partition"；另外还支持添加多个变量，比如考虑到批次效应 model_formula_str = "~embryo.time + batch"
然后看看哪些基因是与时间因素相关的：
# coefficient_table()默认使用 Benjamini and Hochberg（BH）方法进行了p值的校正，得到了q值
fit_coefs = coefficient_table(gene_fits)
# 挑出时间相关的组分
emb_time_terms = fit_coefs %>% filter(term == "embryo.time")
#fit_coefs 中筛选出了 term 列中值为 "embryo.time" 的行
# coefficient_table()默认使用 Benjamini and Hochberg（BH）方法进行了p值的校正，得到了q值
#emb_time_terms %>% filter (q_value < 0.05) %>%
# select(gene_short_name, term, q_value, estimate)
画图展示：
plot_genes_violin(cds_subset, group_cells_by="embryo.time.bin", ncol=2) +
theme(axis.text.x=element_text(angle=45, hjust=1))
方法二：Graph-autocorrelation analysis
# 按照“集群和分类单元”一节中的描述重新加载和重新处理数据
expression_matrix <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("https://depts.washington.edu:/trapnell-lab/software/monocle3/celegans/data/cao_l2_rowData.rds"))
# 制作CDS对象
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-5)
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")
神经元的子集：
neurons_cds <- cds[,grepl("neurons", colData(cds)$assigned_cell_type, ignore.case=TRUE)]
plot_cells(neurons_cds, color_cells_by="partition")
从图中可以看到，存在不同类型的神经元细胞，它们也许来自不同亚群。为了探索这些不同类型的细胞之间是哪些基因差异表达，一种方法就是上面的线性回归。不过对于UMAP或tSNE的图，graph_test()这个函数可以进行一个叫做Moran’s I 的spatial autocorrelation分析方法
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))
将共同作用的基因合并成模块:
gene_module_df = find_gene_modules(neurons_cds[pr_deg_ids,], resolution=1e-2)
要看哪个module和哪个cluster/partition相关，有两种可视化方法：
- 第一种：
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)
- 第二种：针对少数量的module可以看得比较清楚，这里选了4个
plot_cells(neurons_cds,
genes=gene_module_df %>% filter(module %in% c(16,38,33,42)),
group_cells_by="partition",
color_cells_by="partition",
show_trajectory_graph=FALSE)
找到影响发育轨迹的基因
ciliated_cds_pr_test_res = graph_test(cds, neighbor_graph="principal_graph", cores=4)
# 使用neighbor_graph="principal_graph"来检验轨迹相邻的细胞的表达是否相关
pr_deg_ids = row.names(subset(ciliated_cds_pr_test_res, q_value < 0.05))
# 从pr_deg_ids中挑选几个基因，然后可视化
plot_cells(cds, genes=c("hlh-4", "gcy-8", "dac-1", "oig-8"),
show_trajectory_graph=FALSE,
label_cell_groups=FALSE,
label_leaves=FALSE)
# 可以将这些在轨迹上变化的pr_deg_ids继续分为小模块
gene_module_df = monocle3:::find_gene_modules(cds[pr_deg_ids,], resolution=c(0,10^seq(-6,-1)))
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")
plot_cells(cds,
genes=gene_module_df %>% filter(module %in% c(29,20, 11,22)),
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
另外一种方法：
先画完整的轨迹图：
plot_cells(cds, show_trajectory_graph=FALSE)
然后挑某一个细胞亚群对应的一段轨迹：
# 假设这里AFD细胞对应一段22、28、35组成的轨迹线
AFD_genes = c("gcy-8", "dac-1", "oig-8")
AFD_lineage_cds = cds[rowData(cds)$gene_short_name %in% AFD_genes,
clusters(cds) %in% c(22, 28, 35)]
plot_genes_in_pseudotime(AFD_lineage_cds,
color_cells_by="embryo.time.bin",
min_expr=0.5)
参考文章
开源鸿蒙跨平台开发社区汇聚开发者与厂商，共建“一次开发，多端部署”的开源生态，致力于降低跨端开发门槛，推动万物智联创新。
更多推荐
14
0- 0
已为社区贡献1条内容
所有评论(0)