胞數(shù)據(jù)分析實(shí)戰(zhàn)指南)
1. 背景與核心概念1.1 Jupyter Notebook 是什么Jupyter Notebook 是一個(gè)基于 Web 的交互式開(kāi)發(fā)環(huán)境允許我們?cè)跒g覽器中編寫(xiě)代碼、運(yùn)行代碼、查看輸出結(jié)果同時(shí)還能在代碼之間插入 Markdown 說(shuō)明文字、公式、圖片和表格。最早它作為 IPython Notebook 項(xiàng)目誕生后來(lái)逐步發(fā)展為支持 Python、R、Julia 等多種語(yǔ)言的項(xiàng)目也因此得名 JupyterJulia Python R 的合稱(chēng)。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域Jupyter Notebook 幾乎是標(biāo)配工具。原因很簡(jiǎn)單單細(xì)胞數(shù)據(jù)的分析過(guò)程是一條長(zhǎng)鏈路包括數(shù)據(jù)讀取、質(zhì)量控制、歸一化、降維、聚類(lèi)、注釋等多個(gè)環(huán)節(jié)每一步之間高度依賴(lài)人工判斷。Notebook 這種“一段代碼 一段結(jié)果 一段說(shuō)明”的組織方式非常適合探索性數(shù)據(jù)分析和流程復(fù)現(xiàn)。你可以在同一個(gè)文件中既完成代碼調(diào)試又留下對(duì)參數(shù)選擇的理解記錄后續(xù)想要回溯分析思路也會(huì)輕松很多。1.2 Jupyter Notebook 與 JupyterLab 的區(qū)別很多初學(xué)者在安裝時(shí)會(huì)碰到一個(gè)搜索熱點(diǎn)Jupyter Notebook 和 JupyterLab 到底有什么區(qū)別簡(jiǎn)單來(lái)說(shuō)JupyterLab 是 Jupyter Notebook 的下一代交互式開(kāi)發(fā)界面相當(dāng)于把 Notebook 的能力放到了一個(gè)更像 IDE 的工作臺(tái)中。JupyterLab 支持多標(biāo)簽頁(yè)、拖拽布局、文件管理器、終端、代碼控制臺(tái)等功能在同一個(gè)窗口中可以同時(shí)打開(kāi)多個(gè) Notebook、查看數(shù)據(jù)文件、運(yùn)行終端命令。而經(jīng)典 Jupyter Notebook 則更聚焦于單個(gè) Notebook 文件的編輯和運(yùn)行界面更加簡(jiǎn)潔適合只專(zhuān)注分析任務(wù)的時(shí)候使用。兩者使用同一個(gè) Notebook 文件格式.ipynb運(yùn)行時(shí)也共享同一套內(nèi)核機(jī)制。也就是說(shuō)你在經(jīng)典 Notebook 里寫(xiě)的代碼在 JupyterLab 里也能正常打開(kāi)運(yùn)行。現(xiàn)在 Anaconda 默認(rèn)安裝的是 JupyterLab 和 Notebook 兩套組件啟動(dòng)命令分別對(duì)應(yīng)jupyter notebook jupyter lab如果只是在學(xué)習(xí)單細(xì)胞分析的基礎(chǔ)流程使用哪個(gè)都可以如果希望邊寫(xiě)代碼邊觀察文件目錄、邊調(diào)試邊查看圖表JupyterLab 的體驗(yàn)會(huì)更好。本文的代碼在兩種界面下均可直接運(yùn)行。1.3 Scanpy 是什么Scanpy 是一個(gè)基于 Python 的高性能單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析工具包全稱(chēng)是 Single-Cell Analysis in Python。它的核心數(shù)據(jù)結(jié)構(gòu)基于 AnnData能夠高效處理大規(guī)模稀疏矩陣支持從數(shù)據(jù)讀取、質(zhì)量控制、標(biāo)準(zhǔn)化、特征選擇、降維、聚類(lèi)到差異分析、可視化以及偽時(shí)間分析等完整流程。在單細(xì)胞數(shù)據(jù)分析領(lǐng)域R 語(yǔ)言生態(tài)中的 Seurat 一直是主流工具而 Scanpy 的出現(xiàn)讓 Python 用戶(hù)也能完成同等級(jí)的分析工作。相比 SeuratScanpy 有以下幾個(gè)特點(diǎn)構(gòu)建在 Anndata、NumPy、SciPy 和 Pandas 之上與 Python 數(shù)據(jù)科學(xué)生態(tài)天然融合內(nèi)存占用控制較好適合處理幾萬(wàn)甚至更多細(xì)胞的單細(xì)胞數(shù)據(jù)集分析流程高度模塊化每個(gè)步驟對(duì)應(yīng)一個(gè)明確的 API支持豐富的可視化圖形輸出可直接嵌入 Notebook 中展示。對(duì)于已經(jīng)熟悉 Python 的數(shù)據(jù)分析人員來(lái)說(shuō)Scanpy 是入門(mén)單細(xì)胞數(shù)據(jù)分析非常合適的切入點(diǎn)。我們?cè)诒疚闹袝?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)聚類(lèi)效果這一過(guò)程天然適合 Jupyter Notebook 這種交互式環(huán)境。用 Jupyter Notebook 跑 Scanpy每個(gè)分析步驟都能立即看到輸出和圖表比如質(zhì)量控制后的細(xì)胞數(shù)變化、UMAP 降維后的聚類(lèi)分布、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)備與版本說(shuō)明2.1 安裝 AnacondaScanpy 依賴(lài)的第三方庫(kù)較多手動(dòng)用 pip 逐個(gè)安裝容易遇到依賴(lài)沖突所以我們首先推薦通過(guò) Anaconda 來(lái)管理 Python 環(huán)境。Anaconda 自帶 Python 解釋器、Jupyter Notebook、常用科學(xué)計(jì)算庫(kù)以及包管理器 conda安裝完成后即可快速使用。到 Anaconda 官網(wǎng)下載對(duì)應(yīng)操作系統(tǒng)的安裝包安裝過(guò)程中保持默認(rèn)選項(xiàng)即可。安裝完成后打開(kāi)命令行窗口輸入conda --version如果能看到 conda 版本號(hào)說(shuō)明安裝成功。為了保證環(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)境可以避免庫(kù)與庫(kù)之間的版本沖突這是單細(xì)胞分析項(xiàng)目里非常重要的一步。2.2 安裝 Scanpy激活虛擬環(huán)境后使用 conda 安裝 Scanpyconda install -c conda-forge scanpy python-igraph leidenalg這里同時(shí)安裝了python-igraph和leidenalg它們用于 Leiden 聚類(lèi)算法。Leiden 是目前單細(xì)胞聚類(lèi)中較推薦的方法相比早期的 Louvain 算法社區(qū)劃分更穩(wěn)定對(duì)分辨率的調(diào)節(jié)也更友好。如果 conda 安裝速度慢也可以使用 pippip install scanpy pip install leidenalg不過(guò) pip 安裝時(shí)需要注意 NumPy、Pandas、SciPy 等底層庫(kù)的版本兼容性如果出現(xiàn)版本報(bào)錯(cuò)建議仍然回到 conda 安裝方式。安裝完成后在命令行輸入python -c import scanpy as sc; print(sc.__version__)能正常打印版本號(hào)就說(shuō)明環(huán)境已經(jīng)準(zhǔn)備好了。2.3 啟動(dòng) Jupyter Notebook在scanpy-env環(huán)境中啟動(dòng) Notebookjupyter notebook啟動(dòng)后終端會(huì)輸出類(lèi)似下面的信息[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)打開(kāi) Notebook 的主頁(yè)面。如果瀏覽器沒(méi)有自動(dòng)打開(kāi)可以把終端里顯示的http://localhost:8888/?tokenxxx地址手動(dòng)復(fù)制到瀏覽器地址欄。在 Notebook 主頁(yè)面的右上角點(diǎn)擊 New - Python 3即可新建一個(gè) Notebook 文件。建議把工作目錄切換到專(zhuān)門(mén)存放單細(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)行過(guò)程中自動(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í)鍵盤(pán)快捷鍵作用于整個(gè) Notebook可以刪除、復(fù)制、移動(dòng) Cell編輯模式按 Enter 進(jìn)入此時(shí)可以編輯當(dāng)前 Cell 中的代碼或文字。新建的 Cell 默認(rèn)是代碼類(lèi)型如果要寫(xiě)說(shuō)明文字需要把 Cell 切換為 Markdown 類(lèi)型。切換方式是在命令模式下按M切回代碼類(lèi)型則按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ù)說(shuō)明不用頻繁去查文檔。3.3 Markdown 與魔法命令Markdown Cell 中可以使用標(biāo)準(zhǔn)的 Markdown 語(yǔ)法包括#標(biāo)題、**加粗**、-列表、行內(nèi)代碼等。此外 Jupyter 還支持 LaTeX 數(shù)學(xué)公式用$$包裹即可。在單細(xì)胞分析筆記中我習(xí)慣在每個(gè)步驟前用 Markdown Cell 寫(xiě)清楚“這一步在做什么、為什么做”分析結(jié)束之后再回頭看整個(gè)文件就像一份可執(zhí)行的分析報(bào)告。魔法命令是 Jupyter 提供的一組擴(kuò)展指令以%開(kāi)頭。常用例子如下%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 如何在其他瀏覽器打開(kāi) Notebook有同學(xué)在 Windows 上會(huì)遇到默認(rèn)瀏覽器打不開(kāi) Notebook 的情況。常見(jiàn)的解決辦法有三種。第一種直接復(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 打開(kāi) c.NotebookApp.browser C:/Program Files/Google/Chrome/Application/chrome.exe修改完成后重啟 Notebook新瀏覽器配置即可生效。4. Scanpy 核心數(shù)據(jù)結(jié)構(gòu)與數(shù)據(jù)讀取4.1 AnnData 對(duì)象Scanpy 中所有分析操作都圍繞一個(gè)叫AnnData的對(duì)象展開(kāi)。AnnData 可以理解為“帶注釋的數(shù)據(jù)”它包含以下幾個(gè)核心部分adata.X表達(dá)量矩陣通常是稀疏矩陣或稠密矩陣行是細(xì)胞列是基因adata.obs細(xì)胞維度的元數(shù)據(jù)比如細(xì)胞 ID、樣本信息、批次信息、聚類(lèi)結(jié)果等adata.var基因維度的元數(shù)據(jù)比如基因名、基因類(lèi)型、是否高變基因等adata.uns非結(jié)構(gòu)化數(shù)據(jù)存放一些分析過(guò)程中的中間結(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é)在聚類(lèi)結(jié)束后找不到聚類(lèi)標(biāo)簽其實(shí)它就保存在adata.obs[leiden]中。4.2 數(shù)據(jù)讀取方式Scanpy 支持讀取多種單細(xì)胞數(shù)據(jù)格式常見(jiàn)的有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 對(duì)象。因?yàn)椴煌臄?shù)據(jù)來(lái)源列名和行名定義不一樣最穩(wěn)妥的方式是先了解自己數(shù)據(jù)的格式再選擇合適的讀取函數(shù)。4.3 常用質(zhì)量控制指標(biāo)單細(xì)胞數(shù)據(jù)在下游分析前必須經(jīng)過(guò)嚴(yán)格的質(zhì)量控制否則聚類(lèi)結(jié)果會(huì)被低質(zhì)量細(xì)胞帶偏。常用的 QC 指標(biāo)有三個(gè)每個(gè)細(xì)胞檢測(cè)到的基因數(shù)n_genes_by_counts每個(gè)細(xì)胞的總 UMI 數(shù)total_counts每個(gè)細(xì)胞的線(xiàn)粒體基因表達(dá)比例pct_counts_mt。線(xiàn)粒體基因比例高通常意味著細(xì)胞質(zhì) RNA 丟失嚴(yán)重細(xì)胞可能已經(jīng)破裂或?yàn)l臨死亡這類(lèi)細(xì)胞應(yīng)該過(guò)濾掉。在 Scanpy 中可以用以下方式添加線(xiàn)粒體基因比例指標(biāo)adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)對(duì)于人類(lèi)數(shù)據(jù)線(xiàn)粒體基因以MT-開(kāi)頭小鼠數(shù)據(jù)則以mt-開(kāi)頭具體前綴要根據(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)典的入門(mén)數(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)過(guò)初篩的 AnnData 對(duì)象包含 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ì)量控制與過(guò)濾先計(jì)算線(xiàn)粒體基因比例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ù)圖形中分布情況確定過(guò)濾閾值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過(guò)濾掉檢測(cè)到的基因數(shù)少于 200 的細(xì)胞這類(lèi)細(xì)胞可能是空液滴或質(zhì)量較差filter_genes過(guò)濾掉在超過(guò) 3 個(gè)細(xì)胞中表達(dá)的基因去掉在數(shù)據(jù)集中幾乎不表達(dá)的基因保留線(xiàn)粒體基因比例小于 5% 的細(xì)胞保留基因數(shù)少于 2500 的細(xì)胞避免雙細(xì)胞或異常高表達(dá)細(xì)胞干擾分析。執(zhí)行完畢后再次查看adata發(fā)現(xiàn)細(xì)胞數(shù)和基因數(shù)都減少了。5.3 歸一化、對(duì)數(shù)化與高變基因表達(dá)量的原始 counts 數(shù)據(jù)受測(cè)序深度影響很大因此需要做歸一化使得不同細(xì)胞之間可比。Scanpy 中的標(biāo)準(zhǔn)做法是先按每個(gè)細(xì)胞的總 counts 縮放再取對(duì)數(shù)sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)normalize_total把每個(gè)細(xì)胞的 counts 總數(shù)縮放到1e4log1p是對(duì)每個(gè)值計(jì)算log(1 x)。對(duì)數(shù)化之后表達(dá)量的分布更接近正態(tài)分布有利于后續(xù) PCA 等算法的計(jì)算。接著篩選高變基因。高變基因是那些在細(xì)胞之間表達(dá)差異顯著的基因它們攜帶了區(qū)分細(xì)胞類(lèi)型的主要信息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)行線(xiàn)性回歸去除 counts 和線(xiàn)粒體比例的影響然后進(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ù)因素帶來(lái)的干擾scale則是將每個(gè)基因的表達(dá)標(biāo)準(zhǔn)化到零均值、單位方差。pca_variance_ratio圖中曲線(xiàn)會(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 聚類(lèi)sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.8)聚類(lèi)結(jié)果會(huì)寫(xiě)入adata.obs[leiden]。繪制 UMAP 圖不同顏色即代表不同的聚類(lèi)群體sc.pl.umap(adata, color[leiden], legend_locon data)5.5 Marker 基因與細(xì)胞類(lèi)型注釋聚類(lèi)完成后我們通常需要找出每個(gè) cluster 的特異性高表達(dá)基因也就是 marker 基因并結(jié)合文獻(xiàn)知識(shí)判斷每個(gè) cluster 屬于什么細(xì)胞類(lèi)型。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把這些注釋信息寫(xiě)入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 與聚類(lèi) 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 寫(xiě)清楚分析說(shuō)明這樣最終的文件就是一份可復(fù)現(xiàn)、可分享的分析報(bào)告。6. 常見(jiàn)問(wèn)題與排查思路在使用 Jupyter Notebook 和 Scanpy 的過(guò)程中有一些高頻問(wèn)題值得單獨(dú)整理。問(wèn)題現(xiàn)象常見(jiàn)原因解決思路Windows 下 notebook 打開(kāi)后頁(yè)面空白瀏覽器兼容問(wèn)題或 Notebook 版本過(guò)舊清理瀏覽器緩存更換 Chrome/Edge升級(jí) notebook 版本瀏覽器沒(méi)有自動(dòng)彈出 Notebook 頁(yè)面瀏覽器配置或防火墻攔截手動(dòng)復(fù)制終端輸出的 URL 到瀏覽器打開(kāi)pip 安裝 scanpy 后 import 報(bào)錯(cuò)依賴(lài)庫(kù)版本沖突使用 conda 安裝或升級(jí) numpy/pandas/scipy運(yùn)行代碼時(shí)內(nèi)核持續(xù)顯示 Busy當(dāng)前 Cell 數(shù)據(jù)計(jì)算量過(guò)大檢查數(shù)據(jù)規(guī)模拆分 Cell必要時(shí)重啟內(nèi)核圖表無(wú)法顯示未啟用 inline 顯示執(zhí)行%matplotlib inline或設(shè)置sc.settings.autoshowTruefilter_cells過(guò)濾后細(xì)胞數(shù)過(guò)少閾值設(shè)置過(guò)于嚴(yán)格根據(jù) violin 圖重新觀察分布調(diào)整閾值Leiden 聚類(lèi)結(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 空白頁(yè)問(wèn)題的詳細(xì)排查在 Windows 上很多同學(xué)執(zhí)行jupyter notebook之后瀏覽器打開(kāi)的是空白頁(yè)面。這個(gè)問(wèn)題通常與瀏覽器緩存、安全軟件攔截或 Notebook 版本過(guò)舊有關(guān)。排查順序建議如下強(qiáng)制刷新頁(yè)面按住Ctrl F5清除當(dāng)前頁(yè)面緩存復(fù)制終端中完整的 URL包括token參數(shù)打開(kāi)新的無(wú)痕窗口訪(fǎng)問(wèn)關(guān)閉殺毒軟件或防火墻的網(wǎng)頁(yè)攔截功能局部端口8888放行升級(jí) Notebook 相關(guān)包pip install --upgrade notebook如果仍然無(wú)法解決可以嘗試在啟動(dòng)時(shí)指定 IP 和端口jupyter notebook --ip127.0.0.1 --port8889 --no-browser然后手動(dòng)訪(fǎng)問(wèn)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 安裝中的依賴(lài)問(wèn)題Scanpy 很依賴(lài)numpy、scipy、pandas、matplotlib、h5py等庫(kù)舊版本可能因?yàn)?API 改變而報(bào)錯(cuò)。安裝 Scanpy 時(shí)盡量直接使用 conda-forge 頻道conda 會(huì)自動(dòng)解析依賴(lài)版本。如果使用 pip 安裝后出現(xiàn)ImportError: cannot import name X的情況常見(jiàn)做法是把這些基礎(chǔ)庫(kù)一起升級(jí)到較新的兼容版本pip install --upgrade numpy pandas scipy h5py7. 最佳實(shí)踐與工程建議7.1 用虛擬環(huán)境隔離項(xiàng)目依賴(lài)單細(xì)胞分析項(xiàng)目周期長(zhǎng)、依賴(lài)多且不同項(xiàng)目之間對(duì) Scanpy、NumPy 等庫(kù)的版本要求可能不一樣。建議每個(gè)項(xiàng)目都使用獨(dú)立的 conda 虛擬環(huán)境并且把環(huán)境依賴(lài)導(dǎo)出保存conda env export environment.yml這樣無(wú)論是換電腦還是項(xiàng)目交接都能通過(guò)一份環(huán)境文件快速還原分析環(huán)境。7.2 Notebook 中記錄分析參數(shù)單細(xì)胞分析非常強(qiáng)調(diào)可復(fù)現(xiàn)性。每個(gè)關(guān)鍵步驟的參數(shù)例如 QC 閾值、聚類(lèi)分辨率、高變基因篩選條件都應(yīng)該在 Notebook 的 Markdown Cell 中記錄下來(lái)。建議每個(gè)步驟采用“目的 參數(shù) 結(jié)果說(shuō)明”的格式。例如## 質(zhì)量控制 - 目的過(guò)濾掉低質(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ù)到最終聚類(lèi)結(jié)果會(huì)經(jīng)過(guò)多個(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ù)想要重新畫(huà)圖不需要重新跑前面的全部計(jì)算。7.4 計(jì)算資源的合理使用pbmc3k 只有 2700 個(gè)細(xì)胞在本地電腦上運(yùn)行毫無(wú)壓力。但如果數(shù)據(jù)量達(dá)到十萬(wàn)甚至百萬(wàn)級(jí)細(xì)胞就需要考慮計(jì)算資源問(wèn)題。以下幾個(gè)方面可以?xún)?yōu)化過(guò)濾步驟盡量靠前減少后續(xù)分析的細(xì)胞數(shù)和基因數(shù)PCA 和 neighbor 計(jì)算時(shí)可指定n_pcs降低維度對(duì)超大數(shù)據(jù)集考慮使用 GPU 加速版本的 Scanpy 或分塊處理策略操作大矩陣時(shí)保存中間結(jié)果避免重復(fù)計(jì)算占用內(nèi)存。7.5 數(shù)據(jù)來(lái)源與生物意義的驗(yàn)證Scanpy 分析結(jié)果最終要回到實(shí)驗(yàn)和生物學(xué)意義中驗(yàn)證。每次聚類(lèi)完成后不能只看 UMAP 圖還應(yīng)該檢查 marker 基因的表達(dá)情況結(jié)合文獻(xiàn)確認(rèn)細(xì)胞類(lèi)型注釋是否合理。如果某個(gè) cluster 的 marker 基因列表混雜了多種已知細(xì)胞類(lèi)型的特征可能說(shuō)明聚類(lèi)分辨率過(guò)高或數(shù)據(jù)中存在批次效應(yīng)需要進(jìn)一步調(diào)整參數(shù)或進(jìn)行批次校正。8. 總結(jié)與下一步學(xué)習(xí)建議到這里我們已經(jīng)完成了從 Jupyter Notebook 基礎(chǔ)操作到 Scanpy 單細(xì)胞分析入門(mén)全流程的梳理。現(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 聚類(lèi)和 marker 基因分析流程對(duì) Windows 下 Notebook 空白頁(yè)、內(nèi)核連接失敗、依賴(lài)沖突等常見(jiàn)問(wèn)題進(jìn)行定位和修復(fù)。下一步可以嘗試把這套流程應(yīng)用到自己的數(shù)據(jù)上。先從 10X 官網(wǎng)下載一個(gè)公開(kāi)數(shù)據(jù)集練習(xí)read_10x_h5和read_10x_mtx兩種讀取方式之后再為分析 Notebook 補(bǔ)全細(xì)胞類(lèi)型注釋嘗試使用sc.tl.diffmap做擴(kuò)散圖或者用sc.tl.score_genes給細(xì)胞打分熟悉更多 Scanpy 分析模塊。單細(xì)胞分析的學(xué)習(xí)曲線(xiàn)雖然比普通數(shù)據(jù)處理陡峭一些但結(jié)合 Jupyter Notebook 的交互式環(huán)境完全可以把整個(gè)路線(xiàn)逐步打通。如果在環(huán)境配置或代碼運(yùn)行中遇到報(bào)錯(cuò)優(yōu)先去官方文檔查對(duì)應(yīng)函數(shù)的參數(shù)說(shuō)明社區(qū)的討論帖也值得多翻一翻。把這套 pbmc3k 流程真正吃透之后再上真實(shí)數(shù)據(jù)你會(huì)順手很多。