monocle多样本拟时序分析

生信技能树

发布于 2024-07-05 14:44:15

25600

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

前面已经是介绍了单个样品的单细胞转录组表达量矩阵的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 <- 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 )