单细胞测序技术的发展日新月异,新的分析工具也层出不穷。每个工具都有它的优势与不足,在没有权威工具和流程的单细胞生信江湖里,多掌握几种分析方法和工具,探索数据时常常会有意想不到的惊喜。
单细胞初级8讲和高级分析8讲 单细胞分析十八般武艺1:harmony 单细胞分析十八般武艺2:LIGER 单细胞分析十八般武艺3:fastMNN
大家第一次听到RNA速率/速度,是不是感觉和我一样懵呢?再看文献提到的时间导数、微分方程和一大堆公式,是不是更搞不清作者在讲什么?我们常用的基因表达矩阵,如velocyto的作者在文献中所讲,反映的是细胞转录组的瞬间快照。其实在生命体中,每时每刻都发生着mRNA的转录、剪接和降解,它们当然都是有速度的,作者在文中用α、β和γ表示这些速度。不仅mRNA的转录、剪接和降解有速度,前体mRNA和成熟mRNA的丰度变化也有速度,作者所讲的RNA速度特指成熟mRNA的丰度变化速度。RNA速度可以通过前体mRNA和成熟mRNA的指标来估计,并且可以用来预测细胞未来数小时内的状态。在velocyto的动力学模型中,作者假定转录速度 α恒定,此时unspliced mRNA(u)和未来的spliced mRNA(s)的丰度高度相关:
系数 λ 综合了降解率和剪接率、基因调控特性、内含子和外显子长度的比值以及内部启动位点的数量等因素。绝大多数的基因在不同细胞类型中,系数 λ 固定不变,但是11%的基因在不同组织中 λ 值不同。为了证明此模型的正确性,作者使用小鼠肝脏的节律基因来验证:
总的来说:unspliced mRNA丰度与spliced mRNA丰度高度相关且具有先导性,因此可以据此推断细胞的分化方向与速率。
本文演示数据来自velocyto官方教程:http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html
library(Seurat)
library(velocyto.R)
library(tidyverse)
library(SeuratWrappers)
##数据基础分析
# 读取loom文件
velo <- read.loom.matrices(url("http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom"))
# 转换为seurat对象
velo <- as.Seurat(x = velo)
# 降维聚类
velo <- velo %>% SCTransform(assay="spliced") %>% RunPCA(verbose=F)
ElbowPlot(velo, ndims = 50)
nPC=1:20
velo <- FindNeighbors(velo, dims = nPC) %>% FindClusters() %>%
RunUMAP(dims = nPC) %>% RunTSNE(dims = nPC)
##给细胞分配颜色
ident.colors <- (scales::hue_pal())(n = length(x = levels(x = velo)))
names(x = ident.colors) <- levels(x = velo)
cell.colors <- ident.colors[Idents(object = velo)]
names(x = cell.colors) <- colnames(x = velo)
##速率分析
velo <- RunVelocity(velo, deltaT = 1, kCells = 25, fit.quantile = 0.02,
spliced.average = 0.2, unspliced.average = 0.05, ncores = 18)
#kCells:用于斜率平滑度计算最近邻细胞数量,越大越平滑,越小越能反映局部特征
#fit.quantile:0.02代表对基因表达量最高2%与最低2%的值执行gamma拟合
#spliced.average:过滤低表达丰度基因的标准,计算的是基因在cluster内的平均counts值
#unspliced.average:同上
##全局速率可视化
emb = Embeddings(velo, reduction = "umap")
vel = Tool(velo, slot = "RunVelocity")
show.velocity.on.embedding.cor(emb = emb, vel = vel, n = 200, scale = "sqrt",
cell.colors = ac(cell.colors, alpha = 0.5), cex = 0.8, arrow.scale = 3,
show.grid.flow = TRUE, min.grid.cell.mass = 0.5, grid.n = 40,
arrow.lwd = 1, do.par = FALSE, cell.border.alpha = 0.1)
##特定基因速率可视化
gene = "Camp"
RunVelocity(velo, deltaT=1, kCells=25, fit.quantile=0.02, old.fit=vel,
cell.emb=emb, cell.colors=cell.colors, show.gene=gene, do.par=T)
交流探讨:如果您阅读此文有所疑惑,或有不同见解,亦或其他问题,可以点击阅读原文联系探讨。
往期回顾
OSCA单细胞数据分析笔记-4 Overview pipeline
如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程