前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >RNAvelocity4:velocyto.R的使用

RNAvelocity4:velocyto.R的使用

作者头像
生信技能树jimmy
发布2021-09-15 15:42:31
2.5K0
发布2021-09-15 15:42:31
举报
文章被收录于专栏:单细胞天地单细胞天地

RNA注释

这里以inDrop实验数据举例,spliced/unspliced的RNA可以通过:

  1. 使用dropEst 输出管道生成 类似10X的 bam 文件:
代码语言:javascript
复制
~/dropEst/build/dropest -m -F -L eiEIBA -o run1 -g cellranger/refdata-cellranger-mm10-1.2.0/genes/genes.gtf -c ~/dropEst/configs/indrop_v3.xml *.bam
  1. 使用velocyto.py来注释spliced/unspliced的reads,写出标准loom文件:
代码语言:javascript
复制
velocyto run -u Gene -o out -e SCG71 -m mm10_rmsk_srt.gtf -v SCG_71_tophat.filtered.sorted.bam UCSC/mm10/Annotation/Genes/genes.gtf

(注意,也可以使用-V参数直接通过 DropEst 注释spliced/unspliced的reads:)

代码语言:javascript
复制
~/dropEst/dropest -V -C 6000 -m -g ucsc_mm10_exons.gtf.gz -c ~/dropEst/configs/indrop_v3.xml *.aligned.bam)

下面的例子从 velocyto.py 生成的loom文件开始,使用 pagoda2[1] 获取细胞clusters/embedding,然后估计/可视化RNA速率。

加载R包:

代码语言:javascript
复制
library(velocyto.R)

Loading required package: Matrix

数据下载

下载预先计算的loom矩阵

代码语言:javascript
复制
wget http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom

读取矩阵

代码语言:javascript
复制
ldat <-read.loom.matrices(url("http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom"))
#Error: is.character(name) is not TRUE

使用pagoda2标准化和聚类细胞

使用spliced 表达矩阵作为pagoda2的输入,并查看分布。

代码语言:javascript
复制
emat <- ldat$spliced
hist(log10(colSums(emat)),col='wheat',xlab='cell size')
代码语言:javascript
复制
# this dataset has already been pre-filtered, but this is where one woudl do some filtering

emat <- emat[,colSums(emat)>=1e3]

pagoda2处理 pagoda2用于生成细胞嵌入、细胞聚类以及更精确的细胞距离矩阵。您也可以使用其他工具生成这些工具,如 Seurat等。

创建pagoda2对象,调整方差:

代码语言:javascript
复制
library(pagoda2)
r <- Pagoda2$new(emat,modelType='plain',trim=10,log.scale=T)
#2600 cells, 7301 genes; normalizing ... using plain model winsorizing ... log scale ... done.
r$adjustVariance(plot=T,do.par=T,gam.k=10)
#calculating variance fit ... using gam 137 overdispersed genes ... 137 persisting ... done.

运行基本分析步骤以生成细胞嵌入和聚类,并可视化:

代码语言:javascript
复制
r$calculatePcaReduction(nPcs=100,n.odgenes=3e3,maxit=300)

running PCA using 3000 OD genes .. Loading required package: irlba .. done

代码语言:javascript
复制
r$makeKnnGraph(k=30,type='PCA',center=T,distance='cosine');
代码语言:javascript
复制
r$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=T)

绘制嵌入、集群标记(左)和"Xist"表达图(将男性和女性分开展示)

代码语言:javascript
复制
par(mfrow=c(1,2))
r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=10,shuffle.colors=F,mark.cluster.cex=1,alpha=0.3,main='cell clusters')
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$depth,main='depth')  

treating colors as a gradient with zlim: 1000.9 2939

速率估计

准备矩阵和聚类数据:

代码语言:javascript
复制
emat <- ldat$spliced; nmat <- ldat$unspliced; # disregarding spanning reads, as there are too few of them
emat <- emat[,rownames(r$counts)]; nmat <- nmat[,rownames(r$counts)]; # restrict to cells that passed p2 filter
# take cluster labels

cluster.label <- r$clusters$PCA$multilevel # take the cluster factor that was calculated by p2
cell.colors <- pagoda2:::fac2col(cluster.label)

# take embedding form p2

emb <- r$embeddings$PCA$tSNE

除了聚类和 t-SNE 嵌入,从 p2 处理,我们将采取细胞距离,优于默认的R 通常会使用全转录体相关距离。

代码语言:javascript
复制
cell.dist <- as.dist(1-armaCor(t(r$reductions$PCA)))

基于最低平均表达量(至少在其中一个cluster中)筛选基因,输出生成的有效基因总数:

代码语言:javascript
复制
emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.2)
nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
length(intersect(rownames(emat),rownames(nmat)))

[1] 2541

估计RNA速率:

代码语言:javascript
复制
fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=fit.quantile)

可视化 t-SNE 嵌入上的速率:

代码语言:javascript
复制
show.velocity.on.embedding.cor(emb,rvel.cd,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=F,cell.border.alpha = 0.1)

特定基因可视化:

代码语言:javascript
复制
gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 25,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,old.fit=rvel.cd,do.par=T)

调整视图

调整kCells,这会给我们一个更理想化的图像视图,差异肉眼可见:

代码语言:javascript
复制
gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 100,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,do.par=T)

文中链接

[1]

pagoda2: https://github.com/hms-dbmi/pagoda2

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-09-03,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 单细胞天地 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • RNA注释
  • 加载R包:
  • 数据下载
  • 使用pagoda2标准化和聚类细胞
  • 绘制嵌入、集群标记(左)和"Xist"表达图(将男性和女性分开展示)
  • 速率估计
    • 准备矩阵和聚类数据:
      • 估计RNA速率:
        • 可视化 t-SNE 嵌入上的速率:
          • 特定基因可视化:
            • 调整视图
              • 文中链接
              领券
              问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档