胞數(shù)據(jù)分析實(shí)戰(zhàn)指南)
1. 背景與核心概念1.1 Jupyter Notebook 是什么Jupyter Notebook 是一個(gè)基于 Web 的交互式開發(fā)環(huán)境允許我們在瀏覽器中編寫代碼、運(yùn)行代碼、查看輸出結(jié)果同時(shí)還能在代碼之間插入 Markdown 說明文字、公式、圖片和表格。最早它作為 IPython Notebook 項(xiàng)目誕生后來逐步發(fā)展為支持 Python、R、Julia 等多種語言的項(xiàng)目也因此得名 JupyterJulia Python R 的合稱。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域Jupyter Notebook 幾乎是標(biāo)配工具。原因很簡單單細(xì)胞數(shù)據(jù)的分析過程是一條長鏈路包括數(shù)據(jù)讀取、質(zhì)量控制、歸一化、降維、聚類、注釋等多個(gè)環(huán)節(jié)每一步之間高度依賴人工判斷。Notebook 這種“一段代碼 一段結(jié)果 一段說明”的組織方式非常適合探索性數(shù)據(jù)分析和流程復(fù)現(xiàn)。你可以在同一個(gè)文件中既完成代碼調(diào)試又留下對參數(shù)選擇的理解記錄后續(xù)想要回溯分析思路也會(huì)輕松很多。1.2 Jupyter Notebook 與 JupyterLab 的區(qū)別很多初學(xué)者在安裝時(shí)會(huì)碰到一個(gè)搜索熱點(diǎn)Jupyter Notebook 和 JupyterLab 到底有什么區(qū)別簡單來說JupyterLab 是 Jupyter Notebook 的下一代交互式開發(fā)界面相當(dāng)于把 Notebook 的能力放到了一個(gè)更像 IDE 的工作臺(tái)中。JupyterLab 支持多標(biāo)簽頁、拖拽布局、文件管理器、終端、代碼控制臺(tái)等功能在同一個(gè)窗口中可以同時(shí)打開多個(gè) Notebook、查看數(shù)據(jù)文件、運(yùn)行終端命令。而經(jīng)典 Jupyter Notebook 則更聚焦于單個(gè) Notebook 文件的編輯和運(yùn)行界面更加簡潔適合只專注分析任務(wù)的時(shí)候使用。兩者使用同一個(gè) Notebook 文件格式.ipynb運(yùn)行時(shí)也共享同一套內(nèi)核機(jī)制。也就是說你在經(jīng)典 Notebook 里寫的代碼在 JupyterLab 里也能正常打開運(yùn)行?，F(xiàn)在 Anaconda 默認(rèn)安裝的是 JupyterLab 和 Notebook 兩套組件啟動(dòng)命令分別對應(yīng)jupyter notebook jupyter lab如果只是在學(xué)習(xí)單細(xì)胞分析的基礎(chǔ)流程使用哪個(gè)都可以如果希望邊寫代碼邊觀察文件目錄、邊調(diào)試邊查看圖表JupyterLab 的體驗(yàn)會(huì)更好。本文的代碼在兩種界面下均可直接運(yùn)行。1.3 Scanpy 是什么Scanpy 是一個(gè)基于 Python 的高性能單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析工具包全稱是 Single-Cell Analysis in Python。它的核心數(shù)據(jù)結(jié)構(gòu)基于 AnnData能夠高效處理大規(guī)模稀疏矩陣支持從數(shù)據(jù)讀取、質(zhì)量控制、標(biāo)準(zhǔn)化、特征選擇、降維、聚類到差異分析、可視化以及偽時(shí)間分析等完整流程。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域R 語言生態(tài)中的 Seurat 一直是主流工具而 Scanpy 的出現(xiàn)讓 Python 用戶也能完成同等級(jí)的分析工作。相比 SeuratScanpy 有以下幾個(gè)特點(diǎn)構(gòu)建在 Anndata、NumPy、SciPy 和 Pandas 之上與 Python 數(shù)據(jù)科學(xué)生態(tài)天然融合內(nèi)存占用控制較好適合處理幾萬甚至更多細(xì)胞的單細(xì)胞數(shù)據(jù)集分析流程高度模塊化每個(gè)步驟對應(yīng)一個(gè)明確的 API支持豐富的可視化圖形輸出可直接嵌入 Notebook 中展示。對于已經(jīng)熟悉 Python 的數(shù)據(jù)分析人員來說Scanpy 是入門單細(xì)胞數(shù)據(jù)分析非常合適的切入點(diǎn)。我們在本文中會(huì)使用一個(gè)經(jīng)典示例數(shù)據(jù)集 pbmc3k完整跑一遍單細(xì)胞轉(zhuǎn)錄組的標(biāo)準(zhǔn)分析流程。1.4 為什么把 Jupyter Notebook 與 Scanpy 放在一起學(xué)習(xí)單細(xì)胞數(shù)據(jù)分析并不是“把數(shù)據(jù)丟進(jìn)一個(gè)函數(shù)就能出結(jié)果”的事。它需要分析者反復(fù)調(diào)整閾值參數(shù)、查看中間圖表、確認(rèn)聚類效果這一過程天然適合 Jupyter Notebook 這種交互式環(huán)境。用 Jupyter Notebook 跑 Scanpy每個(gè)分析步驟都能立即看到輸出和圖表比如質(zhì)量控制后的細(xì)胞數(shù)變化、UMAP 降維后的聚類分布、marker 基因的表達(dá)情況。遇到某個(gè)步驟結(jié)果不理想可以直接修改參數(shù)重新運(yùn)行不需要從頭啟動(dòng)整個(gè)腳本。這種探索式的工作流正是 Scanpy 教程和單細(xì)胞分析論文代碼大量采用 Notebook 形式的原因。所以本文的安排是先掌握 Jupyter Notebook 的基礎(chǔ)操作再用它跑通 Scanpy 的單細(xì)胞分析流程最終形成一個(gè)可以直接照抄、修改、復(fù)用到自己數(shù)據(jù)上的分析模板。2. 環(huán)境準(zhǔn)備與版本說明2.1 安裝 AnacondaScanpy 依賴的第三方庫較多手動(dòng)用 pip 逐個(gè)安裝容易遇到依賴沖突所以我們首先推薦通過 Anaconda 來管理 Python 環(huán)境。Anaconda 自帶 Python 解釋器、Jupyter Notebook、常用科學(xué)計(jì)算庫以及包管理器 conda安裝完成后即可快速使用。到 Anaconda 官網(wǎng)下載對應(yīng)操作系統(tǒng)的安裝包安裝過程中保持默認(rèn)選項(xiàng)即可。安裝完成后打開命令行窗口輸入conda --version如果能看到 conda 版本號(hào)說明安裝成功。為了保證環(huán)境干凈建議單獨(dú)創(chuàng)建一個(gè)用于單細(xì)胞分析的虛擬環(huán)境conda create -n scanpy-env python3.9 -y conda activate scanpy-env關(guān)于 Python 版本需要根據(jù)你的項(xiàng)目實(shí)際情況調(diào)整這里選擇 3.9 是兼容性較好的通用版本。如果你的電腦上已經(jīng)存在其他 Python 項(xiàng)目使用獨(dú)立虛擬環(huán)境可以避免庫與庫之間的版本沖突這是單細(xì)胞分析項(xiàng)目里非常重要的一步。2.2 安裝 Scanpy激活虛擬環(huán)境后使用 conda 安裝 Scanpyconda install -c conda-forge scanpy python-igraph leidenalg這里同時(shí)安裝了python-igraph和leidenalg它們用于 Leiden 聚類算法。Leiden 是目前單細(xì)胞聚類中較推薦的方法相比早期的 Louvain 算法社區(qū)劃分更穩(wěn)定對分辨率的調(diào)節(jié)也更友好。如果 conda 安裝速度慢也可以使用 pippip install scanpy pip install leidenalg不過 pip 安裝時(shí)需要注意 NumPy、Pandas、SciPy 等底層庫的版本兼容性如果出現(xiàn)版本報(bào)錯(cuò)建議仍然回到 conda 安裝方式。安裝完成后在命令行輸入python -c import scanpy as sc; print(sc.__version__)能正常打印版本號(hào)就說明環(huán)境已經(jīng)準(zhǔn)備好了。2.3 啟動(dòng) Jupyter Notebook在scanpy-env環(huán)境中啟動(dòng) Notebookjupyter notebook啟動(dòng)后終端會(huì)輸出類似下面的信息[I 2025-01-01 10:00:00.123 NotebookApp] Serving notebooks from local directory: /Users/xxx/scanpy-tutorial [I 2025-01-01 10:00:00.125 NotebookApp] Jupyter Notebook 6.5.4 is running at: [I 2025-01-01 10:00:00.125 NotebookApp] http://localhost:8888/?tokenxxxxxxxx瀏覽器會(huì)自動(dòng)打開 Notebook 的主頁面。如果瀏覽器沒有自動(dòng)打開可以把終端里顯示的http://localhost:8888/?tokenxxx地址手動(dòng)復(fù)制到瀏覽器地址欄。在 Notebook 主頁面的右上角點(diǎn)擊 New - Python 3即可新建一個(gè) Notebook 文件。建議把工作目錄切換到專門存放單細(xì)胞分析代碼的文件夾例如scanpy-tutorial這樣后續(xù)導(dǎo)入數(shù)據(jù)文件時(shí)路徑管理更方便。2.4 示例項(xiàng)目結(jié)構(gòu)本文的實(shí)戰(zhàn)代碼按以下目錄組織scanpy-tutorial/ ├── pbmc3k.h5ad # 示例數(shù)據(jù)運(yùn)行過程中自動(dòng)下載 └── scanpy_analysis.ipynb # 分析 Notebook 文件實(shí)際操作中你的文件路徑不一定要完全一致但建議保持?jǐn)?shù)據(jù)和 Notebook 在同一級(jí)目錄減少路徑出錯(cuò)的可能。3. Jupyter Notebook 核心操作3.1 Cell 與兩種模式Notebook 文件由一個(gè)個(gè)單元格Cell組成每個(gè) Cell 可以存放代碼也可以存放 Markdown 文檔。運(yùn)行 Cell 時(shí)Jupyter 會(huì)把代碼發(fā)送給內(nèi)核執(zhí)行并把結(jié)果展示在 Cell 下方。Notebook 有兩種模式命令模式按 Esc 進(jìn)入此時(shí)鍵盤快捷鍵作用于整個(gè) Notebook可以刪除、復(fù)制、移動(dòng) Cell編輯模式按 Enter 進(jìn)入此時(shí)可以編輯當(dāng)前 Cell 中的代碼或文字。新建的 Cell 默認(rèn)是代碼類型如果要寫說明文字需要把 Cell 切換為 Markdown 類型。切換方式是在命令模式下按M切回代碼類型則按Y。3.2 常用快捷鍵掌握以下快捷鍵可以明顯提升操作效率快捷鍵作用Shift Enter運(yùn)行當(dāng)前 Cell 并跳轉(zhuǎn)到下一個(gè) CellCtrl Enter運(yùn)行當(dāng)前 Cell不跳轉(zhuǎn)Esc進(jìn)入命令模式A在當(dāng)前 Cell 上方插入新 CellB在當(dāng)前 Cell 下方插入新 CellD D連續(xù)按兩次 D刪除當(dāng)前 CellM將當(dāng)前 Cell 切換為 MarkdownY將當(dāng)前 Cell 切換為代碼Shift Tab查看函數(shù)或?qū)ο蟮奈臋n摘要其中 Shift Tab 在調(diào)用 Scanpy 函數(shù)時(shí)尤其實(shí)用。比如你忘了sc.pp.filter_cells的參數(shù)含義把光標(biāo)放到函數(shù)名上按 Shift Tab就能快速看到函數(shù)簽名和參數(shù)說明不用頻繁去查文檔。3.3 Markdown 與魔法命令Markdown Cell 中可以使用標(biāo)準(zhǔn)的 Markdown 語法包括#標(biāo)題、**加粗**、-列表、行內(nèi)代碼等。此外 Jupyter 還支持 LaTeX 數(shù)學(xué)公式用$$包裹即可。在單細(xì)胞分析筆記中我習(xí)慣在每個(gè)步驟前用 Markdown Cell 寫清楚“這一步在做什么、為什么做”分析結(jié)束之后再回頭看整個(gè)文件就像一份可執(zhí)行的分析報(bào)告。魔法命令是 Jupyter 提供的一組擴(kuò)展指令以%開頭。常用例子如下%matplotlib inline # 讓 matplotlib 圖表直接顯示在 Notebook 中%time sc.pp.pca(adata) # 顯示該行代碼的運(yùn)行耗時(shí)%load_ext autoreload %autoreload 2 # 修改模塊后自動(dòng)重新加載適合調(diào)試自己的 Python 模塊在單細(xì)胞分析中%matplotlib inline是最常用的魔法命令。如果不執(zhí)行這一句部分版本的 matplotlib 可能不會(huì)在 Notebook 中自動(dòng)顯示圖形。3.4 如何在其他瀏覽器打開 Notebook有同學(xué)在 Windows 上會(huì)遇到默認(rèn)瀏覽器打不開 Notebook 的情況。常見的解決辦法有三種。第一種直接復(fù)制終端輸出的 URL 到其他瀏覽器比如 Chrome 或 Edgehttp://localhost:8888/tree?tokenxxxxxxxx第二種啟動(dòng)時(shí)指定瀏覽器jupyter notebook --browser chrome第三種修改 Jupyter 配置文件。首先生成配置文件jupyter notebook --generate-config然后編輯生成的jupyter_notebook_config.py找到下面這一段并修改# 使用 Chrome 打開 c.NotebookApp.browser C:/Program Files/Google/Chrome/Application/chrome.exe修改完成后重啟 Notebook新瀏覽器配置即可生效。4. Scanpy 核心數(shù)據(jù)結(jié)構(gòu)與數(shù)據(jù)讀取4.1 AnnData 對象Scanpy 中所有分析操作都圍繞一個(gè)叫AnnData的對象展開。AnnData 可以理解為“帶注釋的數(shù)據(jù)”它包含以下幾個(gè)核心部分adata.X表達(dá)量矩陣通常是稀疏矩陣或稠密矩陣行是細(xì)胞列是基因adata.obs細(xì)胞維度的元數(shù)據(jù)比如細(xì)胞 ID、樣本信息、批次信息、聚類結(jié)果等adata.var基因維度的元數(shù)據(jù)比如基因名、基因類型、是否高變基因等adata.uns非結(jié)構(gòu)化數(shù)據(jù)存放一些分析過程中的中間結(jié)果adata.obsm細(xì)胞維度的降維結(jié)果比如 PCA 坐標(biāo)、UMAP 坐標(biāo)adata.varm基因維度的降維結(jié)果adata.layers多個(gè)表達(dá)量矩陣層比如原始 counts 和歸一化后的數(shù)據(jù)可以同時(shí)保存。理解 AnnData 的結(jié)構(gòu)是正確調(diào)用 Scanpy 各種函數(shù)的前提。很多報(bào)錯(cuò)都出在“想要的數(shù)據(jù)放錯(cuò)了位置”例如有些同學(xué)在聚類結(jié)束后找不到聚類標(biāo)簽其實(shí)它就保存在adata.obs[leiden]中。4.2 數(shù)據(jù)讀取方式Scanpy 支持讀取多種單細(xì)胞數(shù)據(jù)格式常見的有import scanpy as sc # 讀取 10X 的 h5 文件 adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) # 讀取 10X 的 mtx 目錄 adata sc.read_10x_mtx(filtered_feature_bc_matrix/) # 讀取 h5ad 文件 adata sc.read_h5ad(data.h5ad) # 讀取 CSV / TSV / txt 表達(dá)矩陣 adata sc.read_csv(matrix.csv) adata sc.read_tsv(matrix.tsv)如果是讀取文本格式的矩陣通常還需要手動(dòng)構(gòu)造 AnnData 對象。因?yàn)椴煌臄?shù)據(jù)來源列名和行名定義不一樣最穩(wěn)妥的方式是先了解自己數(shù)據(jù)的格式再選擇合適的讀取函數(shù)。4.3 常用質(zhì)量控制指標(biāo)單細(xì)胞數(shù)據(jù)在下游分析前必須經(jīng)過嚴(yán)格的質(zhì)量控制否則聚類結(jié)果會(huì)被低質(zhì)量細(xì)胞帶偏。常用的 QC 指標(biāo)有三個(gè)每個(gè)細(xì)胞檢測到的基因數(shù)n_genes_by_counts每個(gè)細(xì)胞的總 UMI 數(shù)total_counts每個(gè)細(xì)胞的線粒體基因表達(dá)比例pct_counts_mt。線粒體基因比例高通常意味著細(xì)胞質(zhì) RNA 丟失嚴(yán)重細(xì)胞可能已經(jīng)破裂或?yàn)l臨死亡這類細(xì)胞應(yīng)該過濾掉。在 Scanpy 中可以用以下方式添加線粒體基因比例指標(biāo)adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)對于人類數(shù)據(jù)線粒體基因以MT-開頭小鼠數(shù)據(jù)則以mt-開頭具體前綴要根據(jù)物種確認(rèn)。5. 單細(xì)胞數(shù)據(jù)分析完整實(shí)戰(zhàn)本節(jié)我們使用 Scanpy 官方示例數(shù)據(jù)集pbmc3k即 10X Genomics 提供的 2700 個(gè)外周血單核細(xì)胞數(shù)據(jù)。這個(gè)數(shù)據(jù)集規(guī)模適中是 Scanpy 文檔中最經(jīng)典的入門數(shù)據(jù)。5.1 加載數(shù)據(jù)在 Notebook 的代碼 Cell 中輸入import scanpy as sc sc.settings.set_figure_params(dpi100, facecolorwhite) adata sc.datasets.pbmc3k()第一次運(yùn)行時(shí)Scanpy 會(huì)自動(dòng)下載數(shù)據(jù)并緩存到本地。pbmc3k()返回的是已經(jīng)經(jīng)過初篩的 AnnData 對象包含 2700 個(gè)細(xì)胞和 32738 個(gè)基因。查看數(shù)據(jù)結(jié)構(gòu)adata輸出結(jié)果會(huì)顯示 AnnData 的基本信息例如AnnData object with n_obs × n_vars 2700 × 32738 var: gene_ids這里的n_obs是細(xì)胞數(shù)n_vars是基因數(shù)。5.2 質(zhì)量控制與過濾先計(jì)算線粒體基因比例adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)繪制 QC 指標(biāo)分布圖直觀判斷閾值sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue)運(yùn)行后會(huì)顯示三個(gè)小提琴圖。根據(jù)圖形中分布情況確定過濾閾值sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 5, :] adata adata[adata.obs.n_genes_by_counts 2500, :]這段代碼的含義是filter_cells過濾掉檢測到的基因數(shù)少于 200 的細(xì)胞這類細(xì)胞可能是空液滴或質(zhì)量較差filter_genes過濾掉在超過 3 個(gè)細(xì)胞中表達(dá)的基因去掉在數(shù)據(jù)集中幾乎不表達(dá)的基因保留線粒體基因比例小于 5% 的細(xì)胞保留基因數(shù)少于 2500 的細(xì)胞避免雙細(xì)胞或異常高表達(dá)細(xì)胞干擾分析。執(zhí)行完畢后再次查看adata發(fā)現(xiàn)細(xì)胞數(shù)和基因數(shù)都減少了。5.3 歸一化、對數(shù)化與高變基因表達(dá)量的原始 counts 數(shù)據(jù)受測序深度影響很大因此需要做歸一化使得不同細(xì)胞之間可比。Scanpy 中的標(biāo)準(zhǔn)做法是先按每個(gè)細(xì)胞的總 counts 縮放再取對數(shù)sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)normalize_total把每個(gè)細(xì)胞的 counts 總數(shù)縮放到1e4log1p是對每個(gè)值計(jì)算log(1 x)。對數(shù)化之后表達(dá)量的分布更接近正態(tài)分布有利于后續(xù) PCA 等算法的計(jì)算。接著篩選高變基因。高變基因是那些在細(xì)胞之間表達(dá)差異顯著的基因它們攜帶了區(qū)分細(xì)胞類型的主要信息sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5)篩選完成后數(shù)據(jù)集中會(huì)多出highly_variable這一列基因注釋。我們可以把后續(xù)主成分分析限定在高變基因上以降低計(jì)算量并減少噪聲adata adata[:, adata.var.highly_variable]5.4 主成分分析與鄰接圖首先將高變基因的表達(dá)量進(jìn)行線性回歸去除 counts 和線粒體比例的影響然后進(jìn)行主成分分析PCAsc.pp.regress_out(adata, [total_counts, pct_counts_mt]) sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) sc.pl.pca_variance_ratio(adata, n_pcs50)regress_out的作用是去除技術(shù)因素帶來的干擾scale則是將每個(gè)基因的表達(dá)標(biāo)準(zhǔn)化到零均值、單位方差。pca_variance_ratio圖中曲線會(huì)隨著主成分序號(hào)增加而下降我們通常選取拐點(diǎn)前的主成分?jǐn)?shù)量用于后續(xù)分析這里直接使用默認(rèn)的 50 個(gè)主成分也可以。接著基于 PCA 結(jié)果構(gòu)建鄰接圖計(jì)算細(xì)胞與細(xì)胞之間的相似度sc.pp.neighbors(adata, n_neighbors10, n_pcs40)n_pcs表示使用前多少個(gè)主成分計(jì)算鄰接關(guān)系。之后運(yùn)行 UMAP 降維和 Leiden 聚類sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.8)聚類結(jié)果會(huì)寫入adata.obs[leiden]。繪制 UMAP 圖不同顏色即代表不同的聚類群體sc.pl.umap(adata, color[leiden], legend_locon data)5.5 Marker 基因與細(xì)胞類型注釋聚類完成后我們通常需要找出每個(gè) cluster 的特異性高表達(dá)基因也就是 marker 基因并結(jié)合文獻(xiàn)知識(shí)判斷每個(gè) cluster 屬于什么細(xì)胞類型。sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)運(yùn)行后會(huì)得到一個(gè)多面板圖表每個(gè)面板是一個(gè) cluster縱軸排列的是該 cluster 的 top marker 基因。我們可以進(jìn)一步查看某個(gè) cluster 的 top 基因表result adata.uns[rank_genes_groups] table sc.get.rank_genes_groups_df(adata, group0) table.head(10)在 pbmc3k 數(shù)據(jù)中cluster 0 的 marker 基因通常包括IL7R、S100A8等結(jié)合已知的免疫細(xì)胞 marker可以大致推斷T 細(xì)胞CD3D、IL7RB 細(xì)胞MS4A1、CD79A單核細(xì)胞LYZ、S100A8NK 細(xì)胞NKG7、GNLY把這些注釋信息寫入adata.obsnew_cluster_names { 0: T cells, 1: Monocytes } adata.obs[cell_type] adata.obs[leiden].map(new_cluster_names).astype(category) sc.pl.umap(adata, colorcell_type)5.6 完整分析流程代碼匯總為了方便復(fù)制運(yùn)行這里把整個(gè)流程整合在一個(gè) Notebook Cell 中import scanpy as sc sc.settings.set_figure_params(dpi100, facecolorwhite) # 1. 加載數(shù)據(jù) adata sc.datasets.pbmc3k() # 2. 質(zhì)量控制 adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue) sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue) sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 5, :] adata adata[adata.obs.n_genes_by_counts 2500, :] # 3. 歸一化與高變基因 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5) adata adata[:, adata.var.highly_variable] # 4. 回歸、標(biāo)準(zhǔn)化與 PCA sc.pp.regress_out(adata, [total_counts, pct_counts_mt]) sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) sc.pl.pca_variance_ratio(adata, n_pcs50) # 5. 鄰接圖、UMAP 與聚類 sc.pp.neighbors(adata, n_neighbors10, n_pcs40) sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.8) sc.pl.umap(adata, color[leiden], legend_locon data) # 6. Marker 基因 sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)在 Notebook 中逐段運(yùn)行時(shí)建議每段代碼中間插入 Markdown Cell 寫清楚分析說明這樣最終的文件就是一份可復(fù)現(xiàn)、可分享的分析報(bào)告。6. 常見問題與排查思路在使用 Jupyter Notebook 和 Scanpy 的過程中有一些高頻問題值得單獨(dú)整理。問題現(xiàn)象常見原因解決思路Windows 下 notebook 打開后頁面空白瀏覽器兼容問題或 Notebook 版本過舊清理瀏覽器緩存更換 Chrome/Edge升級(jí) notebook 版本瀏覽器沒有自動(dòng)彈出 Notebook 頁面瀏覽器配置或防火墻攔截手動(dòng)復(fù)制終端輸出的 URL 到瀏覽器打開pip 安裝 scanpy 后 import 報(bào)錯(cuò)依賴庫版本沖突使用 conda 安裝或升級(jí) numpy/pandas/scipy運(yùn)行代碼時(shí)內(nèi)核持續(xù)顯示 Busy當(dāng)前 Cell 數(shù)據(jù)計(jì)算量過大檢查數(shù)據(jù)規(guī)模拆分 Cell必要時(shí)重啟內(nèi)核圖表無法顯示未啟用 inline 顯示執(zhí)行%matplotlib inline或設(shè)置sc.settings.autoshowTruefilter_cells過濾后細(xì)胞數(shù)過少閾值設(shè)置過于嚴(yán)格根據(jù) violin 圖重新觀察分布調(diào)整閾值Leiden 聚類結(jié)果不穩(wěn)定resolution 參數(shù)不合適調(diào)整 resolution 在 0.2 到 1.2 之間多次嘗試讀取自己的 h5 文件報(bào)錯(cuò)文件版本與 read_10x_h5 不兼容確認(rèn)文件是否為 10X 格式必要時(shí)改用 read_10x_mtx6.1 Notebook 空白頁問題的詳細(xì)排查在 Windows 上很多同學(xué)執(zhí)行jupyter notebook之后瀏覽器打開的是空白頁面。這個(gè)問題通常與瀏覽器緩存、安全軟件攔截或 Notebook 版本過舊有關(guān)。排查順序建議如下強(qiáng)制刷新頁面按住Ctrl F5清除當(dāng)前頁面緩存復(fù)制終端中完整的 URL包括token參數(shù)打開新的無痕窗口訪問關(guān)閉殺毒軟件或防火墻的網(wǎng)頁攔截功能局部端口8888放行升級(jí) Notebook 相關(guān)包pip install --upgrade notebook如果仍然無法解決可以嘗試在啟動(dòng)時(shí)指定 IP 和端口jupyter notebook --ip127.0.0.1 --port8889 --no-browser然后手動(dòng)訪問http://127.0.0.1:8889/tree。6.2 內(nèi)核連接失敗的處理如果在 Notebook 中執(zhí)行代碼時(shí)提示A connection to the notebook server could not be established通常是因?yàn)閮?nèi)核進(jìn)程異常退出。解決方法是在 Notebook 的 Kernel 菜單中選擇 Restart Kernel清理該 Notebook 的運(yùn)行時(shí)狀態(tài)如果頻繁出現(xiàn)可以重啟 Jupyter 服務(wù)。6.3 Scanpy 安裝中的依賴問題Scanpy 很依賴numpy、scipy、pandas、matplotlib、h5py等庫舊版本可能因?yàn)?API 改變而報(bào)錯(cuò)。安裝 Scanpy 時(shí)盡量直接使用 conda-forge 頻道conda 會(huì)自動(dòng)解析依賴版本。如果使用 pip 安裝后出現(xiàn)ImportError: cannot import name X的情況常見做法是把這些基礎(chǔ)庫一起升級(jí)到較新的兼容版本pip install --upgrade numpy pandas scipy h5py7. 最佳實(shí)踐與工程建議7.1 用虛擬環(huán)境隔離項(xiàng)目依賴單細(xì)胞分析項(xiàng)目周期長、依賴多且不同項(xiàng)目之間對 Scanpy、NumPy 等庫的版本要求可能不一樣。建議每個(gè)項(xiàng)目都使用獨(dú)立的 conda 虛擬環(huán)境并且把環(huán)境依賴導(dǎo)出保存conda env export environment.yml這樣無論是換電腦還是項(xiàng)目交接都能通過一份環(huán)境文件快速還原分析環(huán)境。7.2 Notebook 中記錄分析參數(shù)單細(xì)胞分析非常強(qiáng)調(diào)可復(fù)現(xiàn)性。每個(gè)關(guān)鍵步驟的參數(shù)例如 QC 閾值、聚類分辨率、高變基因篩選條件都應(yīng)該在 Notebook 的 Markdown Cell 中記錄下來。建議每個(gè)步驟采用“目的 參數(shù) 結(jié)果說明”的格式。例如## 質(zhì)量控制 - 目的過濾掉低質(zhì)量細(xì)胞和雙細(xì)胞。 - 參數(shù)min_genes200pct_counts_mt 5n_genes 2500。 - 結(jié)果細(xì)胞數(shù)從 2700 降到 2638。這樣的筆記在文章發(fā)表、項(xiàng)目匯報(bào)和校友合作時(shí)價(jià)值極高。7.3 保存中間結(jié)果Scanpy 的分析鏈路上從原始數(shù)據(jù)到最終聚類結(jié)果會(huì)經(jīng)過多個(gè)步驟而每一步的參數(shù)調(diào)整可能產(chǎn)生不同的中間結(jié)果。建議在關(guān)鍵節(jié)點(diǎn)使用write_h5ad保存結(jié)果文件adata.write_h5ad(data/qc.h5ad) adata.write_h5ad(data/clustered.h5ad)保存 h5ad 文件的好處是它同時(shí)保留了表達(dá)矩陣、obs 注釋、uns 中間結(jié)果和降維坐標(biāo)后續(xù)想要重新畫圖不需要重新跑前面的全部計(jì)算。7.4 計(jì)算資源的合理使用pbmc3k 只有 2700 個(gè)細(xì)胞在本地電腦上運(yùn)行毫無壓力。但如果數(shù)據(jù)量達(dá)到十萬甚至百萬級(jí)細(xì)胞就需要考慮計(jì)算資源問題。以下幾個(gè)方面可以優(yōu)化過濾步驟盡量靠前減少后續(xù)分析的細(xì)胞數(shù)和基因數(shù)PCA 和 neighbor 計(jì)算時(shí)可指定n_pcs降低維度對超大數(shù)據(jù)集考慮使用 GPU 加速版本的 Scanpy 或分塊處理策略操作大矩陣時(shí)保存中間結(jié)果避免重復(fù)計(jì)算占用內(nèi)存。7.5 數(shù)據(jù)來源與生物意義的驗(yàn)證Scanpy 分析結(jié)果最終要回到實(shí)驗(yàn)和生物學(xué)意義中驗(yàn)證。每次聚類完成后不能只看 UMAP 圖還應(yīng)該檢查 marker 基因的表達(dá)情況結(jié)合文獻(xiàn)確認(rèn)細(xì)胞類型注釋是否合理。如果某個(gè) cluster 的 marker 基因列表混雜了多種已知細(xì)胞類型的特征可能說明聚類分辨率過高或數(shù)據(jù)中存在批次效應(yīng)需要進(jìn)一步調(diào)整參數(shù)或進(jìn)行批次校正。8. 總結(jié)與下一步學(xué)習(xí)建議到這里我們已經(jīng)完成了從 Jupyter Notebook 基礎(chǔ)操作到 Scanpy 單細(xì)胞分析入門全流程的梳理。現(xiàn)在你應(yīng)該能夠獨(dú)立完成以下工作安裝并啟動(dòng) Jupyter Notebook熟練使用 Cell、Markdown 和快捷鍵理解 Notebook 與 JupyterLab 的區(qū)別根據(jù)需求選擇合適工具了解 AnnData 的五大結(jié)構(gòu)字段跑通 pbmc3k 數(shù)據(jù)集的 QC、歸一化、高變基因篩選、PCA、鄰接圖、UMAP、Leiden 聚類和 marker 基因分析流程對 Windows 下 Notebook 空白頁、內(nèi)核連接失敗、依賴沖突等常見問題進(jìn)行定位和修復(fù)。下一步可以嘗試把這套流程應(yīng)用到自己的數(shù)據(jù)上。先從 10X 官網(wǎng)下載一個(gè)公開數(shù)據(jù)集練習(xí)read_10x_h5和read_10x_mtx兩種讀取方式之后再為分析 Notebook 補(bǔ)全細(xì)胞類型注釋嘗試使用sc.tl.diffmap做擴(kuò)散圖或者用sc.tl.score_genes給細(xì)胞打分熟悉更多 Scanpy 分析模塊。單細(xì)胞分析的學(xué)習(xí)曲線雖然比普通數(shù)據(jù)處理陡峭一些但結(jié)合 Jupyter Notebook 的交互式環(huán)境完全可以把整個(gè)路線逐步打通。如果在環(huán)境配置或代碼運(yùn)行中遇到報(bào)錯(cuò)優(yōu)先去官方文檔查對應(yīng)函數(shù)的參數(shù)說明社區(qū)的討論帖也值得多翻一翻。把這套 pbmc3k 流程真正吃透之后再上真實(shí)數(shù)據(jù)你會(huì)順手很多。 Jupyter Notebook與Scanpy單細(xì)胞數(shù)據(jù)分析實(shí)戰(zhàn)指南
1. 背景與核心概念1.1 Jupyter Notebook 是什么Jupyter Notebook 是一個(gè)基于 Web 的交互式開發(fā)環(huán)境允許我們在瀏覽器中編寫代碼、運(yùn)行代碼、查看輸出結(jié)果同時(shí)還能在代碼之間插入 Markdown 說明文字、公式、圖片和表格。最早它作為 IPython Notebook 項(xiàng)目誕生后來逐步發(fā)展為支持 Python、R、Julia 等多種語言的項(xiàng)目也因此得名 JupyterJulia Python R 的合稱。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域Jupyter Notebook 幾乎是標(biāo)配工具。原因很簡單單細(xì)胞數(shù)據(jù)的分析過程是一條長鏈路包括數(shù)據(jù)讀取、質(zhì)量控制、歸一化、降維、聚類、注釋等多個(gè)環(huán)節(jié)每一步之間高度依賴人工判斷。Notebook 這種“一段代碼 一段結(jié)果 一段說明”的組織方式非常適合探索性數(shù)據(jù)分析和流程復(fù)現(xiàn)。你可以在同一個(gè)文件中既完成代碼調(diào)試又留下對參數(shù)選擇的理解記錄后續(xù)想要回溯分析思路也會(huì)輕松很多。1.2 Jupyter Notebook 與 JupyterLab 的區(qū)別很多初學(xué)者在安裝時(shí)會(huì)碰到一個(gè)搜索熱點(diǎn)Jupyter Notebook 和 JupyterLab 到底有什么區(qū)別簡單來說JupyterLab 是 Jupyter Notebook 的下一代交互式開發(fā)界面相當(dāng)于把 Notebook 的能力放到了一個(gè)更像 IDE 的工作臺(tái)中。JupyterLab 支持多標(biāo)簽頁、拖拽布局、文件管理器、終端、代碼控制臺(tái)等功能在同一個(gè)窗口中可以同時(shí)打開多個(gè) Notebook、查看數(shù)據(jù)文件、運(yùn)行終端命令。而經(jīng)典 Jupyter Notebook 則更聚焦于單個(gè) Notebook 文件的編輯和運(yùn)行界面更加簡潔適合只專注分析任務(wù)的時(shí)候使用。兩者使用同一個(gè) Notebook 文件格式.ipynb運(yùn)行時(shí)也共享同一套內(nèi)核機(jī)制。也就是說你在經(jīng)典 Notebook 里寫的代碼在 JupyterLab 里也能正常打開運(yùn)行?，F(xiàn)在 Anaconda 默認(rèn)安裝的是 JupyterLab 和 Notebook 兩套組件啟動(dòng)命令分別對應(yīng)jupyter notebook jupyter lab如果只是在學(xué)習(xí)單細(xì)胞分析的基礎(chǔ)流程使用哪個(gè)都可以如果希望邊寫代碼邊觀察文件目錄、邊調(diào)試邊查看圖表JupyterLab 的體驗(yàn)會(huì)更好。本文的代碼在兩種界面下均可直接運(yùn)行。1.3 Scanpy 是什么Scanpy 是一個(gè)基于 Python 的高性能單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析工具包全稱是 Single-Cell Analysis in Python。它的核心數(shù)據(jù)結(jié)構(gòu)基于 AnnData能夠高效處理大規(guī)模稀疏矩陣支持從數(shù)據(jù)讀取、質(zhì)量控制、標(biāo)準(zhǔn)化、特征選擇、降維、聚類到差異分析、可視化以及偽時(shí)間分析等完整流程。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域R 語言生態(tài)中的 Seurat 一直是主流工具而 Scanpy 的出現(xiàn)讓 Python 用戶也能完成同等級(jí)的分析工作。相比 SeuratScanpy 有以下幾個(gè)特點(diǎn)構(gòu)建在 Anndata、NumPy、SciPy 和 Pandas 之上與 Python 數(shù)據(jù)科學(xué)生態(tài)天然融合內(nèi)存占用控制較好適合處理幾萬甚至更多細(xì)胞的單細(xì)胞數(shù)據(jù)集分析流程高度模塊化每個(gè)步驟對應(yīng)一個(gè)明確的 API支持豐富的可視化圖形輸出可直接嵌入 Notebook 中展示。對于已經(jīng)熟悉 Python 的數(shù)據(jù)分析人員來說Scanpy 是入門單細(xì)胞數(shù)據(jù)分析非常合適的切入點(diǎn)。我們在本文中會(huì)使用一個(gè)經(jīng)典示例數(shù)據(jù)集 pbmc3k完整跑一遍單細(xì)胞轉(zhuǎn)錄組的標(biāo)準(zhǔn)分析流程。1.4 為什么把 Jupyter Notebook 與 Scanpy 放在一起學(xué)習(xí)單細(xì)胞數(shù)據(jù)分析并不是“把數(shù)據(jù)丟進(jìn)一個(gè)函數(shù)就能出結(jié)果”的事。它需要分析者反復(fù)調(diào)整閾值參數(shù)、查看中間圖表、確認(rèn)聚類效果這一過程天然適合 Jupyter Notebook 這種交互式環(huán)境。用 Jupyter Notebook 跑 Scanpy每個(gè)分析步驟都能立即看到輸出和圖表比如質(zhì)量控制后的細(xì)胞數(shù)變化、UMAP 降維后的聚類分布、marker 基因的表達(dá)情況。遇到某個(gè)步驟結(jié)果不理想可以直接修改參數(shù)重新運(yùn)行不需要從頭啟動(dòng)整個(gè)腳本。這種探索式的工作流正是 Scanpy 教程和單細(xì)胞分析論文代碼大量采用 Notebook 形式的原因。所以本文的安排是先掌握 Jupyter Notebook 的基礎(chǔ)操作再用它跑通 Scanpy 的單細(xì)胞分析流程最終形成一個(gè)可以直接照抄、修改、復(fù)用到自己數(shù)據(jù)上的分析模板。2. 環(huán)境準(zhǔn)備與版本說明2.1 安裝 AnacondaScanpy 依賴的第三方庫較多手動(dòng)用 pip 逐個(gè)安裝容易遇到依賴沖突所以我們首先推薦通過 Anaconda 來管理 Python 環(huán)境。Anaconda 自帶 Python 解釋器、Jupyter Notebook、常用科學(xué)計(jì)算庫以及包管理器 conda安裝完成后即可快速使用。到 Anaconda 官網(wǎng)下載對應(yīng)操作系統(tǒng)的安裝包安裝過程中保持默認(rèn)選項(xiàng)即可。安裝完成后打開命令行窗口輸入conda --version如果能看到 conda 版本號(hào)說明安裝成功。為了保證環(huán)境干凈建議單獨(dú)創(chuàng)建一個(gè)用于單細(xì)胞分析的虛擬環(huán)境conda create -n scanpy-env python3.9 -y conda activate scanpy-env關(guān)于 Python 版本需要根據(jù)你的項(xiàng)目實(shí)際情況調(diào)整這里選擇 3.9 是兼容性較好的通用版本。如果你的電腦上已經(jīng)存在其他 Python 項(xiàng)目使用獨(dú)立虛擬環(huán)境可以避免庫與庫之間的版本沖突這是單細(xì)胞分析項(xiàng)目里非常重要的一步。2.2 安裝 Scanpy激活虛擬環(huán)境后使用 conda 安裝 Scanpyconda install -c conda-forge scanpy python-igraph leidenalg這里同時(shí)安裝了python-igraph和leidenalg它們用于 Leiden 聚類算法。Leiden 是目前單細(xì)胞聚類中較推薦的方法相比早期的 Louvain 算法社區(qū)劃分更穩(wěn)定對分辨率的調(diào)節(jié)也更友好。如果 conda 安裝速度慢也可以使用 pippip install scanpy pip install leidenalg不過 pip 安裝時(shí)需要注意 NumPy、Pandas、SciPy 等底層庫的版本兼容性如果出現(xiàn)版本報(bào)錯(cuò)建議仍然回到 conda 安裝方式。安裝完成后在命令行輸入python -c import scanpy as sc; print(sc.__version__)能正常打印版本號(hào)就說明環(huán)境已經(jīng)準(zhǔn)備好了。2.3 啟動(dòng) Jupyter Notebook在scanpy-env環(huán)境中啟動(dòng) Notebookjupyter notebook啟動(dòng)后終端會(huì)輸出類似下面的信息[I 2025-01-01 10:00:00.123 NotebookApp] Serving notebooks from local directory: /Users/xxx/scanpy-tutorial [I 2025-01-01 10:00:00.125 NotebookApp] Jupyter Notebook 6.5.4 is running at: [I 2025-01-01 10:00:00.125 NotebookApp] http://localhost:8888/?tokenxxxxxxxx瀏覽器會(huì)自動(dòng)打開 Notebook 的主頁面。如果瀏覽器沒有自動(dòng)打開可以把終端里顯示的http://localhost:8888/?tokenxxx地址手動(dòng)復(fù)制到瀏覽器地址欄。在 Notebook 主頁面的右上角點(diǎn)擊 New - Python 3即可新建一個(gè) Notebook 文件。建議把工作目錄切換到專門存放單細(xì)胞分析代碼的文件夾例如scanpy-tutorial這樣后續(xù)導(dǎo)入數(shù)據(jù)文件時(shí)路徑管理更方便。2.4 示例項(xiàng)目結(jié)構(gòu)本文的實(shí)戰(zhàn)代碼按以下目錄組織scanpy-tutorial/ ├── pbmc3k.h5ad # 示例數(shù)據(jù)運(yùn)行過程中自動(dòng)下載 └── scanpy_analysis.ipynb # 分析 Notebook 文件實(shí)際操作中你的文件路徑不一定要完全一致但建議保持?jǐn)?shù)據(jù)和 Notebook 在同一級(jí)目錄減少路徑出錯(cuò)的可能。3. Jupyter Notebook 核心操作3.1 Cell 與兩種模式Notebook 文件由一個(gè)個(gè)單元格Cell組成每個(gè) Cell 可以存放代碼也可以存放 Markdown 文檔。運(yùn)行 Cell 時(shí)Jupyter 會(huì)把代碼發(fā)送給內(nèi)核執(zhí)行并把結(jié)果展示在 Cell 下方。Notebook 有兩種模式命令模式按 Esc 進(jìn)入此時(shí)鍵盤快捷鍵作用于整個(gè) Notebook可以刪除、復(fù)制、移動(dòng) Cell編輯模式按 Enter 進(jìn)入此時(shí)可以編輯當(dāng)前 Cell 中的代碼或文字。新建的 Cell 默認(rèn)是代碼類型如果要寫說明文字需要把 Cell 切換為 Markdown 類型。切換方式是在命令模式下按M切回代碼類型則按Y。3.2 常用快捷鍵掌握以下快捷鍵可以明顯提升操作效率快捷鍵作用Shift Enter運(yùn)行當(dāng)前 Cell 并跳轉(zhuǎn)到下一個(gè) CellCtrl Enter運(yùn)行當(dāng)前 Cell不跳轉(zhuǎn)Esc進(jìn)入命令模式A在當(dāng)前 Cell 上方插入新 CellB在當(dāng)前 Cell 下方插入新 CellD D連續(xù)按兩次 D刪除當(dāng)前 CellM將當(dāng)前 Cell 切換為 MarkdownY將當(dāng)前 Cell 切換為代碼Shift Tab查看函數(shù)或?qū)ο蟮奈臋n摘要其中 Shift Tab 在調(diào)用 Scanpy 函數(shù)時(shí)尤其實(shí)用。比如你忘了sc.pp.filter_cells的參數(shù)含義把光標(biāo)放到函數(shù)名上按 Shift Tab就能快速看到函數(shù)簽名和參數(shù)說明不用頻繁去查文檔。3.3 Markdown 與魔法命令Markdown Cell 中可以使用標(biāo)準(zhǔn)的 Markdown 語法包括#標(biāo)題、**加粗**、-列表、行內(nèi)代碼等。此外 Jupyter 還支持 LaTeX 數(shù)學(xué)公式用$$包裹即可。在單細(xì)胞分析筆記中我習(xí)慣在每個(gè)步驟前用 Markdown Cell 寫清楚“這一步在做什么、為什么做”分析結(jié)束之后再回頭看整個(gè)文件就像一份可執(zhí)行的分析報(bào)告。魔法命令是 Jupyter 提供的一組擴(kuò)展指令以%開頭。常用例子如下%matplotlib inline # 讓 matplotlib 圖表直接顯示在 Notebook 中%time sc.pp.pca(adata) # 顯示該行代碼的運(yùn)行耗時(shí)%load_ext autoreload %autoreload 2 # 修改模塊后自動(dòng)重新加載適合調(diào)試自己的 Python 模塊在單細(xì)胞分析中%matplotlib inline是最常用的魔法命令。如果不執(zhí)行這一句部分版本的 matplotlib 可能不會(huì)在 Notebook 中自動(dòng)顯示圖形。3.4 如何在其他瀏覽器打開 Notebook有同學(xué)在 Windows 上會(huì)遇到默認(rèn)瀏覽器打不開 Notebook 的情況。常見的解決辦法有三種。第一種直接復(fù)制終端輸出的 URL 到其他瀏覽器比如 Chrome 或 Edgehttp://localhost:8888/tree?tokenxxxxxxxx第二種啟動(dòng)時(shí)指定瀏覽器jupyter notebook --browser chrome第三種修改 Jupyter 配置文件。首先生成配置文件jupyter notebook --generate-config然后編輯生成的jupyter_notebook_config.py找到下面這一段并修改# 使用 Chrome 打開 c.NotebookApp.browser C:/Program Files/Google/Chrome/Application/chrome.exe修改完成后重啟 Notebook新瀏覽器配置即可生效。4. Scanpy 核心數(shù)據(jù)結(jié)構(gòu)與數(shù)據(jù)讀取4.1 AnnData 對象Scanpy 中所有分析操作都圍繞一個(gè)叫AnnData的對象展開。AnnData 可以理解為“帶注釋的數(shù)據(jù)”它包含以下幾個(gè)核心部分adata.X表達(dá)量矩陣通常是稀疏矩陣或稠密矩陣行是細(xì)胞列是基因adata.obs細(xì)胞維度的元數(shù)據(jù)比如細(xì)胞 ID、樣本信息、批次信息、聚類結(jié)果等adata.var基因維度的元數(shù)據(jù)比如基因名、基因類型、是否高變基因等adata.uns非結(jié)構(gòu)化數(shù)據(jù)存放一些分析過程中的中間結(jié)果adata.obsm細(xì)胞維度的降維結(jié)果比如 PCA 坐標(biāo)、UMAP 坐標(biāo)adata.varm基因維度的降維結(jié)果adata.layers多個(gè)表達(dá)量矩陣層比如原始 counts 和歸一化后的數(shù)據(jù)可以同時(shí)保存。理解 AnnData 的結(jié)構(gòu)是正確調(diào)用 Scanpy 各種函數(shù)的前提。很多報(bào)錯(cuò)都出在“想要的數(shù)據(jù)放錯(cuò)了位置”例如有些同學(xué)在聚類結(jié)束后找不到聚類標(biāo)簽其實(shí)它就保存在adata.obs[leiden]中。4.2 數(shù)據(jù)讀取方式Scanpy 支持讀取多種單細(xì)胞數(shù)據(jù)格式常見的有import scanpy as sc # 讀取 10X 的 h5 文件 adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) # 讀取 10X 的 mtx 目錄 adata sc.read_10x_mtx(filtered_feature_bc_matrix/) # 讀取 h5ad 文件 adata sc.read_h5ad(data.h5ad) # 讀取 CSV / TSV / txt 表達(dá)矩陣 adata sc.read_csv(matrix.csv) adata sc.read_tsv(matrix.tsv)如果是讀取文本格式的矩陣通常還需要手動(dòng)構(gòu)造 AnnData 對象。因?yàn)椴煌臄?shù)據(jù)來源列名和行名定義不一樣最穩(wěn)妥的方式是先了解自己數(shù)據(jù)的格式再選擇合適的讀取函數(shù)。4.3 常用質(zhì)量控制指標(biāo)單細(xì)胞數(shù)據(jù)在下游分析前必須經(jīng)過嚴(yán)格的質(zhì)量控制否則聚類結(jié)果會(huì)被低質(zhì)量細(xì)胞帶偏。常用的 QC 指標(biāo)有三個(gè)每個(gè)細(xì)胞檢測到的基因數(shù)n_genes_by_counts每個(gè)細(xì)胞的總 UMI 數(shù)total_counts每個(gè)細(xì)胞的線粒體基因表達(dá)比例pct_counts_mt。線粒體基因比例高通常意味著細(xì)胞質(zhì) RNA 丟失嚴(yán)重細(xì)胞可能已經(jīng)破裂或?yàn)l臨死亡這類細(xì)胞應(yīng)該過濾掉。在 Scanpy 中可以用以下方式添加線粒體基因比例指標(biāo)adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)對于人類數(shù)據(jù)線粒體基因以MT-開頭小鼠數(shù)據(jù)則以mt-開頭具體前綴要根據(jù)物種確認(rèn)。5. 單細(xì)胞數(shù)據(jù)分析完整實(shí)戰(zhàn)本節(jié)我們使用 Scanpy 官方示例數(shù)據(jù)集pbmc3k即 10X Genomics 提供的 2700 個(gè)外周血單核細(xì)胞數(shù)據(jù)。這個(gè)數(shù)據(jù)集規(guī)模適中是 Scanpy 文檔中最經(jīng)典的入門數(shù)據(jù)。5.1 加載數(shù)據(jù)在 Notebook 的代碼 Cell 中輸入import scanpy as sc sc.settings.set_figure_params(dpi100, facecolorwhite) adata sc.datasets.pbmc3k()第一次運(yùn)行時(shí)Scanpy 會(huì)自動(dòng)下載數(shù)據(jù)并緩存到本地。pbmc3k()返回的是已經(jīng)經(jīng)過初篩的 AnnData 對象包含 2700 個(gè)細(xì)胞和 32738 個(gè)基因。查看數(shù)據(jù)結(jié)構(gòu)adata輸出結(jié)果會(huì)顯示 AnnData 的基本信息例如AnnData object with n_obs × n_vars 2700 × 32738 var: gene_ids這里的n_obs是細(xì)胞數(shù)n_vars是基因數(shù)。5.2 質(zhì)量控制與過濾先計(jì)算線粒體基因比例adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)繪制 QC 指標(biāo)分布圖直觀判斷閾值sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue)運(yùn)行后會(huì)顯示三個(gè)小提琴圖。根據(jù)圖形中分布情況確定過濾閾值sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 5, :] adata adata[adata.obs.n_genes_by_counts 2500, :]這段代碼的含義是filter_cells過濾掉檢測到的基因數(shù)少于 200 的細(xì)胞這類細(xì)胞可能是空液滴或質(zhì)量較差filter_genes過濾掉在超過 3 個(gè)細(xì)胞中表達(dá)的基因去掉在數(shù)據(jù)集中幾乎不表達(dá)的基因保留線粒體基因比例小于 5% 的細(xì)胞保留基因數(shù)少于 2500 的細(xì)胞避免雙細(xì)胞或異常高表達(dá)細(xì)胞干擾分析。執(zhí)行完畢后再次查看adata發(fā)現(xiàn)細(xì)胞數(shù)和基因數(shù)都減少了。5.3 歸一化、對數(shù)化與高變基因表達(dá)量的原始 counts 數(shù)據(jù)受測序深度影響很大因此需要做歸一化使得不同細(xì)胞之間可比。Scanpy 中的標(biāo)準(zhǔn)做法是先按每個(gè)細(xì)胞的總 counts 縮放再取對數(shù)sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)normalize_total把每個(gè)細(xì)胞的 counts 總數(shù)縮放到1e4log1p是對每個(gè)值計(jì)算log(1 x)。對數(shù)化之后表達(dá)量的分布更接近正態(tài)分布有利于后續(xù) PCA 等算法的計(jì)算。接著篩選高變基因。高變基因是那些在細(xì)胞之間表達(dá)差異顯著的基因它們攜帶了區(qū)分細(xì)胞類型的主要信息sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5)篩選完成后數(shù)據(jù)集中會(huì)多出highly_variable這一列基因注釋。我們可以把后續(xù)主成分分析限定在高變基因上以降低計(jì)算量并減少噪聲adata adata[:, adata.var.highly_variable]5.4 主成分分析與鄰接圖首先將高變基因的表達(dá)量進(jìn)行線性回歸去除 counts 和線粒體比例的影響然后進(jìn)行主成分分析PCAsc.pp.regress_out(adata, [total_counts, pct_counts_mt]) sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) sc.pl.pca_variance_ratio(adata, n_pcs50)regress_out的作用是去除技術(shù)因素帶來的干擾scale則是將每個(gè)基因的表達(dá)標(biāo)準(zhǔn)化到零均值、單位方差。pca_variance_ratio圖中曲線會(huì)隨著主成分序號(hào)增加而下降我們通常選取拐點(diǎn)前的主成分?jǐn)?shù)量用于后續(xù)分析這里直接使用默認(rèn)的 50 個(gè)主成分也可以。接著基于 PCA 結(jié)果構(gòu)建鄰接圖計(jì)算細(xì)胞與細(xì)胞之間的相似度sc.pp.neighbors(adata, n_neighbors10, n_pcs40)n_pcs表示使用前多少個(gè)主成分計(jì)算鄰接關(guān)系。之后運(yùn)行 UMAP 降維和 Leiden 聚類sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.8)聚類結(jié)果會(huì)寫入adata.obs[leiden]。繪制 UMAP 圖不同顏色即代表不同的聚類群體sc.pl.umap(adata, color[leiden], legend_locon data)5.5 Marker 基因與細(xì)胞類型注釋聚類完成后我們通常需要找出每個(gè) cluster 的特異性高表達(dá)基因也就是 marker 基因并結(jié)合文獻(xiàn)知識(shí)判斷每個(gè) cluster 屬于什么細(xì)胞類型。sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)運(yùn)行后會(huì)得到一個(gè)多面板圖表每個(gè)面板是一個(gè) cluster縱軸排列的是該 cluster 的 top marker 基因。我們可以進(jìn)一步查看某個(gè) cluster 的 top 基因表result adata.uns[rank_genes_groups] table sc.get.rank_genes_groups_df(adata, group0) table.head(10)在 pbmc3k 數(shù)據(jù)中cluster 0 的 marker 基因通常包括IL7R、S100A8等結(jié)合已知的免疫細(xì)胞 marker可以大致推斷T 細(xì)胞CD3D、IL7RB 細(xì)胞MS4A1、CD79A單核細(xì)胞LYZ、S100A8NK 細(xì)胞NKG7、GNLY把這些注釋信息寫入adata.obsnew_cluster_names { 0: T cells, 1: Monocytes } adata.obs[cell_type] adata.obs[leiden].map(new_cluster_names).astype(category) sc.pl.umap(adata, colorcell_type)5.6 完整分析流程代碼匯總為了方便復(fù)制運(yùn)行這里把整個(gè)流程整合在一個(gè) Notebook Cell 中import scanpy as sc sc.settings.set_figure_params(dpi100, facecolorwhite) # 1. 加載數(shù)據(jù) adata sc.datasets.pbmc3k() # 2. 質(zhì)量控制 adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue) sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue) sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 5, :] adata adata[adata.obs.n_genes_by_counts 2500, :] # 3. 歸一化與高變基因 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5) adata adata[:, adata.var.highly_variable] # 4. 回歸、標(biāo)準(zhǔn)化與 PCA sc.pp.regress_out(adata, [total_counts, pct_counts_mt]) sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) sc.pl.pca_variance_ratio(adata, n_pcs50) # 5. 鄰接圖、UMAP 與聚類 sc.pp.neighbors(adata, n_neighbors10, n_pcs40) sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.8) sc.pl.umap(adata, color[leiden], legend_locon data) # 6. Marker 基因 sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)在 Notebook 中逐段運(yùn)行時(shí)建議每段代碼中間插入 Markdown Cell 寫清楚分析說明這樣最終的文件就是一份可復(fù)現(xiàn)、可分享的分析報(bào)告。6. 常見問題與排查思路在使用 Jupyter Notebook 和 Scanpy 的過程中有一些高頻問題值得單獨(dú)整理。問題現(xiàn)象常見原因解決思路Windows 下 notebook 打開后頁面空白瀏覽器兼容問題或 Notebook 版本過舊清理瀏覽器緩存更換 Chrome/Edge升級(jí) notebook 版本瀏覽器沒有自動(dòng)彈出 Notebook 頁面瀏覽器配置或防火墻攔截手動(dòng)復(fù)制終端輸出的 URL 到瀏覽器打開pip 安裝 scanpy 后 import 報(bào)錯(cuò)依賴庫版本沖突使用 conda 安裝或升級(jí) numpy/pandas/scipy運(yùn)行代碼時(shí)內(nèi)核持續(xù)顯示 Busy當(dāng)前 Cell 數(shù)據(jù)計(jì)算量過大檢查數(shù)據(jù)規(guī)模拆分 Cell必要時(shí)重啟內(nèi)核圖表無法顯示未啟用 inline 顯示執(zhí)行%matplotlib inline或設(shè)置sc.settings.autoshowTruefilter_cells過濾后細(xì)胞數(shù)過少閾值設(shè)置過于嚴(yán)格根據(jù) violin 圖重新觀察分布調(diào)整閾值Leiden 聚類結(jié)果不穩(wěn)定resolution 參數(shù)不合適調(diào)整 resolution 在 0.2 到 1.2 之間多次嘗試讀取自己的 h5 文件報(bào)錯(cuò)文件版本與 read_10x_h5 不兼容確認(rèn)文件是否為 10X 格式必要時(shí)改用 read_10x_mtx6.1 Notebook 空白頁問題的詳細(xì)排查在 Windows 上很多同學(xué)執(zhí)行jupyter notebook之后瀏覽器打開的是空白頁面。這個(gè)問題通常與瀏覽器緩存、安全軟件攔截或 Notebook 版本過舊有關(guān)。排查順序建議如下強(qiáng)制刷新頁面按住Ctrl F5清除當(dāng)前頁面緩存復(fù)制終端中完整的 URL包括token參數(shù)打開新的無痕窗口訪問關(guān)閉殺毒軟件或防火墻的網(wǎng)頁攔截功能局部端口8888放行升級(jí) Notebook 相關(guān)包pip install --upgrade notebook如果仍然無法解決可以嘗試在啟動(dòng)時(shí)指定 IP 和端口jupyter notebook --ip127.0.0.1 --port8889 --no-browser然后手動(dòng)訪問http://127.0.0.1:8889/tree。6.2 內(nèi)核連接失敗的處理如果在 Notebook 中執(zhí)行代碼時(shí)提示A connection to the notebook server could not be established通常是因?yàn)閮?nèi)核進(jìn)程異常退出。解決方法是在 Notebook 的 Kernel 菜單中選擇 Restart Kernel清理該 Notebook 的運(yùn)行時(shí)狀態(tài)如果頻繁出現(xiàn)可以重啟 Jupyter 服務(wù)。6.3 Scanpy 安裝中的依賴問題Scanpy 很依賴numpy、scipy、pandas、matplotlib、h5py等庫舊版本可能因?yàn)?API 改變而報(bào)錯(cuò)。安裝 Scanpy 時(shí)盡量直接使用 conda-forge 頻道conda 會(huì)自動(dòng)解析依賴版本。如果使用 pip 安裝后出現(xiàn)ImportError: cannot import name X的情況常見做法是把這些基礎(chǔ)庫一起升級(jí)到較新的兼容版本pip install --upgrade numpy pandas scipy h5py7. 最佳實(shí)踐與工程建議7.1 用虛擬環(huán)境隔離項(xiàng)目依賴單細(xì)胞分析項(xiàng)目周期長、依賴多且不同項(xiàng)目之間對 Scanpy、NumPy 等庫的版本要求可能不一樣。建議每個(gè)項(xiàng)目都使用獨(dú)立的 conda 虛擬環(huán)境并且把環(huán)境依賴導(dǎo)出保存conda env export environment.yml這樣無論是換電腦還是項(xiàng)目交接都能通過一份環(huán)境文件快速還原分析環(huán)境。7.2 Notebook 中記錄分析參數(shù)單細(xì)胞分析非常強(qiáng)調(diào)可復(fù)現(xiàn)性。每個(gè)關(guān)鍵步驟的參數(shù)例如 QC 閾值、聚類分辨率、高變基因篩選條件都應(yīng)該在 Notebook 的 Markdown Cell 中記錄下來。建議每個(gè)步驟采用“目的 參數(shù) 結(jié)果說明”的格式。例如## 質(zhì)量控制 - 目的過濾掉低質(zhì)量細(xì)胞和雙細(xì)胞。 - 參數(shù)min_genes200pct_counts_mt 5n_genes 2500。 - 結(jié)果細(xì)胞數(shù)從 2700 降到 2638。這樣的筆記在文章發(fā)表、項(xiàng)目匯報(bào)和校友合作時(shí)價(jià)值極高。7.3 保存中間結(jié)果Scanpy 的分析鏈路上從原始數(shù)據(jù)到最終聚類結(jié)果會(huì)經(jīng)過多個(gè)步驟而每一步的參數(shù)調(diào)整可能產(chǎn)生不同的中間結(jié)果。建議在關(guān)鍵節(jié)點(diǎn)使用write_h5ad保存結(jié)果文件adata.write_h5ad(data/qc.h5ad) adata.write_h5ad(data/clustered.h5ad)保存 h5ad 文件的好處是它同時(shí)保留了表達(dá)矩陣、obs 注釋、uns 中間結(jié)果和降維坐標(biāo)后續(xù)想要重新畫圖不需要重新跑前面的全部計(jì)算。7.4 計(jì)算資源的合理使用pbmc3k 只有 2700 個(gè)細(xì)胞在本地電腦上運(yùn)行毫無壓力。但如果數(shù)據(jù)量達(dá)到十萬甚至百萬級(jí)細(xì)胞就需要考慮計(jì)算資源問題。以下幾個(gè)方面可以優(yōu)化過濾步驟盡量靠前減少后續(xù)分析的細(xì)胞數(shù)和基因數(shù)PCA 和 neighbor 計(jì)算時(shí)可指定n_pcs降低維度對超大數(shù)據(jù)集考慮使用 GPU 加速版本的 Scanpy 或分塊處理策略操作大矩陣時(shí)保存中間結(jié)果避免重復(fù)計(jì)算占用內(nèi)存。7.5 數(shù)據(jù)來源與生物意義的驗(yàn)證Scanpy 分析結(jié)果最終要回到實(shí)驗(yàn)和生物學(xué)意義中驗(yàn)證。每次聚類完成后不能只看 UMAP 圖還應(yīng)該檢查 marker 基因的表達(dá)情況結(jié)合文獻(xiàn)確認(rèn)細(xì)胞類型注釋是否合理。如果某個(gè) cluster 的 marker 基因列表混雜了多種已知細(xì)胞類型的特征可能說明聚類分辨率過高或數(shù)據(jù)中存在批次效應(yīng)需要進(jìn)一步調(diào)整參數(shù)或進(jìn)行批次校正。8. 總結(jié)與下一步學(xué)習(xí)建議到這里我們已經(jīng)完成了從 Jupyter Notebook 基礎(chǔ)操作到 Scanpy 單細(xì)胞分析入門全流程的梳理。現(xiàn)在你應(yīng)該能夠獨(dú)立完成以下工作安裝并啟動(dòng) Jupyter Notebook熟練使用 Cell、Markdown 和快捷鍵理解 Notebook 與 JupyterLab 的區(qū)別根據(jù)需求選擇合適工具了解 AnnData 的五大結(jié)構(gòu)字段跑通 pbmc3k 數(shù)據(jù)集的 QC、歸一化、高變基因篩選、PCA、鄰接圖、UMAP、Leiden 聚類和 marker 基因分析流程對 Windows 下 Notebook 空白頁、內(nèi)核連接失敗、依賴沖突等常見問題進(jìn)行定位和修復(fù)。下一步可以嘗試把這套流程應(yīng)用到自己的數(shù)據(jù)上。先從 10X 官網(wǎng)下載一個(gè)公開數(shù)據(jù)集練習(xí)read_10x_h5和read_10x_mtx兩種讀取方式之后再為分析 Notebook 補(bǔ)全細(xì)胞類型注釋嘗試使用sc.tl.diffmap做擴(kuò)散圖或者用sc.tl.score_genes給細(xì)胞打分熟悉更多 Scanpy 分析模塊。單細(xì)胞分析的學(xué)習(xí)曲線雖然比普通數(shù)據(jù)處理陡峭一些但結(jié)合 Jupyter Notebook 的交互式環(huán)境完全可以把整個(gè)路線逐步打通。如果在環(huán)境配置或代碼運(yùn)行中遇到報(bào)錯(cuò)優(yōu)先去官方文檔查對應(yīng)函數(shù)的參數(shù)說明社區(qū)的討論帖也值得多翻一翻。把這套 pbmc3k 流程真正吃透之后再上真實(shí)數(shù)據(jù)你會(huì)順手很多。