如何進行細胞軌跡分析?
目前進行細胞軌跡分析的方法和軟件非常之多，大致可以概括為兩種方法，一種是以monocle軟件為代表的擬時序分析（pseudotime analysis），另一種則是以velocyto /scVelo為主的RNA速度分析（RNA velocity）。
1、擬時序分析
2019年Nature biotechnology的一篇文章A comparison of single-cell trajectory inference methods [1] 就對45種擬時序分析的方法進行了評估，這些軟件在性能、可擴展性、穩(wěn)定性和易用性上存在差異，當然也各有其優(yōu)缺點。軟件推斷的細胞分化軌跡主要有7種類型，包括環(huán)狀（cycle）、線性（linear）、分叉（bifurcation）、多分叉（Multifurcation）、樹型（Tree）、多種軌跡并存的連接圖（Connected graph），及多種不相連軌跡并存的斷開圖（disconnected graph）。（圖1）
![圖1 ?軌跡分析算法評估流程及7種細胞分化軌跡類型[1]](http://specchiomagico.net/wp-content/uploads/2022/02/640.png)
圖1 ?軌跡分析算法評估流程及7種細胞分化軌跡類型[1]
2、RNA速度分析
目前絕大多數(shù)軟件利用的是細胞的基因表達信息，2018年發(fā)表在Nature上的一篇文獻RNA velocity of single cells[3]，提出了一個新的概念，即RNA速度（RNA velocity ），利用的是reads到表達的信息，通過分析新轉(zhuǎn)錄的、未發(fā)生剪接的前體mRNA，與成熟的、發(fā)生剪接的mRNA（reads是否比對到內(nèi)含子）豐度之間的比值，來表征細胞定向的動態(tài)信息。scVelo[4]軟件可以推斷包括轉(zhuǎn)錄、剪接和降解在內(nèi)的”基因特異性速度(gene-specific rates)”，并還原細胞過程中的”潛在時間(latent-time)”，而潛在時間近似為細胞分化的真實時間，同時可以識別調(diào)控機制，并系統(tǒng)性的推定驅(qū)動基因。scVelo分析流程及代碼貼在文章末尾，感興趣的老師不妨實操看看。
圖2 ?scVelo軟件分析原理示意圖[4]
相較于擬時序分析來說，RNA velocity分析的是細胞在某一時刻的潛在分化方向，可以直接從具有時間向量的結(jié)果讀取到分化方向，無需基于生物學背景指定分化起點和終點；同時RNA velocity的數(shù)據(jù)可以與Seurat降維數(shù)據(jù)兼容，也可以嵌入monocle構(gòu)建細胞軌跡，這對于聯(lián)合兩種分析策略十分便利，利用RNA velocity分析揭示所有未知的細胞分化關系，然后篩選關注的細胞分化過程進行擬時序分析，挖掘分化過程中的調(diào)控機制。
如何應用細胞軌跡分析?
1、細胞分化過程重現(xiàn)
基于單細胞測序分析數(shù)據(jù)，對其他技術(shù)手段獲取的已知細胞分化關系進行驗證。
2020年Nature Communication發(fā)表的Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy[5]，利用RNA速度表征了膠質(zhì)母細胞瘤的分化過程，證實了膠質(zhì)母細胞瘤層次結(jié)構(gòu)頂點是腫瘤祖細胞，星形細胞、間充質(zhì)、少突膠質(zhì)細胞和神經(jīng)元癌細胞比腫瘤祖細胞的分化程度更高，并且在每個患者中均可見特定細胞分化方向（圖3）。
圖 ?3 ?膠質(zhì)母細胞瘤RNA速度分析[5]
圖4 ?人睪丸生殖細胞亞群及分化軌跡鑒定[6]
2、潛在的細胞分化路徑挖掘
除了鑒定已知的細胞分化過程，通過細胞軌跡分析也能推斷未被報道的、可能的未知細胞分化路徑，如神經(jīng)干細胞發(fā)育、腫瘤/疾病微環(huán)境中免疫細胞分化過程，甚至挖掘到一些稀少的、常規(guī)技術(shù)手段很難獲得的中間狀態(tài)細胞。
2021年Nature Communication發(fā)表的Single-cell sequencing of immune cells from anticitrullinated peptide antibody positive and negative rheumatoid arthritis[7]，研究首先通過比較兩種亞型類風濕性關節(jié)炎（RA）患者與健康對照外周血免疫細胞圖譜差異，關注到了B細胞類型豐度變化，發(fā)現(xiàn)ACPA-RA組IGHG4+ plasma B細胞明顯升高，HLA-DRB5+ plasma B細胞明顯降低（圖5，左圖）；隨后對B細胞進行擬時序分析探究B細胞亞型分化關系，發(fā)現(xiàn)在健康對照和RA患者之間B細胞亞型存在不同的分化路徑：在健康對照中na?ve B細胞有兩條分化路徑，1條是分化為HLA-DRB5+plasma B細胞再分化為plasma B細胞，另1條分化為Mem-unsw B細胞到Mem-sw B細胞；而RA患者中，na?ve B細胞直接分化為IGHG4+ plasma B細胞（圖5，右圖）。
圖5 ?類風濕關節(jié)炎患者外周血B細胞分析[7]
圖6 ??膠質(zhì)母細胞瘤干細胞樣細胞分化軌跡[8]
3、細胞分化過程調(diào)控機制解析
通過細胞軌跡分析，分析處于重要分化節(jié)點的關鍵細胞亞群，尋找隨著分化時間呈現(xiàn)特定表達模式的關鍵基因，解析驅(qū)動細胞亞群分化的調(diào)控機制。
2021年 Nature發(fā)表的Single-cell transcriptomic characterization of a gastrulating human embryo[9]，為了探索原腸胚形成過程中外胚層的變化，利用RNA速度分析了外胚層（Epiblast）、外胚層（羊膜/胚胎）（Ectoderm (amniotic/embryonic)）、新生中胚層（Nascent mesoderm）、原條（Primitive Streak），由Epiblast起，發(fā)生兩個分支，一側(cè)由Primitive Streak發(fā)育到Nascent mesoderm，一側(cè)發(fā)育到ectoderm（圖7，左圖）；隨后利用擬時分析推斷Epiblast分化過程中的基因表達變化，檢測到ectoderm marker genes（DLX5、TFAP2A 和GATA3）表達上調(diào)，而早期神經(jīng)誘導marker genes（SOX1、SOX3 和PAX6）和分化的神經(jīng)元marker genes（TUBB3、OLIG2和NEUROD1）表達極低甚至檢測不到，這表明CS7中神經(jīng)分化尚未開始（圖7，右圖）。
圖7 ?RNA速度分析解析原腸胚分化軌跡[9]
如何進行RNA速度分析?
1、軟件安裝：
pip?install?-U?scvelo?#需要3.6以上版本python
pip?install?velocyto
2、adata 數(shù)據(jù)結(jié)構(gòu)：
– adata.X：數(shù)據(jù)矩陣
– adata.obs：細胞的注釋
–?adata.var：基因的注釋
– adata.nus：非結(jié)構(gòu)化注釋，例如圖
– adata.layers：附加層數(shù)據(jù)，如spliced/unspliced矩陣
3、加載需要的庫及一些配置
import?seaborn?as?sns
import?scvelo?as?scv
import?pandas?as?pd
import?scanpy?as?sc
import?numpy?as?np
import?loompy
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
4、計算獲取unspliced/spliced矩陣
##?運行run10x
/software/python3.9.6/bin/velocyto?run10x??????????\
#?相應的repeatmasker的重復序列注釋文件（可從USCU中獲取）
#?http://genome.ucsc.edu/cgi-bin/hgTables
-m?RMSK
cellranger_results??#?cellranger的輸出結(jié)果路徑
sample_name???????????#?樣本名
seurat_GTF??????????????#?用于cellranger定量的基因組注釋文件
##?輸出結(jié)果，在cellranger的結(jié)果中會創(chuàng)建單個樣本的velocyto文件夾，其中的loom文件即為最終結(jié)果
data/SC/Normal1-sc/velocyto:
Normal1-sc.loom
data/SC/Normal2-sc/velocyto:
Normal2-sc.loom
data/SC/Normal3-sc/velocyto:
Normal3-sc.loom
data/SC/Normal4-sc/velocyto:
Normal4-sc.loom
5 、loom文件的讀取及合并
#?讀取單個樣本loom文件
Normal_loom?=?sc.read(‘./Normal1.loom’,?sparse=True,?cache=True)
#?合并同組樣本的loom文件
import?loompy
Normal_loom_list?=?[Normal1-sc.loom,
Normal2-sc.loom,
Normal3-sc.loom,
Normal4-sc.loom]
loompy.combine(Normal_loom_list,?output_file?=?‘Normal_combined.loom’)
#?讀取
Normal_loom?=?sc.read(‘./Normal_combined.loom’,?sparse=True,?cache=True)
6 ?導入Seurat對象的降維聚類結(jié)果
使用scVelo進行分析時需要對數(shù)據(jù)進行降維聚類，若上游分析的降維聚類注釋結(jié)果基于Seurat，可以直接將Seurat的降維聚類結(jié)果添加到adata對象中，下面方法僅供參考。此外，如果上游分析基于scanpy，則可以直接調(diào)用scanpy的降維聚類相關的方法。注：此處并沒有進行注釋，在實際操作中將 seurat_cluster替換為注釋后的信息即可，導入seurat之后的所有調(diào)用的方法中的groupby參數(shù)設置為 seurat_cluster。
library(Seurat)
library(tidyverse)
seurat_obj@meta.data?%>%?rownames_to_column(var?=?‘barcodes’)?%>%
write.table(‘seurat_meta_df.txt’,?sep?=?‘\t’,?quote?=?F)
seurat_obj@reductions$umap@cell.embeddings?%>%?rownames_to_column(var?=?‘barcodes’)?%>%
write.table(‘seurat_umap.txt’,?sep?=?‘\t’,?quote?=?F)
def?using_seurat_meta_umap(loom_file,?seurat_meta_df,?seurat_umap):
#?讀入seurat的meta.data以及umap降維結(jié)果
seurat_meta_data_raw?=?pd.read_csv(seurat_meta_df)
seurat_umap_raw?=?pd.read_csv(seurat_umap)
#?loom文件原始的index順序
anndata_index?=?loom_file.obs.index.values
#?合并
meta_data?=?loom_file.obs.reset_index().?\
merge(seurat_meta_data_raw,
left_on=”CellID”,
right_on=”barcodes”,?how=’left’).?\
set_index(‘CellID’).reindex(anndata_index)
meta_data[‘seurat_clusters’]?=?meta_data[‘seurat_clusters’].?\
apply(lambda?x:?“%d”?%?int(x)?if?~np.isnan(x)?else?x).astype(‘category’)
loom_file.obs?=?meta_data
loom_file?=?loom_file[loom_file.obs.dropna().index,?:]
loom_file.obsm[‘umap’]?=?seurat_umap_raw.set_index(‘barcodes’).?\
reindex(loom_file.obs.index).to_numpy().astype(‘float64’)
return?loom_file
Normal_loom?=?using_seurat_meta_umap(Normal_loom,
‘seurat_meta_df.txt’,
‘seurat_umap.txt’)
#?查看是否替換成功，與seurat結(jié)果進行比較
sc.pl.umap(Normal_loom,?color=’seurat_clusters’)
分析結(jié)果展示
拼接/未拼接計數(shù)的比例，通常有 10%-25% 的未拼接分子包含內(nèi)含子序列。
scv.pl.proportions(Normal_loom,?groupby=’seurat_clusters’)
RNA-速度：速度是基因表達空間中的向量，代表單個細胞運動的方向和速度。速度是通過模擬剪接動力學的轉(zhuǎn)錄動力學獲得的，對于每個基因，擬合了成熟（未剪接）和成熟（剪接）mRNA?計數(shù)的穩(wěn)態(tài)比率，這構(gòu)成了恒定的轉(zhuǎn)錄狀態(tài)。正速度表明基因被上調(diào)，這種情況發(fā)生在該基因未剪接?mRNA?豐度高于預期的穩(wěn)定狀態(tài)的細胞中。相反，負速度表明基因被下調(diào)。
#?預處理
scv.pp.filter_genes(Normal_loom,?min_shared_counts=20)
scv.pp.normalize_per_cell(Normal_loom)
scv.pp.filter_genes_dispersion(Normal_loom,?n_top_genes=2000)
scv.pp.log1p(Normal_loom)
scv.pp.filter_and_normalize(Normal_loom,
min_shared_counts=20,
n_top_genes=2000)
#?運行scVelo動態(tài)模型，如果算力不夠可以考慮靜態(tài)或者隨機方法
scv.tl.recover_dynamics(Normal_loom,?n_jobs=30)
scv.tl.velocity(Normal_loom,?mode=’dynamical’)
scv.tl.velocity_graph(Normal_loom,?n_jobs=n_jobs)
#?保存計算結(jié)果
Normal_loom.write(‘Normal_dy_model.h5ad’,?compression=’gzip’)
scv.pl.velocity_embedding_stream(Normal_loom,
basis=’umap’,
color=’seurat_clusters’)
scv.pl.velocity_embedding(Normal_loom,
arrow_length=3,
arrow_size=2,
show=False,
color=’seurat_clusters’)
潛伏時間：動力學模型描述了細胞過程的潛伏時間（這個潛伏時間代表細胞的內(nèi)部時鐘，并近似于細胞在分化時經(jīng)歷的真實時間，僅基于其轉(zhuǎn)錄動力學。）
scv.tl.latent_time(Normal_loom)
scv.pl.scatter(Normal_loom,?color=’latent_time’,
color_map=’gnuplot’,?size=80,?show=False)
驅(qū)動基因：顯示出明顯的動態(tài)行為，并通過動態(tài)模型中的高可能性特征系統(tǒng)地檢測到。
top_genes?=?Normal_loom.var[‘fit_likelihood’].?\
sort_values(ascending=False).index[:300]
scv.pl.heatmap(Normal_loom,
var_names=top_genes,?sortby=’latent_time’,
col_color=’seurat_clusters’,?n_convolve=100)
scv.pl.scatter(Normal_loom,?basis=top_genes[:15],
ncols=5,?frameon=False,
color=’seurat_clusters’)
基于RNA-速度的PAGA分析：PAGA提供了軌跡推斷構(gòu)圖的方法，其中加權(quán)邊對應于兩個集群之間的連接。在這里，PAGA通過速度推斷的方向性進行了擴展。
scv.tl.paga(Normal_loom,?groups=’seurat_clusters’)
scv.pl.paga(Normal_loom,?basis=’umap’,
size=50,?alpha=.1,
min_edge_width=2,?node_size_scale=1.5)
參考文獻：
京公網(wǎng)安備 11011302003368號