3. 分析框架和工具介紹
3.1 單細(xì)胞分析框架和強(qiáng)大系統(tǒng)
如前所述胳赌，獲得計(jì)數(shù)矩陣后萝喘，探索性數(shù)據(jù)分析階段開始酥诽。由于數(shù)據(jù)的大小和復(fù)雜性鞍泉，需要專門的工具。雖然在單細(xì)胞分析的早期肮帐，人們習(xí)慣于使用自定義腳本來分析數(shù)據(jù)咖驮，但現(xiàn)在已經(jīng)存在專門用于此目的的框架。三個最受歡迎的選項(xiàng)是基于R的Bioconductor和Seurat生態(tài)系統(tǒng)以及基于Python的scverse生態(tài)系統(tǒng)训枢。它們不僅在所使用的編程語言方面有所不同托修，而且在底層數(shù)據(jù)結(jié)構(gòu)和可用的專門分析工具方面也有所不同。
Bioconductor是一個開發(fā)恒界、支持和共享免費(fèi)開源軟件的項(xiàng)目睦刃，重點(diǎn)是對包括單細(xì)胞在內(nèi)的許多不同生物測定的數(shù)據(jù)進(jìn)行嚴(yán)格且可重復(fù)的分析。大量的開發(fā)人員和優(yōu)質(zhì)的用戶體驗(yàn)以及帶有用戶友好小插圖的豐富文檔是Bioconductor的最大優(yōu)勢十酣。Seurat是一款備受推崇的R軟件包涩拙，專為分析單細(xì)胞數(shù)據(jù)而設(shè)計(jì)。它為分析的所有步驟提供工具耸采，包括多模式和空間數(shù)據(jù)兴泥。然而，對于極大的數(shù)據(jù)集（超過 50 萬個細(xì)胞）虾宇，這兩種基于R的選擇都會遇到可擴(kuò)展性問題搓彻，這促使基于Python的社區(qū)開發(fā)scverse生態(tài)系統(tǒng)。scverse 是一個致力于生命科學(xué)基礎(chǔ)工具的組織和生態(tài)系統(tǒng)文留，最初重點(diǎn)關(guān)注單細(xì)胞好唯。可擴(kuò)展性燥翅、可擴(kuò)展性以及與現(xiàn)有Python數(shù)據(jù)和機(jī)器學(xué)習(xí)工具的強(qiáng)大互操作性是scverse生態(tài)系統(tǒng)的一些優(yōu)勢骑篙。
在以下部分中，將更詳細(xì)地介紹scverse生態(tài)系統(tǒng)森书，并解釋最重要的概念靶端，重點(diǎn)關(guān)注最重要的數(shù)據(jù)結(jié)構(gòu)谎势。
3.2 使用AnnData存儲單模態(tài)數(shù)據(jù)
如前所述，在比對和基因注釋之后杨名，基因組數(shù)據(jù)通常被總結(jié)為特征矩陣脏榆。該矩陣的形式為number_observations x number_variables——其中scRNA-seq觀察結(jié)果是細(xì)胞條形碼，變量是注釋基因台谍。在分析過程中须喂，該矩陣的觀察結(jié)果和變量用計(jì)算得出的測量值（例如質(zhì)量控制指標(biāo)或潛在空間嵌入）和先驗(yàn)知識（例如源供體或替代基因標(biāo)識符）進(jìn)行注釋。在scverse生態(tài)系統(tǒng)中趁蕊，AnnData用于將數(shù)據(jù)矩陣與這些注釋相關(guān)聯(lián)坞生。為了實(shí)現(xiàn)快速且內(nèi)存高效的轉(zhuǎn)換，AnnData還支持稀疏矩陣和部分讀取掷伙。
雖然AnnData與R生態(tài)系統(tǒng)的數(shù)據(jù)結(jié)構(gòu)（例如Bioconductor的SummarizedExperiment或Seurat的object）大致相似是己，但R包使用轉(zhuǎn)置特征矩陣。
AnnData對象的核心是在X中存儲稀疏或密集矩陣（scRNA-Seq 中的計(jì)數(shù)矩陣）任柜。該矩陣的維度為obs_names x var_names卒废，其中obs（觀察值）對應(yīng)于細(xì)胞的條形碼var（變量）對應(yīng)于基因標(biāo)識符。該矩陣X被Pandas DataFramesobs和var包含宙地，它們分別保存細(xì)胞和基因的注釋摔认。此外，AnnData保存具有相應(yīng)維度的觀測值 (obsm) 或變量 (varm) 的整個計(jì)算矩陣宅粥。將細(xì)胞與細(xì)胞或基因與基因相關(guān)聯(lián)的圖形結(jié)構(gòu)通常保存在obsp和varp中。任何不適合其他slot的非結(jié)構(gòu)化數(shù)據(jù)都將保存在uns中风纠。還可以在layers中存儲更多的X值潜索。例如，在countslayer中存儲原始的拗窃、未標(biāo)準(zhǔn)化的計(jì)數(shù)數(shù)據(jù)九默，在未命名的默認(rèn)層中存儲標(biāo)準(zhǔn)化數(shù)據(jù)诈铛。
AnnData主要針對單模態(tài)（例如scRNA-Seq）數(shù)據(jù)而設(shè)計(jì)。然而，AnnData的擴(kuò)展允許高效存儲和訪問多模態(tài)數(shù)據(jù)。
3.2.1 安裝
AnnData可以在PyPI和Conda上安裝使用：
pip install anndata
conda install -c conda-forge anndata
3.2.2 初始化一個AnnData對象
首先創(chuàng)建一個具有稀疏計(jì)數(shù)信息的簡單AnnData對象罕扎，可以代表基因表達(dá)計(jì)數(shù)杆查。首先，我們導(dǎo)入所需的包。
import numpy as np
import pandas as pd
import anndata as ad
from scipy.sparse import csr_matrix
下一步，我們使用隨機(jī)泊松分布數(shù)據(jù)初始化AnnData對象搏明。將AnnData對象命名為adata。
counts = csr_matrix(np.random.poisson(1, size=(100, 2000)), dtype=np.float32)
adata = ad.AnnData(counts)
adata
輸出結(jié)果：
AnnData object with n_obs × n_vars = 100 × 2000
獲得的AnnData對象有100個觀測值和2000個變量。這相當(dāng)于100個細(xì)胞和2000個基因。我們傳遞的初始數(shù)據(jù)可以使用adata.X作為稀疏矩陣進(jìn)行訪問襟锐。
adata.X
輸出結(jié)果：
<100x2000 sparse matrix of type '<class 'numpy.float32'>'
with 126413 stored elements in Compressed Sparse Row format>
現(xiàn)在莫杈，我們分別使用.obs_names和.var_names為obs和var提供索引腥光。
adata.obs_names = [f"Cell_{i:d}" for i in range(adata.n_obs)]
adata.var_names = [f"Gene_{i:d}" for i in range(adata.n_vars)]
print(adata.obs_names[:10])
輸出結(jié)果：
Index(['Cell_0', 'Cell_1', 'Cell_2', 'Cell_3', 'Cell_4', 'Cell_5', 'Cell_6',
'Cell_7', 'Cell_8', 'Cell_9'],
dtype='object')
3.2.3 添加匹配的元數(shù)據(jù)metadata
3.2.3.1 觀察obs或變量var水平
我們的AnnData對象的核心部分已經(jīng)完整了艘儒。下一步，我們在觀察和變量級別添加元數(shù)據(jù)说铃。此類注釋存儲在AnnData對象的.obs和.var槽中窒篱，分別用于細(xì)胞和基因注釋。
ct = np.random.choice(["B", "T", "Monocyte"], size=(adata.n_obs,))
adata.obs["cell_type"] = pd.Categorical(ct) # Categoricals are preferred for efficiency
adata.obs
輸出結(jié)果
cell_type
Cell_0 B
Cell_1 T
Cell_2 B
Cell_3 T
Cell_4 Monocyte
... ...
Cell_95 T
Cell_96 T
Cell_97 B
Cell_98 T
Cell_99 Monocyte
100 rows × 1 columns
如果我們現(xiàn)在再次檢查AnnData對象墙杯，我們會注意到它已經(jīng)對obs中的cell_type信息進(jìn)行了更新配并。
adata
輸出結(jié)果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
3.2.3.2 使用元數(shù)據(jù)進(jìn)行子集化
我們還可以使用隨機(jī)生成的細(xì)胞類型對AnnData對象進(jìn)行子集化。
bdata = adata[adata.obs.cell_type == "B"]
bdata
輸出結(jié)果：
View of AnnData object with n_obs × n_vars = 26 × 2000
obs: 'cell_type'
3.2.4 觀察/變量水平矩陣
我們還可以擁有具有多個維度的元數(shù)據(jù)高镐，例如數(shù)據(jù)的UMAP嵌入溉旋。對于此類元數(shù)據(jù)，AnnData具有.obsm/.varm屬性嫉髓。我們使用keys來標(biāo)識我們插入的不同矩陣观腊。.obsm/.varm的要求是.obsm矩陣的長度必須等于觀測值的數(shù)量，因?yàn)?code>.n_obs和.varm矩陣的長度必須等于.n_vars岩喷。它們各自可以獨(dú)立地具有不同數(shù)量的維度恕沫。
讓我們從一個隨機(jī)生成的矩陣開始，我們可以將其解釋為我們想要存儲的數(shù)據(jù)的UMAP嵌入纱意，以及一些隨機(jī)的基因級元數(shù)據(jù)婶溯。
adata.obsm["X_umap"] = np.random.normal(0, 1, size=(adata.n_obs, 2))
adata.varm["gene_stuff"] = np.random.normal(0, 1, size=(adata.n_vars, 5))
adata.obsm
輸出結(jié)果：
AxisArrays with keys: X_umap
AnnData的內(nèi)容再次更新。
adata
輸出結(jié)果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
obsm: 'X_umap'
varm: 'gene_stuff'
關(guān)于.obsm/.varm的注意事項(xiàng)：
- 1.“數(shù)組樣”的元數(shù)據(jù)可以來源于Pandas DataFrame偷霉、scipy稀疏矩陣或numpy密集數(shù)組迄委。
- 2.使用scanpy時，它們的值（列）不容易繪制类少，而
.obs中的內(nèi)容很容易繪制在圖上(例如UMAP)叙身。
3.2.5 非結(jié)構(gòu)化元數(shù)據(jù)
如上所述，AnnData有.uns硫狞，它允許任何非結(jié)構(gòu)化元數(shù)據(jù)信轿。這可以是任何東西，例如包含一些對數(shù)據(jù)分析有用的一般信息的列表或字典残吩。最好僅將此slot用于無法有效存儲在其他slots中的數(shù)據(jù)财忽。
adata.uns["random"] = [1, 2, 3]
adata.uns
輸出結(jié)果：
OverloadedDict, wrapping:
OrderedDict([('random', [1, 2, 3])])
With overloaded keys:
['neighbors'].
3.2.6 Layers層
最后，我們可能有不同形式的原始核心數(shù)據(jù)泣侮，可能一種是標(biāo)準(zhǔn)化的即彪，也可能不是標(biāo)準(zhǔn)化的。這些可以存儲在AnnData的不同層中活尊。例如隶校，讓我們對原始數(shù)據(jù)進(jìn)行l(wèi)og轉(zhuǎn)換并將其存儲在一個層中漏益。
adata.layers["log_transformed"] = np.log1p(adata.X)
adata
輸出結(jié)果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
uns: 'random'
obsm: 'X_umap'
varm: 'gene_stuff'
layers: 'log_transformed'
我們的原始矩陣X沒有修改并且仍然可以訪問。我們可以通過將原始X與新層進(jìn)行比較來驗(yàn)證這一點(diǎn)深胳。
(adata.X != adata.layers["log_transformed"]).nnz == 0
輸出結(jié)果：
False
3.2.7 轉(zhuǎn)換為DataFrames
我們可以從其中一層layer獲取Pandas DataFrame绰疤。
adata.to_df(layer="log_transformed")
輸出結(jié)果：
Gene_0 Gene_1 Gene_2 Gene_3 Gene_4 Gene_5 Gene_6 Gene_7 Gene_8 Gene_9 ... Gene_1990 Gene_1991 Gene_1992 Gene_1993 Gene_1994 Gene_1995 Gene_1996 Gene_1997 Gene_1998 Gene_1999
Cell_0 0.693147 0.000000 1.386294 1.386294 0.693147 0.000000 0.693147 1.098612 0.000000 1.098612 ... 0.693147 0.693147 0.000000 0.693147 1.098612 0.000000 0.000000 1.098612 0.693147 0.000000
Cell_1 1.098612 1.098612 0.000000 0.000000 0.000000 0.693147 0.000000 0.693147 1.386294 0.000000 ... 1.098612 1.386294 1.098612 0.000000 0.693147 0.693147 1.386294 0.000000 0.000000 0.693147
Cell_2 0.693147 0.693147 1.098612 0.693147 1.386294 1.386294 0.693147 0.693147 0.000000 1.098612 ... 0.000000 0.693147 0.693147 0.000000 1.386294 0.693147 0.000000 0.000000 0.693147 0.693147
Cell_3 0.000000 0.000000 0.000000 0.693147 0.693147 0.000000 1.386294 1.098612 0.693147 0.693147 ... 0.693147 0.693147 0.000000 0.693147 0.693147 0.693147 0.693147 0.000000 1.386294 0.693147
Cell_4 0.000000 0.693147 0.000000 1.609438 1.098612 0.693147 0.000000 0.000000 0.693147 0.000000 ... 1.098612 1.098612 0.000000 1.386294 1.791759 1.098612 1.609438 0.693147 0.000000 0.693147
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Cell_95 0.693147 1.098612 0.693147 1.386294 1.098612 0.693147 0.000000 0.000000 0.000000 0.693147 ... 0.000000 0.693147 0.000000 0.000000 0.000000 0.693147 0.000000 0.693147 1.386294 0.693147
Cell_96 0.693147 0.000000 0.000000 1.098612 0.000000 0.693147 0.693147 1.098612 0.000000 0.000000 ... 1.386294 0.693147 1.098612 0.693147 0.000000 0.693147 1.098612 1.386294 0.693147 0.000000
Cell_97 0.693147 0.000000 0.000000 0.000000 0.000000 0.000000 1.098612 0.000000 0.000000 0.693147 ... 0.000000 0.693147 1.386294 0.000000 0.000000 0.000000 0.693147 0.693147 1.609438 0.693147
Cell_98 0.693147 0.693147 0.000000 1.098612 0.000000 0.693147 0.693147 1.098612 0.000000 0.000000 ... 0.000000 0.000000 0.693147 0.693147 0.000000 0.693147 1.386294 0.693147 1.098612 1.386294
Cell_99 1.098612 1.791759 0.000000 0.000000 0.000000 0.000000 1.098612 1.609438 0.693147 0.693147 ... 0.693147 0.000000 0.000000 0.693147 0.693147 0.693147 0.693147 0.000000 0.000000 1.386294
100 rows × 2000 columns
3.2.8 AnnData對象的讀取和寫入
AnnData對象可以在磁盤上以分層數(shù)組存儲（如 HDF5 或 Zarr）。AnnData帶有自己的基于HDF5的持久文件格式：h5ad〕硗溃現(xiàn)在峦睡，我們將以h5ad格式保存AnnData對象。
adata.write("my_results.h5ad", compression="gzip")
!h5ls 'my_results.h5ad'
輸出結(jié)果：
X Group
layers Group
obs Group
obsm Group
uns Group
var Group
varm Group
以及將其讀回权埠。
adata_new = ad.read_h5ad("my_results.h5ad")
adata_new
輸出結(jié)果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
uns: 'random'
obsm: 'X_umap'
varm: 'gene_stuff'
layers: 'log_transformed'
3.3 scanpy分析單模態(tài)數(shù)據(jù)
現(xiàn)在我們了解了單模態(tài)單細(xì)胞分析的基本數(shù)據(jù)結(jié)構(gòu)榨了，但是我們?nèi)绾螌?shí)際分析存儲的數(shù)據(jù)呢？scverse生態(tài)系統(tǒng)中存在多種用于分析特定組學(xué)數(shù)據(jù)的工具攘蔽。例如龙屉，scanpy提供用于一般RNA-Seq重點(diǎn)分析的工具，squidpy重點(diǎn)關(guān)注空間轉(zhuǎn)錄組學(xué)满俗，而scirpy提供用于分析T細(xì)胞受體 (TCR) 和B細(xì)胞受體 (BCR) 數(shù)據(jù)的工具转捕。盡管存在許多針對各種數(shù)據(jù)模式的橫向擴(kuò)展，但它們通常使用scanpy的一些預(yù)處理和可視化功能唆垃。
更具體地說五芝，scanpy是一個建立在AnnData之上的Python包，用于促進(jìn)單細(xì)胞基因表達(dá)數(shù)據(jù)的分析辕万∈嗖剑可通過scanpy訪問預(yù)處理、嵌入渐尿、可視化醉途、聚類、差異基因表達(dá)測試砖茸、偽時間和軌跡推斷以及基因調(diào)控網(wǎng)絡(luò)模擬的多種方法隘擎。基于Python數(shù)據(jù)科學(xué)和機(jī)器學(xué)習(xí)庫凉夯，使scanpy能夠擴(kuò)展到數(shù)百萬個細(xì)胞货葬。
3.3.1 安裝
Scanpy可以在PyPI和Conda安裝使用：
pip install scanpy
conda install -c conda-forge scanpy
3.3.2 Scanpy的API設(shè)計(jì)
scanpy的框架設(shè)計(jì)是將屬于同一步驟的函數(shù)分組到相應(yīng)的模塊中。例如劲够，所有預(yù)處理功能都可以在scanpy.preprocessing模塊中使用宝惰，所有非預(yù)處理的數(shù)據(jù)矩陣轉(zhuǎn)換都可以在scanpy.tools中使用，所有可視化都可以在scanpy.plot中使用再沧。這些模塊通常在導(dǎo)入scanpy之后訪問，例如import scanpy as sc以及相應(yīng)的縮寫sc.pp（用于預(yù)處理）尊残、sc.tl（用于工具）和sc.pl（用于繪圖）炒瘸。所有讀或?qū)憯?shù)據(jù)的模塊都是直接訪問的淤堵。此外，用于各種數(shù)據(jù)集的模塊可用作sc.datasets顷扩。所有函數(shù)及其相應(yīng)的參數(shù)和潛在的示例圖都記錄在scanpy API文檔中拐邪。
3.3.3 Scanpy示范
在下面的內(nèi)容中，我們將很快演示使用scanpy進(jìn)行分析的工作流程隘截。
選擇的數(shù)據(jù)集是健康捐獻(xiàn)者的2700個外周血單核細(xì)胞的數(shù)據(jù)集扎阶，這些細(xì)胞在Illumina NextSeq 500上進(jìn)行了測序。
作為第一步婶芭，我們導(dǎo)入scanpy东臀。 我們設(shè)置Matplotlib繪圖作為所有scanpy繪圖的默認(rèn)值，并且最后打印scanpy的header犀农。包含當(dāng)前環(huán)境中所有相關(guān)Python包的版本惰赋，包括Scanpy和AnnData。
import scanpy as sc
sc.settings.set_figure_params(dpi=80, facecolor="white")
sc.logging.print_header()
輸出結(jié)果：
scanpy==1.9.1 anndata==0.7.8 umap==0.5.2 numpy==1.20.3 scipy==1.7.2 pandas==1.4.3 scikit-learn==1.0.1 statsmodels==0.13.1 python-igraph==0.9.10 pynndescent==0.5.5
接下來呵哨，我們使用scanpy作為AnnData對象獲取3k PBMC數(shù)據(jù)集赁濒。數(shù)據(jù)集會自動下載并保存，也可以使用scanpy設(shè)置指定的路徑孟害。默認(rèn)情況下它是在文件夾data拒炎。
adata = sc.datasets.pbmc3k()
adata
輸出結(jié)果：
100%|██████████| 5.58M/5.58M [00:01<00:00, 3.25MB/s]
AnnData object with n_obs × n_vars = 2700 × 32738
var: 'gene_ids'
返回的AnnData對象有2700個細(xì)胞和32738個基因，var中還包含基因ID挨务。
adata.var
輸出結(jié)果：
gene_ids
index
MIR1302-10 ENSG00000243485
FAM138A ENSG00000237613
OR4F5 ENSG00000186092
RP11-34P13.7 ENSG00000238009
RP11-34P13.8 ENSG00000239945
... ...
AC145205.1 ENSG00000215635
BAGE5 ENSG00000268590
CU459201.1 ENSG00000251180
AC002321.2 ENSG00000215616
AC002321.1 ENSG00000215611
32738 rows × 1 columns
如上所述击你，scanpy的所有分析函數(shù)都可以通過sc.[pp, tl, pl]訪問。作為第一步耘子，我們使用scanpy來顯示在所有細(xì)胞中每個單細(xì)胞中產(chǎn)生最高計(jì)數(shù)分?jǐn)?shù)的基因果漾。我們只需調(diào)用sc.pl.highest_expr_genes函數(shù)，傳遞AnnData對象谷誓，該對象在幾乎所有情況下都是任何scanpy函數(shù)的第一個參數(shù)绒障，并指定我們希望顯示前20個表達(dá)基因。
sc.pl.highest_expr_genes(adata, n_top=20)
顯然捍歪，MALAT1是表達(dá)最多的基因户辱，經(jīng)常在Poly-A捕獲的scRNA-Seq數(shù)據(jù)中檢測到，與實(shí)驗(yàn)方案無關(guān)糙臼。該基因已被證明與細(xì)胞健康呈負(fù)相關(guān)庐镐。特別是死亡/垂死的細(xì)胞具有更高的MALAT1表達(dá)。
現(xiàn)在变逃，我們使用scanpy的預(yù)處理模塊過濾檢測到的基因少于200個的細(xì)胞以及在少于3個細(xì)胞中發(fā)現(xiàn)的基因必逆，以獲得粗略的質(zhì)量閾值。
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
單細(xì)胞RNA-Seq分析的一個常見步驟是通過主成分分析 (PCA) 等進(jìn)行降維，以揭示主要的差異變量名眉。這也會對數(shù)據(jù)進(jìn)行去噪粟矿。scanpy的preprocessing或者tools函數(shù)提供PCA方法。
sc.tl.pca(adata, svd_solver="arpack")
相應(yīng)的繪圖函數(shù)允許我們將基因傳遞給顏色參數(shù)损拢。相應(yīng)的值會自動從AnnData對象中提取陌粹。
sc.pl.pca(adata, color="CST3")
任何高級嵌入和下游計(jì)算的基本步驟是使用數(shù)據(jù)矩陣的PCA表示來計(jì)算鄰域圖。它會自動用于需要它的其他工具福压，例如UMAP的計(jì)算掏秩。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
我們現(xiàn)在使用計(jì)算鄰域圖將單元格嵌入到UMAP中，UMAP是在scanpy中實(shí)現(xiàn)的許多高級降維算法之一荆姆。
sc.tl.umap(adata)
sc.pl.umap(adata, color=["CST3", "NKG7", "PPBP"])
scanpy的文檔提供了教程蒙幻，需要的可以認(rèn)真閱讀。
3.4 后記
其實(shí)除了處理單模態(tài)數(shù)據(jù)外胞枕，AnnData還可以存儲和操作多模態(tài)數(shù)據(jù)杆煞，如CITE-seq等同時測量RNA和表面蛋白數(shù)據(jù)，但這一部分是比較進(jìn)階的內(nèi)容腐泻，目前先只介紹簡單的普通單細(xì)胞RNA-seq分析的內(nèi)容决乎，多模態(tài)的數(shù)據(jù)，之后再詳細(xì)分享給大家派桩。
