十八般武艺 5 的单细胞分析:monocle3
单细胞测序技术的发展日新月异,新的分析工具也层出不穷。每个工具都有它的优势与不足,在没有权威工具和流程的单细胞生信江湖里,多掌握几种分析方法和工具,探索数据时常常会有意想不到的惊喜。
往期专题
单细胞初级8讲和高级分析8讲 单细胞分析十八般武艺1:harmony 单细胞分析十八般武艺2:LIGER 单细胞分析十八般武艺3:fastMNN 单细胞分析十八般武艺4:velocyto
monocle3简介
monocel3的优势
- 从UMAP图识别发育轨迹,可以继承Seurat的质控、批次校正和降维分析结果,实现“一张图”展现细胞的聚类、鉴定和轨迹分析结果。
- 自动对UMAP图分区(partition),可以选择多个起点,轨迹分析算法的逻辑更符合生物学现实。
除了轨迹分析的主要功能,monocle3差异分析方法也有其独到之处,可以做一些与seurat不好实现的分析。
monocel3的安装
先安装一些依赖包,大家安装前可以查看一下这些包是否已经安装过了。
BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
'limma', 'S4Vectors', 'SingleCellExperiment',
'SummarizedExperiment', 'batchelor', 'Matrix.utils'))
然后安装实现umap图分区的包leidenbase,最后安装monole3
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')
安装有困难的朋友可以使用我的镜像kinesin/rstudio:1.2,下载链接见《kinesin_rstudio的日常升级二》,使用方法见《华为云配置单细胞分析环境及报错处理》。
monocle3分析实践
数据来源
数据来自Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer,数据集GEO编号:GSE139324。建议大家自己下载,也可以加Kinesin微信获取数据的百度云链接。原数据集有63个scRNA的数据,都是分选的CD45+免疫细胞。考虑到计算资源问题,挑选了10个样本用于此次演示。
最后用于轨迹分析的细胞是提取的T细胞亚群。
创建seurat对象
library(Seurat)
library(monocle3)
library(tidyverse)
library(patchwork)
rm(list=ls())
##==准备seurat对象列表==##
dir <- dir("GSE139324/") #GSE139324是存放数据的目录
dir <- paste0("GSE139324/",dir)
sample_name <- c('HNC01PBMC', 'HNC01TIL', 'HNC10PBMC', 'HNC10TIL', 'HNC20PBMC',
'HNC20TIL', 'PBMC1', 'PBMC2', 'Tonsil1', 'Tonsil2')
scRNAlist <- list()
for(i in 1:length(dir)){
counts <- Read10X(data.dir = dir[i])
scRNAlist[[i]] <- CreateSeuratObject(counts, project=sample_name[i], min.cells=3, min.features = 200)
scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-")
}
scRNA <- merge(scRNAlist[[1]], scRNAlist[2:10])
##提取表达矩阵用于细胞类型鉴定
count <- GetAssayData(scRNA, assay = "RNA", slot = "counts")
saveRDS(count, "count.rds")
利用azimuth鉴定细胞类型
Seurat官网近期推出了在线细胞类型鉴定服务,可以准确鉴定pbmc细胞,建议大家尝试一下。
点击箭头所指的Browse按钮,将上一步保存的count.rds
文件上传到网站,网址:http://azimuth.satijalab.org/app/azimuth
可以按自己的需要设置质控标准,网页会同步展示小提琴图
质控参数确定后,点击Map cells to reference就可以鉴定细胞类型了。
比对成功后,先点击左边红框标注的Download Results,然后下载预测结果
提取T细胞亚群
## 读取azimuth鉴定的结果
pred <- read.delim("azimuth_pred.tsv")
head(pred, 2)
# cell predicted.id predicted.score mapping.score
#1 HNC01PBMC_AAACCTGAGGAGCGTT-1 CD14 Mono 1.0000000 0.8842671
#2 HNC01PBMC_AAACCTGCACGGACAA-1 CD8 TEM 0.3651781 0.8838794
## 提取只含T细胞的子集
pred.T <- subset(pred, pred$predicted.id %in% c('CD4 CTL',
'CD4 Naive',
'CD4 Proliferating',
'CD4 TCM',
'CD4 TEM',
'CD8 Naive',
'CD8 TCM',
'CD8 Proliferating',
'CD8 TEM',
'MAIT',
'Treg',
'dnT',
'gdT'))
scRNAsub <- scRNA[,as.character(pred.T$cell)]
## T细胞子集去除批次效应
scRNAsub <- SCTransform(scRNAsub) %>% RunPCA()
scRNAsub <- RunHarmony(scRNAsub, group.by.vars = "orig.ident",
assay.use = "SCT", max.iter.harmony = 10)
## 降维聚类
ElbowPlot(scRNA, ndims = 50)
pc.num=1:30
scRNAsub <- RunUMAP(scRNAsub, reduction="harmony", dims=pc.num) %>%
FindNeighbors(reduction="harmony", dims=pc.num) %>%
FindClusters(resolution=0.8)
pred.T <- data.frame(pred.T[, c(2,3,4)], row.names = pred.T$cell)
scRNAsub <- AddMetaData(scRNAsub, metadata = pred.T)
DimPlot(scRNAsub, reduction = "umap", group.by = "predicted.id", label = T) + NoLegend()
monocle3轨迹分析
##创建CDS对象并预处理数据
data <- GetAssayData(scRNAsub, assay = 'RNA', slot = 'counts')
cell_metadata <- scRNAsub@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
#preprocess_cds函数相当于seurat中NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds, num_dim = 50)
#umap降维
cds <- reduce_dimension(cds, preprocess_method = "PCA")
p1 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="celltype") + ggtitle('cds.umap')
##从seurat导入整合过的umap坐标
cds.embed <- cds@int_colData$reducedDims$UMAP
int.embed <- Embeddings(scRNAsub, reduction = "umap")
int.embed <- int.embed[rownames(cds.embed),]
cds@int_colData$reducedDims$UMAP <- int.embed
p2 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="celltype") + ggtitle('int.umap')
monocle降维的结果
导入整合过的umap坐标后作图
## Monocle3聚类分区
cds <- cluster_cells(cds)
p1 <- plot_cells(cds, show_trajectory_graph = FALSE) + ggtitle("label by clusterID")
p2 <- plot_cells(cds, color_cells_by = "partition", show_trajectory_graph = FALSE) +
ggtitle("label by partitionID")
p = wrap_plots(p1, p2)
## 识别轨迹
cds <- learn_graph(cds)
p = plot_cells(cds, label_groups_by_cluster = FALSE, label_leaves = FALSE,
label_branch_points = FALSE)
##细胞按拟时排序
# cds <- order_cells(cds) 存在bug,使用辅助线选择root细胞
p + geom_vline(xintercept = seq(-7,-6,0.25)) + geom_hline(yintercept = seq(0,1,0.25))
embed <- data.frame(Embeddings(scRNAsub, reduction = "umap"))
embed <- subset(embed, UMAP_1 > -6.75 & UMAP_1 < -6.5 & UMAP_2 > 0.24 & UMAP_2 < 0.25)
root.cell <- rownames(embed)
cds <- order_cells(cds, root_cells = root.cell)
plot_cells(cds, color_cells_by = "pseudotime", label_cell_groups = FALSE,
label_leaves = FALSE, label_branch_points = FALSE)
画几条辅助线用于找到root细胞
完成拟时分析的细胞排序结果
monocle3差异分析
##寻找拟时轨迹差异基因
#graph_test分析最重要的结果是莫兰指数(morans_I),其值在-1至1之间,0代表此基因没有
#空间共表达效应,1代表此基因在空间距离相近的细胞中表达值高度相似。
Track_genes <- graph_test(cds, neighbor_graph="principal_graph", cores=10)
#挑选top10画图展示
Track_genes_sig <- Track_genes %>% top_n(n=10, morans_I) %>%
pull(gene_short_name) %>% as.character()
#基因表达趋势图
plot_genes_in_pseudotime(cds[Track_genes_sig,], color_cells_by="predicted.id",
min_expr=0.5, ncol = 2)
#FeaturePlot图
plot_cells(cds, genes=Track_genes_sig, show_trajectory_graph=FALSE,
label_cell_groups=FALSE, label_leaves=FALSE)
##寻找共表达模块
genelist <- pull(Track_genes, gene_short_name) %>% as.character()
gene_module <- find_gene_modules(cds[genelist,], resolution=1e-2, cores = 10)
cell_group <- tibble::tibble(cell=row.names(colData(cds)),
cell_group=colData(cds)$predicted.id)
agg_mat <- aggregate_gene_expression(cds, gene_module, cell_group)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
pheatmap::pheatmap(agg_mat, scale="column", clustering_method="ward.D2")
基因表达趋势图
FeaturePlot图
共表达模块热图
交流探讨:如果您阅读此文有所疑惑,或有不同见解,亦或其他问题,可以点击阅读原文联系探讨。
推荐阅读
-
2021-05-26 HARMONY 和 Seurat 的单细胞分析
-
正负偏差变量 即 d2+、d2- 分别表示决策值中超出和未达到目标值的部分。而 di+、di- 均大于 0 刚性约束和目标约束(柔性目标约束有偏差) 在多目标规划中,>=/<= 在刚性约束中保持不变。当需要将约束条件转换为柔性约束条件时,需要将 >=/<= 更改为 =(因为已经有 d2+、d2- 用来表示正负偏差),并附加上 (+dii-di+) 注意这里是 +di、-di+!之所以是 +di,-di+,是因为需要将目标还原为最接近的原始刚性约束条件 优先级因素和权重因素 对多个目标进行优先排序和优先排序 目标规划的目标函数 是所有偏差变量的加权和。值得注意的是,这个加权和都取最小值。而 di+ 和 dii- 并不一定要出现在每个不同的需求层次中。具体分析需要具体问题具体分析 下面是一个例子: 题目中说设备 B 既要求充分利用,又要求尽可能不加班,那么列出的时间计量表达式即为:min z = P3 (d3- + d3 +) 使用 + 而不是 -d3 + 的原因是:正负偏差不可能同时存在,必须有 di+di=0 (因为判定值不可能同时大于目标值和小于目标值),而前面是 min,所以只要取 + 并让 di+ 和 dii- 都为正值即可。因此,得出以下规则: 最后,给出示例和相应的解法: 问题:某企业生产 A 和 B 两种产品,需要使用 A、B、C 三种设备。下表显示了与工时和设备使用限制有关的产品利润率。问该企业应如何组织生产以实现下列目标? (1) 力争利润目标不低于 1 500 美元; (2) 考虑到市场需求,A、B 两种产品的生产比例应尽量保持在 1:2; (3)设备 A 是贵重设备,严禁超时使用; (4)设备 C 可以适当加班,但要控制;设备 B 要求充分利用,但尽量不加班。 从重要性来看,设备 B 的重要性是设备 C 的三倍。 建立相应的目标规划模型并求解。 解:设企业生产 A、B 两种产品的件数分别为 x1、x2,并建立相应的目标计划模型: 以下为顺序求解法,利用 LINGO 求解: 1 级目标: 模型。 设置。 variable/1..2/:x;! s_con_num/1...4/:g,dplus,dminus;!所需软约束数量(g=dplus=dminus 数量)及相关参数; s_con(s_con_num);! s_con(s_con_num,variable):c;!软约束系数; 结束集 数据。 g=1500 0 16 15. c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(1);!第一个目标函数;!对应于 min=z 的第一小部分;! 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); !使用设置完成的数据构建软约束表达式; ! !软约束表达式 @for(variable:@gin(x)); !将变量约束为整数; ! 结束 此时,第一级目标的最优值为 0,第一级偏差为 0: 第二级目标: !求 dminus(1)=0,然后求解第二级目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(2)+dplus(2);!第二个目标函数 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标结果 @for(variable:@gin(x)); ! 结束 此时,第二个目标的最优值为 0,偏差为 0: 第三目标 !求 dminus(2)=0,然后求解第三个目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=3*dminus(3)+3*dplus(3)+dminus(4);!第三个目标函数。 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标约束条件; ! dminus(2)+dplus(2)=0; !第二个目标约束条件 @for(variable:@gin(x));! 结束 最终结果为 x1=2,x2=4,dplus(1)=100,最优利润为
-
网络在单细胞转录组数据分析中的应用
-
5 分钟知识:白人也能理解的数据包络分析算法
-
最新版 cellphoneDB(v4)解读的 10 倍单细胞空间通信分析
-
实验 5 连续信号的频谱分析
-
5/6 线性系统的频域分析(下)
-
原生 APP 与 H5 的定位偏差?地图定位原理分析
-
scRepertoire||单细胞免疫组化库分析:R (I) 中的应用
-
分析来自 TMDB 的 5,000 条电影数据