0.背景知识
做拟时序分析是为了探索自己感兴趣的几种细胞之间的发育关系，一般不是用全部类型的细胞来做的。例如本例中选择了CD14和CD16单核细胞。
如果让ai来说拟时序的目的那就比我说的多:
拟时序分析（Pseudo-time analysis）是一种用于理解细胞状态和细胞命运变化的计算生物学方法。这种分析通常用于单细胞测序数据，其目的是重建细胞发育或疾病进程中的时间序列，即使实际的时间信息不可用。以下是进行拟时序分析的几个主要目的：
- 细胞状态推断：通过分析单个细胞的基因表达模式，推断细胞在生物学过程中所处的状态。
- 细胞轨迹追踪：重建细胞从一种状态转变到另一种状态的路径，这有助于理解细胞分化、发育或疾病进展的过程。
- 细胞命运预测：预测细胞随时间可能采取的发育轨迹，尤其是在干细胞分化或癌症发展的研究中。
- 动态过程建模：创建细胞状态变化的动态模型，这有助于揭示细胞行为的内在规律和调控机制。
- 细胞亚群识别：在复杂的细胞群体中识别不同的细胞亚群，并理解它们在生物学过程中的作用。
- 基因调控网络推断：通过分析基因表达随“拟时间”的变化，推断基因调控网络和信号传导途径。
- 疾病机理探索：在疾病研究中，拟时序分析有助于揭示疾病发生和发展的分子机制。
- 药物作用机制研究：通过观察药物处理前后细胞状态的变化，研究药物的作用机制和效果。 拟时序分析是一种强大的工具，它可以帮助研究者在没有直接时间标记的情况下，通过基因表达数据来探索细胞状态的变化和动态过程。这种方法在单细胞生物学、发育生物学、癌症生物学和神经科学等领域有着广泛的应用。
今天的代码是处理多样本数据的，单样本的在前一篇。
使用的数据sce.all.Rdata是Seurat降维聚类分群注释的结果，有两个样本。在生信星球聊天框回复“959monocle”可以拿到我的示例数据
rm(list = ls())
library(Seurat)
library(monocle)
library(dplyr)
load("../monocle_one_sample/sce.all.Rdata")
table(sce.all$celltype)
##
## B B Activated CD14 Mono CD16 Mono CD4 Memory T CD4 Naive T
## 975 386 4355 1044 1762 2501
## CD8 T DC Mk NK pDC T activated
## 814 472 236 618 132 631
scRNA = subset(sce.all,idents = c("CD14 Mono","CD16 Mono"))
table(scRNA$celltype)
##
## CD14 Mono CD16 Mono
## 4355 1044
table(scRNA$orig.ident)
##
## IMMUNE_CTRL IMMUNE_STIM
## 2718 2681
本文的输入数据是seurat做完降维聚类分群注释的数据，并将注释的结果添加到了meta表格里面成为了celltype列。
因为只想要CD14和CD16单核细胞，所以提取出来相应的子对象。
因为monocle和seurat是两个不同的体系，所以需要将seurat对象转换为monocle可以接受的CellDataSet对象。虽然monocle3已经出来很久了，但大家都不约而同的选择monocle2，大概就是习惯了吧。。
ct <- scRNA@assays$RNA$counts
gene_ann <- data.frame(
gene_short_name = row.names(ct),
row.names = row.names(ct)
)
pd <- new("AnnotatedDataFrame",
data=scRNA@meta.data)
fd <- new("AnnotatedDataFrame",
data=gene_ann)
sc_cds <- newCellDataSet(
ct,
phenoData = pd,
featureData =fd,
expressionFamily = negbinomial.size(),
lowerDetectionLimit=1)
sc_cds
## CellDataSet (storageMode: environment)
## assayData: 14053 features, 5399 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: AAACATACATTTCC.1 AAACATACCAGAAA.1 ... TTTGACTGCCCTAC.1
## (5399 total)
## varLabels: orig.ident nCount_RNA ... Size_Factor (19 total)
## varMetadata: labelDescription
## featureData
## featureNames: AL627309.1 RP11-206L10.2 ... LRRC3DN (14053 total)
## fvarLabels: gene_short_name
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
1.构建细胞发育轨迹
选择基因有不同的策略，比如 1.使用seurat给出的高变化基因 2.按照平均表达量大于某个数字(比如0.1，官网用的是这个)的基因 3.使用不同细胞类型之间的差异基因，differentialGeneTest计算。 我们默认使用的是最后一个策略。
sc_cds <- estimateSizeFactors(sc_cds)
sc_cds <- estimateDispersions(sc_cds)
table(scRNA@meta.data$celltype)
##
## CD14 Mono CD16 Mono
## 4355 1044
fdif = "diff_test_res.Rdata"
if(!file.exists(fdif)){
diff_test_res <- differentialGeneTest(sc_cds,
fullModelFormulaStr = " ~ celltype + orig.ident",
reducedModelFormulaStr = " ~ orig.ident",
cores = 8)
save(diff_test_res,file = fdif)
}
load(fdif)
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
#然后是查看基因，设置为排序要使用的基因
head(ordering_genes)
## [1] "NOC2L" "ISG15" "AGRN" "RP4-758J18.2" "SSU72"
## [6] "GNB1"
sc_cds <- setOrderingFilter(sc_cds, ordering_genes)
plot_ordering_genes(sc_cds)
#降维
sc_cds <- reduceDimension(sc_cds,residualModelFormulaStr = "~orig.ident")
#细胞排序
sc_cds <- orderCells(sc_cds)
2.绘图展示
发育轨迹图
library(ggsci)
p1 = plot_cell_trajectory(sc_cds)+ scale_color_nejm()
p2 = plot_cell_trajectory(sc_cds, color_by = 'Pseudotime')
p3 = plot_cell_trajectory(sc_cds, color_by = 'celltype') + scale_color_npg()
library(patchwork)
p2+p1/p3
这三种着色方式放在一起非常的带劲，很清晰的展示了pseudotime、state和celltype是怎样变化的。
以orig.ident着色可以看出，不同样本中的细胞基本是均匀分布在轨迹上的，说明前面的代码很好的去除了样本间的批次效应
plot_cell_trajectory(sc_cds, color_by = 'orig.ident')
经典的拟时序热图
展示了一些基因是如何随着时间轨迹的变化而渐变的，这个渐变不同于findmarkers，是体现变化过程的，而不是直接给出差异表达的基因。
gene_to_cluster = diff_test_res %>% arrange(qval) %>% head(50) %>% pull(gene_short_name);head(gene_to_cluster)
## [1] "S100A9" "S100A8" "FCGR3A" "IL8" "PLAC8" "TMSB4X"
plot_pseudotime_heatmap(sc_cds[gene_to_cluster,],
num_clusters = nlevels(Idents(scRNA)),
show_rownames = TRUE,
cores = 4,return_heatmap = TRUE,
hmcols = colorRampPalette(c("navy", "white", "firebrick3"))(100))
用感兴趣的基因给轨迹图着色，gs可以换成你想换的基因
gs = head(gene_to_cluster)
plot_cell_trajectory(sc_cds,markers=gs,
use_color_gradient=T)
也可以是jitter
plot_genes_jitter(sc_cds[gs,],
grouping = "celltype",
color_by = "celltype",
nrow= 3,
ncol = NULL )
