拟时分析

生信喵

构建细胞分化轨迹
一、Trajectory analysis 轨迹推断
       轨迹推断（Trajectory Inference，TI），也称为细胞分析轨迹（differentiation trajectories）基于单细胞转录组数据，利用模型预测细胞分化过程。
       由于分析结果是基于某个组织一瞬间的群体单细胞数据预测得到的，并非真实时间维度上的分化结果，因此也被称为伪时序分析 (pseudotime analysis)，也可以叫做拟时分析。
       拟时分析通过分析关键基因的表达模式，将所有细胞按照发育时间的先后排布在拟时间轴上，模拟发育过程中的细胞分化过程。通过对细胞轨迹的分析，我们可以挖掘出细胞分化过程中经历的细胞类型变化、伴随发育过程变化的动态变化基因、祖细胞不同的分化命运等与生命发育息息相关的信息。

二、monocle3 进行细胞轨迹推断
文档：https://cole-trapnell-lab.github.io/monocle3/docs/trajectories/

1. #1 读取数据
   
2. getwd()
   
3. setwd('/home/xhs/jyxy/16.rscrna/')
   
4. expression_matrix <- readRDS("./monocle3/packer/packer_embryo_expression.rds")
   
5. cell_metadata <- readRDS("./monocle3/packer/packer_embryo_colData.rds")
   
6. gene_annotation <- readRDS("./monocle3/packer/packer_embryo_rowData.rds")
   
7. 
   
8. cds <- new_cell_data_set(expression_matrix,cell_metadata = cell_metadata,gene_metadata = gene_annotation)
   
9. 
   
10. #2 数据处理以及去除批次效应
    
11. cds <- preprocess_cds(cds, num_dim = 50)
    
12. 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")
    
13. 
    
14. 
    
15. #3 数据降维以及结果可视化
    
16. cds <- reduce_dimension(cds)
    
17. #设置group_label_size = 3,cell_size = 0.7调整文本和点的大小
    
18. plot_cells(cds, label_groups_by_cluster=FALSE,color_cells_by = "cell.type",group_label_size = 3,cell_size = 0.7)
    
19. 
    
20. #基因沿着轨迹变化
    
21. ciliated_genes <- c("che-1","hlh-17","nhr-6","dmd-6","ceh-36","ham-1")
    
22. plot_cells(cds,genes=ciliated_genes,label_cell_groups=FALSE,show_trajectory_graph=FALSE)
    
23. 
    
24. #4 聚类
    
25. cds <- cluster_cells(cds)
    
26. plot_cells(cds, color_cells_by = "partition")
    
27. 
    
28. #5 构建Learn the trajectory graph
    
29. cds <- learn_graph(cds)
    
30. plot_cells(cds,color_cells_by = "cell.type",label_groups_by_cluster=FALSE,label_leaves=FALSE,label_branch_points=FALSE)
    
31. 
    
32. #6排序 Order the cells in pseudotime
    
33. 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)
    
34. 
    
35. #7 选取root跟节点
    
36. cds <- order_cells(cds)
    
37. 
    
38. plot_cells(cds,
    
39.            color_cells_by = "pseudotime",
    
40.            label_cell_groups=FALSE,
    
41.            label_leaves=FALSE,
    
42.            label_branch_points=FALSE,
    
43.            graph_label_size=1.5)
    
44. 
    
45. #8 自定义函数，确定根节点
    
46. get_earliest_principal_node <- function(cds, time_bin="130-170"){
    
47.   cell_ids <- which(colData(cds)[, "embryo.time.bin"] == time_bin)
    
48.   closest_vertex <-
    
49.     cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
    
50.   closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
    
51.   root_pr_nodes <-
    
52.     igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
    
53.                                                               (which.max(table(closest_vertex[cell_ids,]))))]
    
54.   root_pr_nodes
    
55. }
    
56. cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
    
57. plot_cells(cds,color_cells_by = "pseudotime",label_cell_groups=FALSE,label_leaves=FALSE,label_branch_points=FALSE,graph_label_size=1.5)
    
58. 
    
59. #9 选取树枝Subset cells by branch
    
60. cds_sub <- choose_graph_segments(cds)
    
61. 
    
62. 
    
63. #10 绘制3D 轨迹图
    
64. cds_3d <- reduce_dimension(cds, max_components = 3)
    
65. cds_3d <- cluster_cells(cds_3d)
    
66. cds_3d <- learn_graph(cds_3d)
    
67. cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))
    
68. 
    
69. cds_3d_plot_obj <- plot_cells_3d(cds_3d, color_cells_by="partition")

复制代码