通過動力學建模將RNA速率推廣到瞬時細胞狀態(tài)
RNA速率開辟了利用scRNA測序數據研究細胞分化的新途徑惭婿。RNA速率基于剪接和非剪接mRNA的比率描述了在一個給定時間點上巷蚪，一個基因表達變化的速率帮孔。然而克胳，如果共同剪接率的核心假設和對穩(wěn)定狀態(tài)mRNA水平的完整剪接動力學的觀察被違反德挣，速率估計就會出現誤差恭垦。在這里，作者提出了scVelo方法格嗅，通過使用基于可能性的動力學模型求解剪接動力學的完整轉錄動力學來克服這些限制番挺。這將RNA速率推廣到具有瞬時細胞狀態(tài)的系統(tǒng)，這種狀態(tài)在發(fā)育和對擾動的反應中是常見的屯掖。
背景簡介
單細胞轉錄學使得在單細胞分辨率下對細胞分化和譜系選擇等生物過程的無偏見研究成為可能玄柏。由此產生的計算問題被稱為軌跡推斷。軌跡推斷算法從處于發(fā)育過程不同階段的一組細胞開始贴铜，旨在重建導致潛在細胞命運的轉錄變化的發(fā)育序列粪摘。軌跡推斷的一個中心挑戰(zhàn)是單細胞RNA測序(scRNA-seq)的破壞性，它只揭示細胞的 某一瞬時靜止狀態(tài)绍坝。為了從描述性軌跡模型轉向預測性軌跡模型徘意，就需要用額外的信息來約束可能導致相同軌跡的可能動力的空間。因此轩褐，譬如血統(tǒng)追蹤分析便可以通過修改基因添加信息椎咧，以重建血統(tǒng)關系。然而把介，這些分析并不容易建立勤讽，而且在許多系統(tǒng)中受到技術限制，如人類組織拗踢。
RNA速率的概念可以通過新轉錄出來的未剪接的前mRNA和成熟的剪接后的mRNA來恢復定向動力信息得以確定脚牍。
假設一個簡單的每個基因的反應模型來將未剪接和剪接的mRNA的豐度聯(lián)系起來，就可以推斷出mRNA豐度的變化巢墅，稱為RNA速率诸狭。正的RNA速率表明一個基因表達上調，這種情況發(fā)生在細胞中砂缩，該基因的未剪接mRNA的豐度比在穩(wěn)定狀態(tài)下預期的更高作谚。相反三娩，負的RNA速率表明一個基因的表達下調庵芭。然后，跨基因的速率組合起來就可以用來估計單個細胞的未來狀態(tài)雀监。
原始的模型評估RNA速率是在假設誘導和抑制基因表達的轉錄階段持續(xù)足夠長時間使其達到活躍轉錄和不活躍轉錄平衡狀態(tài)的前提下双吆。在推斷處于恒定轉錄穩(wěn)定狀態(tài)的未剪接和剪接mRNA豐度的比率后眨唬，速度被確定為觀察到的比率與其穩(wěn)定狀態(tài)比率的偏差。推測穩(wěn)態(tài)比率有兩個基本假設好乐，即(1)在基因水平上匾竿，捕捉到轉錄誘導、抑制和穩(wěn)態(tài)mRNA水平的完整剪接動態(tài)蔚万；(2)在細胞水平上岭妖，所有基因都有共同的剪接率。這些假設經常被違反反璃，特別是當一個種群包括多個具有不同動力學的異質性亞群時昵慌。我們將這種原始的建模方法稱為“穩(wěn)態(tài)模型”。
scVelo淮蜈，基于可能性的動力學模型斋攀，可以解決完整的基因轉錄動力學。從而將RNA速率預測推廣到瞬時系統(tǒng)和具有異質性亞群動力學的系統(tǒng)梧田。我們在一個有效的期望最大化(EM)框架中推導出轉錄淳蔼、剪接和降解的基因特異性反應速度和潛在的基因共享潛伏期。推斷的潛伏時間代表了細胞的內部時鐘裁眯，它準確地描述了細胞在潛在生物過程中的位置鹉梨。與現有的基于相似性的偽時間方法不同，這種潛在時間僅基于轉錄動力學未状，并考慮了運動的速度和方向俯画。推斷的潛伏時間能夠重建轉錄事件和細胞命運的時間序列。與穩(wěn)態(tài)模型相比司草，動態(tài)模型通常產生更一致的跨相鄰細胞的速度估計艰垂，并準確地識別轉錄狀態(tài)。此外埋虹，scVelo還能確定調控例如過渡狀態(tài)和細胞命運轉變的階段變化的機制猜憎。scVelo能確定這些轉錄態(tài)變化的假定驅動基因。驅動基因表現出明顯的動態(tài)行為搔课，并通過動力模型中的高概率特征被系統(tǒng)地檢測到胰柑。這一過程提供了一種基于動力學的標準差異表達圖譜的替代。
如何在單細胞分辨率下求解完整的基因轉錄動力學爬泥？
scVelo提出的基于可能性的模型柬讨，由兩組參數控制：
(1)轉錄
αk(t)、剪接β和降解 γ的反應速度(2)細胞特定的轉錄狀態(tài)
ki和時間ti的潛在變量袍啡。通過最小化細胞到當前相軌跡的距離來為每個細胞分配時間點ti踩官。轉錄狀態(tài)ki是通過將可能性與軌跡上的各個片段相關聯(lián)來分配的--即【誘導】、【抑制】和【活躍和不活躍的穩(wěn)定狀態(tài)】通過EM算法境输，極大似然估計蔗牡，反復迭代來推斷參數颖系。
scVelo計算RNA速率——基本分析
*引自https://scvelo.readthedocs.io/VelocityBasics/
import scvelo as scv
scv.logging.print_version()
scv.settings.verbosity = 3 # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True # set max width size for presenter view
scv.set_figure_params('scvelo') # for beautified visualization
加載數據
scVelo包的示例數據是胰腺癌數據。
To run velocity analysis on your own data, read your file (loom, h5ad, csv …) to an AnnData object with adata = scv.read('path/file.loom', cache=True). If you want to merge your loom file into an already existing AnnData object, use scv.utils.merge(adata, adata_loom).
adata = scv.datasets.pancreas()
adata
scVelo is based on
adata:an object that stores a data matrix
adata.X,annotation of observations
adata.obs,variables
adata.var,and unstructured annotations
adata.uns.Names of observations and variables can be accessed via
adata.obs_names and adata.var_names, respectively.AnnData objects can be sliced like dataframes, for example,
adata_subset = adata[:, list_of_gene_names].For more details, see the anndata docs.
scv.pl.proportions(adata,save=True)
數據預處理
預處理包括以下幾個步驟：
gene selection by detection (with a minimum number of counts) and high variability (dispersion), normalizing every cell by its total size and logarithmizing X. Filtering and normalization is applied in the same vein to spliced/unspliced counts and X. Logarithmizing is only applied to X. If X is already preprocessed from former analysis, it will not be touched.
預處理的過程使用的函數為scv.pp.filter_and_normalize辩越，其執(zhí)行以下四個操作
#這四行代碼無需單獨執(zhí)行嘁扼，只是展示函數scv.pp.filter_and_normalize（）具體做了那些操作
scv.pp.filter_genes(adata, min_shared_counts=20)
scv.pp.normalize_per_cell(adata)
scv.pp.filter_genes_dispersion(adata, n_top_genes=2000)
scv.pp.log1p(adata)
Further, we need the first and second order moments (means and uncentered variances) computed among nearest neighbors in PCA space, summarized in scv.pp.moments, which internally computes scv.pp.pca and scv.pp.neighbors. First order is needed for deterministic velocity estimation, while stochastic estimation also requires second order moments.
#預處理過程執(zhí)行這兩行
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
Further preprocessing (such as batch effect correction) may be used to remove unwanted sources of variability. See the best practices for further details. Note, that any additional preprocessing step only affects X and is not applied to spliced/unspliced counts.
估計RNA速率
Velocities are vectors in gene expression space and represent the direction and speed of movement of the individual cells. The velocities are obtained by modeling transcriptional dynamics of splicing kinetics, either stochastically (default) or deterministically (by setting mode='deterministic').
For each gene, a steady-state-ratio of pre-mature (unspliced) and mature (spliced) mRNA counts is fitted, which constitutes a constant transcriptional state. Velocities are then obtained as residuals from this ratio. Positive velocity indicates that a gene is up-regulated, which occurs for cells that show higher abundance of unspliced mRNA for that gene than expected in steady state. Conversely, negative velocity indicates that a gene is down-regulated.
The solution to the full dynamical model is obtained by setting mode='dynamical', which requires to run scv.tl.recover_dynamics(adata) beforehand. We will elaborate more on the dynamical model in the next tutorial.
scv.tl.velocity(adata)
The computed velocities are stored in
adata.layers just like the count matrices.
The combination of velocities across genes can then be used to estimate the future state of an individual cell. In order to project the velocities into a lower-dimensional embedding, transition probabilities of cell-to-cell transitions are estimated. That is, for each velocity vector we find the likely cell transitions that are accordance with that direction. The transition probabilities are computed using cosine correlation between the potential cell-to-cell transitions and the velocity vector, and are stored in a matrix denoted as velocity graph. The resulting velocity graph has dimension nobs×nobs and summarizes the possible cell state changes that are well explained through the velocity vectors (for runtime speedup it can also be computed on reduced PCA space by setting approx=True).
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.
As mentioned, it is internally used to project the velocities into a low-dimensional embedding by applying the mean transition with respect to the transition probabilities, obtained with scv.tl.velocity_embedding. Further, we can trace cells along the Markov chain to their origins and potential fates, thereby getting root cells and end points within a trajectory, obtained via scv.tl.terminal_states.
可視化結果
Finally, the velocities are projected onto any embedding, specified by basis, and visualized in one of these ways: - on cellular level with scv.pl.velocity_embedding, - as gridlines with scv.pl.velocity_embedding_grid, - or as streamlines with scv.pl.velocity_embedding_stream.
Note, that the data has an already pre-computed UMAP embedding, and annotated clusters. When applying to your own data, these can be obtained with scv.tl.umap and scv.tl.louvain. For more details, see the scanpy tutorial. Further, all plotting functions are defaulted to using basis='umap' and color='clusters', which you can set accordingly.
scv.pl.velocity_embedding_stream(adata, basis='umap',save=f"./stream.pdf")
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120,save=f"./embedding.pdf")
解釋RNA速率
This is perhaps the most important part as we advise the user not to limit biological conclusions to the projected velocities, but to examine individual gene dynamics via phase portraits to understand how inferred directions are supported by particular genes.
See the gif here to get an idea of how to interpret a spliced vs. unspliced phase portrait. Gene activity is orchestrated by transcriptional regulation. Transcriptional induction for a particular gene results in an increase of (newly transcribed) precursor unspliced mRNAs while, conversely, repression or absence of transcription results in a decrease of unspliced mRNAs. Spliced mRNAs is produced from unspliced mRNA and follows the same trend with a time lag. Time is a hidden/latent variable. Thus, the dynamics needs to be inferred from what is actually measured: spliced and unspliced mRNAs as displayed in the phase portrait.
Now, let us examine the phase portraits of some marker genes, visualized with scv.pl.velocity(adata, gene_names) or scv.pl.scatter(adata, gene_names).
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2,save=f"./Interprete_the_velocities.pdf")
黑色虛線代表穩(wěn)態(tài)，黑線上表明基因上調黔攒，黑線下表明基因下調
For instance Cpe explains the directionality in the up-regulated Ngn3 (yellow) to Pre-endocrine (orange) to β-cells (green), while Adk explains the directionality in the down-regulated Ductal (dark green) to Ngn3 (yellow) to the remaining endocrine cells.
scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'],
add_outline='Ngn3 high EP, Pre-endocrine, Beta',save=f"./Cpe.pdf")
scv.pl.scatter(adata, 'Adk', color=['clusters', 'velocity'],save=f"./Adk.pdf")
確定一些重要的驅動基因
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,save=f"./Ngn3_high_EP.pdf")
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs,save=f"./Pre-endocrine.pdf")
The genes Ptprs, Pclo, Pam, Abcc8, Gnas, for instance, support the directionality from Ngn3 high EP (yellow) to Pre-endocrine (orange) to Beta (green).
基因Ptprs, Pclo, Pam, Abcc8支持黃->橙->綠的方向
分析細胞周期中的RNA速率
The cell cycle detected by RNA velocity, is biologically affirmed by cell cycle scores (standardized scores of mean expression levels of phase marker genes).
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95],save=f"./cell_cycle.pdf")
For the cycling Ductal cells, we may screen through S and G2M phase markers. The previous module also computed a spearmans correlation score, which we can use to rank/sort the phase marker genes to then display their phase portraits.
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,save=f"./g2m_genes.pdf")
Particularly Hells and Top2a are well-suited to explain the vector field in the cycling progenitors. Top2a gets assigned a high velocity shortly before it actually peaks in the G2M phase. There, the negative velocity then perfectly matches the immediately following down-regulation.
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True,save=f"./Hells-Top2a.pdf")
分析速度和一致性
Two more useful stats: - The speed or rate of differentiation is given by the length of the velocity vector. - The coherence of the vector field (i.e., how a velocity vector correlates with its neighboring velocities) provides a measure of confidence.
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95],save=f"./coherence.pdf")
These provide insights where cells differentiate at a slower/faster pace, and where the direction is un-/determined.
On cluster-level, we find that differentiation substantially speeds up after cell cycle exit (Ngn3 low EP), keeping the pace during Beta cell production while slowing down during Alpha cell production.
速度圖和偽時間
We can visualize the velocity graph to portray all velocity-inferred cell-to-cell connections/transitions. It can be confined to high-probability transitions by setting a threshold. The graph, for instance, indicates two phases of Epsilon cell production, coming from early and late Pre-endocrine cells.
scv.pl.velocity_graph(adata, threshold=.1,save=f"./velocity_graph.pdf")
Further, the graph can be used to draw descendents/anscestors coming from a specified cell. Here, a pre-endocrine cell is traced to its potential fate.
x, y = scv.utils.get_cell_transitions(adata, basis='umap', starting_cell=70)
ax = scv.pl.velocity_graph(adata, c='lightgrey', edge_width=.05, show=False)
ax = scv.pl.scatter(adata, x=x, y=y, s=120, c='ascending', cmap='gnuplot', ax=ax,save=f"./one_cell_velocity_graph.pdf")
Finally, based on the velocity graph, a velocity pseudotime can be computed. After inferring a distribution over root cells from the graph, it measures the average number of steps it takes to reach a cell after walking along the graph starting from the root cells.
Contrarily to diffusion pseudotime, it implicitly infers the root cells and is based on the directed velocity graph instead of the similarity-based diffusion kernel.
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot',save=f"./pseudotime.pdf")
PAGA速度圖
PAGA graph abstraction has benchmarked as top-performing method for trajectory inference. It provides a graph-like map of the data topology with weighted edges corresponding to the connectivity between two clusters. Here, PAGA is extended by velocity-inferred directionality.
#需要python-igraph包趁啸，安裝命令為：conda install python-igraph
# this is needed due to a current bug - bugfix is coming soon.
adata.uns['neighbors']['distances'] = adata.obsp['distances']
adata.uns['neighbors']['connectivities'] = adata.obsp['connectivities']
scv.tl.paga(adata, groups='clusters')
scv.pl.paga(adata, basis='umap', size=50, alpha=.1,
min_edge_width=2, node_size_scale=1.5,save=f"./PAGA.pdf")
scVelo計算RNA速率——動態(tài)模型（最新）
We run the dynamical model to learn the full transcriptional dynamics of splicing kinetics.
It is solved in a likelihood-based expectation-maximization framework, by iteratively estimating the parameters of reaction rates and latent cell-specific variables, i.e. transcriptional state and cell-internal latent time. It thereby aims to learn the unspliced/spliced phase trajectory for each gene.
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
Running the dynamical model can take a while. Hence, you may want to store the results for re-use, with adata.write('data/pancreas.h5ad'), which can later be read with data = scv.read('data/pancreas.h5ad').
#adata.write('data/pancreas.h5ad', compression='gzip')
#adata = scv.read('data/pancreas.h5ad')
scv.pl.velocity_embedding_stream(adata, basis='umap',save=f"./Umap.pdf")
計算動力學參數
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)
plt.savefig("Kinetic_rate_paramters.pdf",bbox_inches = 'tight')
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).
計算潛伏時間
The dynamical model recovers the latent time of the underlying cellular processes. This latent time represents the cell’s internal clock and approximates the real time experienced by cells as they differentiate, based only on its transcriptional dynamics.
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=80,save=f"./latent_time.pdf")
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,save=f"./heatmap.pdf")
找Top-likelihood 基因
Driver genes display pronounced dynamic behavior and are systematically detected via their characterization by high likelihoods in the dynamic model.
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, frameon=False,save=f"./top_genes.pdf")
var_names = ['Actn4', 'Ppp3ca', 'Cpe', 'Nnat']
scv.pl.scatter(adata, var_names, frameon=False,save=f"./top_genes2.pdf")
scv.pl.scatter(adata, x='latent_time', y=var_names, frameon=False,save=f"./top_genes3.pdf")
找cluster特異性的Top-likelihood 基因
Moreover, partial gene likelihoods can be computed for a each cluster of cells to enable cluster-specific identification of potential drivers.
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)
plt.savefig(cluster+"_top_genes.pdf",bbox_inches = 'tight')
