最近學(xué)習(xí)畫(huà)單細(xì)胞的擬時(shí)序軌跡圖有梆,有一些收獲音五。大家經(jīng)常使用的一般是monocle亿汞,DPT等;RNA velocity是一個(gè)比較復(fù)雜的方法杆兵,它通過(guò)細(xì)胞轉(zhuǎn)錄時(shí)細(xì)胞內(nèi)成熟mRNA和mRNA前體的含量來(lái)計(jì)算細(xì)胞接下來(lái)生長(zhǎng)分化的趨勢(shì)。也就是說(shuō)仔夺,它在計(jì)算細(xì)胞擬時(shí)序時(shí)考慮了細(xì)胞內(nèi)在的結(jié)構(gòu)琐脏,而非通過(guò)細(xì)胞中基因的表達(dá)量。所以它需要用到原始的bam文件缸兔。我們今天來(lái)解析一下這個(gè)方法的安裝和使用日裙。
一、安裝指導(dǎo)
http://velocyto.org/velocyto.py/install/index.html
1惰蜜、python環(huán)境下 velocyto.py
包括:
- 一個(gè)命令行工具昂拂,用于從mapping的bam文件中得到的loom文件【主要是得到spliced mRNA count矩陣和unspliced mRNA count矩陣】
- 一個(gè)分析工具,用于對(duì)loom文件做后續(xù)的RNA velocity分析
代碼:
通過(guò)pip或者conda安裝(如果出現(xiàn)錯(cuò)誤類型:ImportError: dlopen: cannot load any more object with static TLS抛猖;有可能是導(dǎo)入包的順序的問(wèn)題格侯,嘗試退出,重新進(jìn)入财著,然后首先導(dǎo)入sklearn包):
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install -U --no-deps velocyto
#或者:
git clone https://github.com/velocyto-team/velocyto.py.git
cd velocyto.py
pip install -e . #注意后面的這個(gè)點(diǎn)
veocyto的安裝成功與否和python版本有關(guān)联四,必須3.6.0以上的版本。R語(yǔ)言的velocyto不能獨(dú)立使用撑教,必須在python版本下生成loom文件朝墩,才可以使用。
2伟姐、R環(huán)境下 velocyto.R
包括:
R版本的velocyto沒(méi)有辦法根據(jù)bam文件計(jì)算得到成熟mRNA的spliced count矩陣和mRNA前體的unspliced count矩陣的loom文件收苏。但是據(jù)說(shuō)可以使用dropEst計(jì)算得到loom文件亿卤。
[不過(guò)本人并沒(méi)有嘗試過(guò),如果今后有機(jī)會(huì)嘗試成功了鹿霸,一定回來(lái)告訴大家排吴。]
- 一個(gè)分析工具,用于對(duì)loom文件做后續(xù)的RNA velocity分析
代碼:
library(devtools)
install_github("velocyto-team/velocyto.R")
且需要安裝其他的包作為輔助:
library(Seurat)
library(tidyverse)
library(SeuratWrappers)
3杜跷、scvelo
scVelo傍念,是Volker Bergen團(tuán)隊(duì)2020年發(fā)表在NBT上的一個(gè)計(jì)算RNA velocity的工具,和velocyto相比葛闷,它的可視化結(jié)果更易解讀且美觀——scVelo安裝指南地址憋槐。它除了基于steady-state的模型計(jì)算RNA velocity,還提供了基于dynamic模型計(jì)算RNA velocity淑趾。
但是阳仔,它也只能提供loom文件后續(xù)的計(jì)算,所以扣泊,還是需要先安裝velocyto.py
安裝代碼:
#盡量先使用第一個(gè)方法安裝
pip install -U scvelo
#開(kāi)發(fā)者版本
#or
pip install git+https://github.com/theislab/scvelo@develop
#or
git clone https://github.com/theislab/scvelo && cd scvelo
git checkout --track origin/develop
pip install -e .
二近范、使用指南
(1)通過(guò)bam文件生成loom文件
1、10X
Usage: velocyto run10x [OPTIONS] SAMPLEFOLDER GTFFILE
Runs the velocity analysis for a Chromium 10X Sample
10XSAMPLEFOLDER specifies the cellranger sample folder
GTFFILE genome annotation file
Options:
-s, --metadatatable FILE Table containing metadata of the various samples (csv fortmated rows are samples and cols are entries)
-m, --mask FILE .gtf file containing intervals to mask
-l, --logic TEXT The logic to use for the filtering (default: Default)
-M, --multimap Consider not unique mappings (not reccomended)
-@, --samtools-threads INTEGER The number of threads to use to sort the bam by cellID file using samtools
--samtools-memory INTEGER The number of MB used for every thread by samtools to sort the bam file
-t, --dtype TEXT The dtype of the loom file layers - if more than 6000 molecules/reads per gene per cell are expected set uint32 to avoid truncation (default run_10x: uint16)
-d, --dump TEXT For debugging purposes only: it will dump a molecular mapping report to hdf5. --dump N, saves a cell every N cells. If p is prepended a more complete (but huge) pickle report is printed (default: 0)
-v, --verbose Set the vebosity level: -v (only warinings) -vv (warinings and info) -vvv (warinings, info and debug)
--help Show this message and exit.
例子:
velocyto run10x -m mypath/GRCh38_rmsk.gtf mypath2/Sample_01 somepath/refdata-gex-GRCh38-2020-A/genes/genes.gtf
其中:
-
-m:是遮蔽指令延蟹,GRCh38_rmsk.gtf 是需要mark注釋文件评矩。
可以從UCSA網(wǎng)站下載得到:hg38_rmsk.gtf(人類);
可以直接點(diǎn)擊鏈接:UCSC_hg38阱飘;然后點(diǎn)擊get_output
mypath2/Sample_01是cellranger 輸出文件斥杜,Sample_01文件夾下包含一個(gè)outs文件夾,里面存儲(chǔ)了bam文件沥匈。
genes.gtf是cellranger用的基因組注釋gtf文件蔗喂,位于reference目錄下的gene文件夾
注:
- 需要1.6以上版本的samtools,在命令行輸入samtools --version 可以獲取當(dāng)前環(huán)境的版本信息高帖。否則會(huì)出現(xiàn)以下錯(cuò)誤提醒缰儿。
MemoryError: bam file #0 could not be sorted by cells.
This is probably related to an old version of samtools, please install samtools >= 1.6.
In alternative this could be a memory error, try to set the - your system. Otherwise sort manually by samtools ``sort -l [compression] -m [mb_to_use]M -t [tagname] -O BAM -@ [threads_to_use] -o cellsorted_[bamfile] [bamfile]
輸出:
在cellranger輸出文件夾中,會(huì)出現(xiàn)一個(gè)velocyto文件夾散址,里面存儲(chǔ)了該樣本的loom文件:Sample_01.loom
(2)loom文件讀取乖阵、輸出、整合
1预麸、python
#合并loom文件;
#somepath代表文件所存儲(chǔ)的路徑
import loompy as loompy
files = ["somepath/Sample_01.loom","somepath2/Sample_02.loom"]
loompy.combine(files, 'somepath/merge.loom',key="Accession")
#讀取loom文件
import velocyto as vcy
vlm = vcy.VelocytoLoom("somepath/merge.loom")
2义起、R
讀取:
library(velocyto.R)
library(Seurat)
library(SeuratWrappers)
#讀取loom文件
input_loom <- "/velocyto/merge.loom"
sample.loom <- read.loom.matrices(input_loom,engine = "hdf5r")
可以將loom文件整合到seurat文件中师崎,接下來(lái)直接在R中利用seurat對(duì)象分析
#'一定要注意默终,seurat對(duì)象中的細(xì)胞和基因一定要和loom文件中一致
seurat.obj[["spliced"]] <- CreateAssayObject(counts = sample.loom$spliced)
seurat.obj[["unspliced"]] <- CreateAssayObject(counts = sample.loom$unspliced)
seurat.obj[["ambiguous"]] <- CreateAssayObject(counts = sample.loom$ambiguous)
輸出:轉(zhuǎn)化為H5格式輸出
#
library(SeuratDisk)
SaveH5Seurat(seurat.obj, filename = "/velocyto/seurat_loom.h5Seurat")
Convert("/velocyto/seurat_loom.h5Seurat", dest = "h5ad")
(3)計(jì)算RNA速率
1、python
這里展示直接通過(guò)scvelo軟件的分析
scvelo
import sklearn; sklearn.show_versions()
import scvelo as scv
adata = scv.read("/velocyto/seurat_loom.h5ad")
adata
#數(shù)據(jù)統(tǒng)計(jì)
scv.pl.proportions(adata, groupby='cellType', highlight='unspliced', show=True, save="/merge_data_statistics.pdf")
#對(duì)數(shù)據(jù)進(jìn)行過(guò)濾
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)
#流形
scv.pl.velocity_embedding_stream(adata, basis="umap",figsize = (14,10),color="cellType",show = True,save = "/merge_data_velo_stream.pdf")
#arrow
scv.pl.velocity_embedding(adata,arrow_length = 3,color="cellType",arrow_size = 2,dpi = 120,save = "/merge_data_velo_arrow.pdf")
2、R
#輸入,前面生成的seurat_loom文件
#速率分析
velo <- RunVelocity(seurat.obj, deltaT = 1, kCells = 25, fit.quantile = 0.02,
spliced.average = 0.2, unspliced.average = 0.05, ncores = 18)
#全局速率可視化
emb = Embeddings(velo, reduction = "umap")
vel = Tool(velo, slot = "RunVelocity")
#可視化結(jié)果
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)