hello，大家好，今天我們來(lái)分享多組學(xué)數(shù)據(jù)推斷RNA Velocyto撒踪，大家做velocyto的時(shí)候估計(jì)都是用RNA數(shù)據(jù)來(lái)分析的塞颁，其實(shí)浦箱，對(duì)應(yīng)的ATAC數(shù)據(jù)也含有豐富的動(dòng)力學(xué)信息，多組學(xué)推斷更加準(zhǔn)確一點(diǎn)祠锣，參考的文獻(xiàn)在Single-cell multi-omic velocity infers dynamic and decoupled gene regulation酷窥，非常棒的文章，我們先來(lái)分享文章伴网，最后看看示例代碼蓬推。
Abstract
單細(xì)胞多組學(xué)數(shù)據(jù)集，其中在同一細(xì)胞內(nèi)分析多種分子模式澡腾，提供了一個(gè)獨(dú)特的機(jī)會(huì)來(lái)發(fā)現(xiàn)細(xì)胞表觀基因組和轉(zhuǎn)錄組變化之間的關(guān)系（這里我們就是指 RNA + ATAC）沸伏。為了實(shí)現(xiàn)這一目標(biāo)，開(kāi)發(fā)了 MultiVelo动分，這是一種基因表達(dá)機(jī)制模型毅糟，可擴(kuò)展 RNA 速度框架以納入表觀基因組數(shù)據(jù)。 MultiVelo 使用概率潛在變量模型從單細(xì)胞數(shù)據(jù)中估計(jì)染色質(zhì)可及性和基因表達(dá)的轉(zhuǎn)換時(shí)間和速率參數(shù)刺啦，提供表觀基因組和轉(zhuǎn)錄組變化之間時(shí)間關(guān)系的定量總結(jié)留特。與僅從 RNA 估計(jì)的速度相比，結(jié)合染色質(zhì)可及性數(shù)據(jù)顯著提高了細(xì)胞命運(yùn)預(yù)測(cè)的準(zhǔn)確性。在來(lái)自大腦蜕青、皮膚和血細(xì)胞的單細(xì)胞多組學(xué)數(shù)據(jù)集上擬合 MultiVelo 揭示了兩類不同的基因苟蹈，其區(qū)別在于染色質(zhì)是在轉(zhuǎn)錄停止之前還是之后關(guān)閉。MultiVelo模型還識(shí)別了四種類型的細(xì)胞狀態(tài)——兩種表觀基因組和轉(zhuǎn)錄組耦合的狀態(tài)和兩種不同的解耦狀態(tài)右核。 MultiVelo 推斷的參數(shù)量化了基因占據(jù)四種狀態(tài)中每一種狀態(tài)的時(shí)間長(zhǎng)度慧脱，通過(guò)轉(zhuǎn)錄組和表觀基因組之間的耦合程度對(duì)基因進(jìn)行排序。最后贺喝，分析確定了轉(zhuǎn)錄因子表達(dá)和結(jié)合位點(diǎn)可及性之間以及疾病相關(guān) SNP 可及性和相關(guān)基因表達(dá)之間的時(shí)間滯后菱鸥。
Introduction
從 DNA 到 RNA 再到蛋白質(zhì)的基因表達(dá)調(diào)控是控制細(xì)胞命運(yùn)的關(guān)鍵過(guò)程。協(xié)調(diào)的躏鱼、逐步的基因表達(dá)變化——基因按特定順序打開(kāi)和關(guān)閉——是細(xì)胞特化發(fā)育過(guò)程的基礎(chǔ)氮采。 越來(lái)越多的高通量單細(xì)胞測(cè)序技術(shù)被用于揭示這些逐步的基因表達(dá)變化。 但是染苛，由于實(shí)驗(yàn)測(cè)量會(huì)破壞細(xì)胞鹊漠，因此只能進(jìn)行時(shí)間快照測(cè)量，并且不可能觀察到相同的單個(gè)細(xì)胞隨時(shí)間的變化茶行。
計(jì)算方法可以利用單細(xì)胞快照來(lái)推斷發(fā)育過(guò)程中的順序基因表達(dá)變化躯概。例如，細(xì)胞軌跡推斷算法使用成對(duì)的細(xì)胞相似性將細(xì)胞映射到與預(yù)測(cè)的發(fā)育進(jìn)展相對(duì)應(yīng)的“偽時(shí)間”軸上畔师。然而娶靡，基于相似性的軌跡推斷無(wú)法預(yù)測(cè)細(xì)胞轉(zhuǎn)變的方向或相對(duì)速率。推斷 RNA Velocyto的方法通過(guò)擬合微分方程系統(tǒng)來(lái)解決這些限制看锉，該系統(tǒng)使用剪接和未剪接的轉(zhuǎn)錄物計(jì)數(shù)來(lái)描述轉(zhuǎn)錄變化的方向和速率姿锭。最初的 RNA Velocyto方法依賴于穩(wěn)態(tài)假設(shè)來(lái)擬合模型參數(shù)，但后來(lái)的工作開(kāi)發(fā)了一個(gè)動(dòng)力學(xué)模型（就是我們常用的Scvelo） 度陆，除了穩(wěn)態(tài)外艾凯，該模型明確適合基因表達(dá)的誘導(dǎo)和抑制階段献幔。至關(guān)重要的是懂傀，這種 RNA Velocyto的動(dòng)力學(xué)模型還推斷出每個(gè)細(xì)胞的潛在時(shí)間值，提供了一種在細(xì)胞分化過(guò)程中重建基因表達(dá)變化順序的機(jī)械手段蜡感。最近的一篇論文進(jìn)一步擴(kuò)展了 RNA 速度框架蹬蚁，以包括來(lái)自相同細(xì)胞的基因表達(dá)和蛋白質(zhì)測(cè)量，但使用穩(wěn)態(tài)假設(shè)來(lái)估計(jì)參數(shù)郑兴，因此沒(méi)有估計(jì)每個(gè)細(xì)胞的潛伏時(shí)間值犀斋。單細(xì)胞表觀基因組值也被單獨(dú)用于推斷細(xì)胞分化的未來(lái)方向，但這些方法沒(méi)有結(jié)合基因表達(dá)情连。
單細(xì)胞多組學(xué)測(cè)量提供了將表觀基因組數(shù)據(jù)納入轉(zhuǎn)錄機(jī)制模型的機(jī)會(huì)叽粹。 例如，SNARE-seq、SHARE-seq 和 10X Genomics Multiome 等新技術(shù)可以量化同一細(xì)胞中 RNA 和染色質(zhì)的可及性虫几。 表觀基因組和轉(zhuǎn)錄組在細(xì)胞分化過(guò)程中都會(huì)發(fā)生變化锤灿，因此單細(xì)胞多組學(xué)數(shù)據(jù)集中的時(shí)間快照可能揭示這些分子層之間的相互作用。 例如辆脸，如果表觀基因組譜系啟動(dòng)發(fā)生在特定的基因組位點(diǎn)但校，單細(xì)胞多組學(xué)數(shù)據(jù)可能會(huì)揭示基因的染色質(zhì)重塑與其轉(zhuǎn)錄之間存在顯著的時(shí)間滯后。 同樣啡氢，觀察轉(zhuǎn)錄因子表達(dá)和假定結(jié)合位點(diǎn)染色質(zhì)可及性的動(dòng)態(tài)變化可以揭示它們的時(shí)間關(guān)系状囱。
現(xiàn)有的 RNA 速度模型假設(shè)基因的轉(zhuǎn)錄率在基因表達(dá)的整個(gè)誘導(dǎo)階段是一致的。 然而倘是，表觀基因組的變化在調(diào)節(jié)基因表達(dá)方面起著關(guān)鍵作用亭枷，例如收緊或放松啟動(dòng)子和增強(qiáng)子區(qū)域的染色質(zhì)compaction。 例如搀崭，從常染色質(zhì)到異染色質(zhì)的轉(zhuǎn)變會(huì)顯著降低該位點(diǎn)的轉(zhuǎn)錄率奶栖，因?yàn)檗D(zhuǎn)錄機(jī)制無(wú)法訪問(wèn) DNA。 因此门坷，一個(gè)更現(xiàn)實(shí)的模型將反映增強(qiáng)子和啟動(dòng)子染色質(zhì)可及性對(duì)轉(zhuǎn)錄率的影響 宣鄙。
開(kāi)發(fā)了 MultiVelo，這是一種計(jì)算方法默蚌，用于從單細(xì)胞多組學(xué)數(shù)據(jù)集中推斷基因表達(dá)的表觀基因組調(diào)控冻晤。 模型擴(kuò)展了動(dòng)態(tài) RNA 速度模型以結(jié)合多組學(xué)測(cè)量，以更準(zhǔn)確地預(yù)測(cè)每個(gè)細(xì)胞的過(guò)去和未來(lái)狀態(tài)绸吸，共同推斷每種模式的瞬時(shí)誘導(dǎo)或抑制率鼻弧，并確定模式之間的耦合程度或時(shí)間滯后 。MultiVelo 使用概率潛在變量模型來(lái)估計(jì)基因調(diào)控的轉(zhuǎn)換時(shí)間和速率參數(shù)锦茁，提供表觀基因組和轉(zhuǎn)錄組變化之間時(shí)間關(guān)系的定量總結(jié)攘轩。
示例證明 MultiVelo 準(zhǔn)確地恢復(fù)了細(xì)胞譜系，并量化了染色質(zhì)可及性和基因表達(dá)暫時(shí)不同步的引發(fā)和去耦間隔的長(zhǎng)度码俩。 微分方程模型準(zhǔn)確地?cái)M合了來(lái)自胚胎小鼠大腦度帮、胚胎人類大腦的單細(xì)胞多組學(xué)數(shù)據(jù)集，以及來(lái)自人類造血干細(xì)胞和祖細(xì)胞的新生成的數(shù)據(jù)集。 此外，模型通過(guò)染色質(zhì)可及性預(yù)測(cè)了兩種不同的基因表達(dá)調(diào)控機(jī)制精钮，并且在研究的所有組織中確定了兩種機(jī)制的明確例子鼓拧。 最后，使用 MultiVelo 來(lái)推斷轉(zhuǎn)錄因子 (TF) 與其結(jié)合位點(diǎn)之間以及 GWAS SNP 與其連鎖基因之間的時(shí)間關(guān)系。 總之，MultiVelo 提供了對(duì)表觀基因組變化在細(xì)胞命運(yùn)轉(zhuǎn)變過(guò)程中調(diào)節(jié)基因表達(dá)的機(jī)制的基本見(jiàn)解。
Results
MultiVelo: A Mechanistic Model of Gene Expression Incorporating Chromatin Accessibility（這部分內(nèi)容涉及到很多數(shù)學(xué)符號(hào)冕臭，就不給大家翻譯了）
也就是說(shuō)估計(jì)三部分腺晾，ATAC，pre-RNA辜贵，RNA丘喻。
MultiVelo 模型的數(shù)學(xué)公式立即導(dǎo)致了關(guān)于基因表達(dá)過(guò)程中染色質(zhì)可及性與轉(zhuǎn)錄之間關(guān)系的兩個(gè)重要見(jiàn)解。 首先念颈，染色質(zhì)可及性和 RNA 轉(zhuǎn)錄狀態(tài)有多種數(shù)學(xué)上可行的組合泉粉。 也就是說(shuō)，當(dāng)轉(zhuǎn)錄被誘導(dǎo)或抑制時(shí)榴芳，染色質(zhì)可以打開(kāi)或關(guān)閉嗡靡。 這意味著多個(gè)事件順序是可能的：染色質(zhì)關(guān)閉可以在轉(zhuǎn)錄抑制開(kāi)始之前或之后發(fā)生。 將第一個(gè)排序（染色質(zhì)關(guān)閉在轉(zhuǎn)錄抑制之前開(kāi)始）稱為模型 1窟感，將第二個(gè)排序稱為模型 2讨彼。Note that there are other mathematically possible orderings where transcription occurs before chromatin opening, but these are not biologically plausible, and we do not find convincing evidence that they occur in the datasets we analyzed。
MultiVelo Accurately Fits Simulated Data
進(jìn)行了模擬柿祈，以確定 MultiVelo 是否可以恢復(fù)速率參數(shù)和切換時(shí)間哈误，并在存在噪聲的情況下區(qū)分模型 1 和模型 2。 結(jié)果表明躏嚎，MultiVelo 可以準(zhǔn)確擬合噪聲數(shù)據(jù)蜜自，并且可以恢復(fù)底層參數(shù)。 此外卢佣，發(fā)現(xiàn) MultiVelo 以高精度區(qū)分模型 1 和模型 2（98.5% 的模擬基因根據(jù)模型似然被正確分配）重荠。 分析還證實(shí)，在擬合 ODE 參數(shù)之前虚茶，可以通過(guò)簡(jiǎn)單地比較穩(wěn)態(tài)線上方和下方頂部分位數(shù)中的細(xì)胞數(shù)量來(lái)區(qū)分模型 1 和模型 2 基因（95.8% 的模擬基因被正確分配 ) 戈鲁。
MultiVelo Distinguishes Two Models of Gene Expression Regulation in Embryonic Mouse Brain
首先將 MultiVelo 應(yīng)用于來(lái)自胚胎小鼠大腦 (E18) 的 10X Multiome 數(shù)據(jù)。 MultiVelo 準(zhǔn)確擬合了觀察到的染色質(zhì)可及性嘹叫、未剪接的前 mRNA 和剪接的 mRNA 計(jì)數(shù)在整個(gè)腦細(xì)胞群中婆殿，識(shí)別出 426 個(gè)基因，其模式極有可能符合模型罩扇。 MultiVelo 推斷的所得速度向量和潛在時(shí)間值準(zhǔn)確地恢復(fù)了哺乳動(dòng)物皮層發(fā)育的已知軌跡婆芦。 具體來(lái)說(shuō)，外腦室下區(qū) (OSVZ) 中的徑向神經(jīng)膠質(zhì) (RG) 細(xì)胞會(huì)產(chǎn)生神經(jīng)元暮蹂、星形膠質(zhì)細(xì)胞和少突膠質(zhì)細(xì)胞寞缝。 皮質(zhì)層在神經(jīng)元遷移過(guò)程中以由內(nèi)而外的方式形成，新生細(xì)胞移動(dòng)到上層仰泻，老細(xì)胞留在更深的層。 RG細(xì)胞可以分裂成作為神經(jīng)干細(xì)胞的中間祖細(xì)胞（IPC）滩届，并在不同層進(jìn)一步產(chǎn)生各種成熟的興奮性神經(jīng)元集侯。
與 scVelo 等僅使用 RNA 的模型相比被啼，結(jié)合染色質(zhì)可及性和基因表達(dá)提高了velocyto估計(jì)的準(zhǔn)確性。 特別是棠枉，僅 RNA 模型預(yù)測(cè)了上層神經(jīng)元內(nèi)在生物學(xué)上不可信的回流浓体。 細(xì)胞周期評(píng)分表明發(fā)育過(guò)程始于 RG 附近的循環(huán) population，證實(shí)了 MultiVelo 推斷的延遲時(shí)間辈讶。 MultiVelo 和 scVelo 使用類似的參數(shù)設(shè)置和估計(jì)算法命浴，表明表觀基因組數(shù)據(jù)提供了有關(guān)細(xì)胞過(guò)去和未來(lái)狀態(tài)的重要附加信息，而不僅僅是轉(zhuǎn)錄組數(shù)據(jù)提供的信息贱除。
分析預(yù)計(jì)染色質(zhì)可及性的添加最有助于區(qū)分染色質(zhì)重塑和基因表達(dá)不同步的細(xì)胞狀態(tài)生闲，例如當(dāng)基因的啟動(dòng)子和增強(qiáng)子開(kāi)始打開(kāi)但幾乎沒(méi)有發(fā)生轉(zhuǎn)錄時(shí)。兩個(gè)明顯的例子是 Eomes 和 Tle4月幌，它們是 IPC 和深層神經(jīng)元的典型標(biāo)記碍讯。這些基因的 RNA 轉(zhuǎn)錄物僅在一種或兩種特定細(xì)胞類型中高度表達(dá)。剩余的細(xì)胞密集地聚集在 (u, s) 相圖的原點(diǎn)附近扯躺，使得 RNA 速度方法難以區(qū)分它們的相對(duì)順序捉兴。然而，這些基因的染色質(zhì)可及性在基因表達(dá)之前開(kāi)始上升录语，揭示了僅從基因表達(dá)中看不到的逐漸變化倍啥。換句話說(shuō)，結(jié)合染色質(zhì)使我們能夠推斷出 3D 速度向量澎埠，指示每個(gè)細(xì)胞對(duì)每個(gè)基因的預(yù)測(cè)分化逗栽，比單獨(dú)來(lái)自 RNA 的 2D 相圖更好地解決細(xì)胞差異。
MultiVelo 確定了在此數(shù)據(jù)集中由模型 1 和模型 2 最好地描述的基因的清晰示例失暂。 比較分配給模型 1 和模型 2 的基因的相圖顯示最大染色質(zhì)可及性的時(shí)間存在明顯差異彼宠，與模型預(yù)測(cè)一致。 模型 1 基因（如 Satb2）在轉(zhuǎn)錄誘導(dǎo)階段（相圖上的對(duì)角穩(wěn)態(tài)線以上）達(dá)到最大染色質(zhì)可及性弟塞，而模型 2 基因（如 Gria2）在轉(zhuǎn)錄抑制階段（低于對(duì)角穩(wěn)態(tài)線）的可及性最高 -狀態(tài)線）凭峡。 在檢查 c, u 和 c, s 的成對(duì)相圖時(shí)，模型 1 和模型 2 之間的區(qū)別也很明顯决记。 但是摧冀，不能通過(guò)單獨(dú)檢查 u, s 的相圖中的 RNA 信息來(lái)區(qū)分模型； 區(qū)分需要來(lái)自染色質(zhì)的額外信息系宫。
進(jìn)一步研究了模型 1 和模型 2 基因索昂，看看它們是否具有任何特征。 Gene ontology（GO）分析表明扩借，M2基因顯著富集了與細(xì)胞周期相關(guān)的術(shù)語(yǔ)椒惨，如“細(xì)胞周期的正調(diào)控”、“有絲分裂細(xì)胞周期”和“細(xì)胞周期相變的調(diào)控”潮罪。 此外康谆，模型 2 基因往往比模型 1 基因更早地在潛伏時(shí)間內(nèi)達(dá)到最高剪接表達(dá)（p = 9 × 10?7领斥，Wilcoxon 秩和單邊檢驗(yàn)）。 假設(shè)細(xì)胞可以使用模型 2 來(lái)快速沃暗、瞬時(shí)激活不需要維持表達(dá)的基因月洛，而模型 1 可能對(duì)需要穩(wěn)定表達(dá)的基因有用。
接下來(lái)研究了每種類型的基因表達(dá)動(dòng)力學(xué)（僅誘導(dǎo)孽锥、僅抑制嚼黔、模型 1 或模型 2）發(fā)生的頻率。 大多數(shù)高度可變的基因同時(shí)顯示誘導(dǎo)和抑制階段（完整軌跡）惜辑，對(duì)于只有部分軌跡的基因唬涧，僅誘導(dǎo)相圖出現(xiàn)的頻率高于僅抑制相圖（29.5% 對(duì) 2.4% 的可變基因。Note韵丑，因?yàn)槟Ｐ?1 和模型 2 在誘導(dǎo)階段做出了相同的預(yù)測(cè)爵卒，無(wú)法區(qū)分模型 1 和模型 2 對(duì)于僅誘導(dǎo)基因。 1（41.4%的可變基因）撵彻，而其余的最適合模型2（26.7%的可變基因）模型1更常見(jiàn)的事實(shí)與染色質(zhì)狀態(tài)變化通常先于mRNA表達(dá)變化的預(yù)期一致钓株。
無(wú)論基因是否具有完整或部分動(dòng)力學(xué)，MultiVelo 都適合描述其染色質(zhì)可及性和基因表達(dá)動(dòng)力學(xué)的三維軌跡的 ODE 參數(shù)陌僵。 通過(guò)對(duì)隨時(shí)間變化的轉(zhuǎn)錄率進(jìn)行建模轴合，MultiVelo 能夠更好地捕捉 RNA 相圖中不同類型的曲率，而純 RNA 模型無(wú)法捕捉到這種曲率差異碗短。 具有不同模型分配和動(dòng)力學(xué)的基因在可能性或總計(jì)數(shù)方面沒(méi)有表現(xiàn)出顯著差異受葛，表明technical artifacts無(wú)法解釋這些現(xiàn)象。
MultiVelo Identifies Epigenomic Priming and Decoupling in Embryonic Mouse Brain
MultiVelo 的一個(gè)重要的特性是它能夠量化染色質(zhì)可及性和分化細(xì)胞內(nèi)基因表達(dá)之間的不一致和一致偎谁。 具體來(lái)說(shuō)总滩，MultiVelo 推斷開(kāi)關(guān)時(shí)間參數(shù)，這些參數(shù)確定每個(gè)基因處于四種可能狀態(tài)（啟動(dòng)巡雨、耦合開(kāi)啟闰渔、解耦和耦合關(guān)閉）之一的間隔。 接下來(lái)研究了這些推斷的狀態(tài)和時(shí)間間隔是否可以準(zhǔn)確地捕捉胚胎小鼠腦細(xì)胞中表觀基因組和轉(zhuǎn)錄組變化之間的相互作用铐望。
MultiVelo 識(shí)別了 10X Multiome 數(shù)據(jù)中四種狀態(tài)中每一種的清晰示例冈涧。 例如，Grin2b 是一種僅誘導(dǎo)基因正蛙，其表達(dá)朝著神經(jīng)元命運(yùn)增加督弓，因此僅預(yù)測(cè)了該基因的誘導(dǎo)狀態(tài)——啟動(dòng)和耦合。 模型 1 基因 Nfix 的相圖具有完整的軌跡形狀乒验，并標(biāo)有所有四種狀態(tài)愚隧。 相反，Epha5 是 Model 2 基因徊件，它的可及性在整個(gè)時(shí)間范圍內(nèi)持續(xù)上升奸攻，沒(méi)有觀察到關(guān)閉階段蒜危，所以它只占據(jù)耦合開(kāi)啟和解耦狀態(tài) 虱痕。
通過(guò)在 UMAP 坐標(biāo)上繪制可訪問(wèn)性 (c) 和表達(dá)式（u 和 s）并并排檢查它們睹耐，可以定性地確認(rèn)狀態(tài)分配。在視覺(jué)上部翘，觀察到 c 和 u UMAP 圖的顏色在狀態(tài)分配耦合打開(kāi)或耦合關(guān)閉時(shí)匹配硝训，并且在分配狀態(tài)啟動(dòng)或解耦時(shí)出現(xiàn)顏色差異。例如新思，Robo2 RNA 表達(dá)和染色質(zhì)可及性之間的最大差異發(fā)生在圓圈區(qū)域窖梁，該區(qū)域預(yù)計(jì)處于解耦狀態(tài)。 Robo2 是一個(gè) Model 1 基因夹囚；染色質(zhì)閉合開(kāi)始后纵刘，表達(dá)保持在相對(duì)較高的水平，即使它的可及性已經(jīng)經(jīng)歷了成熟神經(jīng)元的下降荸哟。同樣假哎，Gria2 的可訪問(wèn)性與解耦狀態(tài)的 RNA 不同。模型 2 基因 Gria2 的染色質(zhì)可及性在轉(zhuǎn)錄誘導(dǎo)期之后繼續(xù)增加鞍历。此外舵抹，基因 Grin2b 顯示了染色質(zhì)引發(fā)階段的一個(gè)明顯例子，在此期間染色質(zhì)在 RNA 產(chǎn)生之前打開(kāi)劣砍。
沿著每個(gè)基因的推斷時(shí)間 t 繪制 c惧蛹、u 和 s 能夠詳細(xì)檢查狀態(tài)轉(zhuǎn)換。首先刑枝，Robo2 的 u(t) 和 s(t) 值在轉(zhuǎn)錄抑制階段顯示出兩個(gè)拐點(diǎn)香嗓，對(duì)應(yīng)于從耦合到解耦狀態(tài)和從解耦到耦合關(guān)閉狀態(tài)的轉(zhuǎn)變。這種模式表明染色質(zhì)關(guān)閉和轉(zhuǎn)錄抑制的不同影響在 u(t) 和 s(t) 中可見(jiàn)装畅。換句話說(shuō)靠娱，MultiVelo 預(yù)測(cè)對(duì)于 Robo2，染色質(zhì)關(guān)閉會(huì)降低整體轉(zhuǎn)錄率洁灵，因?yàn)?RNA 水平在染色質(zhì)轉(zhuǎn)換后立即開(kāi)始下降饱岸。隨后轉(zhuǎn)錄率從正值到零的轉(zhuǎn)換導(dǎo)致第二個(gè)拐點(diǎn)，導(dǎo)致 RNA 表達(dá)的更快速下調(diào)徽千。 Gria2 的 c(t)苫费、u(t) 和 s(t) 的圖顯示出相反的趨勢(shì)：即使在切換到轉(zhuǎn)錄抑制后，c 仍繼續(xù)上升双抽，導(dǎo)致 c 和 u 在解耦狀態(tài)期間向相反的方向移動(dòng).在 Grin2b 的長(zhǎng)啟動(dòng)階段百框，c(t) 開(kāi)始上升，而 u(t) 和 s(t) 保持為零 牍汹。
由于 MultiVelo 適合每個(gè)基因的速率和轉(zhuǎn)換時(shí)間參數(shù)铐维，分析提供了觀察基因調(diào)控總體趨勢(shì)的機(jī)會(huì)柬泽。首先，為了確定不同基因的狀態(tài)是否在時(shí)間上是協(xié)調(diào)的嫁蛇，計(jì)算了每個(gè)細(xì)胞每個(gè)狀態(tài)中高似然基因的數(shù)量锨并。確實(shí)存在通過(guò)神經(jīng)元cluster的級(jí)聯(lián)狀態(tài)轉(zhuǎn)換；每個(gè)細(xì)胞的多個(gè)基因通常同時(shí)處于啟動(dòng)或去耦狀態(tài)睬棚。其次第煮，尋找轉(zhuǎn)換時(shí)間和速率參數(shù)的趨勢(shì)。將每個(gè)基因的誘導(dǎo)/抑制周期置于 0 到 1 之間的時(shí)間尺度上抑党，發(fā)現(xiàn)耦合開(kāi)啟和耦合關(guān)閉狀態(tài)比啟動(dòng)和解耦狀態(tài)占基因表達(dá)過(guò)程的更大比例包警。這是有道理的，因?yàn)榧词够蛟趦煞N模式之間經(jīng)歷某種程度的解耦和時(shí)間滯后底靠，染色質(zhì)可及性和基因表達(dá)仍應(yīng)普遍相關(guān)害晦。中位數(shù)啟動(dòng)間隔長(zhǎng)度為總時(shí)間的 21%，中位數(shù)解耦間隔長(zhǎng)度為總時(shí)間的 19%暑中。此外壹瘟，可以通過(guò)啟動(dòng)和解耦間隔的時(shí)間對(duì)基因進(jìn)行排序，以找到可訪問(wèn)性和表達(dá)之間不一致的例子痒芝。此外俐筋，發(fā)現(xiàn)染色質(zhì)通常以相似的速度打開(kāi)和關(guān)閉：推斷的染色質(zhì)關(guān)閉率 (αcc) 和染色質(zhì)打開(kāi)率 (αco) 之間的中值比率幾乎正好是 1。
MultiVelo Quantifies Epigenomic Priming in SHARE-seq Data from Mouse Hair Follicle
最近的一項(xiàng)研究使用 SHARE-seq 來(lái)研究毛囊組織中轉(zhuǎn)運(yùn)放大細(xì)胞 (TAC) 的快速增殖严衬，這些細(xì)胞會(huì)產(chǎn)生幾種成熟的效應(yīng)細(xì)胞澄者，包括內(nèi)根鞘 (IRS) 和毛干層：角質(zhì)層、皮質(zhì) 層和髓質(zhì)请琳。 當(dāng)應(yīng)用于該數(shù)據(jù)集時(shí)粱挡，MultiVelo 正確識(shí)別了從 TAC 到 IRS 和毛干細(xì)胞的分化方向，與最初論文中報(bào)告的擴(kuò)散圖分析一致俄精。 潛伏時(shí)間預(yù)測(cè) TAC 是根細(xì)胞——與生物學(xué)預(yù)期一致——而單獨(dú)使用 RNA 的速度分析未能捕捉到毛干分化方向询筏。 與 mouse brai相比，在該數(shù)據(jù)集中觀察到明顯更多的僅誘導(dǎo)基因和更少的模型 2 基因竖慧。
最初的 SHARE-seq 論文的主要結(jié)果之一是鑒定了啟動(dòng)子和增強(qiáng)子染色質(zhì)可及性預(yù)示基因表達(dá)的基因嫌套，作者將這種現(xiàn)象稱為“染色質(zhì)潛力”。這種現(xiàn)象最明顯的例子是 Wnt3圾旨，它編碼旁分泌信號(hào)分子踱讨，在控制頭發(fā)生長(zhǎng)方面很重要。事實(shí)上砍的，UMAP 圖按可及性著色痹筛，未剪接和剪接的 mRNA 表達(dá)顯示出明顯的時(shí)間延遲。接下來(lái)檢查了 SHARE-seq 論文中確定的其他基因。擬合模型顯示 MultiVelo 忠實(shí)地捕捉了每個(gè)基因的動(dòng)態(tài)帚稠，并清楚地說(shuō)明了啟動(dòng)和去耦區(qū)域谣旁。例如，Wnt3 和 Dsc1 顯示僅誘導(dǎo)模式和開(kāi)始時(shí)的啟動(dòng)狀態(tài)滋早，而 Cux1榄审、Dlx3 和 Cobll1 具有誘導(dǎo)和抑制狀態(tài) 中間有一個(gè)短的解耦期。
為了進(jìn)一步量化可及性馆衔、未拼接表達(dá)和拼接表達(dá)之間的時(shí)間關(guān)系瘟判，使用動(dòng)態(tài)時(shí)間扭曲 (DTW) 來(lái)對(duì)齊每個(gè)分子層的時(shí)間序列值怨绣。 DTW 非線性扭曲兩個(gè)時(shí)間序列以最大化它們的相似性并識(shí)別可能的滯后相關(guān)性角溃。 Wnt3 上的 DTW 結(jié)果表明，最佳扭曲函數(shù)在時(shí)間上向前映射 c 時(shí)間序列上的每個(gè)點(diǎn)篮撑，與基因表達(dá)之前的染色質(zhì)可及性一致减细。未拼接和拼接的表達(dá)顯示出相似的模式，但延遲時(shí)間較短赢笨。由于 DTW 將前面曲線上的每個(gè)時(shí)間點(diǎn)映射到后面曲線上的時(shí)間點(diǎn)未蝌，因此可以通過(guò)減去匹配點(diǎn)的時(shí)間來(lái)計(jì)算每個(gè)時(shí)間點(diǎn)的時(shí)滯。該分析表明茧妒，c 和 s 之間的延遲以及 u 和 s 之間的延遲在整個(gè)觀察時(shí)間內(nèi)都保持為正值萧吠。此外，在整個(gè)觀測(cè)范圍內(nèi)桐筏，c 和 s 之間的延遲比 u 和 s 之間的延遲長(zhǎng)纸型，最大 c 和 s 延遲達(dá)到 0.6（在總時(shí)間范圍 1 之外）。
MultiVelo Reveals Early Epigenomic and Transcriptomic Changes in Human Hematopoietic Stem and Progenitor Cells
造血祖細(xì)胞由干細(xì)胞樣細(xì)胞群組成梅忌，這些細(xì)胞群在進(jìn)入更多譜系受限狀態(tài)時(shí)狰腌，可快速且持續(xù)地分化為各種中等和成熟血細(xì)胞類型，自我更新能力逐漸降低
將純化的人類 CD34+ 細(xì)胞培養(yǎng)了 7 天牧氮，然后使用 10X Multiome 平臺(tái)對(duì)其進(jìn)行測(cè)序琼腔。 使用單核 RNA-seq 和 ATAC-seq 數(shù)據(jù)過(guò)濾后獲得了 11,605 個(gè)高質(zhì)量細(xì)胞。 使用先前描述的標(biāo)記基因踱葛，確定了類似于許多早期血液發(fā)育群體的cluster丹莲，包括 HSC、多能祖細(xì)胞 (MPP)尸诽、淋巴引發(fā)的多能祖細(xì)胞 (LMPP)甥材、粒細(xì)胞-巨噬細(xì)胞祖細(xì)胞 (GMP) 和巨核細(xì)胞- 紅細(xì)胞祖細(xì)胞 (MEP)。 還確定了類似于早期粒細(xì)胞逊谋、紅細(xì)胞擂达、樹(shù)突細(xì)胞 (DC) 和血小板的cluster。
血細(xì)胞分化是用 RNA 速度建模的一個(gè)具有挑戰(zhàn)性的系統(tǒng)，但發(fā)現(xiàn)結(jié)合染色質(zhì)信息顯著提高了預(yù)測(cè)細(xì)胞方向的局部一致性和生物學(xué)準(zhǔn)確性板鬓。 相比之下悲敷，僅從 RNA 推斷的速度向量并不能準(zhǔn)確反映 HSPC 的已知分化層次。 與小鼠大腦一樣俭令，MultiVelo 預(yù)測(cè)模型 1 在該數(shù)據(jù)集中比模型 2 更常見(jiàn)后德； 僅歸納是第三個(gè)最常見(jiàn)的基因類。 觀察到的啟動(dòng)和解耦區(qū)間的中值長(zhǎng)度比耦合階段的短抄腔。 這些模式與我們?cè)谛∈蟠竽X數(shù)據(jù)集中觀察到的一致瓢湃，表明可能存在共同的潛在生物學(xué)機(jī)制 。
與小鼠大腦數(shù)據(jù)集一樣赫蛇，HSPC 數(shù)據(jù)集中的模型 2 基因顯著豐富了與細(xì)胞周期相關(guān)的 GO terms绵患。 terms“有絲分裂細(xì)胞周期的調(diào)節(jié)”、“有絲分裂中期/后期轉(zhuǎn)變的調(diào)節(jié)”和“有絲分裂姐妹染色單體分離的調(diào)節(jié)”都富集在模型2基因中悟耘，F(xiàn)DR < 0.002落蝙。 如果我們檢查髓系、紅細(xì)胞系和血小板譜系的不同軌跡暂幼，許多 G2/M 期標(biāo)記基因 18 顯示清晰的模型 2 模式筏勒，在表達(dá)開(kāi)始下降后染色質(zhì)可及性最高 。
進(jìn)一步研究了模型 1 和模型 2 基因的組蛋白修飾譜是否不同旺嬉。 由于可以使用 FACS 對(duì)經(jīng)典定義的 HSPC 亞群進(jìn)行分類管行，因此我們分析中的某些細(xì)胞亞群可以獲得大量 ChIP-seq 數(shù)據(jù)。 使用這些大量數(shù)據(jù)集邪媳，我們比較了 FACS 純化的 HSC 中 H3K4me3捐顷、H3K4me1 和 H3K27ac 的水平，染色質(zhì)可及性峰值與模型 1 和模型 2 基因相關(guān)悲酷。 我們發(fā)現(xiàn)模型 2 基因顯示出顯著更高的 H3K4me3（p = 0.016套菜，單側(cè) Wilcoxon 秩和檢驗(yàn)），這是活躍啟動(dòng)子的標(biāo)志设易。 相比之下逗柴，模型 1 基因顯示出稍高的 H3K4me1 (p = 0.097)，這是一種啟動(dòng)的增強(qiáng)子標(biāo)記顿肺。 兩種模型在 HSC 中都顯示出相似的 H3K27ac（一種活性增強(qiáng)子標(biāo)記物）（p = 0.48）戏溺。
MultiVelo 擬合的基因模型揭示了許多啟動(dòng)的例子。 幾種終末細(xì)胞類型特異性標(biāo)記物顯示出僅誘導(dǎo)動(dòng)力學(xué)屠尊，染色質(zhì)可及性增加旷祸，隨后基因表達(dá)增加（GMP 中的 AZU1、紅細(xì)胞中的 HBD讼昆、粒細(xì)胞中的 HDC托享、DC 祖細(xì)胞中的 LYZ 和巨核細(xì)胞 (MK) 祖細(xì)胞中的 PF4） 方向）。 在 HSPC 中，我們?cè)俅慰吹揭恍┟黠@的長(zhǎng)啟動(dòng)期示例闰围，例如 LYZ 和 PF4
繪制Velocyto使我們能夠更詳細(xì)地檢查局部染色質(zhì)和 RNA 趨勢(shì)赃绊。 雖然染色質(zhì)在這些基因開(kāi)始時(shí)顯示出最大的潛力（最高速度），但對(duì)于 RNA羡榴，干細(xì)胞群（如 HSC碧查、MPP、MEP 和 GMP）在向一個(gè)譜系的分化過(guò)程中顯示出增加的潛力校仑。 更多分化的細(xì)胞類型失去了保持這種潛力的能力忠售，并逐漸接近平衡（零速度），即使表達(dá)仍在有所增加迄沫。 請(qǐng)注意稻扬，即使整體表達(dá)式升高且速度保持為正，局部加速度仍可以切換符號(hào)邢滑。 由于染色質(zhì)和 mRNA 的聯(lián)合數(shù)學(xué)建模腐螟，MultiVelo 能夠捕獲關(guān)于分化方向和速率的如此豐富的信息。 添加染色質(zhì)顯著豐富了從 RNA 中可用的信息困后，這可以通過(guò)檢查僅 RNA 相圖看出。
MultiVelo Relates Transcription Factors, Polymorphic Sites, and Gene Expression in Developing Human Brain
接下來(lái)將 MultiVelo 應(yīng)用到最近發(fā)布的來(lái)自開(kāi)發(fā)人類皮層的 10X Multiome 數(shù)據(jù)集衬廷。 與胚胎小鼠大腦數(shù)據(jù)集一樣摇予，MultiVelo 推斷出的速度向量與已知的腦細(xì)胞發(fā)育模式一致。 MultiVelo 正確地推斷出放射狀膠質(zhì)細(xì)胞附近的循環(huán)細(xì)胞群是潛伏時(shí)間最早的細(xì)胞類型吗跋。 相比之下侧戴，在沒(méi)有染色質(zhì)信息的情況下推斷的速度向量預(yù)測(cè)了中間祖細(xì)胞和上層興奮性神經(jīng)元的不協(xié)調(diào)回流。
與小鼠大腦數(shù)據(jù)集一樣跌宛，確定了模型 1 和模型 2 基因的清晰示例酗宋，盡管在人類數(shù)據(jù)集中預(yù)測(cè)遵循模型 2 的基因較少。 有趣的是疆拘，模型 2 基因 MEF2C 被僅 RNA 模型預(yù)測(cè)為具有主要的抑制相蜕猫，可能是因?yàn)?u - s 相圖的“寬度”很窄。然而哎迄，染色質(zhì)信息的添加允許正確的預(yù)測(cè)該基因同時(shí)具有誘導(dǎo)期和抑制期回右。
MultiVelo 的一個(gè)關(guān)鍵優(yōu)勢(shì)是它能夠?qū)⒓?xì)胞置于從染色質(zhì)和表達(dá)數(shù)據(jù)推斷出的潛伏時(shí)間尺度上。我們推斷漱挚，潛伏時(shí)間可以識(shí)別基因座的表達(dá)和可及性之間的時(shí)間滯后翔烁，而不僅僅是那些緊鄰基因的基因座。例如旨涝，潛伏時(shí)間可用于計(jì)算轉(zhuǎn)錄因子 (TF) 的表達(dá)與其結(jié)合位點(diǎn)的可及性之間的時(shí)間長(zhǎng)度蹬屹。為此，我們使用 chromVar40 計(jì)算每個(gè)細(xì)胞的峰與每個(gè) TF 的結(jié)合位點(diǎn)的總可及性，子集僅針對(duì)數(shù)據(jù)集中可變表達(dá)的 TF慨默。然后我們使用動(dòng)態(tài)時(shí)間扭曲 (DTW) 將每個(gè) TF 的時(shí)間序列表達(dá)與其綁定站點(diǎn)的可訪問(wèn)性對(duì)齊秃踩。這揭示了一個(gè)一致的模式，其中轉(zhuǎn)錄因子最高 RNA 表達(dá)的時(shí)間先于相應(yīng)的下游靶標(biāo)的高可及性時(shí)間业筏。由 TF 表達(dá)和結(jié)合位點(diǎn)可訪問(wèn)性著色的 UMAP 圖在視覺(jué)上證實(shí)了這種模式憔杨。所有表達(dá)的 TF 的中位時(shí)間滯后為正，表明在大多數(shù)情況下 TF 表達(dá)先于結(jié)合位點(diǎn)可訪問(wèn)性蒜胖。如果沒(méi)有額外的數(shù)據(jù)消别，我們無(wú)法最終確定這些時(shí)間滯后背后的機(jī)制。然而台谢，轉(zhuǎn)錄后和翻譯后調(diào)節(jié)寻狂、影響染色質(zhì)重塑復(fù)合物活性的因素以及細(xì)胞間信號(hào)傳導(dǎo)都可能導(dǎo)致這種現(xiàn)象。
MultiVelo 推斷的潛伏時(shí)間也可用于將疾病相關(guān)變異基因座的染色質(zhì)可及性與附近基因的表達(dá)聯(lián)系起來(lái)朋沮。我們收集了 6968 個(gè)單核苷酸多態(tài)性 (SNP) 及其關(guān)聯(lián)基因的列表蛇券，這些基因與精神疾病（包括雙相情感障礙和精神分裂癥）的全基因組關(guān)聯(lián)研究有關(guān)樊拓。我們進(jìn)一步將這些 SNP 子集到與我們模型擬合的基因相關(guān)的重疊染色質(zhì)可及性峰纠亚，總共 757 個(gè) SNP。許多這些變異發(fā)生在神經(jīng)元轉(zhuǎn)錄因子和其他發(fā)育重要基因附近筋夏。然后我們計(jì)算了 400 b.p. 每個(gè)細(xì)胞的染色質(zhì)可及性蒂胞。窗口以每個(gè) SNP 為中心。使用 MultiVelo 的潛伏時(shí)間条篷，我們確定了每個(gè) SNP 的最大可達(dá)性時(shí)間以及 SNP 可達(dá)性與其連鎖基因的最大表達(dá)之間的時(shí)間滯后骗随。該分析揭示了 3 個(gè)主要的 SNP 組，根據(jù)它們的最大可接近性是在潛伏時(shí)間的早期還是晚期以及相關(guān)基因表達(dá)之前或之后進(jìn)行區(qū)分赴叹。 SNP 可及性和相關(guān)基因表達(dá)的 UMAP 圖證實(shí)這些 SNP 組具有不同的性質(zhì)鸿染。這些分組對(duì)于理解 SNP 的功能很重要；例如乞巧，與完全分化的細(xì)胞相比涨椒，僅在潛伏期早期才可接近的 SNP 在發(fā)育細(xì)胞中可能發(fā)揮更大的作用。同樣摊欠，可接近性先于基因表達(dá)的 SNP 比可接近性落后的 SNP 更有可能參與調(diào)節(jié)其表達(dá)丢烘。
Discussion
總之，MultiVelo 準(zhǔn)確地恢復(fù)了細(xì)胞譜系些椒，并量化了染色質(zhì)可及性和基因表達(dá)暫時(shí)不同步的引發(fā)和解偶聯(lián)間隔的長(zhǎng)度播瞳。 我們的模型準(zhǔn)確地?cái)M合了來(lái)自胚胎小鼠大腦、小鼠背部皮膚免糕、胚胎人類大腦和人類造血干細(xì)胞的單細(xì)胞多組學(xué)數(shù)據(jù)集赢乓。 此外忧侧，我們的模型確定了兩類基因，它們?cè)谌旧|(zhì)關(guān)閉和轉(zhuǎn)錄抑制的相對(duì)順序上有所不同牌芋，我們?cè)谖覀冄芯康乃薪M織中找到了這兩種機(jī)制的清晰例子蚓炬。 我們預(yù)計(jì) MultiVelo 將深入了解一系列生物環(huán)境中基因表達(dá)的表觀基因組調(diào)控，包括正常細(xì)胞分化躺屁、重編程和疾病肯夏。
Method（又到了數(shù)學(xué)的時(shí)間）
Previous Approaches: RNA velocity
Differential Equation Model of Gene Expression Incorporating Chromatin Accessibility
Model Likelihood
Parameter Estimation and Latent Time Inference by Expectation Maximization
Both the cell times t and the ODE parameters are unknown, so we perform expectation-maximization to simultaneously infer them. The E-step involves determining the expected value of latent time for each cell given the current best estimate of the ODE parameters. Because inverting the three-dimensional ODEs analytically is not straightforward, we perform this time estimation by finding the time whose ODE prediction is nearest each data point, selecting the time from a vector of uniformly spaced time points (see Implementation Detail section). In the M-step, we find the ODE parameters that maximize the data likelihood (equivalent to minimizing MSE) given the current time estimates for each cell. We use the Nelder-Mead simplex algorithm to minimize MSE.
Model Pre-Determination and Distinguishing Genes with Partial and Complete Dynamics
Parameter Initialization（最難的部分）
示例代碼
加載
import os
import scipy
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import multivelo as mv
import matplotlib.pyplot as plt
from IPython.core.display import display, HTML
display(HTML("<style>div.output_scroll { height: 50em; }</style>"))
display(HTML("<style>.container { width:80% !important; }</style>"))
scv.settings.verbosity = 3
scv.settings.presenter_view = True
scv.set_figure_params('scvelo')
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 200)
np.set_printoptions(suppress=True)
Reading in unspliced and spliced counts
adata_rna = scv.read("velocyto/10X_multiome_mouse_brain.loom", cache=True)
adata_rna.obs_names = [x.split(':')[1][:-1] + '-1' for x in adata_rna.obs_names]
adata_rna.var_names_make_unique()
sc.pp.filter_cells(adata_rna, min_counts=1000)
sc.pp.filter_cells(adata_rna, max_counts=20000)
# Top 1000 variable genes are used for downstream analyses.
scv.pp.filter_genes(adata_rna, min_shared_counts=10)
scv.pp.normalize_per_cell(adata_rna)
scv.pp.filter_genes_dispersion(adata_rna, n_top_genes=1000)
# Load cell annotations
cell_annot = pd.read_csv('cell_annotations.tsv', sep='\t', index_col=0)
adata_rna = adata_rna[cell_annot.index,:]
adata_rna.obs['celltype'] = cell_annot['celltype']
# We subset for lineages towards neurons and astrocytes.
adata_rna = adata_rna[adata_rna.obs['celltype'].isin(['RG, Astro, OPC',
'IPC',
'V-SVZ',
'Upper Layer',
'Deeper Layer',
'Ependymal cells',
'Subplate'])]
adata_rna
View of AnnData object with n_obs × n_vars = 3653 × 1000
obs: 'n_counts', 'initial_size_spliced', 'initial_size_unspliced', 'initial_size', 'celltype'
var: 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'gene_count_corr', 'means', 'dispersions', 'dispersions_norm', 'highly_variable'
layers: 'ambiguous', 'matrix', 'spliced', 'unspliced'
Preprocessing the ATAC counts
adata_atac = sc.read_10x_mtx('outs/filtered_feature_bc_matrix/', var_names='gene_symbols', cache=True, gex_only=False)
adata_atac = adata_atac[:,adata_atac.var['feature_types'] == "Peaks"]
# We aggregate peaks around each gene as well as those that have high correlations with promoter peak or gene expression.
# Peak annotation contains the metadata for all peaks.
# Feature linkage contains pairs of correlated genomic features.
adata_atac = mv.aggregate_peaks_10x(adata_atac,
'outs/peak_annotation.tsv',
'outs/analysis/feature_linkage/feature_linkage.bedpe',
verbose=True)
# Let's examine the total count distribution and remove outliers.
plt.hist(adata_atac.X.sum(1), bins=100, range=(0, 100000));
sc.pp.filter_cells(adata_atac, min_counts=2000)
sc.pp.filter_cells(adata_atac, max_counts=60000)
# We normalize aggregated peaks with TF-IDF.
mv.tfidf_norm(adata_atac)
Finding shared barcodes and features between RNA and ATAC
shared_cells = pd.Index(np.intersect1d(adata_rna.obs_names, adata_atac.obs_names))
shared_genes = pd.Index(np.intersect1d(adata_rna.var_names, adata_atac.var_names))
len(shared_cells), len(shared_genes)
###(3365, 936)
# We reload in the raw data and continue with a subset of cells.
adata_rna = scv.read("velocyto/10X_multiome_mouse_brain.loom", cache=True)
adata_rna.obs_names = [x.split(':')[1][:-1] + '-1' for x in adata_rna.obs_names]
adata_rna.var_names_make_unique()
adata_rna = adata_rna[shared_cells, shared_genes]
adata_atac = adata_atac[shared_cells, shared_genes]
scv.pp.normalize_per_cell(adata_rna)
scv.pp.log1p(adata_rna)
scv.pp.moments(adata_rna, n_pcs=30, n_neighbors=50)
adata_rna.obs['celltype'] = cell_annot.loc[adata_rna.obs_names, 'celltype']
adata_rna.obs['celltype'] = adata_rna.obs['celltype'].astype('category')
# Reorder the categories for color consistency with the manuscript.
all_clusters = ['Upper Layer',
'Deeper Layer',
'V-SVZ',
'RG, Astro, OPC',
'Ependymal cells',
'IPC',
'Subplate']
adata_rna.obs['celltype'] = adata_rna.obs['celltype'].cat.reorder_categories(all_clusters)
scv.tl.umap(adata_rna)
scv.pl.umap(adata_rna, color='celltype')
Smoothing gene aggregagted peaks by neighbors
# Read in Seurat WNN neighbors.
nn_idx = np.loadtxt("seurat_wnn/nn_idx.txt", delimiter=',')
nn_dist = np.loadtxt("seurat_wnn/nn_dist.txt", delimiter=',')
nn_cells = pd.Index(pd.read_csv("seurat_wnn/nn_cells.txt", header=None)[0])
# Make sure cell names match.
np.all(nn_cells == adata_atac.obs_names)
mv.knn_smooth_chrom(adata_atac, nn_idx, nn_dist)
adata_atac
AnnData object with n_obs × n_vars = 3365 × 936
obs: 'n_counts'
layers: 'Mc'
obsp: 'connectivities'
Running multi-omic dynamical model
MultiVelo incorporates chromatin accessibility information into RNA velocity and achieves better lineage predictions.
The detailed argument list can be shown with "help(mv.recover_dynamics_chrom)".
# This will take a while.
adata_result = mv.recover_dynamics_chrom(adata_rna,
adata_atac,
max_iter=5,
init_mode="invert",
verbose=False,
plot=False,
parallel=True,
save_plot=False,
rna_only=False,
fit=True,
extra_color_key='celltype'
)
# Save the result for use later on
adata_result.write("multivelo_result.h5ad")
adata_result = sc.read_h5ad("multivelo_result.h5ad")
mv.pie_summary(adata_result)
mv.switch_time_summary(adata_result)
mv.likelihood_plot(adata_result)
By default, the velocity genes used for velocity graph is determined as those whose likelihoods are above 0.05. They can be reset with "mv.set_velocity_genes" function upon examining the distributions of variables above if needed.
Computing velocity stream and latent time
mv.velocity_graph(adata_result)
mv.latent_time(adata_result)
mv.velocity_embedding_stream(adata_result, basis='umap', color='celltype')
scv.pl.scatter(adata_result, color='latent_time', color_map='gnuplot', size=80)
Let's examine some example genes
gene_list = ['Grin2b', 'Tle4', 'Robo2', 'Satb2', 'Gria2']
We can plot accessbility and expression against gene time
# Accessibility/expression by gene time, colored by the four potential states.
# The solid black curve indicates anchors.
mv.dynamic_plot(adata_result, gene_list, color_by='state', axis_on=False, frame_on=False)
We can plot velocity against gene time
# Velocity by gene time, colored by the four potential states.
# The solid black curve indicates anchors.
mv.dynamic_plot(adata_result, gene_list, color_by='state', by='velocity', axis_on=False, frame_on=False)
We can examine accessibility and expression against globally shared latent time
# Accessibility/expression by global latent time, colored by cell type assignments.
# The solid black curve indicates the mean.
mv.dynamic_plot(adata_result, gene_list, color_by='celltype', gene_time=False, axis_on=False, frame_on=False)
Phase portraits on the unspliced-spliced, chromatin-unspliced, or 3-dimensional planes can be plotted
# Unspliced-spliced phase portraits, colored by celltype.
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='us', axis_on=False, frame_on=False)
# Unspliced-spliced phase portraits, colored by log chromatin accessibility.
# title_more_info shows more information in each subplot title: model, direction, and likelihood.
mv.scatter_plot(adata_result, gene_list, color_by='c', by='us', cmap='coolwarm', title_more_info=True, axis_on=False, frame_on=False)
# Chromatin-unspliced phase portraits, colored by celltype.
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='cu', axis_on=False, frame_on=False)
# 3D phase portraits, colored by celltype.
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='cus', axis_on=False, downsample=2)
####### We can also add velocity arrows to phase portraits to show the predicted directions
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='us', axis_on=False, frame_on=False, downsample=2, velocity_arrows=True)
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='cu', axis_on=False, frame_on=False, downsample=2, velocity_arrows=True)
mv.scatter_plot(adata_result, gene_list, color_by='celltype', by='cus', downsample=3, velocity_arrows=True)
Github在MultiVelo
示例，小鼠腦犀暑、文章數(shù)據(jù)
生活很好驯击，有你更好
