monocle单样本拟时序分析

0.背景知识

做拟时序分析是为了探索自己感兴趣的几种细胞之间的发育关系，一般不是用全部类型的细胞来做的。例如本例中选择了CD14和CD16单核细胞。

如果让ai来说拟时序的目的那就比我说的多:

拟时序分析（Pseudo-time analysis）是一种用于理解细胞状态和细胞命运变化的计算生物学方法。这种分析通常用于单细胞测序数据，其目的是重建细胞发育或疾病进程中的时间序列，即使实际的时间信息不可用。以下是进行拟时序分析的几个主要目的：

1. 细胞状态推断：通过分析单个细胞的基因表达模式，推断细胞在生物学过程中所处的状态。
2. 细胞轨迹追踪：重建细胞从一种状态转变到另一种状态的路径，这有助于理解细胞分化、发育或疾病进展的过程。
3. 细胞命运预测：预测细胞随时间可能采取的发育轨迹，尤其是在干细胞分化或癌症发展的研究中。
4. 动态过程建模：创建细胞状态变化的动态模型，这有助于揭示细胞行为的内在规律和调控机制。
5. 细胞亚群识别：在复杂的细胞群体中识别不同的细胞亚群，并理解它们在生物学过程中的作用。
6. 基因调控网络推断：通过分析基因表达随“拟时间”的变化，推断基因调控网络和信号传导途径。
7. 疾病机理探索：在疾病研究中，拟时序分析有助于揭示疾病发生和发展的分子机制。
8. 药物作用机制研究：通过观察药物处理前后细胞状态的变化，研究药物的作用机制和效果。 拟时序分析是一种强大的工具，它可以帮助研究者在没有直接时间标记的情况下，通过基因表达数据来探索细胞状态的变化和动态过程。这种方法在单细胞生物学、发育生物学、癌症生物学和神经科学等领域有着广泛的应用。

今天的代码是处理单样本数据的，明天再整一个多样本的。

这里的示例数据scRNA.Rdata是seurat降维聚类分群注释的结果，在生信星球聊天框回复“958monocle”获取。

rm(list = ls())
library(Seurat)
library(monocle)
library(dplyr)
load("scRNA.Rdata")
table(scRNA$celltype)

## 
##            B  B Activated    CD14 Mono    CD16 Mono CD4 Memory T  CD4 Naive T 
##          568          201         2144          537          903         1524 
##        CD8 T           DC           Mk           NK          pDC  T activated 
##          462          214          121          320           81          332

scRNA = subset(scRNA,idents = c("CD14 Mono","CD16 Mono"))
table(scRNA$celltype)

## 
## CD14 Mono CD16 Mono 
##      2144       537

本文的输入数据是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",
          [email protected])
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, 2681 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: AAACATACGCGAAG.1 AAACATACTCAGGT.1 ... TTTGACTGCCCTAC.1
##     (2681 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([email protected]$celltype)

## 
## CD14 Mono CD16 Mono 
##      2144       537

fdif = "diff_test_res.Rdata"
if(!file.exists(fdif)){
  diff_test_res <- differentialGeneTest(sc_cds,
                                        fullModelFormulaStr = "~celltype",
                                        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] "ISG15"    "SSU72"    "GNB1"     "TNFRSF14" "WRAP73"   "C1orf174"

sc_cds <- setOrderingFilter(sc_cds, ordering_genes)
plot_ordering_genes(sc_cds)

#降维
sc_cds <- reduceDimension(sc_cds)
#细胞排序
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是怎样变化的。

经典的拟时序热图

展示了一些基因是如何随着时间轨迹的变化而渐变的，这个渐变不同于findmarkers，是体现变化过程的，而不是直接给出差异表达的基因。

gene_to_cluster = diff_test_res %>% arrange(qval) %>% head(50) %>% pull(gene_short_name);head(gene_to_cluster)

## [1] "FCGR3A" "PLAC8"  "MS4A4A" "MS4A7"  "LYZ"    "VMO1"

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 )