軌跡分析系列：
速率分析系列：
官網(wǎng)： https://scvelo.readthedocs.io/getting_started/
1. 簡(jiǎn)介
測(cè)量單個(gè)細(xì)胞中的基因活性需要破壞這些細(xì)胞以讀取其內(nèi)容，這使得研究動(dòng)態(tài)過(guò)程和了解細(xì)胞命運(yùn)決定具有挑戰(zhàn)性赂乐。 La Manno et al. (Nature, 2018)引入了 RNA velocity 的概念聊疲，利用新轉(zhuǎn)錄的未剪接的前體 mRNA 和成熟的剪接 mRNA 可以在常見(jiàn)的單細(xì)胞 RNA-seq 流程中區(qū)分的事實(shí)办悟，可以恢復(fù)定向動(dòng)態(tài)信息，前者可通過(guò)內(nèi)含子的存在檢測(cè)簸呈。這種不僅測(cè)量基因活性落午，而且測(cè)量它們?cè)趩蝹€(gè)細(xì)胞中的變化（RNA 速率）的概念，開(kāi)辟了研究細(xì)胞分化的新方法鲫咽。最初提出的框架將速率作為觀察到的剪接和未剪接 mRNA 的比率與推斷的穩(wěn)態(tài)的偏差签赃。如果違反了共同剪接速率的中心假設(shè)和對(duì)具有穩(wěn)態(tài) mRNA 水平的完整剪接動(dòng)力學(xué)的觀察，則會(huì)出現(xiàn)速率估計(jì)錯(cuò)誤分尸。
Bergen et al. (Nature Biotechnology, 2020)開(kāi)發(fā)了 scVelo锦聊，通過(guò)使用基于似然的動(dòng)力學(xué)模型求解剪接動(dòng)力學(xué)的完整轉(zhuǎn)錄動(dòng)力解決了這些局限。這將 RNA 速率推廣到包括瞬態(tài)細(xì)胞狀態(tài)的各種系統(tǒng)寓落，這些系統(tǒng)在發(fā)育和對(duì)擾動(dòng)的響應(yīng)中很常見(jiàn)括丁。此外，scVelo 推斷轉(zhuǎn)錄伶选、剪接和降解的基因特異性速率，并恢復(fù)在細(xì)胞過(guò)程的潛在時(shí)間尖昏。這個(gè)僅基于其轉(zhuǎn)錄動(dòng)力學(xué)的潛在時(shí)間代表細(xì)胞的內(nèi)部時(shí)鐘仰税，并近似于細(xì)胞在分化時(shí)經(jīng)歷的真實(shí)時(shí)間。此外抽诉，scVelo 識(shí)別調(diào)節(jié)變化的機(jī)制陨簇，例如細(xì)胞命運(yùn)決定的階段，并在其中系統(tǒng)地檢測(cè)推定的驅(qū)動(dòng)基因迹淌。（參考：單細(xì)胞測(cè)序的軌跡推斷）
目前的scVelo的主要應(yīng)用：
計(jì)算RNA velocity
確定可能和細(xì)胞動(dòng)態(tài)變化過(guò)程相關(guān)的驅(qū)動(dòng)基因
推定轉(zhuǎn)錄事件的時(shí)序
估算RNA相關(guān)的轉(zhuǎn)錄剪切降解速率
做統(tǒng)計(jì)檢驗(yàn)
核心步驟：
scv.pp.filter_and_normalize(adata)選擇基因河绽，然后標(biāo)準(zhǔn)化
scv.pp.moments(adata)計(jì)算一階和二階矩
scv.tl.velocity(adata)計(jì)算velocity，默認(rèn)mode='stochastic'唉窃，也可以設(shè)置成'dynamics' (Dynamical Modeling)
scv.tl.velocity_graph(adata)將velocity投影到降維之后的封裝(t-SNE/umap)之中
可視化：
scv.pl.velocity_embedding(adata, basis='umap')
scv.pl.velocity_embedding_grid(adata, basis='umap')
scv.pl.velocity_embedding_stream(adata, basis='umap')
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'],add_outline='Ngn3 high EP, Pre-endocrine, Beta')
scv.pl.velocity_graph(adata, threshold=.1)
2. 安裝
scVelo 需要 Python 3.6 或更高版本耙饰。建議使用Miniconda。
- PyPI
使用以下命令從PyPI安裝 scVelo ：
pip install -U scvelo
-U 是--upgrade 的縮寫(xiě)纹份。如果出現(xiàn)Permission denied錯(cuò)誤苟跪，請(qǐng)改用pip install -U scvelo --user
- 開(kāi)發(fā)版
要使用最新的開(kāi)發(fā)版本廷痘，請(qǐng)使用以下命令從GitHub安裝：
pip install git+https://github.com/theislab/scvelo
或者
git clone https://github.com/theislab/scvelo
pip install -e scvelo
-e是--editable的縮寫(xiě)，將包鏈接到原始克隆位置件已，這樣拉取的更改也會(huì)反映在環(huán)境中笋额。
- 依賴包
anndata - 帶注釋的數(shù)據(jù)對(duì)象。
scanpy - 用于單細(xì)胞分析的工具包篷扩。
numpy兄猩、scipy、pandas鉴未、scikit-learn枢冤、matplotlib。
部分scVelo（定向 PAGA 和 Louvain 模塊化）需要安裝（可選）：
pip install python-igraph louvain
通過(guò)hnswlib使用快速鄰近搜索進(jìn)一步需要安裝（可選）：
pip install pybind11 hnswlib
3. Getting Started
3.1 導(dǎo)入數(shù)據(jù)
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
分析基于內(nèi)置胰腺數(shù)據(jù)歼狼。
如果要對(duì)自己的數(shù)據(jù)進(jìn)行速率分析掏导，可以通過(guò)adata = scv.read('path/file.loom', cache=True)將文件（loom, h5ad, csv …）讀取到AnnData數(shù)據(jù)對(duì)象。
如果想將loom文件合并到已存在的 AnnData 對(duì)象中羽峰，可以使用scv.utils.merge(adata, adata_loom)趟咆。
adata = scv.datasets.pancreas()
adata
# 查看細(xì)胞類型
adata.obs['clusters'].unique()
scVelo 是基于adata。這個(gè)對(duì)象存儲(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)坯汤。AnnData可以像數(shù)據(jù)框一樣進(jìn)行取切片操作虐唠。例如：adata_subset = adata[:, list_of_gene_names]。更多細(xì)節(jié)可參考anndata docs惰聂。
# 8類細(xì)胞剪切和未剪切的比例
scv.pl.proportions(adata)
這個(gè)圖展示了剪切/未剪切count的比例疆偿。根據(jù)使用單細(xì)胞技術(shù)的不同（Drop-Seq，Smart-Seq）搓幌，數(shù)據(jù)中通常有10%-25%的未剪切分子含有內(nèi)含子序列杆故。作者還建議檢查cluster水平的變化，以驗(yàn)證剪切效率的一致性溉愁。這個(gè)圖顯示处铛，如預(yù)期的變化那樣，在循環(huán)導(dǎo)管細(xì)胞中未剪切的比例略低拐揭，在許多基因開(kāi)始轉(zhuǎn)錄的Ngn3高表達(dá)的和內(nèi)分泌前細(xì)胞中的比例更高撤蟆。
3.2 數(shù)據(jù)預(yù)處理
預(yù)處理包括：
- 通過(guò)檢測(cè)（以最少計(jì)數(shù)檢測(cè)）和高變異性（離散度）進(jìn)行基因選擇，也就是過(guò)濾堂污。
- 按每個(gè)細(xì)胞的初始大小和logX使每個(gè)細(xì)胞標(biāo)準(zhǔn)化家肯。
過(guò)濾是同時(shí)對(duì)spliced/unspliced counts和X進(jìn)行的，Logarithmizing則僅僅對(duì)X進(jìn)行敷鸦。如果X在先前的操作中已經(jīng)被標(biāo)準(zhǔn)化息楔，這里則不會(huì)再次標(biāo)準(zhǔn)化寝贡。
上述的操作在scVelo中都被包裝在一個(gè)函數(shù)中：scv.pp.filter_and_normalize，但實(shí)際上它包含了如下四步：
# 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ì)算的一階矩和二階矩 (first and second order moments) (平均值和去中心化方差)圃泡，這些內(nèi)容封裝在scv.pp.moments。這個(gè)函數(shù)實(shí)際上計(jì)算了scv.pp.pca和scv.pp.neighbors。確定速率估計(jì)需要一階矩，而隨機(jī)估計(jì)需要二階矩撑蒜。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
進(jìn)一步的預(yù)處理（如批次校正）可用于消除不必要的變異源。更多詳細(xì)信息參考最佳實(shí)踐风秤。注意，任何額外的預(yù)處理步驟僅影響count data X扮叨，不應(yīng)用于spliced/unspliced counts缤弦。
3.3 估計(jì)RNA速率
速率是基因表達(dá)空間中的載體，代表單個(gè)細(xì)胞運(yùn)動(dòng)的方向和速率彻磁。無(wú)論是隨機(jī)模型（默認(rèn)）還是確定性模型（通過(guò)設(shè)置mode='deterministic'）碍沐，速率都是通過(guò)對(duì)拼接動(dòng)力學(xué)的轉(zhuǎn)錄動(dòng)力學(xué)進(jìn)行建模獲得。對(duì)于每個(gè)基因衷蜓，預(yù)成熟（未剪切）和成熟（剪切）mRNA計(jì)數(shù)的穩(wěn)定狀態(tài)比累提，構(gòu)成一個(gè)恒定的轉(zhuǎn)錄狀態(tài)。然后磁浇，從此比率中獲取速率作為殘差斋陪。正速率表示基因被上調(diào)，這發(fā)生在細(xì)胞顯示該基因的未剪切mRNA的豐度高于預(yù)期的穩(wěn)定狀態(tài)置吓。相反无虚，負(fù)速表示基因被下調(diào)。
完整的動(dòng)力學(xué)模型的解決方案是通過(guò)設(shè)置mode='dynamical'獲得的衍锚，這需要事先運(yùn)行scv.tl.recover_dynamics(adata)骑科。我們將在下一個(gè)教程中詳細(xì)闡述動(dòng)態(tài)模型。
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中构拳。
查看一下這個(gè)函數(shù)
help(scv.tl.velocity)
adata
然后，基因間速率的組合可用于估計(jì)單個(gè)細(xì)胞的未來(lái)狀態(tài)梁棠。為了將速率投射到低維嵌入中置森，我們估計(jì)了細(xì)胞-細(xì)胞間的轉(zhuǎn)換概率 (transition probabilities) 。也就是說(shuō)符糊，對(duì)于每個(gè)速率矢量凫海，我們尋找了和矢量方向一致的可能的細(xì)胞轉(zhuǎn)換。轉(zhuǎn)換概率是使用潛在細(xì)胞到細(xì)胞的過(guò)渡和速率矢量之間的余弦相關(guān)性 (cosine correlation) 計(jì)算的男娄，并存儲(chǔ)在表示速率圖的矩陣中行贪。由此產(chǎn)生的速率圖具有維度nobs×nobs漾稀，并總結(jié)了通過(guò)速率矢量可以很好地解釋的可能的細(xì)胞狀態(tài)變化（對(duì)于運(yùn)行時(shí)速率，也可以通過(guò)設(shè)置approx=True在減少的 PCA 空間上計(jì)算）建瘫。
scv.tl.velocity_graph(adata)
對(duì)于多種應(yīng)用來(lái)說(shuō)崭捍，通過(guò)應(yīng)用高斯函數(shù)將余弦值相關(guān)性轉(zhuǎn)換為實(shí)際的過(guò)渡概率，速率圖可以轉(zhuǎn)換為transition matrix啰脚。我們可以通過(guò)scv.utils.get_transition_matrix訪問(wèn)這個(gè)馬爾可夫轉(zhuǎn)移矩陣殷蛇。
如前所述，這個(gè)函數(shù)內(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)的軌跡匀们。
3.4 繪制結(jié)果圖
最后，速率被投影到any embedding中 (通過(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 嵌入，并注釋了細(xì)胞群圆存。在應(yīng)用到自己的數(shù)據(jù)時(shí)叼旋，可以使用scv.tl.umap和scv.tl.louvain。詳細(xì)信息參考scanpy tutorial沦辙。此外夫植，所有繪圖功能都默認(rèn)使用basis='umap'和color='clusters'，這些設(shè)置可以根據(jù)需要進(jìn)行更改油讯。
scv.pl.velocity_embedding_stream(adata, basis='umap')
流線顯示的速率矢量可對(duì)發(fā)育過(guò)程進(jìn)行詳細(xì)分析详民。它準(zhǔn)確地勾畫(huà)了導(dǎo)管細(xì)胞和內(nèi)分泌祖細(xì)胞的循環(huán)。此外陌兑，它闡明了細(xì)胞的譜系命運(yùn)沈跨、細(xì)胞周期回歸和內(nèi)分泌細(xì)胞分化狀態(tài)。
圖中的箭頭顯示了單個(gè)細(xì)胞運(yùn)動(dòng)的方向和速率兔综，從未剪切指向剪切饿凛。這個(gè)結(jié)果提示了 Ngn3 細(xì)胞（黃色）的早期內(nèi)分泌命運(yùn)，以及近端α細(xì)胞（藍(lán)色）和瞬態(tài)β細(xì)胞（綠色）之間的明顯分化差異软驰。
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
scv.pl.velocity_embedding_grid(adata, basis='umap')
3.5 解釋速率
這可能是最重要的部分涧窒，因?yàn)槲覀兘ㄗh用戶不要將生物結(jié)論限制在預(yù)測(cè)速率上，而是通過(guò)圖像來(lái)檢查個(gè)體基因動(dòng)力學(xué)锭亏，以了解特定基因如何支持推斷方向纠吴。
查看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)：phase portrait中展示的剪切和未剪切的 mRNA永部。
現(xiàn)在，讓我們來(lái)檢查一些標(biāo)記基因的圖像呐矾，通過(guò)scv.pl.velocity(adata, gene_names)或scv.pl.scatter(adata, gene_names)可視化
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
對(duì)角黑線對(duì)應(yīng)著估計(jì)的"穩(wěn)定狀態(tài)"的未剪切/剪切比率苔埋，即unspliced和spliced的 mRNA 豐度比值。spliced的 mRNA 豐度處于恒定的轉(zhuǎn)錄狀態(tài)蜒犯。特定基因的RNA速率被確定為殘差(residual)组橄，也就是觀察值與穩(wěn)定狀態(tài)線的偏差程度 (見(jiàn)上面鏈接的動(dòng)圖) 。正速率表示基因被向上調(diào)節(jié)罚随，這發(fā)生在細(xì)胞顯示該基因的未剪切mRNA的豐度高于預(yù)期的穩(wěn)定狀態(tài)玉工。相反，負(fù)速率表示基因表達(dá)降低淘菩。
例如遵班，Cpe可以解釋上調(diào)的Ngn3 high EP(黃色)--Pre endocrine -(橙色)--β細(xì)胞(綠色)的分化方向，而Adk則解釋了下調(diào)的Ductal(深綠色)-- Ngn3 high EP(黃色)到remaining endocrine細(xì)胞的方向潮改。
scv.pl.scatter(adata, 'Cpe', color=['clusters', 'velocity'],add_outline='Ngn3 high EP, Pre-endocrine, Beta')
3.6 識(shí)別重要基因
我們需要一種系統(tǒng)的方法來(lái)識(shí)別基因狭郑，這可能有助于解釋由此產(chǎn)生的向量場(chǎng)和推斷的譜系。為此汇在，我們可以測(cè)試哪些基因具有群體特異性微分速率表達(dá)翰萨，與剩余群相比，其比例要高得多/更低糕殉。模塊scv.tl.rank_velocity_genes運(yùn)行差異速率 t-test亩鬼，并生成每個(gè)cluster的基因排名“⒌可以設(shè)置閾值（例如min_corr）辛孵，以限制對(duì)基因候選者選擇的test。
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 高表達(dá)的 EP（黃色）到內(nèi)分泌前（橙色）到Beta（綠色）的方向性。
3.7 循環(huán)祖細(xì)胞的速率（可選宝与，了解一下）
RNA速率檢測(cè)到的細(xì)胞周期焚廊，通過(guò)細(xì)胞周期評(píng)分（相標(biāo)記基因平均表達(dá)水平的標(biāo)準(zhǔn)化分?jǐn)?shù)）在生物學(xué)上得到了證實(shí)冶匹。
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
對(duì)于循環(huán)導(dǎo)管細(xì)胞，我們可以篩選S和G2M期的marker咆瘟。前一個(gè)模塊還計(jì)算了一個(gè)spearmans 相關(guān)分?jǐn)?shù)嚼隘，我們可以用它來(lái)對(duì)phase marker基因進(jìn)行排名/排序，然后顯示它們的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)
特別是Hells 和Top2a非常適合解釋循環(huán)祖細(xì)胞的矢量場(chǎng)飞蛹。Top2a 在 G2M期實(shí)際達(dá)到峰值前不久被分配了高速。在這里灸眼，負(fù)的速率和隨后的下調(diào)完美匹配卧檐。
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
3.8 速率和連貫性
還有兩個(gè)有用的統(tǒng)計(jì)數(shù)據(jù)：
- 分化的速度/速率由速率矢量的長(zhǎng)度給出。
- 矢量場(chǎng)的一致性（即速率矢量如何與其鄰近速率相關(guān)）提供了置信度的衡量標(biāo)準(zhǔn)焰宣。
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
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.
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
3.9 速率圖和擬時(shí)間
我們可以可視化速率圖霉囚，以描繪所有速率推斷的細(xì)胞-細(xì)胞連接/過(guò)渡。它可以通過(guò)設(shè)置threshold限制在高概率轉(zhuǎn)換匕积。例如盈罐，該圖指示來(lái)自早期和晚期內(nèi)分泌細(xì)胞的 Epsilon 細(xì)胞產(chǎn)生的兩個(gè)階段。
scv.pl.velocity_graph(adata, threshold=.1)
此外闪唆，該圖可用于繪制來(lái)自特定細(xì)胞的后代/祖先盅粪。在這里，內(nèi)分泌前細(xì)胞可以追溯到其潛在的命運(yùn)悄蕾。
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í)間笼吟。在從圖形中推斷出 root cells 的分布后库物，它可以計(jì)算從 root cells 開(kāi)始沿著圖形到達(dá)一個(gè)細(xì)胞所需的平均步數(shù)。
與diffusion pseudotime相反贷帮， it implicitly infers the root cells戚揭。而且它是基于directed velocity graph，而不是基于相似性的diffusion kernel撵枢。
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
3.10 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.
# PAGA requires to install igraph, if not done yet.
!pip install python-igraph --upgrade --quiet
# 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).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
reads從左/行到右/列讀取民晒，例如分配了從Ductal到Ngn3 low EP的可信過(guò)渡。
此表可以通過(guò)疊加在 UMAP 嵌入上的定向圖進(jìn)行匯總展示锄禽。
scv.pl.paga(adata, basis='umap', size=50, alpha=.1,
min_edge_width=2, node_size_scale=1.5)
4. Dynamical Modeling
4.1 數(shù)據(jù)準(zhǔn)備
和前面一樣
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.settings.set_figure_params('scvelo') # for beautified visualization
adata = scv.datasets.pancreas()
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
4.2 動(dòng)力學(xué)模型
我們使用動(dòng)力學(xué)模型來(lái)學(xué)習(xí)splicing kinetics的全部轉(zhuǎn)錄動(dòng)態(tài)
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)
# 運(yùn)行動(dòng)力學(xué)模型可能需要一段時(shí)間潜必。因此可以存儲(chǔ)結(jié)果以供重復(fù)使用
#adata.write('data/pancreas.h5ad', compression='gzip')
#adata = scv.read('data/pancreas.h5ad')
scv.pl.velocity_embedding_stream(adata, basis='umap')
4.3 Kinetic rate paramters
無(wú)需實(shí)驗(yàn)即可估算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()
4.4 Latent time（相當(dāng)于軌跡分析）
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)
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)
4.5 Top-likelihood genes
Driver genes display pronounced dynamic behavior and are systematically detected via their characterization by high likelihoods in the dynamic model.
# 查看top15的基因
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)
4.6 Cluster-specific top-likelihood genes
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)
5. Differential Kinetics
One important concern is dealing with systems that represent multiple lineages and processes, where genes are likely to show different kinetic regimes across subpopulations. Distinct cell states and lineages are typically governed by different variations in the gene regulatory networks and may hence exhibit different splicing kinetics. This gives rise to genes that display multiple trajectories in phase space.（這個(gè)是說(shuō)有多個(gè)譜系存在的情況，也是我們分析自己數(shù)據(jù)的時(shí)候最常見(jiàn)到的情況）
為了解決這個(gè)問(wèn)題，可以使用動(dòng)力學(xué)模型來(lái)執(zhí)行差分動(dòng)力學(xué)的似然比測(cè)試垂攘。通過(guò)這種方式维雇，我們可以檢測(cè)到無(wú)法通過(guò)整體動(dòng)力學(xué)的單個(gè)模型很好解釋的動(dòng)力學(xué)行為的集群。將細(xì)胞聚類到其不同的動(dòng)力學(xué)模型中晒他，然后允許分別擬合每個(gè)機(jī)制吱型。
這一部分的演示數(shù)據(jù)集使用的是包含了多種異質(zhì)性亞群的dentate gyrus neurogenesis。
5.1 預(yù)處理和數(shù)據(jù)準(zhǔn)備
和前面一樣
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.settings.set_figure_params('scvelo') # for beautified visualization
adata = scv.datasets.dentategyrus()
scv.pp.filter_and_normalize(adata, min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
5.2 Basic Velocity Estimation
scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis='umap')
5.3 Differential Kinetic Test
不同的細(xì)胞類型和譜系可能表現(xiàn)出不同的動(dòng)力學(xué)機(jī)制陨仅，因?yàn)樗鼈兛梢杂刹煌木W(wǎng)絡(luò)結(jié)構(gòu)控制津滞。即使細(xì)胞類型與譜系相關(guān)，由于可變剪切灼伤、可變多聚腺苷酸環(huán)化和降解調(diào)節(jié)触徐，動(dòng)力也會(huì)有所不同。
動(dòng)力學(xué)模型允許我們通過(guò)差分動(dòng)力學(xué)的似然比測(cè)試來(lái)解決這個(gè)問(wèn)題饺蔑，以檢測(cè)顯示動(dòng)力學(xué)行為的集群/譜系锌介，而整體動(dòng)力學(xué)的的單個(gè)模型無(wú)法充分解釋這些行為。測(cè)試每種細(xì)胞類型是否獨(dú)立擬合產(chǎn)生顯著改善的可能性猾警。
遵循漸進(jìn)卡方分布的似然比可以檢驗(yàn)顯著性孔祸。請(qǐng)注意，出于效率原因发皿，默認(rèn)情況下使用正交回歸而不是完整的相軌跡來(lái)測(cè)試集群是否能被整體動(dòng)力學(xué)很好的解釋或表現(xiàn)出不同的動(dòng)力學(xué)崔慧。
var_names = ['Tmsb10', 'Fam155a', 'Hn1', 'Rpl6']
scv.tl.differential_kinetic_test(adata, var_names=var_names, groupby='clusters')
scv.get_df(adata[:, var_names], ['fit_diff_kinetics', 'fit_pval_kinetics'], precision=2)
kwargs = dict(linewidth=2, add_linfit=True, frameon=False)
scv.pl.scatter(adata, basis=var_names, add_outline='fit_diff_kinetics', **kwargs)
diff_clusters=list(adata[:, var_names].var['fit_diff_kinetics'])
scv.pl.scatter(adata, legend_loc='right', size=60, title='diff kinetics',
add_outline=diff_clusters, outline_width=(.8, .2))
5.4 Testing top-likelihood genes
通過(guò)最高似然基因進(jìn)行篩選，我們發(fā)現(xiàn)了一些顯示多種動(dòng)力學(xué)機(jī)制的基因動(dòng)態(tài)穴墅。
scv.tl.recover_dynamics(adata)
#adata.write('data/pancreas.h5ad', compression='gzip')
#adata = scv.read('data/pancreas.h5ad')
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:100]
scv.tl.differential_kinetic_test(adata, var_names=top_genes, groupby='clusters')
scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, add_outline='fit_diff_kinetics', **kwargs)
scv.pl.scatter(adata, basis=top_genes[15:30], ncols=5, add_outline='fit_diff_kinetics', **kwargs)
5.5 Recompute velocities
最后惶室，可以利用多個(gè)相互競(jìng)爭(zhēng)的動(dòng)力學(xué)機(jī)制的信息重新計(jì)算速度。
scv.tl.velocity(adata, diff_kinetics=True)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding(adata, dpi=120, arrow_size=2, arrow_length=2)
