參考資料：
參考生物技能樹
Velocyto官網(wǎng)Tutorial
scvelo實(shí)戰(zhàn)教程
Seurat-to-RNA-Velocity
分析步驟：
本教程是velocyto基于Seurat對(duì)象中UMAP和細(xì)胞類型進(jìn)行RNA速率分析树瞭，推斷細(xì)胞Cluster命運(yùn)的狀態(tài)（過渡與穩(wěn)定）和分化方向性（軌跡）馆截。
第一步：在linux系統(tǒng)中用velocyto將bam文件轉(zhuǎn)換成loom文件审磁；
第二步：按照scvelo需要的格式獲取seurat中的UMAP坐標(biāo)和Celltype信息戚揭；
第三步：在 python程序中整合loom文件和UMAP坐標(biāo)和Celltype數(shù)據(jù)，并進(jìn)行下游速率分析。
1、將bam文件轉(zhuǎn)換成loom文件（Linux系統(tǒng)操作）
1.1 conda安裝velocyto的一些依賴
需要一些依賴
conda create -n velocyto
conda activate velocyto
conda install samtools
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install pysam
pip install velocyto
Successfully installed loompy-3.0.6 numpy-groupies-0.9.13 pandas-1.2.5 pytz-2021.1 velocyto-0.17.17
1.2 下載特定物種的特殊gtf文件（這個(gè)步驟是可選的）
你如果要屏蔽repetetive elements，可以使用這一步驟
我們這個(gè)單細(xì)胞轉(zhuǎn)錄組使用cellranger流程的話掠手，需要重復(fù)數(shù)據(jù)的gtf文件，rmsk贱傀。
進(jìn)入U(xiǎn)CSC官網(wǎng)惨撇，在Tools菜單中打開Table Brower，選擇目的物種以及個(gè)性化的需求
1.3 bam文件循環(huán)進(jìn)行sort
理論上府寒，這下一步命令魁衙，就可以完成bam轉(zhuǎn)換成loom文件，但是因?yàn)?samtools問題株搔，上面的流程可能是會(huì)失敗剖淀。這個(gè)時(shí)候可以自行先運(yùn)行 samtools sort 命令處理得到 cellsorted_possorted_genome_bam.bam 文件。
ls */outs/possorted_genome_bam.bam|while read id
do new=${id/possorted_genome_bam.bam/cellsorted_possorted_genome_bam.bam}
echo $new
samtools sort -@ 4 -t CB -O BAM -o $new $id
done
1.4 bam文件循環(huán)運(yùn)行velocyto
#! /bin/sh
rmsk_gtf=/lustre/home/acct-medxh/medxh/sccdata/Rawdata/reference/hg38_repeat_rmsk.gtf # 從genome.ucsc.edu下載
cellranger_outDir=/lustre/home/acct-medxh/medxh/sccdata/Rawdata/output_EC3 # 這個(gè)輸出目錄可以隨便選擇文件夾纤房，但是最后一個(gè)文件夾決定了loom的細(xì)胞名字如”output_EU1:TTTCCTCAGTCCGCGTx“纵隔，所以output_EU1盡量為樣品名
cellranger_gtf=/lustre/home/acct-medxh/medxh/yard/reference/refdata-gex-GRCh38-2020-A/genes/genes.gtf # 這個(gè)是cellranger官網(wǎng)提供的
ls -lh $rmsk_gtf $cellranger_outDir $cellranger_gtf
ls -d *-*|while read cellranger_outDir
do
velocyto run10x -m $rmsk_gtf $cellranger_outDir $cellranger_gtf
done
2、Seurat UMAP和Celltype提扰谝獭（R和python）
參考官網(wǎng)https://scvelo.readthedocs.io/getting_started/
2.1 提取seurat中的cells捌刮、UMAP、celltype (R studio)
setwd("~/Desktop/bioinformation/End_analysis/velocyto")
library(dplyr)
library(ggplot2)
library(cowplot)
library(patchwork)
suppressMessages(library(Seurat))
End.combined <- readRDS("~/Desktop/bioinformation/End_analysis/End.combined.rds")
2.2 提取對(duì)象
table(End.combined$group)
Idents(object=End.combined)<-End.combined$group
Eutopic<-subset(End.combined, idents = c("Eutopic"), invert = F)
2.3 修改Cells名字（”ATCCCTGTCACTGAAC-1_1"）保持和loom中cell ID（"output_EU3:TTTGACTGTTGCCGACx"）一致
df<-data.frame(Cells=Cells(Eutopic),sample=paste("output_EU",Eutopic$sampleID,sep = ""))#制做data.frame舒岸，含有cells和output_EU3
df$id<-sapply(df$Cells,function(x)paste(unlist(strsplit(x, "-"))[1],"x",sep = ""))#將“-”之前的字符和“x”連接
df$Cells<-paste(df$sample,df$id,sep = ":")#將“output_EU3”和Cells通過“：”連接
write.csv(df$Cells, file = "cellID_obs.csv", row.names = FALSE) #提取修改后的cell ID
cell_embeddings<-Embeddings(Eutopic, reduction = "umap")#使用Embeddings功能提取seurat UMAP或者TSNE
rownames(cell_embeddings)<-df$Cells #修改rowname的名字使其和cell ID一致
write.csv(cell_embeddings, file = "cell_embeddings.csv")
clusters_obs<-Eutopic$celltype.1 #提取celltype
names(clusters_obs)<-df$Cells #修改名字和保持cell ID一致
write.csv(clusters_obs, file = "clusters_obs.csv")
3绅作、anndata導(dǎo)入loom文件和Seurat中meta-data并進(jìn)行下游速率分析（多樣本整合）
3.1加載依賴包
cona activate velocyto
python
import anndata
import scvelo as scv
import pandas as pd
import numpy as np
import matplotlib as plt
3.2 讀取每個(gè)樣品中的loom文件
sample_one = anndata.read_loom("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/output_EU1.loom")
##Variable names are not unique. To make them unique, call `.var_names_make_unique`.
sample_one.var_names_make_unique()
sample_one
##AnnData object with n_obs × n_vars = 5992 × 36601
# obs: 'Clusters', '_X', '_Y'
# var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand'
# layers: 'matrix', 'ambiguous', 'spliced', 'unspliced'
sample_two = anndata.read_loom("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/output_EU2.loom")
sample_two.var_names_make_unique()
sample_three = anndata.read_loom("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/output_EU3.loom")
sample_three.var_names_make_unique()
3.3 讀取seurat中cells、UMAP和celltype的信息
sample_obs = pd.read_csv("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/cellID_obs.csv")
umap = pd.read_csv("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/cell_embeddings.csv")
cell_clusters = pd.read_csv("/Users/linyu1/Desktop/bioinformation/End_analysis/velocyto/clusters_obs.csv")
3.4 在每個(gè)樣品loom中過濾出seurat中選擇的Cell ID
sample_one = sample_one[np.isin(sample_one.obs.index,sample_obs["x"])]
sample_two = sample_two[np.isin(sample_two.obs.index,sample_obs["x"])]
sample_three = sample_three[np.isin(sample_three.obs.index,sample_obs["x"])]
3.5 將過濾后的文件整合成一個(gè)文件
adata = sample_one.concatenate(sample_two, sample_three)
3.6 將umap坐標(biāo)和celltype添加到anndata 對(duì)象中
需要確保我們添加umapCell ID和anndata對(duì)象中的Cell ID的順序完全匹配蛾派，Cell ID是對(duì)象的觀察層中的行名俄认，因此我們可以使用以下方法更改
一步、我們把index提取成表達(dá)矩陣并更改列名洪乍；
二步眯杏、由于多個(gè)樣品合并成adata時(shí)Cell ID后面多了“-0”（如下圖）,為了完全匹配Cell ID，通過apply和split去除“-0”（保證adata和seurat中barcode名字相同）壳澳；
三步岂贩、umap坐標(biāo)和celltype添加到anndata 對(duì)象。
adata_index = pd.DataFrame(adata.obs.index)
adata_index = adata_index.rename(columns = {0:'Cell ID'})
adata_index = adata_index.rename(columns = {"CellID":'Cell ID'})
rep=lambda x : x.split("-")[0]
adata_index["Cell ID"]=adata_index["Cell ID"].apply(rep)
umap = umap.rename(columns = {'Unnamed: 0':'Cell ID'})#更改umap的列名統(tǒng)一相同的列名Cell ID
umap = umap[np.isin(umap["Cell ID"],adata_index["Cell ID"])] #過濾adata_index在umap中的cell ID
umap=umap.drop_duplicates(subset=["Cell ID"]) #去除重復(fù)值
umap_ordered = adata_index.merge(umap, on = "Cell ID")#依據(jù)adata_index Cell ID順序與umap的數(shù)據(jù)進(jìn)行合并
umap_ordered = umap_ordered.iloc[:,1:] #去除umap_ordered中的第一列
adata.obsm['X_umap'] = umap_ordered.values #將UMAP并入到adata對(duì)象中
cell_clusters=cell_clusters.iloc[:,1:]
cell_clusters = cell_clusters.rename(columns = {'cell':'Cell ID'})
cell_clusters = cell_clusters[np.isin(cell_clusters["Cell ID"],adata_index["Cell ID"])]
cell_clusters=cell_clusters.drop_duplicates(subset=["Cell ID"]) #去除重復(fù)值
cell_clusters_ordered = adata_index.merge(cell_clusters, on = "Cell ID")
cell_clusters_ordered = cell_clusters_ordered.iloc[:,1:]
adata.obs['clusters']=cell_clusters_ordered.values
3.8 運(yùn)行RNA Velocity
我們現(xiàn)在可以運(yùn)行scVelo命令巷波，并根據(jù)Seurat UMAP坐標(biāo)生成RNA速度圖
scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)
scv.tl.velocity(adata, mode = "stochastic")
scv.tl.velocity_graph(adata)
3.9 spliced/unspliced的比例
scv.pl.proportions(adata)
3.10 可視化
scv.pl.velocity_embedding(adata, basis = 'umap')
scv.pl.velocity_embedding_stream(adata, basis = 'umap')
3.11 解釋velocities
這可能是最重要的部分河闰，因?yàn)槲覀兘ㄗh用戶不要將生物學(xué)結(jié)論局限于預(yù)測(cè)的速度科平，而是通過階段畫像來檢查個(gè)體基因動(dòng)力學(xué)褥紫，以了解推斷的方向是如何由特定基因支持的姜性。
參考下圖思考如何解釋spliced vs. unspliced 時(shí)相特征[圖1.6a]∷杩迹基因活性是由轉(zhuǎn)錄調(diào)控調(diào)控的部念。對(duì)特定基因的轉(zhuǎn)錄的誘導(dǎo)表達(dá)，會(huì)導(dǎo)致前體unspliced mRNAs 的增加氨菇，而相反地儡炼，轉(zhuǎn)錄抑制或缺失會(huì)導(dǎo)致unspliced mRNAs 的減少。spliced mRNA由未剪接的unspliced的mRNA產(chǎn)生查蓉，具有相同的趨勢(shì)乌询，但存在時(shí)間滯后。時(shí)間是一個(gè)隱藏/潛在的變量豌研。因此妹田，動(dòng)力學(xué)需要從實(shí)際測(cè)量的東西來推斷: 在相位圖中顯示的剪接和未剪接的mrna。
用scv.pl.velocity(adata, gene_names檢查某些特定基因的時(shí)相特征
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
黑線對(duì)應(yīng)于估計(jì)的“穩(wěn)態(tài)”比率鹃共，即處于恒定轉(zhuǎn)錄狀態(tài)的spliced 和 unspliced mRNA豐度之比鬼佣。一個(gè)特定基因的RNA速度是由殘差決定的，即觀察值偏離穩(wěn)態(tài)線的程度霜浴。正速度表明一個(gè)基因被上調(diào)晶衷，這種現(xiàn)象發(fā)生在細(xì)胞中，在穩(wěn)定狀態(tài)下阴孟，該基因的unspliced mRNA的豐度高于預(yù)期晌纫。相反，負(fù)速度表明基因被下調(diào)永丝。
例如锹漱，Cpe解釋了上調(diào)的Ngn3(黃色)向Pre-endocrine(橙色)向β-cells(綠色)的方向，而Adk解釋了下調(diào)的Ductal(深綠色)向Ngn3(黃色)向其余內(nèi)分泌細(xì)胞的方向类溢。
scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'],
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
3.12 鑒定重要基因
我們需要一個(gè)系統(tǒng)的方法來識(shí)別基因凌蔬，這可能有助于解釋最終的向量場(chǎng)和推斷的譜系。為了做到這一點(diǎn)闯冷，我們可以測(cè)試哪些基因具有集群特有的差異速度表達(dá)砂心，與其他種群相比顯著地更高或更低。模塊scv.tl.rank_velocity_genes運(yùn)行一個(gè)差分速度t檢驗(yàn)蛇耀，并為每個(gè)集群輸出一個(gè)基因排名辩诞。可以設(shè)置閾值(例如min_corr)來限制對(duì)候選基因的測(cè)試纺涤。
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)
例如译暂，基因Ptprs抠忘、Pclo、Pam外永、Abcc8崎脉、Gnas支持從Ngn3高EP(黃色)到Pre-endocrine(橙色)再到Beta(綠色)的方向性。
3.13 細(xì)胞周期的Velocities分析
由RNA速度檢測(cè)的細(xì)胞周期伯顶，在生物學(xué)上由細(xì)胞周期分?jǐn)?shù)(階段標(biāo)志基因平均表達(dá)水平的標(biāo)準(zhǔn)化分?jǐn)?shù)) 確定囚灼。
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)
s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).index
g2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).index
kwargs = dict(frameon=False, ylabel='cell cycle genes')
scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)
特別是hell和Top2a非常適合解釋周期祖先中的向量場(chǎng)。Top2a在G2M階段達(dá)到峰值前不久被賦予了一個(gè)高速祭衩。在這里灶体，負(fù)的速度與緊隨其后的下調(diào)完全匹配。
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
3.14 速度和一致性
兩個(gè)更有用的stats:-速度或分化率是由速度矢量的長(zhǎng)度給出的掐暮。-矢量場(chǎng)的一致性(即蝎抽，速度矢量如何與其鄰近速度相關(guān)聯(lián))提供了一個(gè)可信度的度量。
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
這些提供了細(xì)胞在哪里以較慢/較快的速度分化
在集群水平上路克，我們發(fā)現(xiàn)分化在細(xì)胞周期退出(Ngn3 low EP)后顯著加快樟结，在Beta細(xì)胞生產(chǎn)期間保持速度，而在Alpha細(xì)胞生產(chǎn)期間減慢速度衷戈。
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
3.15 動(dòng)力學(xué)模型
它在一個(gè)基于可能性的期望最大化框架中解決狭吼，通過迭代估計(jì)反應(yīng)速率和潛在細(xì)胞特異性變量的參數(shù)，即轉(zhuǎn)錄狀態(tài)和細(xì)胞內(nèi)部潛伏時(shí)間殖妇。因此刁笙，它的目的是了解每個(gè)基因未拼接/已拼接階段的軌跡。
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
#adata.write('data/pancreas.h5ad', compression='gzip')
#adata = scv.read('data/pancreas.h5ad')
scv.pl.velocity_embedding_stream(adata, basis='umap')
3.16 動(dòng)率參數(shù)
在不需要任何實(shí)驗(yàn)數(shù)據(jù)的情況下谦趣，可以估計(jì)RNA轉(zhuǎn)錄疲吸、剪接和降解的速率。
它們有助于更好地理解細(xì)胞的特性和表型異質(zhì)性前鹅。
df = adata.var
df = df[(df['fit_likelihood'] > .1) & df['velocity_genes'] == True]
kwargs = dict(xscale='log', fontsize=16)
with scv.GridSpec(ncols=3) as pl:
pl.hist(df['fit_alpha'], xlabel='transcription rate', **kwargs)
pl.hist(df['fit_beta'] * df['fit_scaling'], xlabel='splicing rate', xticks=[.1, .4, 1], **kwargs)
pl.hist(df['fit_gamma'], xlabel='degradation rate', xticks=[.1, .4, 1], **kwargs)
scv.get_df(adata, 'fit*', dropna=True).head()
The estimated gene-specific parameters comprise rates of transription (fit_alpha), splicing (fit_beta), degradation (fit_gamma), switching time point (fit_t_), a scaling parameter to adjust for under-represented unspliced reads (fit_scaling), standard deviation of unspliced and spliced reads (fit_std_u, fit_std_s), the gene likelihood (fit_likelihood), inferred steady-state levels (fit_steady_u, fit_steady_s) with their corresponding p-values (fit_pval_steady_u, fit_pval_steady_s), the overall model variance (fit_variance), and a scaling factor to align the gene-wise latent times to a universal, gene-shared latent time (fit_alignment_scaling).
3.17 逆時(shí)分析
動(dòng)態(tài)模型恢復(fù)了潛在細(xì)胞過程的潛在時(shí)間摘悴。這一潛伏時(shí)間代表了細(xì)胞的內(nèi)部時(shí)鐘，僅根據(jù)其轉(zhuǎn)錄動(dòng)力學(xué)舰绘，與細(xì)胞分化時(shí)所經(jīng)歷的實(shí)時(shí)時(shí)間近似蹂喻。
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=80)
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
scv.pl.heatmap(adata, var_names=top_genes, sortby='latent_time', col_color='clusters', n_convolve=100)
2.3 Top-likelihood基因
驅(qū)動(dòng)基因表現(xiàn)出明顯的動(dòng)態(tài)行為，并通過動(dòng)態(tài)模型中的高可能性對(duì)其特征進(jìn)行系統(tǒng)地檢測(cè)捂寿。
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, frameon=False)
var_names = ['Actn4', 'Ppp3ca', 'Cpe', 'Nnat']
scv.pl.scatter(adata, var_names, frameon=False)
scv.pl.scatter(adata, x='latent_time', y=var_names, frameon=False)
3.18 Cluster-specific top-likelihood 基因
此外口四，可以計(jì)算每個(gè)細(xì)胞集群的部分基因可能性，以實(shí)現(xiàn)集群特異性識(shí)別潛在驅(qū)動(dòng)因素秦陋。
scv.tl.rank_dynamical_genes(adata, groupby='clusters')
df = scv.get_df(adata, 'rank_dynamical_genes/names')
df.head(5)
for cluster in ['Ductal', 'Ngn3 high EP', 'Pre-endocrine', 'Beta']:
scv.pl.scatter(adata, df[cluster][:5], ylabel=cluster, frameon=False)
