RNA 速率學(xué)習(xí)筆記

最近學(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

包括:

  1. 一個(gè)命令行工具昂拂,用于從mapping的bam文件中得到的loom文件【主要是得到spliced mRNA count矩陣和unspliced mRNA count矩陣】
  2. 一個(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)告訴大家排吴。]

  1. 一個(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

其中:

  1. -m:是遮蔽指令延蟹,GRCh38_rmsk.gtf 是需要mark注釋文件评矩。
    可以從UCSA網(wǎng)站下載得到:hg38_rmsk.gtf(人類);
    可以直接點(diǎn)擊鏈接:UCSC_hg38阱飘;然后點(diǎn)擊get_output

    UCSC文件下載格式

  2. mypath2/Sample_01是cellranger 輸出文件斥杜,Sample_01文件夾下包含一個(gè)outs文件夾,里面存儲(chǔ)了bam文件沥匈。

  3. genes.gtf是cellranger用的基因組注釋gtf文件蔗喂,位于reference目錄下的gene文件夾

注:

  1. 需要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)

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末齐蔽,一起剝皮案震驚了整個(gè)濱河市两疚,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌含滴,老刑警劉巖诱渤,帶你破解...
    沈念sama閱讀 206,013評(píng)論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場(chǎng)離奇詭異谈况,居然都是意外死亡勺美,警方通過(guò)查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評(píng)論 2 382
  • 文/潘曉璐 我一進(jìn)店門(mén)碑韵,熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)赡茸,“玉大人,你說(shuō)我怎么就攤上這事祝闻≌嘉裕” “怎么了?”我有些...
    開(kāi)封第一講書(shū)人閱讀 152,370評(píng)論 0 342
  • 文/不壞的土叔 我叫張陵联喘,是天一觀的道長(zhǎng)华蜒。 經(jīng)常有香客問(wèn)我,道長(zhǎng)豁遭,這世上最難降的妖魔是什么叭喜? 我笑而不...
    開(kāi)封第一講書(shū)人閱讀 55,168評(píng)論 1 278
  • 正文 為了忘掉前任,我火速辦了婚禮蓖谢,結(jié)果婚禮上域滥,老公的妹妹穿的比我還像新娘。我一直安慰自己蜈抓,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 64,153評(píng)論 5 371
  • 文/花漫 我一把揭開(kāi)白布昂儒。 她就那樣靜靜地躺著沟使,像睡著了一般。 火紅的嫁衣襯著肌膚如雪渊跋。 梳的紋絲不亂的頭發(fā)上腊嗡,一...
    開(kāi)封第一講書(shū)人閱讀 48,954評(píng)論 1 283
  • 那天,我揣著相機(jī)與錄音拾酝,去河邊找鬼燕少。 笑死,一個(gè)胖子當(dāng)著我的面吹牛蒿囤,可吹牛的內(nèi)容都是我干的客们。 我是一名探鬼主播,決...
    沈念sama閱讀 38,271評(píng)論 3 399
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼底挫!你這毒婦竟也來(lái)了恒傻?” 一聲冷哼從身側(cè)響起,我...
    開(kāi)封第一講書(shū)人閱讀 36,916評(píng)論 0 259
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤建邓,失蹤者是張志新(化名)和其女友劉穎盈厘,沒(méi)想到半個(gè)月后,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體官边,經(jīng)...
    沈念sama閱讀 43,382評(píng)論 1 300
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡沸手,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,877評(píng)論 2 323
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了注簿。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片契吉。...
    茶點(diǎn)故事閱讀 37,989評(píng)論 1 333
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖滩援,靈堂內(nèi)的尸體忽然破棺而出栅隐,到底是詐尸還是另有隱情,我是刑警寧澤玩徊,帶...
    沈念sama閱讀 33,624評(píng)論 4 322
  • 正文 年R本政府宣布租悄,位于F島的核電站,受9級(jí)特大地震影響恩袱,放射性物質(zhì)發(fā)生泄漏泣棋。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,209評(píng)論 3 307
  • 文/蒙蒙 一畔塔、第九天 我趴在偏房一處隱蔽的房頂上張望潭辈。 院中可真熱鬧,春花似錦澈吨、人聲如沸把敢。這莊子的主人今日做“春日...
    開(kāi)封第一講書(shū)人閱讀 30,199評(píng)論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)修赞。三九已至,卻和暖如春桑阶,著一層夾襖步出監(jiān)牢的瞬間柏副,已是汗流浹背。 一陣腳步聲響...
    開(kāi)封第一講書(shū)人閱讀 31,418評(píng)論 1 260
  • 我被黑心中介騙來(lái)泰國(guó)打工蚣录, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留割择,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 45,401評(píng)論 2 352
  • 正文 我出身青樓萎河,卻偏偏與公主長(zhǎng)得像荔泳,于是被迫代替她去往敵國(guó)和親蕉饼。 傳聞我的和親對(duì)象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,700評(píng)論 2 345

推薦閱讀更多精彩內(nèi)容