一旦找到了感兴趣的差异基因，Monocle可以将它们分组为module。通过调用find_gene_modules()函数，在基因（而不是细胞）上运行UMAP，然后使用Louvain community分析将它们分组到module中。
gene_module_df <- find_gene_modules(neurons_cds[pr_deg_ids,], resolution=1e-2)
数据框gene_module_df的行是基因，并标识其所属的module。要查看哪些module在各个簇或分区中被表达，可以使用两种不同的可视化方法。
第一种方法是制作一个简单的表格，显示每个module中所有基因在所有簇中的汇总表达。Monocle提供了一个名为aggregate_gene_expression的简单工具函数来实现这一目的。
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)
第二种方法是将gene_module_df直接传递给plot_cells()。如果模块数量较多，可能很难看出每个module的表达位置，因此我们将只查看其中的一部分：
plot_cells(neurons_cds,
genes=gene_module_df %>% filter(module %in% c(8, 28, 33, 37)),
group_cells_by="partition",
color_cells_by="partition",
show_trajectory_graph=FALSE)
寻找随pseudotime变化的基因
我们如何找到在轨迹不同路径上差异表达的基因？我们如何找到那些仅限于轨迹开始部分的基因？或者被排除在外的基因？
我们再次使用graph_test()，使用参数neighbor_graph="principal_graph"，这告诉它测试轨迹上相似位置的细胞是否具有相关的表达。
ciliated_cds_pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=4)
pr_deg_ids <- row.names(subset(ciliated_cds_pr_test_res, q_value < 0.05))
plot_cells(cds, genes=c("hlh-4", "gcy-8", "dac-1", "oig-8"),
show_trajectory_graph=FALSE,
label_cell_groups=FALSE,
label_leaves=FALSE)
gene_module_df <- find_gene_modules(cds[pr_deg_ids,], resolution=c(10^seq(-6,-1)))
如果想看具体的某几个基因随pseudotime变化的趋势，可以用函数plot_genes_in_pseudotime()：
AFD_genes <- c("gcy-8", "dac-1", "oig-8")
AFD_lineage_cds <- cds[rowData(cds)$gene_short_name %in% AFD_genes,
colData(cds)$cell.type %in% c("AFD")]
AFD_lineage_cds <- order_cells(AFD_lineage_cds)
plot_genes_in_pseudotime(AFD_lineage_cds,
color_cells_by="embryo.time.bin",
min_expr=0.5)
