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