monocle多样本拟时序分析

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 
gene_ann   gene_short_name = row.names(ct), 
  row.names = row.names(ct)
)

pd           [email protected])
fd           data=gene_ann)

sc_cds   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 sc_cds table([email protected]$celltype)

## 
## CD14 Mono CD16 Mono 
##      4355      1044

fdif = "diff_test_res.Rdata"
if(!file.exists(fdif)){
  diff_test_res                                         fullModelFormulaStr = " ~ celltype + orig.ident", 
                                           reducedModelFormulaStr = " ~ orig.ident", 

                                        cores = 8)
  save(diff_test_res,file = fdif)
}
load(fdif)
ordering_genes 0.01))
#然后是查看基因，设置为排序要使用的基因
head(ordering_genes)

## [1] "NOC2L"        "ISG15"        "AGRN"         "RP4-758J18.2" "SSU72"       
## [6] "GNB1"

sc_cds plot_ordering_genes(sc_cds)

#降维
sc_cds "~orig.ident")
#细胞排序
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