scVelo用于RNA 速率基本流程
在這里，您將學(xué)習(xí)RNA速率分析的基本流程。
示例數(shù)據(jù)菜用胰腺發(fā)育數(shù)據(jù)集，其譜系包含四個(gè)主要細(xì)胞命運(yùn)：α、β、δ和ε細(xì)胞。有關(guān)詳細(xì)信息，請(qǐng)參閱此處?？梢园凑障嗤乃悸窇?yīng)用于您自己的數(shù)據(jù)。
[ ]:
# update to the latest version, if not done yet.
!pip install scvelo --upgrade --quiet
[1]:
import scvelo as scv
scv.logging.print_version()
Running scvelo 0.2.0 (python 3.8.2) on 2020-05-16 12:55.
[2]:
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
加載數(shù)據(jù)
分析基于內(nèi)置胰腺數(shù)據(jù)。
要對(duì)自己的數(shù)據(jù)進(jìn)行速率分析，請(qǐng)通過(guò)adata = scv.read('path/file.loom', cache=True)將您的文件（loom, h5ad, csv …）讀取到AnnData數(shù)據(jù)對(duì)象。
如果您想將loom文件合并到已存在的 AnnData 對(duì)象中，請(qǐng)使用scv.utils.merge(adata, adata_loom)
[3]:
adata = scv.datasets.pancreas()
adata
[3]:
AnnData object with n_obs × n_vars = 3696 × 27998
obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
var: 'highly_variable_genes'
uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
obsm: 'X_pca', 'X_umap'
layers: 'spliced', 'unspliced'
scVelo 是基于adata，存儲(chǔ)了數(shù)據(jù)矩陣adata.X、觀測(cè)注釋adata.obs、變量adata.var和非結(jié)構(gòu)化注釋的對(duì)象adata.uns。觀測(cè)和變量的名稱可分別通過(guò)adata.obs_names和adata.var_names訪問(wèn)。例如，adata對(duì)象可以像adata_subset = adata[:, list_of_gene_names]數(shù)據(jù)幀一樣切片。有關(guān)詳細(xì)信息，請(qǐng)參閱[anndata docs]。
[4]:
scv.pl.proportions(adata)
這里，展示了剪切/未剪切count的比例。根據(jù)使用的協(xié)議（Drop-Seq，Smart-Seq），我們通常有10%-25%的未剪切分子含有內(nèi)含子序列。我們還建議您檢查cluster水平的變化，以驗(yàn)證剪切效率的一致性。在這里，我們發(fā)現(xiàn)變化如預(yù)期的那樣，在循環(huán)導(dǎo)管細(xì)胞中未剪切的比例略低，在許多基因開始轉(zhuǎn)錄的Ngn3高表達(dá)的和內(nèi)分泌前細(xì)胞中的比例更高。
預(yù)處理數(shù)據(jù)
預(yù)處理包括通過(guò)檢測(cè)（count數(shù)最少）和高變異性（分散）進(jìn)行基因選擇，按其總大小使每個(gè)細(xì)胞標(biāo)準(zhǔn)化。
所有這些都?xì)w入一個(gè)函數(shù)scv.pp.filter_and_normalize，主要運(yùn)行如下：
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)
此外，我們需要在PCA空間中最近的鄰居之間計(jì)算的第一和第二順序時(shí)刻，匯總在scv.pp.moments，它計(jì)算了scv.pp.pca和scv.pp.neighbors。確定速率估計(jì)需要第一個(gè)順序，而隨機(jī)估計(jì)也需要第二個(gè)順序時(shí)刻。
[5]:
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
Filtered out 20801 genes that are detected in less than 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Logarithmized X.
computing neighbors
finished (0:00:03) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:00) --> added
'Ms' and 'Mu', moments of spliced/unspliced abundances (adata.layers)
進(jìn)一步的預(yù)處理（如批次校正）可用于消除不必要的變異源。有關(guān)更多詳細(xì)信息，請(qǐng)參閱最佳實(shí)踐。請(qǐng)注意，任何額外的預(yù)處理步驟僅影響count data X，不應(yīng)用于剪切/未剪切計(jì)數(shù)。
估計(jì)RNA速率
速率是基因表達(dá)空間中的載體，代表單個(gè)細(xì)胞運(yùn)動(dòng)的方向和速率。速率通過(guò)對(duì)拼接動(dòng)力學(xué)的轉(zhuǎn)錄動(dòng)力學(xué)進(jìn)行建模獲得，無(wú)論是隨機(jī)模型（默認(rèn)）還是確定性模型（通過(guò)設(shè)置mode='deterministic'）。對(duì)于每個(gè)基因，預(yù)成熟（未剪切）和成熟（剪切）mRNA計(jì)數(shù)的穩(wěn)定狀態(tài)比，構(gòu)成一個(gè)恒定的轉(zhuǎn)錄狀態(tài)。然后，從此比率中獲取速率作為殘差。正速率表示基因被向上調(diào)節(jié)，這發(fā)生在細(xì)胞顯示該基因的未剪切mRNA的豐度高于預(yù)期的穩(wěn)定狀態(tài)。相反，負(fù)速表示基因被降低調(diào)節(jié)。
完整的動(dòng)力學(xué)模型的解決方案是通過(guò)設(shè)置mode='dynamical'獲得的，這需要事先運(yùn)行scv.tl.recover_dynamics(adata)。我們將在下一個(gè)教程中詳細(xì)闡述動(dòng)態(tài)模型。
[6]:
scv.tl.velocity(adata)
computing velocities
finished (0:00:01) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
計(jì)算的速率存儲(chǔ)在計(jì)數(shù)矩陣adata.layers中。
然后，基因間速率的組合可用于估計(jì)單個(gè)細(xì)胞的未來(lái)狀態(tài)。為了將速率投射到低維嵌入中，估計(jì)了細(xì)胞到細(xì)胞過(guò)渡的概率。即，對(duì)于每個(gè)速率矢量，我們發(fā)現(xiàn)可能的細(xì)胞過(guò)渡方向。過(guò)渡概率是使用潛在細(xì)胞到細(xì)胞的過(guò)渡和速率矢量之間的共生相關(guān)性計(jì)算的，并存儲(chǔ)在表示速率圖的矩陣中。由此產(chǎn)生的速率圖具有維度n obs×n obs，并總結(jié)了通過(guò)速率矢量很好地解釋的可能的細(xì)胞狀態(tài)變化（對(duì)于運(yùn)行時(shí)速率，也可以通過(guò)設(shè)置approx=True在減少的 PCA 空間上計(jì)算）。
[7]:scv.tl.velocity_graph(adata)
computing velocity graph finished (0:00:10) --> added 'velocity_graph', sparse matrix with cosine correlations (adata.uns)
對(duì)于各種應(yīng)用，速率圖可以轉(zhuǎn)換為過(guò)渡矩陣，通過(guò)應(yīng)用高斯函數(shù)將余弦值相關(guān)性轉(zhuǎn)換為實(shí)際的過(guò)渡概率。這需要設(shè)置scv.utils.get_transition_matrix
如前所述，它內(nèi)部用于將速率投射到低維嵌入中，通過(guò)對(duì)與 scv.tl.velocity_embedding獲得的平均概率過(guò)渡概率進(jìn)行平均轉(zhuǎn)換。此外，我們可以通過(guò)scv.tl.terminal_states沿著馬爾代夫轉(zhuǎn)移矩陣追蹤細(xì)胞的起源和潛在命運(yùn)，從而獲得根細(xì)胞和終點(diǎn)的軌跡。
預(yù)測(cè)速率
最后，速率被投影到任何嵌入，通過(guò)basis指定，并以以下方式之一可視化：
細(xì)胞水平 scv.pl.velocity_embedding，
網(wǎng)格線scv.pl.velocity_embedding_grid
流線型scv.pl.velocity_embedding_stream
請(qǐng)注意，數(shù)據(jù)已預(yù)計(jì)算 UMAP 嵌入，并注釋了群集。在應(yīng)用到您自己的數(shù)據(jù)時(shí)，可以使用scv.tl.umap和scv.tl.louvain。有關(guān)詳細(xì)信息，請(qǐng)參閱[scanpy tutorial]。此外，所有繪圖功能都默認(rèn)使用，并且，您可以相應(yīng)地設(shè)置basis='umap'和color='clusters'
[8]:scv.pl.velocity_embedding_stream(adata, basis='umap')
computing velocity embedding finished (0:00:00) --> added 'velocity_umap', embedded velocity vectors (adata.obsm)
流線顯示的速率矢量可對(duì)發(fā)育過(guò)程進(jìn)行詳細(xì)分析。它準(zhǔn)確地勾畫了導(dǎo)管細(xì)胞和內(nèi)分泌祖細(xì)胞的循環(huán)。此外，它闡明了細(xì)胞的譜系命運(yùn)、細(xì)胞周期回歸和內(nèi)分泌細(xì)胞分化狀態(tài)。
我們?cè)趩渭?xì)胞水平上獲得的速率矢量的最詳細(xì)分辨率，每個(gè)箭頭顯示單個(gè)細(xì)胞運(yùn)動(dòng)的方向和速率。這揭示了 Ngn3 細(xì)胞（黃色）的早期內(nèi)分泌命運(yùn)，以及近端α細(xì)胞（藍(lán)色）和瞬態(tài)β細(xì)胞（綠色）之間的明顯差異。
[9]:scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
解釋速率
這可能是最重要的部分，因?yàn)槲覀兘ㄗh用戶不要將生物結(jié)論限制在預(yù)測(cè)速率上，而是通過(guò)圖像來(lái)檢查個(gè)體基因動(dòng)力學(xué)，以了解特定基因如何支持推斷方向。
請(qǐng)參閱此處的 GIF，了解如何解釋剪切與未剪切的圖像?；蚧顒?dòng)是由轉(zhuǎn)錄調(diào)節(jié)的。特定基因的轉(zhuǎn)錄誘導(dǎo)導(dǎo)致（新轉(zhuǎn)錄的）前體未剪切 mRNA 增加，而相反，抑制或沒(méi)有轉(zhuǎn)錄會(huì)導(dǎo)致未轉(zhuǎn)錄 mRNA 的減少。拼接的 mRNA 由未剪切的 mRNA 生成，并遵循相同的趨勢(shì)，并具有時(shí)滯。時(shí)間是一個(gè)隱藏/潛在的變量。因此，需要從實(shí)際測(cè)量中推斷出動(dòng)態(tài)：相像中顯示的剪切和未剪切的 mRNA。
現(xiàn)在，讓我們來(lái)檢查一些標(biāo)記基因的圖像，通過(guò)scv.pl.velocity(adata, gene_names)或scv.pl.scatter(adata, gene_names)`可視化
[10]:scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
黑線對(duì)應(yīng)于估計(jì)的"穩(wěn)定狀態(tài)"比率，即未剪切與剪切的 mRNA 豐度之比，后者處于恒定的轉(zhuǎn)錄狀態(tài)。特定基因的RNA速率被確定為殘留的，即觀察與穩(wěn)定狀態(tài)線的偏差程度。正速率表示基因被向上調(diào)節(jié)，這發(fā)生在細(xì)胞顯示該基因的未剪切mRNA的豐度高于預(yù)期的穩(wěn)定狀態(tài)。相反，負(fù)速表示基因被降低調(diào)節(jié)。
例如，Cpe解釋了上調(diào)節(jié) Ngn3（黃色）到前內(nèi)分泌（橙色）到β細(xì)胞（綠色）的方向，而Adk則解釋了向下調(diào)節(jié)的Ductal（深綠色）到 Ngn3（黃色）到剩余內(nèi)分泌細(xì)胞的方向。
[11]:scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'], add_outline='Ngn3 high EP, Pre-endocrine, Beta')
識(shí)別重要基因
我們需要一種系統(tǒng)的方法來(lái)識(shí)別基因，這可能有助于解釋由此產(chǎn)生的向量場(chǎng)和推斷的譜系。為此，我們可以測(cè)試哪些基因具有群體特異性微分速率表達(dá)，與剩余群相比，其比例要高得多/更低。模塊scv.tl.rank_velocity_genes運(yùn)行差異速率 t 測(cè)試，并生成每個(gè)集群的基因排名?？梢栽O(shè)置閾值（例如`min_corr），以限制對(duì)基因候選者選擇的測(cè)試。
[12]:scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])df.head()#ranking velocity genes# finished (0:00:02) --> added# 'rank_velocity_genes', sorted scores by group ids (adata.uns)# 'spearmans_score', spearmans correlation scores (adata.var)[12]:
[13]: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 高表達(dá)的 EP（黃色）到內(nèi)分泌前（橙色）到Beta（綠色）的方向性。
循環(huán)祖細(xì)胞的速率
RNA速率檢測(cè)到的細(xì)胞周期，通過(guò)細(xì)胞周期分?jǐn)?shù)（相標(biāo)記基因平均表達(dá)水平的標(biāo)準(zhǔn)化分?jǐn)?shù)）在生物學(xué)上得到了確認(rèn)。
[14]:scv.tl.score_genes_cell_cycle(adata)scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
calculating cell cycle phase--> 'S_score' and 'G2M_score', scores of cell cycle phases (adata.obs)
對(duì)于循環(huán)導(dǎo)管細(xì)胞，我們可以篩選通過(guò)S和G2M相標(biāo)記。前一個(gè)模塊還計(jì)算了一個(gè)spearmans 相關(guān)分?jǐn)?shù)，我們可以用它來(lái)對(duì)相標(biāo)記基因進(jìn)行排名/排序，然后顯示其相像。
[15]:s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).indexg2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).indexkwargs = dict(frameon=False, ylabel='cell cycle genes')scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)
特別是Hells 和Top2a非常適合解釋循環(huán)祖細(xì)胞的矢量場(chǎng)。Top2a 在 G2M*階段實(shí)際達(dá)到峰值前不久被分配了高速。在那里，負(fù)速然后完美匹配了接下來(lái)的下調(diào)。
[16]:scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
速率和連貫性
還有兩個(gè)有用的統(tǒng)計(jì)數(shù)據(jù)：
- 分化的速率或速率由速率矢量的長(zhǎng)度給出。
- 矢量場(chǎng)的一致性（即速率矢量如何與其鄰近速率相關(guān)）提供了置信度的衡量標(biāo)準(zhǔn)。
[17]:scv.tl.velocity_confidence(adata)keys = 'velocity_length', 'velocity_confidence'scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
--> added 'velocity_length' (adata.obs)--> added 'velocity_confidence' (adata.obs)--> added 'velocity_confidence_transition' (adata.obs)
這些提供了細(xì)胞在哪里以較慢/更快的速率分化以及方向未確定的見(jiàn)解。
在集群層面，我們發(fā)現(xiàn)，細(xì)胞周期回歸（Ngn3低表達(dá)的EP）后分化速率顯著加快，在β細(xì)胞生產(chǎn)過(guò)程中保持步調(diào)，同時(shí)在阿爾法細(xì)胞生產(chǎn)過(guò)程中減速。
[18]:df = adata.obs.groupby('clusters')[keys].mean().Tdf.style.background_gradient(cmap='coolwarm', axis=1)[18]:
速率圖和偽時(shí)間
我們可以可視化速率圖，以描繪所有速率推斷的細(xì)胞到細(xì)胞連接/過(guò)渡。它可以通過(guò)設(shè)置 threshold.限制在高概率轉(zhuǎn)換。例如，該圖指示來(lái)自早期和晚期內(nèi)分泌細(xì)胞的 Epsilon 細(xì)胞產(chǎn)生的兩個(gè)階段。
[19]:scv.pl.velocity_graph(adata, threshold=.1)
此外，該圖可用于繪制來(lái)自特定細(xì)胞的后代/祖先。在這里，內(nèi)分泌前細(xì)胞可以追溯到其潛在的命運(yùn)。
[20]: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)
最后，根據(jù)速率圖，可以計(jì)算速率偽時(shí)間。在從圖形中推斷出在根細(xì)胞上的分布后，它測(cè)量從根細(xì)胞開始沿著圖形行走后到達(dá)細(xì)胞所需的平均步數(shù)。
與diffusion 偽時(shí)間分析相反，它含蓄地推斷根細(xì)胞，并基于定向速率圖，而不是基于相似性的diffusion kernel。
[21]:scv.tl.velocity_pseudotime(adata)scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
computing terminal states identified 2 regions of root cells and 1 region of end points finished (0:00:00) --> added 'root_cells', root cells of Markov diffusion process (adata.obs) 'end_points', end points of Markov diffusion process (adata.obs)
PAGA速率圖
PAGA圖形概念已作為軌跡推理的性能最高的方法進(jìn)行基準(zhǔn)測(cè)試。它提供了數(shù)據(jù)拓?fù)鋱D的圖形狀地圖，其加權(quán)邊緣與兩個(gè)群集之間的連接相對(duì)應(yīng)。在這里，PAGA 通過(guò)速率推斷方向性擴(kuò)展。
[ ]:# PAGA requires to install igraph, if not done yet.!pip install python-igraph --upgrade --quiet[22]:# 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')df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).Tdf.style.background_gradient(cmap='Blues').format('{:.2g}')
running PAGA finished (0:00:01) --> added 'paga/transitions_confidence', connectivities adjacency (adata.uns) 'paga/connectivities', connectivities adjacency (adata.uns) 'paga/connectivities_tree', connectivities subtree (adata.uns)[22]:
reads從左/行到右/列讀取，例如分配了從Ductal到Ngn3 low EP的可信過(guò)渡。
此表可以通過(guò)疊加在 UMAP 嵌入上的定向圖進(jìn)行匯總展示。
[23]:scv.pl.paga(adata, basis='umap', size=50, alpha=.1, min_edge_width=2, node_size_scale=1.5)
