velocyto簡單使用

軟件介紹

軟件安裝

• R版本

library(devtools)
install_github("velocyto-team/velocyto.R")

需要安裝hdf5，詳細見hdf5的安裝
然后報錯，需要安裝pcaMethods

BiocManager::install("pcaMethods")

繼續安裝報錯

yum install boost-devel

安裝完成

安裝pagoda2

BiocManager::install(c("AnnotationDbi", "BiocGenerics", "GO.db", "pcaMethods"))
devtools::install_github("hms-dbmi/pagoda2")

安裝依賴

yum install openssl-devel libcurl-devel

• python版本
  python 版本大於3.6

conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install pysam
pip install velocyto

軟件的使用

這里是python版本的使用方法

• 從cellranger得到loom文件，命令如下

velocyto run10x -m repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf

• note
  • Repeat_msk.gft 需要從UCSA網站下載得到：hg38_rmsk.gtf
  • samtools 需要是1.6版本
  • genes.gtf是cellranger用的gtf文件
  • mypath/sample01是10x輸出的文件夾，是cellranger運行后得到的樣品的文件夾，不是使filtered_gene_bc_matrices中的文件夾（只包含barcodes.tsv.gz, features.tsv.gz, matrix.mtx.gz）
    這樣就可以得到loom文件，loom文件被保存在cellranger的sample文件下velocyto文件夾
    官網說這一步大約需要3個小時

合並多個樣品的loom文件

files = ["file1.loom","file2.loom","file3.loom","file4.loom"]
# on the command line do: cp file1.loom merged.loom
ds = loompy.connect("merged.loom")
for fn in files[1:]:
    ds.add_loom(fn, batch_size=1000)

然后使用得到的loom文件進行下游分析
參考博客：http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html
官網也推薦使用的流程
http://velocyto.org/

velocyto代碼

#導入包
library(velocyto.R)

input_loom <- "SCG71.loom"
ldat <- read.loom.matrices(input_loom)

#使用剪切位點的表達量作為輸入
emat <- ldat$spliced
#做直方圖查看數據分步
hist(log10(colSums(emat)),col='wheat',xlab='cell size')

#對數據進行過濾，如果過濾了可以不進行
emat <- emat[,colSums(emat)>=1e3]

# 使用Pagoda2 processing, 導入pagoda2包
# pagoda2用來生成細胞矩陣，對細胞聚類
library(pagoda2)
#讀入數據進行標准化
r <- Pagoda2$new(emat,modelType='plain',trim=10,log.scale=T)
#查看作圖結果
r$adjustVariance(plot=T,do.par=T,gam.k=10)

對細胞進行聚類和細胞嵌合分析

r$calculatePcaReduction(nPcs=100,n.odgenes=3e3,maxit=300)

r$makeKnnGraph(k=30,type='PCA',center=T,distance='cosine')

r$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=T)

作圖展示聚類結果

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') 

准備矩陣和聚類結果

#忽略跨剪切位點的數據
emat <- ldat$spliced
nmat <- ldat$unspliced
#通過p2對細胞進行過濾
emat <- emat[,rownames(r$counts)]
nmat <- nmat[,rownames(r$counts)]
#對分類數據進行標記
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

計算細胞間的距離

cell.dist <- as.dist(1-armaCor(t(r$reductions$PCA)))

基於最小平均表達量篩選基因（至少在一個簇中），輸出產生的有效基因數

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)))

計算RNA速率

fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=fit.quantile)

可視化RNA速率結果

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)

可視化特定的基因

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)

增加k的數目

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)