單細(xì)胞分析跑到后半程最容易出現(xiàn)兩種極端一種是聚類圖很漂亮但一問(wèn)你“這堆細(xì)胞到底是什么”答不上來(lái)另一種是差異表達(dá)跑完基因列表拉出來(lái)幾百行挑來(lái)挑去挑不出幾個(gè)能站得住腳的 marker。scanpy 的那套示例流程走到第三塊拼的其實(shí)不是代碼量而是對(duì)矩陣狀態(tài)、統(tǒng)計(jì)口徑和注釋邏輯的把控。這一篇接著往下講差異表達(dá)、細(xì)胞類型注釋和可視化這套組合拳順帶把軌跡推斷和細(xì)胞周期這些進(jìn)階動(dòng)作串起來(lái)。適合已經(jīng)跑通過(guò)質(zhì)控、歸一化、PCA、鄰居圖和 Leiden 聚類手里有一個(gè)帶leiden列的AnnData但面對(duì)rank_genes_groups的輸出不知道怎么下手的同學(xué)。全文的代碼都基于常見的公開外周血單核細(xì)胞示例數(shù)據(jù)集參數(shù)和坑點(diǎn)按我自己的實(shí)際操作習(xí)慣給。1. 走到這一步先別急著跑差異表達(dá)差異表達(dá)這步最冤的地方在于代碼只有一行跑完也能出結(jié)果但結(jié)果對(duì)不對(duì)沒人告訴你。我見過(guò)太多人聚類剛出圖順手就sc.tl.rank_genes_groups出來(lái)的 top 基因全是線粒體基因和核糖體基因然后開始懷疑人生。問(wèn)題通常不在統(tǒng)計(jì)方法而在喂進(jìn)去的矩陣根本不是你以為的那個(gè)矩陣。1.1 三種矩陣狀態(tài)決定你的結(jié)果可不可信單細(xì)胞流程里adata.X會(huì)經(jīng)歷至少三次變身原始 counts、歸一化加 log1p 之后的對(duì)數(shù)表達(dá)、以及sc.pp.scale之后的 z-score。這三種狀態(tài)在內(nèi)存里長(zhǎng)得一模一樣都是稀疏或稠密矩陣肉眼看不出來(lái)。差異表達(dá)要的是第二種也就是 log 歸一化后的數(shù)據(jù)但很多人做完 scale 直接往下跑X里已經(jīng)是 z-score 了。z-score 有正有負(fù)負(fù)值在做倍數(shù)變化的時(shí)候會(huì)被壓成很小甚至出錯(cuò)最后算出來(lái)的logfoldchanges完全沒有生物學(xué)含義。更麻煩的是這種錯(cuò)誤不會(huì)報(bào)錯(cuò)你只會(huì)看到一堆奇怪的正負(fù)號(hào)。所以我自己的習(xí)慣是在sc.pp.scale之前先把 log 歸一化后的數(shù)據(jù)存進(jìn)adata.raw后續(xù)所有差異表達(dá)和 marker 可視化都走use_rawTrue。# 歸一化與 log 轉(zhuǎn)換 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) # 關(guān)鍵一步把 log 數(shù)據(jù)鎖進(jìn) raw再去做 scale adata.raw adata sc.pp.scale(adata, max_value10)這么做的代價(jià)是內(nèi)存翻倍但換來(lái)的是后面所有分析都能隨時(shí)回到干凈的 log 矩陣。做完 scale 再 PCA、鄰居圖、UMAP、Leiden全部用 scale 后的X這部分沒問(wèn)題。但只要是算表達(dá)量差異、畫小提琴圖、畫點(diǎn)圖就統(tǒng)一use_rawTrue口徑就一致了。注意如果你的流程里壓根沒做 scale那adata.raw可能不存在這時(shí)候跑use_rawTrue會(huì)直接報(bào)錯(cuò)。判斷方法很簡(jiǎn)單print(adata.raw)看一眼是不是None。還有一個(gè)隱蔽的坑normalize_total之后再log1p和先log1p再歸一化結(jié)果不一樣。scanpy 官方流程是前者你照著做就行但如果你中間插了一步sc.pp.filter_genes把低表達(dá)基因過(guò)濾掉要注意過(guò)濾發(fā)生在歸一化之前還是之后。歸一化之后再過(guò)濾基因會(huì)導(dǎo)致每個(gè)細(xì)胞的總 counts 基數(shù)變化嚴(yán)格來(lái)說(shuō)需要重新歸一化。實(shí)際影響不大但如果你對(duì)結(jié)果特別較真就按“過(guò)濾基因 → 過(guò)濾細(xì)胞 → 歸一化 → log1p”的順序來(lái)。1.2 簇編號(hào)是隨機(jī)的別把0號(hào)簇寫死在代碼里L(fēng)eiden 聚類的輸出是一列0, 1, 2, ...的整數(shù)標(biāo)簽但這些數(shù)字沒有任何順序含義純粹是算法內(nèi)部給社區(qū)編的號(hào)。同一個(gè)數(shù)據(jù)集換個(gè)隨機(jī)種子、換個(gè) scanpy 版本甚至只是調(diào)了一下resolution簇的編號(hào)就會(huì)重排。我踩過(guò)一次很典型的坑腳本里寫死了“0號(hào)簇是 T 細(xì)胞、1號(hào)簇是單核細(xì)胞”結(jié)果換了臺(tái)機(jī)器重跑簇編號(hào)全亂注釋表直接對(duì)不上。正確的做法是把簇編號(hào)當(dāng)成一個(gè)臨時(shí)索引真正的身份靠 marker 基因判定。為了減少這種混亂我一般會(huì)在聚類之后、差異表達(dá)之前先跑一遍sc.tl.rank_genes_groups看看每個(gè)簇的 top 基因用一個(gè)字典把簇號(hào)和細(xì)胞類型對(duì)應(yīng)起來(lái)再統(tǒng)一寫回adata.obs[cell_type]。這樣即使中間重跑導(dǎo)致編號(hào)變化你只需要改字典不用改后面所有代碼。resolution的選擇也值得說(shuō)兩句。默認(rèn) 1.0 在幾千個(gè)細(xì)胞的 PBMC 數(shù)據(jù)上通常給出 8 到 12 個(gè)簇粒度差不多。如果你的目標(biāo)是找稀有亞群可以往上調(diào)到 1.5 甚至 2.0如果只是做粗分群0.5 就夠。我不建議一開始就上高分辨率因?yàn)榇匾欢嗖町惐磉_(dá)的比較次數(shù)會(huì)暴漲多重檢驗(yàn)校正之后的pvals_adj會(huì)變得很難看很多真實(shí)的差異會(huì)被壓下去。1.3 多批次數(shù)據(jù)的差異表達(dá)要怎么處理單樣本數(shù)據(jù)沒什么好糾結(jié)的直接按簇分組比較就行。但如果你是多個(gè)樣本合并分析的情況就復(fù)雜了。假設(shè)你做的是合并后的聚類leiden是所有細(xì)胞混在一起分的那么某個(gè)簇里可能同時(shí)包含來(lái)自多個(gè)樣本的細(xì)胞。這時(shí)候直接跑差異表達(dá)得到的基因里會(huì)混進(jìn)大量批次效應(yīng)帶來(lái)的差異而不是真正的細(xì)胞類型差異。我的處理順序是先做批次校正具體用什么方法看數(shù)據(jù)量樣本量少用 Harmony 這類基于嵌入的校正樣本量大可以考慮基于圖的整合在校正后的嵌入空間上做鄰居圖和 Leiden然后用校正前的 log 表達(dá)矩陣做差異表達(dá)。這個(gè)組合聽起來(lái)矛盾其實(shí)很合理聚類需要在去除批次效應(yīng)的空間里做否則同一類型的細(xì)胞會(huì)被批次拆開但差異表達(dá)要的是真實(shí)的表達(dá)量對(duì)比用校正后的矩陣反而會(huì)把生物學(xué)信號(hào)也一起抹掉。如果樣本量足夠更穩(wěn)妥的做法是在差異表達(dá)里顯式加入?yún)f(xié)變量或者干脆按樣本分開跑差異表達(dá)再取交集。rank_genes_groups本身沒有內(nèi)置的協(xié)變量參數(shù)需要在調(diào)用前把數(shù)據(jù)按樣本拆分。這個(gè)操作成本高一些但對(duì)于結(jié)論的可信度提升很明顯。2. 差異表達(dá)的參數(shù)選擇與結(jié)果解讀矩陣狀態(tài)確認(rèn)干凈、簇編號(hào)心里有數(shù)之后就可以正式跑了。這一行代碼sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon)誰(shuí)都寫得出來(lái)但method選什么、結(jié)果怎么篩、閾值卡多少這三個(gè)問(wèn)題決定了你后面注釋的成敗。2.1 四種method到底選哪個(gè)scanpy 提供了幾個(gè)統(tǒng)計(jì)檢驗(yàn)方法差別不小我把常用的幾個(gè)列成表method原理適用場(chǎng)景注意事項(xiàng)wilcoxon秩和檢驗(yàn)非參數(shù)默認(rèn)首選對(duì)表達(dá)分布不正態(tài)的數(shù)據(jù)穩(wěn)健細(xì)胞數(shù)少時(shí) p 值偏保守t-test參數(shù)檢驗(yàn)假設(shè)正態(tài)數(shù)據(jù)量大、分布接近正態(tài)對(duì)離群值敏感稀疏基因容易失真t-test_overestim_vart 檢驗(yàn)但方差不合并簇間細(xì)胞數(shù)差異懸殊時(shí)更敏感容易放出假陽(yáng)性需要更嚴(yán)的校正logreg邏輯回歸想看基因的判別能力而非單純表達(dá)差異結(jié)果解釋方式和前三個(gè)不同logfoldchanges含義要小心我自己的默認(rèn)選擇是wilcoxon。單細(xì)胞數(shù)據(jù)零膨脹嚴(yán)重表達(dá)分布根本不是正態(tài)的t 檢驗(yàn)的前提不成立用非參方法更穩(wěn)妥。t-test_overestim_var在簇大小極不均衡時(shí)確實(shí)更容易撈出差異基因比如一個(gè) 2000 細(xì)胞的大簇和一個(gè) 30 細(xì)胞的小簇比較但這個(gè)“更容易”里包含了不少假陽(yáng)性我一般只在探索性分析時(shí)用正式結(jié)論還是回到 wilcoxon。還有一個(gè)容易被忽略的參數(shù)是pts。默認(rèn)不輸出表達(dá)比例加上ptsTrue之后結(jié)果里會(huì)多一列pts表示該基因在這個(gè)簇里表達(dá)值大于 0的細(xì)胞比例以及pts_rest表示在其他簇里的比例。這兩個(gè)數(shù)字特別有用能幫你過(guò)濾掉那種“統(tǒng)計(jì)顯著但只有三五個(gè)細(xì)胞表達(dá)”的基因。sc.tl.rank_genes_groups( adata, groupbyleiden, methodwilcoxon, use_rawTrue, ptsTrue )2.2 結(jié)果提取與閾值到底卡多少跑完之后結(jié)果存在adata.uns[rank_genes_groups]里是個(gè)結(jié)構(gòu)化的嵌套字典直接看很費(fèi)勁。用sc.get.rank_genes_groups_df轉(zhuǎn)成 DataFrame 才是正常人的操作。df sc.get.rank_genes_groups_df(adata, group0)得到的表里有這么幾列names、scores、logfoldchanges、pvals、pvals_adj、pct_nz_group、pct_nz_reference開了ptsTrue才有后兩列。這里有個(gè)必須說(shuō)清楚的細(xì)節(jié)logfoldchanges的底數(shù)是 2也就是說(shuō)logfoldchanges1代表兩倍差異2代表四倍。很多人誤以為是自然對(duì)數(shù)把閾值卡在 0.25 就以為是 1.28 倍實(shí)際上那是 1.19 倍尺度差了一截。另外 scanpy 在計(jì)算倍數(shù)變化時(shí)如果檢測(cè)到你的數(shù)據(jù)是 log1p 過(guò)的內(nèi)部會(huì)把值還原到線性尺度再算比值所以你不用自己手動(dòng)expm1但也正因?yàn)槿绱藆se_raw一定要指向正確的 log 矩陣否則還原出來(lái)的東西就離譜了。閾值怎么卡我給一個(gè)實(shí)操建議pvals_adj 0.05這是底線多重檢驗(yàn)校正后的顯著性logfoldchanges 0.5約等于 1.4 倍差異能過(guò)濾掉大量微小但顯著的變化pct_nz_group 0.25這個(gè)簇里至少 25% 的細(xì)胞表達(dá)該基因避免被極少數(shù)細(xì)胞帶偏三個(gè)條件同時(shí)卡通常每個(gè)簇能留下二三十個(gè)比較靠譜的候選基因。如果你只是想做初步注釋其實(shí)放寬到pvals_adj 0.05加logfoldchanges 0.25也就夠看了反正最后還要靠 marker 知識(shí)庫(kù)交叉驗(yàn)證。提示scores這一列是檢驗(yàn)統(tǒng)計(jì)量的標(biāo)準(zhǔn)化值可以粗略理解為“差異的強(qiáng)度”但不同 method 之間不可比。用它排序沒問(wèn)題別拿它當(dāng)效應(yīng)量。2.3 一次拿到一整套marker基因表格逐個(gè)簇提 DataFrame 太累我一般寫個(gè)循環(huán)把所有簇的 top 基因匯總成一張表順便把過(guò)濾條件也帶上方便后面查。import pandas as pd results {} for group in adata.obs[leiden].cat.categories: df sc.get.rank_genes_groups_df(adata, groupgroup) df df[(df[pvals_adj] 0.05) (df[logfoldchanges] 0.5)] df[cluster] group results[group] df.head(30) marker_table pd.concat(results.values(), ignore_indexTrue)這張表后面的用法很多可以直接拿 top 基因去搜可以畫熱圖看整體格局也可以作為注釋時(shí)的參考清單。我習(xí)慣再補(bǔ)一步把所有簇的 top 基因取并集去重之后作為后續(xù)點(diǎn)圖的var_names這樣一張圖就能看到全局的 marker 分布比逐簇看 violin 效率高得多。這里再補(bǔ)充一個(gè)經(jīng)驗(yàn)不要只看 top 10。有些細(xì)胞類型的特征基因在差異表達(dá)里排名并不靠前因?yàn)樗谶@個(gè)簇里表達(dá)量不高但特異性很強(qiáng)被那些高表達(dá)但有噪聲的基因擠下去了。我一般會(huì)把每個(gè)簇的前 50 個(gè)基因都導(dǎo)出來(lái)掃一眼找那些“知道是 marker 但沒進(jìn)前十”的基因往往能糾正一些誤判。3. 細(xì)胞類型注釋從看marker到寫標(biāo)簽差異表達(dá)給你的是“哪些基因不一樣”注釋解決的是“這些不一樣的基因說(shuō)明它是什么細(xì)胞”。這兩步之間沒有自動(dòng)化的橋梁必須靠 marker 知識(shí)加人工判斷。這是整個(gè)流程里最像“手藝活”的部分也是最容易糊弄過(guò)去的部分。3.1 一份能直接抄的PBMC marker清單外周血單核細(xì)胞是最經(jīng)典的示例數(shù)據(jù)細(xì)胞類型相對(duì)標(biāo)準(zhǔn)marker 也比較明確。我整理了一份常用的對(duì)照表遇到新數(shù)據(jù)可以先拿這套去套套不上的再單獨(dú)找細(xì)胞類型核心marker輔助marker備注CD4 T 細(xì)胞CD3D, CD3E, IL7RLTB, MALAT1與 CD8 共享 CD3 系列CD8 T 細(xì)胞CD3D, CD8A, CD8BGZMK, CCL5需與 NK 區(qū)分NK 細(xì)胞NKG7, GNLY, KLRD1NCAM1, PRF1不表達(dá) CD3B 細(xì)胞CD79A, MS4A1, CD19CD37, BANK1MS4A1 即 CD20經(jīng)典單核細(xì)胞LYZ, S100A8, S100A9CD14, VCAN占比通常最高非經(jīng)典單核細(xì)胞FCGR3A, MS4A7CDKN1C, LST1FCGR3A 即 CD16樹突狀細(xì)胞FCER1A, CST3CLEC10A, CD1C數(shù)量少容易被并入單核漿細(xì)胞樣樹突IL3RA, CLEC4CGZMB, TCF4IL3RA 即 CD123血小板PPBP, GP9, PF4ITGA2B常被視為雜質(zhì)紅系前體HBB, HBA1, HBA2—高表達(dá)時(shí)注意是否污染用這張表的時(shí)候有個(gè)細(xì)節(jié)要注意marker 是組合看的不是單基因定生死。比如 LYZ 高表達(dá)可能是單核細(xì)胞也可能是中性粒細(xì)胞或者某些 DC 亞群CD3D 陽(yáng)性但 CD8A 陰性大概率是 CD4 T。我一般會(huì)先看每個(gè)簇的 top 基因里有沒有兩個(gè)以上的該類型 marker再?zèng)Q定標(biāo)簽。另外血小板的 PPBP 和紅系的 HBB 如果出現(xiàn)在某個(gè)“看起來(lái)像 T 細(xì)胞”的簇里很可能是雙細(xì)胞或者環(huán)境 RNA 污染。處理方式取決于比例如果只是零星幾個(gè)細(xì)胞直接忽略如果整個(gè)簇都這樣那這個(gè)簇可能就是個(gè)偽簇考慮在聚類前就把這些基因高的細(xì)胞濾掉。3.2 打分法與閾值判定光看 top 基因列表有時(shí)候不夠直觀尤其是當(dāng)一個(gè)簇的 marker 信號(hào)比較弱、top 基因里混著一堆核糖體基因的時(shí)候。這時(shí)候sc.tl.score_genes就派上用場(chǎng)了它能給每個(gè)細(xì)胞算一個(gè)基因集的綜合得分然后你按簇看這個(gè)得分的分布。sc.tl.score_genes( adata, gene_list[CD3D, CD3E, CD2, IL7R], score_nameT_cell_score, ctrl_size50, n_bins25, use_rawTrue )幾個(gè)參數(shù)的含義值得說(shuō)清楚。ctrl_size是從表達(dá)量分箱里隨機(jī)抽多少個(gè)對(duì)照基因默認(rèn) 50調(diào)大一點(diǎn)結(jié)果更穩(wěn)但計(jì)算慢n_bins是按平均表達(dá)量把基因分成多少個(gè)箱子默認(rèn) 25主要目的是讓對(duì)照基因的表達(dá)水平和你目標(biāo)基因集匹配避免因?yàn)槟繕?biāo)基因本身表達(dá)高就天然得分高。use_rawTrue這里同樣重要用 scale 后的數(shù)據(jù)算出來(lái)的得分沒法解釋。得分算完怎么看我一般用sc.pl.violin按簇畫出來(lái)看哪個(gè)簇的得分分布整體上移。判定標(biāo)準(zhǔn)不是“得分大于 0 就是陽(yáng)性”而是看分布的形狀如果某個(gè)簇的得分分布明顯雙峰說(shuō)明這個(gè)簇內(nèi)部可能混了兩類細(xì)胞值得提高分辨率再分一次如果整個(gè)簇的分布都很平那說(shuō)明這套 marker 對(duì)這個(gè)簇沒有區(qū)分度。打分法還有一個(gè)經(jīng)典用途是細(xì)胞周期評(píng)分這個(gè)放到后面第 5 節(jié)講。3.3 自動(dòng)注釋與人工復(fù)核的配合現(xiàn)在基于參考數(shù)據(jù)集的自動(dòng)注釋工具不少思路基本一致拿一個(gè)已經(jīng)注釋好的參考數(shù)據(jù)集把你的數(shù)據(jù)映射過(guò)去用相似度或者分類器給出候選標(biāo)簽。這類方法的好處是快壞處是它一定會(huì)給你一個(gè)答案哪怕你的數(shù)據(jù)里存在參考集里沒有的細(xì)胞類型。我的用法是自動(dòng)注釋出初稿人工復(fù)核定終稿。具體操作上先用自動(dòng)方法給每個(gè)簇一個(gè)候選標(biāo)簽然后拿這個(gè)標(biāo)簽去反查差異表達(dá)結(jié)果如果自動(dòng)注釋說(shuō)是 B 細(xì)胞那這個(gè)簇的 top 基因里應(yīng)該出現(xiàn) CD79A、MS4A1 這些對(duì)得上就通過(guò)對(duì)不上就手動(dòng)改。這個(gè)交叉驗(yàn)證的步驟花不了多少時(shí)間但能擋掉大部分明顯的錯(cuò)誤。還有一種情況是自動(dòng)注釋給出的標(biāo)簽特別細(xì)比如把 T 細(xì)胞分成了七八個(gè)亞群但你的數(shù)據(jù)里每個(gè)亞群只有幾十個(gè)細(xì)胞。這種細(xì)粒度標(biāo)簽在統(tǒng)計(jì)上站不住腳我一般會(huì)往上合并一層統(tǒng)一叫 T 細(xì)胞然后在討論里注明“未進(jìn)一步細(xì)分亞群”。別為了顯得分析深入就硬保留細(xì)分標(biāo)簽細(xì)胞數(shù)不夠的亞群差異表達(dá)結(jié)果基本是噪聲。3.4 把注釋寫回obs并保持可追溯注釋結(jié)果最終要落到adata.obs里通常是新增一列cell_type。這里我強(qiáng)烈建議同時(shí)保留原始的leiden列不要覆蓋。后續(xù)做任何分組可視化、差異表達(dá)你都需要能在“聚類編號(hào)”和“細(xì)胞類型”之間來(lái)回切換。比如你想看某個(gè)細(xì)胞類型內(nèi)部是否還有異質(zhì)性就得回到 leiden 層面看。cluster2celltype { 0: CD4 T, 1: CD14 Monocyte, 2: B cell, 3: CD8 T, 4: NK cell, 5: FCGR3A Monocyte, 6: Dendritic cell, 7: Platelet } adata.obs[cell_type] adata.obs[leiden].map(cluster2celltype).astype(category)用map而不用循環(huán)賦值好處是映射關(guān)系一目了然而且沒映射到的簇會(huì)自動(dòng)變成缺失值你能立刻發(fā)現(xiàn)漏了哪個(gè)。映射完記得看一眼adata.obs[cell_type].value_counts()確認(rèn)沒有大量 NaN。另外把這份字典單獨(dú)存成一個(gè)文件比如 json 或者 csv跟分析腳本放一起下次別人接手或者你自己三個(gè)月后回看能立刻搞清楚當(dāng)時(shí)的判斷依據(jù)。4. 可視化細(xì)節(jié)從能看到能放進(jìn)匯報(bào)圖和圖之間的差距在單細(xì)胞分析里體現(xiàn)得特別明顯。同一份數(shù)據(jù)有人畫出來(lái)的點(diǎn)圖一眼就能看出細(xì)胞類型分布有人畫出來(lái)糊成一團(tuán)什么都看不清。這一節(jié)專講參數(shù)。4.1 dotplot的關(guān)鍵參數(shù)點(diǎn)圖是 marker 展示的主力橫軸基因、縱軸簇點(diǎn)的大小表示表達(dá)比例顏色表示平均表達(dá)量。核心參數(shù)就那么幾個(gè)但每個(gè)都能顯著改變可讀性。sc.pl.dotplot( adata, var_namesmarker_dict, groupbycell_type, standard_scalevar, cmapRdBu_r, dot_max0.8, dot_min0.1, figsize(12, 5), use_rawTrue )standard_scalevar是按基因做標(biāo)準(zhǔn)化把每個(gè)基因的表達(dá)量縮放到 0 到 1 之間。這個(gè)參數(shù)非常關(guān)鍵不同基因的絕對(duì)表達(dá)量差異巨大CD3D 可能均值 3PPBP 可能均值 8不標(biāo)準(zhǔn)化的話顏色分布完全被高表達(dá)基因主導(dǎo)低表達(dá)但特異的基因就看不出來(lái)了。按基因標(biāo)準(zhǔn)化之后每個(gè)基因的顏色都是相對(duì)自己的尺度特異性一目了然。dot_max和dot_min控制點(diǎn)大小的映射范圍。默認(rèn)dot_max是 1.0也就是表達(dá)比例最高的那個(gè)點(diǎn)會(huì)撐滿如果你的數(shù)據(jù)里有個(gè)基因在所有簇里都接近 100% 表達(dá)那其他點(diǎn)就全被壓得很小。把dot_max設(shè)成 0.8 或者 0.6能讓中等表達(dá)的點(diǎn)的差異更明顯。cmapRdBu_r是紅藍(lán)配色配合standard_scalevar的時(shí)候紅色代表高表達(dá)、藍(lán)色代表低表達(dá)中間是白色。我不建議用默認(rèn)的 viridis因?yàn)樗诒硎尽坝姓胸?fù)的相對(duì)高低”時(shí)不如發(fā)散型配色直觀。當(dāng)然如果你要投的期刊有配色要求另說(shuō)。還有一個(gè)實(shí)用技巧基因分組。var_names可以直接傳一個(gè)字典鍵是細(xì)胞類型值是對(duì)應(yīng)的基因列表這樣點(diǎn)圖會(huì)自動(dòng)按組排列基因并加上分組標(biāo)簽比自己手動(dòng)排順序省事得多。marker_dict { T cell: [CD3D, CD3E, IL7R, CD8A], NK: [NKG7, GNLY, KLRD1], B cell: [CD79A, MS4A1, CD19], Monocyte: [LYZ, S100A8, CD14, FCGR3A], DC: [FCER1A, CST3], Platelet: [PPBP, GP9] }4.2 stacked violin / matrixplot / heatmap 的取舍這三種圖都能展示多基因多組的表達(dá)格局但側(cè)重不同選錯(cuò)了會(huì)事倍功半。sc.pl.stacked_violin展示的是每個(gè)基因在每個(gè)組里的完整表達(dá)分布能看出雙峰、零膨脹這些結(jié)構(gòu)信息。缺點(diǎn)是基因一多圖會(huì)變得極長(zhǎng)。我一般只在需要展示兩三個(gè)關(guān)鍵基因的分布細(xì)節(jié)時(shí)用它比如確認(rèn)某個(gè) marker 是不是真的在目標(biāo)簇里特異性高表達(dá)還是所有簇都有一小撮陽(yáng)性細(xì)胞。sc.pl.matrixplot是點(diǎn)圖的簡(jiǎn)化版只保留顏色不保留點(diǎn)的大小適合臨床或者跨學(xué)科讀者信息密度低但易懂。如果匯報(bào)對(duì)象不熟悉單細(xì)胞用這個(gè)比點(diǎn)圖友好。sc.pl.rank_genes_groups_heatmap是直接基于差異表達(dá)結(jié)果畫熱圖每個(gè)簇取 top N 個(gè)基因。這張圖的信息量最大能同時(shí)看到多個(gè)簇的 marker 格局但基因太多時(shí)行標(biāo)簽會(huì)糊成一片。我一般取每個(gè)簇前 5 到 8 個(gè)基因n_genes5配合swap_axes調(diào)整方向出來(lái)的圖比較清爽。sc.pl.rank_genes_groups_heatmap( adata, n_genes5, groupbycell_type, show_gene_labelsTrue, use_rawTrue, cmapRdBu_r )提示熱圖默認(rèn)會(huì)對(duì)表達(dá)量做標(biāo)準(zhǔn)化如果圖上的顏色和你預(yù)期的表達(dá)高低對(duì)不上先檢查standard_scale參數(shù)和use_raw設(shè)置。這兩個(gè)是熱圖翻車的高頻原因。4.3 出圖保存與配色sc.pl系列函數(shù)都支持save參數(shù)傳個(gè)文件名就會(huì)自動(dòng)存到figures/目錄下格式默認(rèn)是 png可以通過(guò)sc.settings.figdir改路徑、sc.settings.file_format_figs改格式。我一般設(shè)成 pdf矢量圖放進(jìn)匯報(bào)里放大不失真。配色上連續(xù)變量用發(fā)散型或者順序型配色分類變量比如 UMAP 上的細(xì)胞類型多用定性配色。scanpy 默認(rèn)的分類配色在類別超過(guò) 10 種之后會(huì)開始重復(fù)或者相近容易混淆。我習(xí)慣手動(dòng)指定把主要的細(xì)胞類型分配差異大的顏色稀有的類型用灰色系這樣 UMAP 上主次分明。sc.settings.figdir ./figures sc.settings.file_format_figs pdf sc.pl.umap( adata, colorcell_type, palette[#4E79A7, #F28E2B, #E15759, #76B7B2, #59A14F, #EDC948, #B07AA1, #9C755F], frameonFalse, legend_locon data, legend_fontsize10, save_celltype.pdf )frameonFalse去掉邊框legend_locon data把標(biāo)簽直接標(biāo)在簇中心比側(cè)邊圖例清晰。但如果你的簇有重疊on data會(huì)擠在一起那就退回默認(rèn)的right margin。5. 往前一步細(xì)胞周期、軌跡與通路差異表達(dá)和注釋做完一個(gè)完整的描述性分析就閉環(huán)了。但如果想再往前走一步看看細(xì)胞狀態(tài)和分化關(guān)系還有幾件事值得做。5.1 細(xì)胞周期回歸增殖活躍的細(xì)胞在轉(zhuǎn)錄組上會(huì)和同類型的靜息細(xì)胞拉開很大距離聚類的時(shí)候經(jīng)常出現(xiàn)一個(gè)“細(xì)胞周期簇”把本該在一起的細(xì)胞拆開。處理方式是給每個(gè)細(xì)胞算 S 期和 G2M 期得分然后在降維之前把這兩個(gè)得分的影響回歸掉。# 使用公開的 S 期與 G2M 期基因列表 sc.tl.score_genes_cell_cycle(adata, s_geness_genes, g2m_genesg2m_genes) sc.pp.regress_out(adata, [S_score, G2M_score]) sc.pp.scale(adata, max_value10)注意順序score_genes_cell_cycle必須在scale之前調(diào)用而且要在 log 歸一化之后。regress_out這一步在細(xì)胞數(shù)上萬(wàn)的時(shí)候會(huì)很慢因?yàn)樗侵鸹蜃鼍€性回歸。如果數(shù)據(jù)量大可以考慮用sc.pp.scale之前先做sc.pp.highly_variable_genes只對(duì)高變基因做回歸能省不少時(shí)間。還有一個(gè)反直覺的點(diǎn)細(xì)胞周期回歸不是必須的。如果你的研究對(duì)象本身就是增殖過(guò)程比如發(fā)育或者腫瘤那細(xì)胞周期差異就是生物學(xué)信號(hào)不該回歸掉。只有在研究細(xì)胞類型構(gòu)成、且確認(rèn)細(xì)胞周期主導(dǎo)了聚類結(jié)果的時(shí)候才需要這一步。5.2 PAGA與DPT想看細(xì)胞類型之間的連接關(guān)系PAGA 是個(gè)好選擇。它把每個(gè)簇壓縮成一個(gè)節(jié)點(diǎn)用簇間的連接強(qiáng)度作為邊的權(quán)重能直觀看出哪些細(xì)胞類型在轉(zhuǎn)錄組上更接近。sc.tl.paga(adata, groupscell_type) sc.pl.paga(adata, plotFalse) sc.tl.umap(adata, init_pospaga) sc.pl.paga_compare(adata, threshold0.03)sc.pl.paga_compare會(huì)同時(shí)畫出 PAGA 圖和 UMAP兩張圖并排看能驗(yàn)證 UMAP 上的鄰近關(guān)系是否和 PAGA 的連接強(qiáng)度一致。如果不一致往往說(shuō)明 UMAP 的局部結(jié)構(gòu)有失真這時(shí)候 PAGA 的結(jié)果更可信。DPT 是沿著 PAGA 的圖結(jié)構(gòu)算偽時(shí)間的需要一個(gè)明確的起點(diǎn)細(xì)胞。這個(gè)起點(diǎn)怎么選取決于你的生物學(xué)問(wèn)題通常選一個(gè)處于分化起始狀態(tài)的細(xì)胞類型里的代表性細(xì)胞。adata.uns[iroot] 0 # 替換成你選定的起始細(xì)胞的索引 sc.tl.dpt(adata)DPT 的結(jié)果對(duì)起點(diǎn)極其敏感換個(gè)起點(diǎn)偽時(shí)間軸就完全反過(guò)來(lái)。所以我一般會(huì)試幾個(gè)候選起點(diǎn)看哪個(gè)得到的偽時(shí)間排序符合已知的生物學(xué)常識(shí)再定下來(lái)。5.3 基因集打分看通路除了細(xì)胞周期score_genes可以用來(lái)評(píng)估任意基因集的活性比如炎癥相關(guān)、氧化磷酸化、細(xì)胞毒性等等。做法和前面一樣關(guān)鍵是基因集來(lái)源要可靠別自己隨手拼。打分結(jié)果可以按細(xì)胞類型畫成熱圖看哪類細(xì)胞在哪個(gè)通路上更活躍。這一步的價(jià)值在于它能把差異表達(dá)得到的“一堆基因”提升到“功能層面”讓結(jié)論更容易解釋。比如你發(fā)現(xiàn)某個(gè)單核細(xì)胞亞群在炎癥基因集上得分顯著更高這比列出一百個(gè)上調(diào)基因有說(shuō)服力得多。需要注意的是通路打分和差異表達(dá)是兩套邏輯前者是基因集整體水平后者是單基因水平結(jié)論有可能不一致。不一致的時(shí)候優(yōu)先相信差異表達(dá)的單基因結(jié)果因?yàn)樗苯油返梅质芑蚣|(zhì)量和大小影響很大適合作為輔助證據(jù)。6. 踩坑實(shí)錄與速查表寫了一堆方法最后還是落到實(shí)際會(huì)遇到什么問(wèn)題。下面這張表是我自己和我?guī)н^(guò)的同學(xué)最常撞的幾類按出現(xiàn)頻率排序。6.1 典型問(wèn)題排查表現(xiàn)象可能原因排查方法處理方式所有簇的 top 基因都是 MT- 或 RPS 開頭質(zhì)控沒做好或差異表達(dá)用了未過(guò)濾數(shù)據(jù)檢查adata.var_names里線粒體基因比例回到質(zhì)控步驟重新過(guò)濾logfoldchanges 出現(xiàn)大量負(fù)值或異常大X是 scale 后的 z-scoreprint(adata.X.max())看是否超過(guò) 10改用use_rawTrue某個(gè)簇沒有任何顯著差異基因簇太小或該簇是雙細(xì)胞/偽簇看該簇細(xì)胞數(shù)和 marker 分布合并到鄰近簇或?yàn)V除點(diǎn)圖上所有點(diǎn)顏色接近沒做standard_scale檢查參數(shù)加standard_scalevarUMAP 上同類細(xì)胞被拆成幾塊批次效應(yīng)未校正按樣本著色看分布做批次校正后重聚類use_rawTrue報(bào)錯(cuò)adata.raw為 Noneprint(adata.raw)在 scale 前補(bǔ)上adata.raw adataleiden 報(bào)錯(cuò)找不到模塊缺少 igraph / leidenalgpip list檢查安裝對(duì)應(yīng)依賴或改用flavorigraph6.2 幾個(gè)容易忽略的實(shí)操細(xì)節(jié)第一個(gè)是隨機(jī)種子。Leiden 聚類、UMAP、t-SNE 都有隨機(jī)性不設(shè)種子的話每次重跑圖都不一樣。在腳本開頭加一句sc.settings.seed 0能保證結(jié)果可復(fù)現(xiàn)。我做過(guò)的項(xiàng)目里有一半以上的“圖怎么變了”問(wèn)題都是沒設(shè)種子導(dǎo)致的。第二個(gè)是數(shù)據(jù)類型的坑。adata.obs里的分組列必須是 categorical 或者至少是字符串如果是整數(shù)類型rank_genes_groups在某些版本里會(huì)把它們當(dāng)連續(xù)變量處理結(jié)果完全不對(duì)。用astype(category)顯式轉(zhuǎn)換一下能省掉很多莫名其妙的報(bào)錯(cuò)。第三個(gè)是基因名的大小寫和別名。同一個(gè)基因在不同注釋版本里可能叫CD3D也可能叫CD3dmarker 列表直接拷過(guò)來(lái)匹配不上是常事。判斷方法是用adata.var_names.isin(gene_list).sum()看能匹配上幾個(gè)如果明顯少于預(yù)期就得檢查命名規(guī)范必要時(shí)做大寫轉(zhuǎn)換或者別名映射。第四個(gè)是保存中間結(jié)果。單細(xì)胞流程跑一遍動(dòng)輒十幾分鐘做完注釋之后的AnnData一定要adata.write(annotated.h5ad)存下來(lái)。我見過(guò)太多人注釋完直接關(guān)了 notebook第二天想畫個(gè)圖發(fā)現(xiàn)得從頭跑。存的時(shí)候注意write默認(rèn)會(huì)保留raw文件會(huì)比較大可以用adata.write(x.h5ad, compressiongzip)壓縮一下。最后一個(gè)也是我覺得最重要的注釋結(jié)論要留證據(jù)。不要只留一個(gè)cell_type列把當(dāng)時(shí)的點(diǎn)圖、差異表達(dá)表格、映射字典都存到一個(gè)文件夾里。半年后你或者別人回頭看能立刻知道每個(gè)標(biāo)簽是怎么來(lái)的。這件事在單獨(dú)分析的時(shí)候感覺多余但在需要復(fù)現(xiàn)或者答辯的時(shí)候能救命。我個(gè)人在實(shí)際操作中的體會(huì)是單細(xì)胞分析真正花時(shí)間的從來(lái)不是寫代碼而是反復(fù)在“這個(gè)簇到底是什么”上糾結(jié)。差異表達(dá)給你的是線索不是答案marker 知識(shí)庫(kù)給你的是參照也不是答案。真正的答案往往來(lái)自多個(gè)證據(jù)的交叉top 基因?qū)Φ蒙稀⒋蚍址植贾С?、自?dòng)注釋也這么認(rèn)為、文獻(xiàn)里這個(gè)組織確實(shí)有這類細(xì)胞。四條里對(duì)上三條就可以定下來(lái)了剩下的那一條對(duì)不上記下來(lái)別硬圓。