RNA速率分析是一個非常重要的分析內(nèi)容，我們這里來深入解析一下這個分析內(nèi)容（python版本，個人喜歡python），參考網(wǎng)址在Velocyto
安裝velocyto其實很簡單，直接pip安裝
pip install velocyto
至于用法：
velocyto --help
Usage: velocyto [OPTIONS] COMMAND [ARGS]...
Options:
--version Show the version and exit.
--help Show this message and exit.
Commands:
run Runs the velocity analysis outputting a loom file
run10x Runs the velocity analysis for a Chromium Sample
run-dropest Runs the velocity analysis on DropEst preprocessed data
run-smartseq2 Runs the velocity analysis on SmartSeq2 data (independent bam file per cell)
tools helper tools for velocyto
其實做velocyto分析使用的是轉(zhuǎn)錄組的可變剪切，那么對于轉(zhuǎn)錄組的要求當(dāng)然是越長越好，最好是全長，也就是說用smartseq2的數(shù)據(jù)做這個分析最為合適，但是10X數(shù)據(jù)最多最火，也可以嘗試。
velocyto run10x --help
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 repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf
產(chǎn)生的loom文件將用于下游的分析。
我們這里詳細看看這個文件及分析是怎樣的。（這里還是python為例，R也可以R讀loom需要hdf5r，非常難裝，裝成功了一定要記錄，這也是我討厭R的原因之一）
首先第一步，我們來看看生成的loom文件里面都有什么
import scvelo as scv
scv.set_figure_params()##設(shè)置可視化的一些參數(shù)
data = scv.read('TH_possorted_genome_bam_YHVE5.loom', cache=True) ###讀取loom文件
當(dāng)然我們這里也可以讀取h5ad格式的文件（由python分析軟件scanpy生成）。
data = scv.read(filename, cache=True)
我們來看看loom文件里面有什么：
data
AnnData object with n_obs × n_vars = 70668 × 27998
var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
layers: 'matrix', 'ambiguous', 'spliced', 'unspliced'
data.var
Accession Chromosome End Start Strand
Xkr4 ENSMUSG00000051951 1 3671498 3205901 -
Gm37381 ENSMUSG00000102343 1 3986215 3905739 -
Rp1 ENSMUSG00000025900 1 4409241 3999557 -
Rp1 ENSMUSG00000109048 1 4409187 4292981 -
Sox17 ENSMUSG00000025902 1 4497354 4490931 -
... ... ... ... ... ...
Gm28672 ENSMUSG00000100492 Y 89225296 89222067 +
Gm28670 ENSMUSG00000099982 Y 89394761 89391528 +
Gm29504 ENSMUSG00000100533 Y 90277501 90275224 +
Gm20837 ENSMUSG00000096178 Y 90433263 90401248 +
Erdr1 ENSMUSG00000096768 Y 90816464 90784738 +
###var里面是一些基因的信息
ldata.layers
Layers with keys: matrix, ambiguous, spliced, unspliced
###layer下面的東西有點多，包括矩陣信息，以及基因是否是可變剪切等信息矩陣。
可以看出這個地方只有關(guān)于可變剪切的信息，沒有我們單細胞聚類得到的信息，所有我們需要將聚類結(jié)果與該loom文件的內(nèi)容進行整合（merge），注意這里是python，需要讀取scanpy的分析結(jié)果h5ad文件.但這里我們沒有h5ad文件，只好讀取Seurat的單獨結(jié)果進行賦值，這里我們只需要聚類信息和二維坐標(biāo)信息
import pandas as pd
import numpy as np
cluster = pd.read_csv(file)
adata.obs['clusters'] = cluster['Cluster']
tsne_axis=open(file,'r')
for i in all_tsne[1:]:
tsne_axis_dict={}
all_tsne=tsne_axis.readlines()
... key=i.strip().split(',')[0]
... tsne_axis_dict[key]=[]
... tsne_axis_dict[key].append(i.strip().split(',')[1])
... tsne_axis_dict[key].append(i.strip().split(',')[2])
... tsne_axis_dict[key]=np.array(tsne_axis_dict[key],dtype='float32')
tsne_axis_list=[]
for i in Barcode:
tsne_axis_list.append(tsne_axis_dict[i])
adata.obsm['X_tsne']=np.array(tsne_axis_list)
###UMAP的賦值也是一樣的，這樣我們就得到了可以直接分析的無縫銜接的文件。
接下來就是可視化了，就相對簡單了
scv.pl.proportions(adata)
這個圖展示的是每個cluster可變切剪的比例值，這個地方我們需要注意，跟聚類結(jié)果嚴格相關(guān)，也就是說我們必須確保聚類結(jié)果相對完美的情況下做這個分析。
接下來我們要進行RNA速率分析，這里我們要注意，我們嫁接了Seurat 的結(jié)果，不要再進行過濾等操作
scv.tl.velocity(adata)
###The computed velocities are stored in adata.layers just like the count matrices.
scv.tl.velocity_graph(adata)
###For a variety of applications, the velocity graph can be converted to a transition matrix by applying a Gaussian kernel to transform the cosine correlations into actual transition probabilities. You can access the Markov transition matrix via scv.utils.get_transition_matrix.
接下來就是一些展示了，很簡單
scv.pl.velocity_embedding_stream(adata, basis='umap')
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
marker gene的展示
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
Identify important genes
We need a systematic way to identify genes that may help explain the resulting vector field and inferred lineages. To do so, we can test which genes have cluster-specific differential velocity expression, being siginificantly higher/lower compared to the remaining population. The module scv.tl.rank_velocity_genes runs a differential velocity t-test and outpus a gene ranking for each cluster. Thresholds can be set (e.g. min_corr) to restrict the test on a selection of gene candidates.
scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
df.head()
kwargs = dict(frameon=False, size=10, linewidth=1.5,
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
scv.pl.scatter(adata, df['Ngn3 high EP'][:5], ylabel='Ngn3 high EP', **kwargs)
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs)
接下來的可視化分析都很簡單，我就不多說了，大家趕緊去試試吧
請保持憤怒，讓王多魚傾家蕩產(chǎn)
