单细胞转录组基础分析六:伪时间分析
单细胞测序技术是近年最大的生命科学突破之一,相关文章频繁发表于各大顶级期刊,然而单细胞数据的分析依然是大家普遍面临的障碍。本专题将针对10X Genomics单细胞转录组数据演示各种主流分析,包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。不足之处请大家批评指正,欢迎添加Kinesin微信交流探讨!
Monocle进行伪时间分析的核心技术是一种机器学习算法——反向图形嵌入 (Reversed Graph Embedding)。它分析的前提需要一张展现细胞转录特征相似性关系的图,Monocle2使用DDTree降维图,Monocle3使用UMAP降维图。Monocle的机器学习算法可以依据上述降维图形,学习描述细胞如何从一种状态过渡到另一种状态的轨迹。Monocle假设轨迹是树状结构,一端是“根”,另一端是“叶”。一个细胞在生物过程的开始,从根开始沿着主干进行,直到它到达第一个分支。然后,该细胞必须选择一条路径,并沿着树移动越来越远,直到它到达一片叶子。一个细胞的假时间值是它返回根所需的距离。降维方面monocle与seurat的过程大同小异,首先进行数据标准化,其次选择部分基因代表细胞转录特征 ,最后选用适当的算法降维。对Monocle原理感兴趣的同学可以登录官网查看:
http://cole-trapnell-lab.github.io/monocle-release/
数据导入与处理
轨迹分析的前提是待分析的细胞有紧密的发育关系,PBMC细胞不是很好的的示例数据,我们选择T细胞群体演示一下。Monocle建议导入原始表达矩阵,由它完成数据标准化和其他预处理。
dir.create("pseudotime")
scRNAsub <- readRDS("scRNAsub.rds") #scRNAsub是上一节保存的T细胞子集seurat对象
data <- as(as.matrix(scRNAsub@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = scRNAsub@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
mycds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
expressionFamily = negbinomial.size())
expressionFamily参数用于指定表达矩阵的数据类型,有几个选项可以选择:
- 稀疏矩阵用negbinomial.size(),
- FPKM值用tobit(),
- logFPKM值用gaussianff()
mycds是Monocle为我们的数据生成的对象,相当于我们在seurat使用的scRNA对象。数据导入后需要进行标准化和其他预处理:
mycds <- estimateSizeFactors(mycds)
mycds <- estimateDispersions(mycds, cores=4, relative_expr = TRUE)
#mycds <- detectGenes(mycds, min_expr = 2) #很多教程不用
与seurat把标准化后的表达矩阵保存在对象中不同,monocle只保存一些中间结果在对象中,需要用时再用这些中间结果转化。经过上面三个函数的计算,mycds对象中多了SizeFactors、Dipersions、num_cells_expressed和num_genes_expressed等信息。
提选择代表性基因
完成数据导入和预处理后,就可以考虑选择哪些基因代表细胞的发育特征,Monocle官网教程提供了4个选择方法:
- 选择发育差异表达基因
- 选择clusters差异表达基因
- 选择离散程度高的基因
- 自定义发育marker基因
前三种都是无监督分析方法,细胞发育轨迹生成完全不受人工干预;最后一种是半监督分析方法,可以使用先验知识辅助分析。第一种方法要求实验设计有不同的时间点,对起点和终点的样本做基因表达差异分析,挑选显著差异的基因进行后续分析。对于没有时序设计的实验样本,可以使用第2、3种方法挑选基因。第2种方法要先对细胞降维聚类,然后用clusters之间差异表达的基因开展后续分析。Monocle有一套自己的降维聚类方法,与seurat的方法大同小异,很多教程直接使用seurat的差异分析结果。第3种方法使用离散程度高的基因开展分析,seurat有挑选高变基因的方法,monocle也有自己选择的算法。本案例数据不具备使用第1、4种方法的条件,因此这里只演示2、3种方法的使用。
##使用clusters差异表达基因
diff.genes <- read.csv('subcluster/diff_genes_wilcox.csv')
diff.genes <- subset(diff.genes,p_val_adj<0.01)$gene
mycds <- setOrderingFilter(mycds, diff.genes)
p1 <- plot_ordering_genes(mycds)
##使用seurat选择的高变基因
var.genes <- VariableFeatures(scRNAsub)
mycds <- setOrderingFilter(mycds, var.genes)
p2 <- plot_ordering_genes(mycds)
##使用monocle选择的高变基因
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
mycds <- setOrderingFilter(mycds, disp.genes)
p3 <- plot_ordering_genes(mycds)
##结果对比
p1|p2|p3
选择不同的基因集,拟时分析的结果不同,实践中可以几种方法都试一下。
降维及细胞排序
使用disp.genes开展后续分析
#降维
mycds <- reduceDimension(mycds, max_components = 2, method = 'DDRTree')
#排序
mycds <- orderCells(mycds)
#State轨迹分布图
plot1 <- plot_cell_trajectory(mycds, color_by = "State")
ggsave("pseudotime/State.pdf", plot = plot1, width = 6, height = 5)
ggsave("pseudotime/State.png", plot = plot1, width = 6, height = 5)
##Cluster轨迹分布图
plot2 <- plot_cell_trajectory(mycds, color_by = "seurat_clusters")
ggsave("pseudotime/Cluster.pdf", plot = plot2, width = 6, height = 5)
ggsave("pseudotime/Cluster.png", plot = plot2, width = 6, height = 5)
##Pseudotime轨迹图
plot3 <- plot_cell_trajectory(mycds, color_by = "Pseudotime")
ggsave("pseudotime/Pseudotime.pdf", plot = plot3, width = 6, height = 5)
ggsave("pseudotime/Pseudotime.png", plot = plot3, width = 6, height = 5)
##合并作图
plotc <- plot1|plot2|plot3
ggsave("pseudotime/Combination.pdf", plot = plotc, width = 10, height = 3.5)
ggsave("pseudotime/Combination.png", plot = plotc, width = 10, height = 3.5)
##保存结果
write.csv(pData(mycds), "pseudotime/pseudotime.csv")
使用diff.genes分析的结果
轨迹图分面显示
p1 <- plot_cell_trajectory(mycds, color_by = "State") + facet_wrap(~State, nrow = 1)
p2 <- plot_cell_trajectory(mycds, color_by = "seurat_clusters") + facet_wrap(~seurat_clusters, nrow = 1)
plotc <- p1/p2
ggsave("pseudotime/trajectory_facet.png", plot = plotc, width = 6, height = 5)
Monocle基因可视化
s.genes <- c("ITGB1","CCR7","KLRB1","GNLY")
p1 <- plot_genes_jitter(mycds[s.genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(mycds[s.genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(mycds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
ggsave("pseudotime/genes_visual.png", plot = plotc, width = 8, height = 4.5)
拟时相关基因聚类热图
Monocle中differentialGeneTest()函数可以按条件进行差异分析,将相关参数设为fullModelFormulaStr = "~sm.ns(Pseudotime)"时,可以找到与拟时相关的差异基因。我们可以按一定的条件筛选基因后进行差异分析,全部基因都输入会耗费比较长的时间。建议使用cluster差异基因或高变基因输入函数计算。分析结果主要依据qval区分差异的显著性,筛选之后可以用plot_pseudotime_heatmap函数绘制成热图。
#cluster差异基因
diff.genes <- read.csv('subcluster/diff_genes_wilcox.csv')
sig_diff.genes <- subset(diff.genes,p_val_adj<0.0001&abs(avg_logFC)>0.75)$gene
sig_diff.genes <- unique(as.character(sig_diff.genes))
diff_test <- differentialGeneTest(mycds[sig_diff.genes,], cores = 1,
fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 0.01))
p1 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=3,
show_rownames=T, return_heatmap=T)
ggsave("pseudotime/pseudotime_heatmap1.png", plot = p1, width = 5, height = 8)
#高变基因
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
diff_test <- differentialGeneTest(mycds[disp.genes,], cores = 4,
fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 1e-04))
p2 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=5,
show_rownames=T, return_heatmap=T)
ggsave("pseudotime/pseudotime_heatmap2.png", plot = p2, width = 5, height = 10)
BEAM分析
单细胞轨迹中通常包括分支,它们的出现是因为细胞的表达模式不同。当细胞做出命运选择时,或者遗传、化学或环境扰动时,就会表现出不同的基因表达模式。BEAM(Branched expression analysis modeling)是一种统计方法,用于寻找以依赖于分支的方式调控的基因。
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
mycds_sub <- mycds[disp.genes,]
plot_cell_trajectory(mycds_sub, color_by = "State")
beam_res <- BEAM(mycds_sub, branch_point = 1, cores = 8)
beam_res <- beam_res[order(beam_res$qval),]
beam_res <- beam_res[,c("gene_short_name", "pval", "qval")]
mycds_sub_beam <- mycds_sub[row.names(subset(beam_res, qval < 1e-4)),]
plot_genes_branched_heatmap(mycds_sub_beam, branch_point = 1, num_clusters = 3, show_rownames = T)
获取帮助
本教程的目的在于把常用的单细胞分析流程串起来,适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释,官网有详细的教程和权威的解释,翻译好的中文教程也有多个版本,有兴趣的可以搜索一下。如果您学习本教程有一定困难,可以点击 “阅读原文” 找到作者联系方式,获取帮助。
往期回顾
- 使用mitmproxy嗅探双向认证ssl链接——嗅探AWS IoT SDK的mqtts
- Django中ORM介绍和字段及其参数
- 几个实现分页的方法
- 如何识别IDA反汇编中动态链接库中的函数
- Kali Linux安装TL-WN821N USB无线网卡驱动(make失败)
- Cookie、Session登陆验证相关介绍和用法
- 题型分析
- .NET CORE 框架ABP的代码生成器(ABP Code Power Tools )使用说明文档
- Lua编写wireshark插件初探——解析Websocket上的MQTT协议
- 在Ubuntu 16.04环境下安装Docker-CE(附视频教程)
- 修改HTML5 input placeholder 颜色及修改失效的解决办法
- 设置同样字体大小,chrome浏览器有时字体偏大的解决办法(转)
- 手机端调用系统相册并上传图片
- select自定义小三角样式
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- leetcode栈之用两个栈实现队列
- Ui Automator 框架和Ui Automator Viewer你会用吗?附送「必备adb命令」拿走不谢 !
- 3分钟短文:Laravel模型读数据的那个“障眼法”
- 你有一份面试题要查收
- @RequestParam,@RequestBody,@PathVariable注解区别
- linux centos 查看cpu是否开启超线程虚拟化和配置信息
- DDD领域驱动设计实战(三)- 理解实体
- 如何理解 Go 中的反射
- Kubernetes 1.19.0——网络
- 如何用VSCode实现一个vue.js项目?
- 网页轮播图案例
- leetcode栈之用队列实现栈
- 3分钟短文:无私分享!Laravel模型使用2个小技巧
- 抖音小程序基础之 如何绑定点击事件(教程含源码)
- 抖音小程序基础之 如何使用内联样式 style和class(教程含源码)