monocle中有两种差异分析模型：
-
Regression analysis: using
fit_models(), you can evaluate whether each gene depends on variables such as time, treatments, etc. -
Graph-autocorrelation analysis: using
graph_test(), you can find genes that vary over a trajectory or between clusters.
1) Regression analysis
Monocle通过对每个基因拟合回归模型来进行分析。指定该模型以考虑实验中的各种因素（如时间、处理等）。
gene_fits <- fit_models(cds_subset, model_formula_str = "~embryo.time")
gene_fits是一个包含每个基因一行的表格。模型列包含广义线性模型对象，每个模型旨在使用上述方程解释细胞中基因的表达。参数model_formula_str应为一个字符串，指定模型公式。使用的模型公式可以包含colData中存在的任何列，包括在Monocle其他分析步骤中添加的列。例如，可以通过使用~cluster或~partition作为模型公式来测试不同cluster和partition之间的基因。也可以包含多个变量，例如~embryo.time + batch。
我们从每个模型中提取系数表。
fit_coefs <- coefficient_table(gene_fits)
| id | gene_short_name | num_cells_expressed | status | term | estimate | std_err | test_val | p_value | normalized_effect | model_component | q_value |
|---|---|---|---|---|---|---|---|---|---|---|---|
| WBGene00001820 | ham-1 | 1618 | OK | (Intercept) | 1.09 | 0.0665 | 16.4 | 4.58e-59 | 0 | count | 1.83e-58 |
| WBGene00001820 | ham-1 | 1618 | OK | embryo.time | -0.00388 | 0.000211 | -18.4 | 9.18e-74 | -0.00558 | count | 4.59e-73 |
| WBGene00007058 | dmd-6 | 873 | OK | (Intercept) | -5.07 | 0.123 | -41.3 | 0 | 0 | count | 0 |
| WBGene00007058 | dmd-6 | 873 | OK | embryo.time | 0.0075 | 0.000209 | 35.9 | 3.39e-256 | 0.00418 | count | 2.03e-255 |
| WBGene00001961 | hlh-17 | 51 | OK | (Intercept) | -6.04 | 0.501 | -12.1 | 4.26e-33 | 0 | count | 8.52e-33 |
| WBGene00001961 | hlh-17 | 51 | OK | embryo.time | 0.00334 | 0.00103 | 3.23 | 0.00123 | 0.00093 | count | 0.00247 |
我们通常不关心截距项β₀，因此可以轻松提取时间项：
emb_time_terms <- fit_coefs %>% filter(term == "embryo.time")
## 提取具有显著时间成分的基因
emb_time_terms %>% filter (q_value < 0.05) %>%
select(gene_short_name, term, q_value, estimate)
coefficient_table()函数测试每个系数在Wald检验下是否显著不同于零。默认情况下，coefficient_table()会使用Benjamini和Hochberg的方法对这些p值进行多重假设检验的调整，调整后的值可以在q_value列中找到。
| gene_short_name | term | q_value | estimate |
|---|---|---|---|
| ham-1 | embryo.time | 4.59e-73 | -0.00388 |
| dmd-6 | embryo.time | 2.03e-255 | 0.0075 |
| hlh-17 | embryo.time | 0.00247 | 0.00334 |
| nhr-6 | embryo.time | 8.41e-28 | 0.00579 |
| ceh-36 | embryo.time | 1.45e-34 | 0.00279 |
差异基因可视化
plot_genes_violin(cds_subset, group_cells_by="embryo.time.bin", ncol=2) +
theme(axis.text.x=element_text(angle=45, hjust=1))
plot_genes_hybrid(cds_subset, group_cells_by="embryo.time.bin", ncol=2) +
theme(axis.text.x=element_text(angle=45, hjust=1))
2) Graph-autocorrelation analysis
数据框pr_graph_test_res包含了cell_data_set中每个基因的Moran's I检验结果，可以根据morans_I列对该表进行排序，该列的范围从-1到+1。值为0表示没有效应，而+1表示完美的正自相关，表明附近细胞在基因表达值上非常相似。显著的值小于零通常很少见。
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))
