万字长文介绍单细胞高级分析!学会这个分析,搞定单细胞套路80%的难题!
大家好,我是晨曦,有了前面知识的铺垫,相信各位小伙伴对于scRNA-seq已经不是懵懂无知的小白啦毕竟现在都在提倡走出舒适圈,那么这次我们就开始我们scRNA-seq的第一个高阶分析的讲解说是高阶分析,其实在scRNA-seq不断普及的今天,在高级的分析,放在全国乃至全球,也会逐渐沦为普通,但是虽然使用的是越来越普遍,但是其难度却没有一丝一毫的降低,我们仍然要有着一颗学徒的心,踏踏实实的完成我们scRNA-seq的学习
晨曦讲解:本教程参考了网上周老师的教程+monocle3官网教程+网上各种教程
高能万字警告!!!!!!
当我们要进行一种分析时,首先我们必须要知道这个分析是干什么的,接下来将通过几个问题,为我们进行拟时序分析前的一些基本知识的讲解
机体在外界调节的“刺激”下,会发生一系列的变化,这些变化归根到底其实都是细胞功能和结构的变化,为了响应这些外界刺激,其细胞会从一种功能“状态”转变为另一种功能“状态”,而这些功能和结构的转变,往往会经历基因表达的升高或沉默,毕竟分子表达差异决定了表型的差异
但是这类细微的变化,往往很难在体外捕捉到,这时候我们就引入了拟时序分析,或者可以称为细胞轨迹分析,通过细胞轨迹分析可以推断出细胞的分化轨迹或细胞亚型的演化过程,主要原理是基于关键基因的表达模式,在拟时间中对单个细胞进行排序,模拟出时间发育过程的动态变化
晨曦解读
原理层面的东西了解个大概即可, 毕竟我们不是专门从事算法研究的,概括一下 ,拟时序分析可以帮助我们构建细胞的发育轨迹,研究在外界刺激下,细胞的动态演变过程,这里的伪时间是一个抽象的分化单位:它只是一个细胞沿着轨迹起点到终点的最短距离,轨迹的总长度是由细胞从起始状态移动到终点状态所经历的总转录变化量来定义的
需要注意几个关键点:
- 只有细胞本身存在随时间变化的特性,才可以做拟时序分析,所以一般用在发育组织的分析中
- 细胞轨迹分析其实也可以看作是差异分析的细化,差异分析只是告诉我们差异,但是不会告诉我们变化,或者可以说是如何变化,而这些细胞轨迹分析会告诉我们
- 这里借鉴周老师对于拟时序分析的看法(个人觉得很有共鸣)
- 拟时序分析主要基于关键基因的表达模式,在拟时间中对单个细胞进行排序,模拟出时间发育过程的动态变化
这里涉及到拟时序分析本身就是一种排序而已,既然是排序,就要知道排序的基本要素
1.对什么排序→对象
2.如何判断先后顺序→顺序
3.如何寻找分支点(如果存在分支的话)→分支
在这里首先我们要知道,排序技术是一种在低维空间排布高维数据的降维技术,所以排序离不开降维;降维离不开特征提取(或者选择)
至此可以概括出拟时序分析的三个主要步骤
1.基因筛选(质控)(如果从Seurat对象导入则可以省略这一步)
2.数据降维(monocle3)
3.细胞排序
这里就不得不提及一篇文献:Reversed graph embedding resolves complex single-cell developmental trajectories
通过这篇文献,我们可以引入我们的工具,也是我们接下来的主角——momocle3
Monocle引入了在伪时间(拟时间)内对单个细胞排序的策略,利用单个细胞的非同步进程,将它们置于与细胞分化等生物学过程相对应的轨迹上,然后利用机器学习技术(反向图嵌入)进行排序
目前已经经历了三代的更新,我们后续的讲解都是以最新的版本来进行讲解
晨曦解读
如果进行了拟时序分析,希望引用下面的文献,当然,这些文献也是我们最好的学习资料
1.The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells
2.Single-cell mRNA quantification and differential analysis with Census
3.Reversed graph embedding resolves complex single-cell trajectories
4.The single-cell transcriptional landscape of mammalian organogenesis
Monocle3功能介绍
R包:Monocle3的功能介绍如下
晨曦解读
功能一:执行常规的scRNA-seq的标准流程(Seurat包)
功能二:拟时序分析(将细胞按照基因表达相似水平排序)(我们重点关注的功能)
功能三:差异表达分析(分析哪些基因会随着拟时序变化,这些基因可能就是我们研究的重点)
Monocle3代码操作
晨曦解读
针对拟时序分析,首先是因为其原理确实是十分的复杂,而且学习完原理后,你走的依旧是较为标准化的流程,所以这里我们在介绍完一些基本的背景知识后,就开始直接进行代码操作的部分,也是为了让大家可以更快的掌握这一门高阶分析的技术
#准备工作
library(monocle3)
#这个包有很多依赖包,下载也不是十分的容易,如果下载不下来建议参考下面的代码
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
'limma', 'S4Vectors', 'SingleCellExperiment',
'SummarizedExperiment', 'batchelor', 'Matrix.utils'))
install.packages("devtools")
devtools::install_github('cole-trapnell-lab/leidenbase')
devtools::install_github('cole-trapnell-lab/monocle3')
晨曦解读
这里粘贴一下我下载这个包的错误,当然这些错误在monocle3官网都有解决办法~
#接下来是数据导入
#monocle3在官网提供了一些读取数据的方式,但是因为我们在scRNA-seq的时候最常使用Seurat包,所以接下来我们将重点展示如何从Seurat对象中提取特定数据,构建我们进行细胞轨迹的CDS对象
#下面我们还是先从官网标准代码入门(该步骤很慢)
expression_matrix<- readRDS(url("http://staff.washington.edu/hpliner/data/packer_embryo_expression.rds"))
cell_metadata<- readRDS(url("http://staff.washington.edu/hpliner/data/packer_embryo_colData.rds"))
//staff.washington.edu/hpliner/data/packer_embryo_rowData.rds")) :
cds<- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
晨曦解读
可以看出来,我们进行细胞轨迹分析的时候,是需要构建CDS对象的,而这个对象构建,需要三个文件即表达矩阵、细胞信息、基因信息
这里是不是特别像我们构建Seurat对象的步骤,也就是说如果我们只有counts矩阵的时候,就可以采用官网这种方式来构建CDS对象
#从Seurat对象中提取数据构建CDS对象
#首先加载完成基本流程的Seurat对象(这里自己走常规标准流程即可)
scRNA <- readRDS("~/scRNA.rds")
#这里我们首先提取我们感兴趣的细胞亚群,根据Seurat的降维聚类结果,我们选择T细胞,来作为我们后续分析的重点,所以我们需要把T细胞亚群从Seurat种提取出来
T <- subset(scRNA_final, scRNA_final@active.ident %in% c('NK',
'CD8 T',
'Naive CD4 T',
'Memory CD4 T'))
晨曦解读
如果各位小伙伴忘记scRNA-seq标准流程,建议阅读下面推文后再阅读本期推文教程
#提取表达矩阵
data <- as(as.matrix(T@assays$RNA@counts), 'sparseMatrix')
#sparseMatrix代表稀疏矩阵
#提取细胞信息
pd <- T@meta.data
#构建基因信息
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
晨曦解读
首先我们来解释第一个可能存在疑惑的问题
1.为什么选择原始表达矩阵?
Monocle内部有内置的标准化步骤,所以这里的表达矩阵推荐使用原始数据,如果自己进行对数据的标准化等操作,可能会影响Monocle的分析
2.从Seurat对象中提取的这三个数据都是什么样的样子?
data数据
pd数据
fData
晨曦解读
从这里我们可以看出来,我们需要的三个数据都是什么样的格式,其实概括来讲还是表达矩阵、细胞信息、基因信息
#构建CDS对象
cds<- new_cell_data_set(data,
cell_metadata = pd,
gene_metadata = fData)
晨曦解读
然后我们来看一下这个CDS对象,下面的这些信息不需要全部了解,将挑选重点的再后面逐步进行介绍
#预处理数据
#monocle3中的轨迹分析之前都需要各种规范化的预处理步骤
#如果数据之间彼此有批次效应,你需要自己进行矫正,也就是单细胞多数据集的整合
cds <- preprocess_cds(cds, num_dim = 50)
#这里的preprocess_cds函数其实本质上相当于Seurat中的NormalizeData+RunPCA+ScaleData
#这里针对RNA-seq数据使用的是PCA
plot_pc_variance_explained(cds)
晨曦解读
上面介绍了preprocess_cds函数相当于Seurat中的三步骤结合,其中这个函数中的方法有三种,分别是PCA or LSI. For LS,RNA-seq数据选择PCA,如果是ATAC-seq则使用LSI
#降维可视化
cds <- reduce_dimension(cds)#默认为UMAP
plot_cells(cds)
晨曦解读
这个图的颜色肯定是不符合我们的要求的,我们现在来探索一下这个图的颜色,既然要为这个图添加颜色,就需要找到可以映射颜色的数据
#通过阅读帮助文档得color_cells_by参数可以为UMAP图添加颜色,添加颜色的数据可以是以下任意
colnames(colData(cds))
#那么我们接下来选择通过用Seurat内的细胞聚类来进行颜色的标注
p1 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters",group_label_size=4)
晨曦解读
上图是通过monocle3包进行降维聚类的结果,对比一下seurat对象的结果,我们可以看出来,虽然图像不一样,但是大致的分群还是类似的
#从Seurat对象导入UMAP坐标,并通过monocle3绘图函数进行绘制
cds.embed<- cds@int_colData$reducedDims$UMAP#获取CDS对象中umap的坐标
int.embed <- Embeddings(scRNA, reduction = "umap")#获取Seurat对象中umap的坐标
#这里之所以不可以拿Seurat对象中的UMAP是因为细胞ID的数量不同,毕竟是两种降维聚类方法,所以我们要提取坐标中细胞ID一致的,然后绘图,才可以更好比较结果
int.embed <- int.embed[rownames(cds.embed),]#取交集
cds@int_colData$reducedDims$UMAP <- int.embed#赋值
p2 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters")
#最好别忘了切换回去,也就是运行下面这个代码
#cds@int_colData$reducedDims$UMAP <- cds.embed
#还有就是后续我们也可以使用monocle3进行细胞注释,但是这一步我们建议是在Seurat中完成,在monocle3完成聚类就可以了~
head(rownames(cds))
plot_cells(cds, genes=c("AL627309.1", "NOC2L", "RP11-206L10.2", "LINC00115"))
#当然也支持tSNE
cds = reduce_dimension(cds, reduction_method="tSNE")
plot_cells(cds, reduction_method="tSNE", color_cells_by="seurat_clusters"
晨曦解读
但是看起来针对monocle3来说,还是UMAP看起来好一些~
#后续的细胞注释包括marker基因的选择+细胞注释,这个流程就在Seurat中完成即可,我们只需要调取Seurat对象中的结果
ccds@colData@listData[["seurat_clusters"]] <- [email protected]
p4 <- plot_cells(cds, reduction_method="UMAP", color_cells_by="seurat_clusters",group_label_size=4)
晨曦解读
对于多数据集联合主张还是先进行批次效应的去除,整合完数据集后再进行细胞轨迹分析,而且根据作者建议,尽管在Seurat包中已经进行标准化,在进行monocle3分析的时候仍然需要再次进行标准化,所以得到一个初步的结论,如果单个数据集使用counts数据,多数据集则使用整合后的数据,然后后续就都是标准流程
plot_cells(cds, color_cells_by="seurat_clusters", label_cell_groups=FALSE)
晨曦解读
一般来说,对于多数据集,如果亚群的分组存在明显的分离,则说明存在批次效应,这个图看上去没有,毕竟只有一个数据集
#cluster your cells
#这里的聚类已经不再是单单的聚类,按照作者的说法,是把细胞进行分区,后去的轨迹分析都会针对这个分区来进行
#我们其实是跳过了monocle3中的寻找marker基因-细胞注释的环节,因为我们认为在Seurat中也可以完成同样的操作
cds <- cluster_cells(cds)
#轨迹学习
cds<- learn_graph(cds)
p = plot_cells(cds, label_groups_by_cluster = FALSE, label_leaves = FALSE,
label_branch_points = FALSE)
p
晨曦解读
接下来顺步进行即可~
#可视化修饰
p5 <-plot_cells(cds,
color_cells_by = "seurat_clusters",
label_groups_by_cluster=FALSE,
label_branch_points=TRUE,
group_label_size=4)
p5
#保存一下
ggsave(" trajectorygraph.pdf", plot = p,width = 8,height = 6)
晨曦解读
调节的具体参数还可以看plot_cells函数的帮助文档,讲解的很全面~
作者在官网建议:上面这个图将被用于许多下游分析,比如分支分析和差异表达分析
#这里插入一个官网的轨迹学习的可视化,目的是解释一下图上的内容
plot_cells(cds,
color_cells_by = "embryo.time.bin",
label_cell_groups=FALSE,
label_leaves=TRUE,
label_branch_points=TRUE,
graph_label_size=1.5)
晨曦解读
数字黑色圈代表节点即分叉点,黑色白色圈代表各种可能的结局,黑色的线则代表轨迹学习的路线,这些数据的可以通过参数进行设置(label_leaves和label_branch_points)
进行完轨迹学习后,根据官网的指导,下面就需要进行细胞的排序,也就是Order the cells in pseudotime这一个步骤
#Order the cells in pseudotime
#其实概括来讲就是我们需要告诉monocle3,哪里是轨迹分析的起点
#原理别进行深究,因为......你懂的
cds<- order_cells(cds)
#弹出一个交互式界面,选择轨迹分析的起点,可以选择多个点,然后进行可视化
plot_cells(cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5)
#如果想要显示细胞亚群标签,调节里面的参数即可~
晨曦解读
我们刚才选择的点以及和其有关被标注了出来,灰色部分代表和我们所选择的点没有发育或者分化关系,所以这些点用灰色表示
#当然我们也可以通过编程来帮助我们进行roots的选择
get_earliest_principal_node <- function(cds, time_bin="CD8 T"){
cell_ids <- which(colData(cds)[, "seurat_clusters"] == time_bin)
closest_vertex <-cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
cds = order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))
#上述代码我们需要改动的地方有两处,一处是time_bin参数,我们需要指定我们感兴趣的亚群名称,第二处我们需要指定细胞亚群的分组名称,这里我们输入的是seurat_clusters
#用手点可能更直观一点
晨曦解读
重点就是细胞轨迹的可视化没有标注出真正的起始点在哪里,也就是说,算法无法明确你的组织究竟是从A到B,还是B到A,需要凭借个人的生物学背景来判断,比如常规来说都是干细胞分化成各种细胞反过来是不行的
#当然也可以选择分支
cds_sub<- choose_graph_segments(cds)
#还可以做酷炫3D图
cds_3d<- reduce_dimension(cds, max_components = 3)
cds_3d<- cluster_cells(cds_3d)
#cds_3d <- learn_graph(cds_3d)
#cds_3d <- order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))
cds_3d_plot_obj<- plot_cells_3d(cds_3d, color_cells_by="seurat_clusters")
晨曦解读
实际上我只想得到上面这个3D版本的聚类图,所以我注释掉了官网中间的两行代码~
#接下来开始探讨差异表达即细胞轨迹中细胞之间的差异基因
晨曦解读
如果只是单纯针对轨迹分析,想要寻找到细胞与细胞之间分化和发育的变化究竟是由哪些基因的表达改变造成的,选择第二种即可~
第一种方法主要探讨的是不同刺激下,比如时间以及治疗措施等等~
#这里我们主要探讨第二种方法
#首先官网教程给的是选择一个分支即可,这样既可以缩短运行的时间,同时也可以提高效率,但是这里我们针对整个CDS对象来使用
subset_pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=4)
#运行完以后得到这下面这个数据框
晨曦解读
这里面比较有意思的数值就是morans_I(空间共表达),其数值越靠近1代表这个基因在空间距离相近的细胞中表达值越相似,0则代表没有空间共表达效应
#下面我们进行绘制单个基因在不同亚群之间的表达量变化
#官网给出的是指定基因
#我们这里就按照morans_I的大小来选择前十个基因
top10<- subset_pr_test_res %>% top_n(n=10, morans_I) %>%
%>% as.character()
#绘图
p<- plot_genes_in_pseudotime(cds[top10,], color_cells_by="seurat_clusters",
min_expr=0.5, ncol = 2)
晨曦解读
这个图具体看其实就是看连线,这几个基因会发现大部分都是直线,说明这些基因的表达在伪时间中其实并没有太多的变化,CCL5基因变化的比较有意思,呈现波动感~
当然还有一些基因在一开始表达很强烈,但是随着时间推移表达量逐渐下降
plot_cells(cds, genes=top10, show_trajectory_graph=FALSE,
label_cell_groups=FALSE, label_leaves=FALSE)
#这里官网还有一个筛选共表达模块,但是我觉得用处有限,这里就不展示了,感兴趣的同学可以访问下面官方教程
#https://cole-trapnell-lab.github.io/monocle3
#把我们细胞轨迹分析的结果返回Seurat对象中
pseudotime <- pseudotime(cds, reduction_method = 'UMAP')
pseudotime <- pseudotime[rownames(T@meta.data)]
T$pseudotime <- pseudotime
saveRDS(T, file = "T_pseudotime.rds")
晨曦解读
细胞轨迹分析就结束啦~
整体来说,这套分析流程的难点我们主要可以分成两大部分,首先光是安装R包其实就劝退了很多同学,所以基本上来说,各位同学需要保证自己的R包可以安装成功
其次就是针对代码中的每一个步骤进行理解,当然可能还会有同学问,图中的可视化为什么是这个样子,为什么没有像文献那样的可视化展示,其实首先如果我们获得了细胞轨迹分析的数据后,只需要运用我们前面所教的“扒图”技能,就可以充分get网上只要有的细胞轨迹可视化的结果,其本质上依旧是ggplot2以及monocle3的运用
那么我们这期推文到这里就结束啦,近期收到了很多小伙伴的提问,也是陆陆续续的在开始进行相关推文的书写,也欢迎各位小伙伴在评论区留言,说出你的困惑,互相讨论,共同提高
我是晨曦,我们下次再见~
PS:回复“细胞轨迹”可以获得推文中的示例数据和代码哦~-
推文内的代码更多的是教学,所以扩展了很多,大家后台获取的代码是纯分析代码,简化了大家的copy难度~
晨曦单细胞文献阅读系列
非肿瘤单细胞分析模板已到位!眼馋单细胞的小伙伴快来看!手把手教你产出第一篇单细胞SCI!
晨曦碎碎念系列传送门(未完待续...)
晨曦单细胞笔记系列传送门
晨曦从零开始学画图系列传送门
晨曦单细胞数据库系列传送门
—
END—
撰文丨晨 曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
关键词
数据集
生信
差异
过程
R包
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。