跟着官网学习单细胞Monocle V3包（下）
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆继续上次的介绍 跟着官网学习单细胞Monocle V3包（上）
Monocle示例2-构建发育轨迹
意义： 在发育生物学中，细胞整个发育阶段中对刺激做出的不同相应导致了细胞功能状态的转化。不同状态的细胞表达不同的基因，从而产生不同的蛋白、代谢物参与不同的生命过程。伴随着细胞状态的转化，转录信息也在不断更换，有的基因起初沉默后来激活。但是这个事情很难去探索，因为要纯化出不同状态并且稳定的细胞基本是不现实的，而scRNA可以不用通过实验方法的纯化细胞，就能研究细胞生命的各个状态。
这也是Monocle最大的亮点—利用算法去学习每个细胞基因表达量的变化，从而对细胞状态进行推断，推断出整体的一个表达轨迹，那么就可以将每个细胞放在特定位置上；另外有可能一个过程有多个结果，那么在轨迹图上就会有多个分支，代表了细胞分化
上面使用的数据集到此结束，下面发育轨迹会用到一套新的数据集（来自Packer & Zhu et al. 的线虫整个胚胎发育的时间线数据集，这里主要节选了神经元的数据）
需要注意的是，下面做出的图和官网给出的图是不同的
载入数据，构建对象 => new_cell_data_set()
expression_matrix = readRDS(url("http://staff.washington.edu/hpliner/data/packer_embryo_expression.rds"))
cell_metadata = readRDS(url("http://staff.washington.edu/hpliner/data/packer_embryo_colData.rds"))
gene_annotation = readRDS(url("http://staff.washington.edu/hpliner/data/packer_embryo_rowData.rds"))
# 文件大小(同样可以公众号后台回复"monocle-2"获取)
#|-- [497K] packer_embryo_colData.rds
#|-- [10.0M] packer_embryo_expression.rds
#`-- [226K] packer_embryo_rowData.rds
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
预处理 => preprocess_cds()
和示例一的操作差不多，只是这里根据原作者（Packer & Zhu et al ）的方法处理了一下批次效应
cds <- preprocess_cds(cds, num_dim = 100, 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")
降维可视化 => reduce_dimension()
强烈建议用UMAP
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)
细胞聚类 => cluster_cells()
Monocle不认为一组数据中的所有细胞都来自同一个”祖先“，很多实验中，它们会有多个发育轨迹。Monocle会通过聚类来判断细胞是否应该归属同一个发育轨迹。之前介绍的cluster_cells()中，细胞可以按cluster细分，还可以按partition归为大类。
# 这里就用partition聚类
cds <- cluster_cells(cds)
plot_cells(cds, color_cells_by = "partition")
利用learn graph在每个partition中寻找主路径 => learn_graph()
理论上，任何的UMAP或者tSNE降维结果都能找到这样的主路径，选用UMAP是考虑了它能更多展示细胞之间的关联，另外最重要是对数据的了解
cds <- learn_graph(cds)
plot_cells(cds,
color_cells_by = "cell.type",
label_groups_by_cluster=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE)
有了初步的轨迹线，接下来就是对细胞出现的先后进行排序 => order_cells()
也就是pseudotime（拟时序分析），推断哪个细胞先出来，哪个细胞后出来。它会对每个细胞在特定的过程（比如分化过程）中的参与度进行评估
很多生命活动中，细胞并非同步发育的，例如在细胞分化阶段，实验捕获到的细胞可能会分布到各个发育过程。因此做单细胞测序时得到的那一群细胞，其中可能包括已经完成某个生命过程的细胞，还可能包括没有开始这个过程的细胞。这种非同步的特性会阻碍我们理解细胞过渡状态究竟发生了怎样的调控。
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分配一个根节点。
不过，怎么简单高效挑选根节点呢？=> get_earliest_principal_node()
首先根据轨迹图中的节点将与它们最相近的细胞分成组，然后计算来自最早时间点的细胞所占组分，最后挑出包含早期细胞数量最多的节点，认为它是就是根节点
# 官方给出了一个函数，这里定义了一个time_bin，选择了最早的时间点区间。
get_earliest_principal_node <- function(cds, time_bin="130-170"){
# 首先找到出现在最早时间区间的细胞ID
cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
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,]))))]
root_pr_nodes
}
cds = order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
然后进行可视化：
plot_cells(cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5)
利用3D的发育轨迹对上述内容做个概述
3D轨迹实际上就是降维时选前3个主成分 => max_components = 3，后续都和2D保持类似
需要注意的是，这个3D图需要等待一段时间加载，加载出来的图可以用鼠标拖动
# 3D trajectories（4步走）
cds_3d = reduce_dimension(cds, max_components = 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差异分析
版本3和版本2的差异分析可以说是完全不同，版本3取代了2中的
differentialGeneTest()andBEAM()
Monocle3提供了不同细胞类型之间寻找差异基因的方法，主要有两种：
Regression analysis：利用
fit_models()，用来评价基因表达是否会受到诸如时间、处理等的影响Graph-autocorrelation analysis：利用
graph_test()，用来寻找一条轨迹或不同cluster中基因的差异
方法一：Regression analysis => fit_models()
这一部分将会根据几种不同的标准进行差异分析，并且不同复杂程度的分析用时也有较大差别（几分钟甚至几小时）。教程只会针对一小部分基因进行最简单的分析，但实际上它可以处理上千基因
测试小基因集为：
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"
然后看看哪些基因是与时间因素相关的：
fit_coefs = coefficient_table(gene_fits)
# 挑出时间相关的组分
emb_time_terms = fit_coefs %>% filter(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)
# 结果数据集中的6个基因有5个是时间相关的差异基因
除此以外，还能画小提琴图：=> plot_genes_violin()
plot_genes_violin(cds_subset, group_cells_by="embryo.time.bin", ncol=2) +
theme(axis.text.x=element_text(angle=45, hjust=1))
如果要考虑批次效应
gene_fits = fit_models(cds_subset, model_formula_str = "~embryo.time + batch")
fit_coefs = coefficient_table(gene_fits)
fit_coefs %>% filter(term != "(Intercept)") %>%
select(gene_short_name, term, q_value, estimate)
模型做出来了，效果到底如何？=> evaluate_fits()
对模型进行评价，使用evaluate_fits()
gene_fits = fit_models(cds_subset, model_formula_str = "~embryo.time")
evaluate_fits(gene_fits)
那么模型中到底要不要考虑批次呢？使用compare_models()比较判断。结果会返回一个likelihood ratio test结果，
time_batch_models = fit_models(cds_subset,
model_formula_str = "~embryo.time + batch",
expression_family="negbinomial")
time_models = fit_models(cds_subset,
model_formula_str = "~embryo.time",
expression_family="negbinomial")
compare_models(time_batch_models, time_models) %>% select(gene_short_name, q_value)
其中第一个time_batch_models叫做full model，这个模型知道细胞的收集时间和批次，来预测其中每个基因的表达；第二个time_models叫做reduced model，这个模型只知道细胞的收集时间。因为第一个模型中信息更丰富，所以它对每个细胞中基因表达预测也更准。而Monocle要回答的问题是：第一个模型相比于第二个，究竟做了多少提升？ 如果提升的部分大多来自第一个模型的批次效应，那么得到的likelihood ratio test结果就越显著
gene_short_name q_value
<chr> <dbl>
1 ham-1 3.74e-20
2 dmd-6 7.11e-37
3 hlh-17 1.04e- 5
4 nhr-6 5.95e- 3
5 ceh-36 6.23e-10
6 che-1 5.06e-6
结果看到，所有基因的likelihood ratio tests结果都是极显著，说明了数据中的确存在显著的批次效应，也验证了添加批次信息的可靠性。
上面所使用的线性模型方法要指定一个基因表达量的分布，大多数研究使用负二项分布，这种分布往往是符合测序reads或UMI count值的，这个方法也被应用在RNA-seq分析（如DESeq2）中
fit_models()函数支持多种分布，默认是：quasipoisson ，它和负二项分布相似，只不过比负二项分布准确性低一点、速度会更快，对于细胞数量多更适用。
关于其他的分布（使用expression_family参数设置）：
方法二：Graph-autocorrelation analysis => graph_test()
注意：下面的代码和图片看看就好（因为不能运行）
# 选一部分神经元细胞
neurons_cds = cds[,colData(cds)$assigned_cell_type == "Neurons"]
plot_cells(neurons_cds, color_cells_by="partition")
从图中可以看到，存在不同类型的神经元细胞，它们也许来自不同亚群。为了探索这些不同类型的细胞之间是哪些基因差异表达，一种方法就是上面的线性回归。不过对于UMAP或tSNE的图，graph_test()这个函数可以进行一个叫做Moran's I 的spatial autocorrelation分析方法
Moran's I这个指标叫：莫兰指数，澳大利亚统计学家莫兰于1950年提出。它是用来度量空间相关性的一个指标，数据经过方差归一化之后，会落在-1.0和1.0之间。这个值大于0表示空间正相关，值越大空间相关性越强；等于0时空间是随机性分布；小于0时表示空间负相关，值越小空间差异越大
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))
pr_graph_test_res结果是一个数据框，存储了cds对象中每个基因的Moran's I检验结果。正值表示基因分布在UMAP空间的一小块”焦点“区域（例如一个或几个cluster）。
现在有了和”焦点“区域和基因，那么怎么把二者结合起来？
将共同作用的基因合并成模块 => find_gene_modules()
根据这些共同作用基因的表达量对所有细胞进行UMAP，然后利用Louvain community analysis将它们聚集成小模块，然后将小模块和大cluster甚至更大的partition联系起来
gene_module_df = find_gene_modules(neurons_cds[pr_deg_ids,], resolution=1e-2)
得到的结果gene_module_df中有一行记录了每个基因和module的对应关系
要看哪个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)
就知道dac-1基因的激活发生在其他两个基因之前
如果使用版本2…
版本2的部分应该是要新写一篇的（后期粘上版本2的链接）
首先创建对象：
需要指定
exprs、phenoData、featureData，后两个都要是AnnotatedDataFrame对象（目的是不让用户随便修改，看来还是版本3比较人性化），然后依然是要与表达矩阵行、列对应
#do not run
HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
HSMM_gene_annotation <- read.delim("gene_annotations.txt")
pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
phenoData = pd, featureData = fd)
体验
一句话就是：满怀期待，收获满满，看得费劲
优点
创建对象变得方便许多，从原来的
new()函数创建对象到现在的数据框，减少了函数的记忆据说是差异分析算法进行了更新，具体新不新有待探索
缺点
Monocle3的示例数据做的不是很走心，数据下载速度受限。这一点不如Seurat做的清晰明了
另外很多地方不能重复出来（比如目前加载的shinny功能还不完善，会出现无响应的状态；文档还有一些小错误）
许多代码操作过程复杂
版本2、3之间代沟很大，像是完全不同的两个版本，而不是像Seurat2、3几个参数、函数的改变
可视化的结果好像比版本2少了一些，不知道这部分是没有改变还是被整合了
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~