Monocle 3 | 太牛了！单细胞必学R包！~（三）（建立单细胞轨迹）

1写在前面

到周末了，今天去骑车的，感觉还不错。🚴🏻

下面就是今天的内容了： 👇

单细胞转录组、蛋白组、表观组学等单细胞技术的发展为研究细胞周期、细胞分化等细胞动态过程提供了新的机会。🤩

使用轨迹推断（TI，trajectory inference）的方法可以根据测序的细胞之间表达模式的相似性对单细胞沿着轨迹进行排序，模拟细胞动态变化的过程。🥸

轨迹推断也称作伪时间分析（pseudotime analysis）。🤩

2用到的包

rm(list = ls())
library(tidyverse)
library(monocle3)

3示例数据

expression_matrix <- readRDS("./packer_embryo_expression.rds")
cell_metadata <- readRDS("./packer_embryo_colData.rds")
gene_annotation <- readRDS("./packer_embryo_rowData.rds")

cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)
cds

4预处理数据

我们先处理一下批次效应，这里有一个原文中提供的方法，我们用一下吧。😏

自己的数据还是要自己的方法处理的。😂

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")

5降维并可视化结果

我们只是用了data的一小部分，但Monocle重建了一个有许多分支的轨迹。🥰

叠加人工注释后，这些分支主要由一种细胞类型占据。😋

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)

6细胞聚类

尽管细胞可以连续地从一个状态过渡到下一个状态，之间没有明显边界，但Monocle并不假设数据集中的所有细胞都是从一个共同的转录"祖先"转化而来。😘

事实上可能有多个不同的轨迹。🧐

比如，对感染的反应时，固有免疫细胞和基质细胞会有非常不同的初始转录组，对感染的反应也很不同，所以它们应该是同一轨迹的一部分。😍

细胞会分成一个个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)

7轨迹排序

首先我们需要找到轨迹的roots，这个可能需要根据你的实验目的来确定。🥸

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)

---

选择root。😀

cds <- order_cells(cds)

---

这里我随便标了2个root，我们看一下效果。🥰

plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           graph_label_size=1.5)

---

我知道你肯定觉得手动挑选root很累。😋

好的，代码方式来了。😂

get_earliest_principal_node <- function(cds, time_bin="130-170"){
  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)

8挑选分支

这里我随便挑选一下分支哦。😜

cds_sub <- choose_graph_segments(cds)

93D展示

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")

---

最后祝大家早日不卷!~

---

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

📍 往期精彩

📍 🤩 LASSO | 不来看看怎么美化你的LASSO结果吗！？
📍 🤣 chatPDF | 别再自己读文献了！让chatGPT来帮你读吧！~
📍 🤩 WGCNA | 值得你深入学习的生信分析方法！~
📍 🤩 ComplexHeatmap | 颜狗写的高颜值热图代码！
📍 🤥 ComplexHeatmap | 你的热图注释还挤在一起看不清吗！？
📍 🤨 Google | 谷歌翻译崩了我们怎么办！？（附完美解决方案）
📍 🤩 scRNA-seq | 吐血整理的单细胞入门教程
📍 🤣 NetworkD3 | 让我们一起画个动态的桑基图吧~
📍 🤩 RColorBrewer | 再多的配色也能轻松搞定！~
📍 🧐 rms | 批量完成你的线性回归
📍 🤩 CMplot | 完美复刻Nature上的曼哈顿图
📍 🤠 Network | 高颜值动态网络可视化工具
📍 🤗 boxjitter | 完美复刻Nature上的高颜值统计图
📍 🤫 linkET | 完美解决ggcor安装失败方案（附教程）
📍 ......

---

---

本文由 mdnice 多平台发布