资讯动态

Scanpy单细胞差异表达、注释与可视化实战指南

发布时间:2026/10/28 3:31:48 来源:尧图企业网站定制
单细胞分析跑到后半程最容易出现两种极端一种是聚类图很漂亮但一问你“这堆细胞到底是什么”答不上来另一种是差异表达跑完基因列表拉出来几百行挑来挑去挑不出几个能站得住脚的 marker。scanpy 的那套示例流程走到第三块拼的其实不是代码量而是对矩阵状态、统计口径和注释逻辑的把控。这一篇接着往下讲差异表达、细胞类型注释和可视化这套组合拳顺带把轨迹推断和细胞周期这些进阶动作串起来。适合已经跑通过质控、归一化、PCA、邻居图和 Leiden 聚类手里有一个带leiden列的AnnData但面对rank_genes_groups的输出不知道怎么下手的同学。全文的代码都基于常见的公开外周血单核细胞示例数据集参数和坑点按我自己的实际操作习惯给。1. 走到这一步先别急着跑差异表达差异表达这步最冤的地方在于代码只有一行跑完也能出结果但结果对不对没人告诉你。我见过太多人聚类刚出图顺手就sc.tl.rank_genes_groups出来的 top 基因全是线粒体基因和核糖体基因然后开始怀疑人生。问题通常不在统计方法而在喂进去的矩阵根本不是你以为的那个矩阵。1.1 三种矩阵状态决定你的结果可不可信单细胞流程里adata.X会经历至少三次变身原始 counts、归一化加 log1p 之后的对数表达、以及sc.pp.scale之后的 z-score。这三种状态在内存里长得一模一样都是稀疏或稠密矩阵肉眼看不出来。差异表达要的是第二种也就是 log 归一化后的数据但很多人做完 scale 直接往下跑X里已经是 z-score 了。z-score 有正有负负值在做倍数变化的时候会被压成很小甚至出错最后算出来的logfoldchanges完全没有生物学含义。更麻烦的是这种错误不会报错你只会看到一堆奇怪的正负号。所以我自己的习惯是在sc.pp.scale之前先把 log 归一化后的数据存进adata.raw后续所有差异表达和 marker 可视化都走use_rawTrue。# 归一化与 log 转换 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) # 关键一步把 log 数据锁进 raw再去做 scale adata.raw adata sc.pp.scale(adata, max_value10)这么做的代价是内存翻倍但换来的是后面所有分析都能随时回到干净的 log 矩阵。做完 scale 再 PCA、邻居图、UMAP、Leiden全部用 scale 后的X这部分没问题。但只要是算表达量差异、画小提琴图、画点图就统一use_rawTrue口径就一致了。注意如果你的流程里压根没做 scale那adata.raw可能不存在这时候跑use_rawTrue会直接报错。判断方法很简单print(adata.raw)看一眼是不是None。还有一个隐蔽的坑normalize_total之后再log1p和先log1p再归一化结果不一样。scanpy 官方流程是前者你照着做就行但如果你中间插了一步sc.pp.filter_genes把低表达基因过滤掉要注意过滤发生在归一化之前还是之后。归一化之后再过滤基因会导致每个细胞的总 counts 基数变化严格来说需要重新归一化。实际影响不大但如果你对结果特别较真就按“过滤基因 → 过滤细胞 → 归一化 → log1p”的顺序来。1.2 簇编号是随机的别把0号簇写死在代码里Leiden 聚类的输出是一列0, 1, 2, ...的整数标签但这些数字没有任何顺序含义纯粹是算法内部给社区编的号。同一个数据集换个随机种子、换个 scanpy 版本甚至只是调了一下resolution簇的编号就会重排。我踩过一次很典型的坑脚本里写死了“0号簇是 T 细胞、1号簇是单核细胞”结果换了台机器重跑簇编号全乱注释表直接对不上。正确的做法是把簇编号当成一个临时索引真正的身份靠 marker 基因判定。为了减少这种混乱我一般会在聚类之后、差异表达之前先跑一遍sc.tl.rank_genes_groups看看每个簇的 top 基因用一个字典把簇号和细胞类型对应起来再统一写回adata.obs[cell_type]。这样即使中间重跑导致编号变化你只需要改字典不用改后面所有代码。resolution的选择也值得说两句。默认 1.0 在几千个细胞的 PBMC 数据上通常给出 8 到 12 个簇粒度差不多。如果你的目标是找稀有亚群可以往上调到 1.5 甚至 2.0如果只是做粗分群0.5 就够。我不建议一开始就上高分辨率因为簇一多差异表达的比较次数会暴涨多重检验校正之后的pvals_adj会变得很难看很多真实的差异会被压下去。1.3 多批次数据的差异表达要怎么处理单样本数据没什么好纠结的直接按簇分组比较就行。但如果你是多个样本合并分析的情况就复杂了。假设你做的是合并后的聚类leiden是所有细胞混在一起分的那么某个簇里可能同时包含来自多个样本的细胞。这时候直接跑差异表达得到的基因里会混进大量批次效应带来的差异而不是真正的细胞类型差异。我的处理顺序是先做批次校正具体用什么方法看数据量样本量少用 Harmony 这类基于嵌入的校正样本量大可以考虑基于图的整合在校正后的嵌入空间上做邻居图和 Leiden然后用校正前的 log 表达矩阵做差异表达。这个组合听起来矛盾其实很合理聚类需要在去除批次效应的空间里做否则同一类型的细胞会被批次拆开但差异表达要的是真实的表达量对比用校正后的矩阵反而会把生物学信号也一起抹掉。如果样本量足够更稳妥的做法是在差异表达里显式加入协变量或者干脆按样本分开跑差异表达再取交集。rank_genes_groups本身没有内置的协变量参数需要在调用前把数据按样本拆分。这个操作成本高一些但对于结论的可信度提升很明显。2. 差异表达的参数选择与结果解读矩阵状态确认干净、簇编号心里有数之后就可以正式跑了。这一行代码sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon)谁都写得出来但method选什么、结果怎么筛、阈值卡多少这三个问题决定了你后面注释的成败。2.1 四种method到底选哪个scanpy 提供了几个统计检验方法差别不小我把常用的几个列成表method原理适用场景注意事项wilcoxon秩和检验非参数默认首选对表达分布不正态的数据稳健细胞数少时 p 值偏保守t-test参数检验假设正态数据量大、分布接近正态对离群值敏感稀疏基因容易失真t-test_overestim_vart 检验但方差不合并簇间细胞数差异悬殊时更敏感容易放出假阳性需要更严的校正logreg逻辑回归想看基因的判别能力而非单纯表达差异结果解释方式和前三个不同logfoldchanges含义要小心我自己的默认选择是wilcoxon。单细胞数据零膨胀严重表达分布根本不是正态的t 检验的前提不成立用非参方法更稳妥。t-test_overestim_var在簇大小极不均衡时确实更容易捞出差异基因比如一个 2000 细胞的大簇和一个 30 细胞的小簇比较但这个“更容易”里包含了不少假阳性我一般只在探索性分析时用正式结论还是回到 wilcoxon。还有一个容易被忽略的参数是pts。默认不输出表达比例加上ptsTrue之后结果里会多一列pts表示该基因在这个簇里表达值大于 0的细胞比例以及pts_rest表示在其他簇里的比例。这两个数字特别有用能帮你过滤掉那种“统计显著但只有三五个细胞表达”的基因。sc.tl.rank_genes_groups( adata, groupbyleiden, methodwilcoxon, use_rawTrue, ptsTrue )2.2 结果提取与阈值到底卡多少跑完之后结果存在adata.uns[rank_genes_groups]里是个结构化的嵌套字典直接看很费劲。用sc.get.rank_genes_groups_df转成 DataFrame 才是正常人的操作。df sc.get.rank_genes_groups_df(adata, group0)得到的表里有这么几列names、scores、logfoldchanges、pvals、pvals_adj、pct_nz_group、pct_nz_reference开了ptsTrue才有后两列。这里有个必须说清楚的细节logfoldchanges的底数是 2也就是说logfoldchanges1代表两倍差异2代表四倍。很多人误以为是自然对数把阈值卡在 0.25 就以为是 1.28 倍实际上那是 1.19 倍尺度差了一截。另外 scanpy 在计算倍数变化时如果检测到你的数据是 log1p 过的内部会把值还原到线性尺度再算比值所以你不用自己手动expm1但也正因为如此use_raw一定要指向正确的 log 矩阵否则还原出来的东西就离谱了。阈值怎么卡我给一个实操建议pvals_adj 0.05这是底线多重检验校正后的显著性logfoldchanges 0.5约等于 1.4 倍差异能过滤掉大量微小但显著的变化pct_nz_group 0.25这个簇里至少 25% 的细胞表达该基因避免被极少数细胞带偏三个条件同时卡通常每个簇能留下二三十个比较靠谱的候选基因。如果你只是想做初步注释其实放宽到pvals_adj 0.05加logfoldchanges 0.25也就够看了反正最后还要靠 marker 知识库交叉验证。提示scores这一列是检验统计量的标准化值可以粗略理解为“差异的强度”但不同 method 之间不可比。用它排序没问题别拿它当效应量。2.3 一次拿到一整套marker基因表格逐个簇提 DataFrame 太累我一般写个循环把所有簇的 top 基因汇总成一张表顺便把过滤条件也带上方便后面查。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 基因去搜可以画热图看整体格局也可以作为注释时的参考清单。我习惯再补一步把所有簇的 top 基因取并集去重之后作为后续点图的var_names这样一张图就能看到全局的 marker 分布比逐簇看 violin 效率高得多。这里再补充一个经验不要只看 top 10。有些细胞类型的特征基因在差异表达里排名并不靠前因为它在这个簇里表达量不高但特异性很强被那些高表达但有噪声的基因挤下去了。我一般会把每个簇的前 50 个基因都导出来扫一眼找那些“知道是 marker 但没进前十”的基因往往能纠正一些误判。3. 细胞类型注释从看marker到写标签差异表达给你的是“哪些基因不一样”注释解决的是“这些不一样的基因说明它是什么细胞”。这两步之间没有自动化的桥梁必须靠 marker 知识加人工判断。这是整个流程里最像“手艺活”的部分也是最容易糊弄过去的部分。3.1 一份能直接抄的PBMC marker清单外周血单核细胞是最经典的示例数据细胞类型相对标准marker 也比较明确。我整理了一份常用的对照表遇到新数据可以先拿这套去套套不上的再单独找细胞类型核心marker辅助marker备注CD4 T 细胞CD3D, CD3E, IL7RLTB, MALAT1与 CD8 共享 CD3 系列CD8 T 细胞CD3D, CD8A, CD8BGZMK, CCL5需与 NK 区分NK 细胞NKG7, GNLY, KLRD1NCAM1, PRF1不表达 CD3B 细胞CD79A, MS4A1, CD19CD37, BANK1MS4A1 即 CD20经典单核细胞LYZ, S100A8, S100A9CD14, VCAN占比通常最高非经典单核细胞FCGR3A, MS4A7CDKN1C, LST1FCGR3A 即 CD16树突状细胞FCER1A, CST3CLEC10A, CD1C数量少容易被并入单核浆细胞样树突IL3RA, CLEC4CGZMB, TCF4IL3RA 即 CD123血小板PPBP, GP9, PF4ITGA2B常被视为杂质红系前体HBB, HBA1, HBA2—高表达时注意是否污染用这张表的时候有个细节要注意marker 是组合看的不是单基因定生死。比如 LYZ 高表达可能是单核细胞也可能是中性粒细胞或者某些 DC 亚群CD3D 阳性但 CD8A 阴性大概率是 CD4 T。我一般会先看每个簇的 top 基因里有没有两个以上的该类型 marker再决定标签。另外血小板的 PPBP 和红系的 HBB 如果出现在某个“看起来像 T 细胞”的簇里很可能是双细胞或者环境 RNA 污染。处理方式取决于比例如果只是零星几个细胞直接忽略如果整个簇都这样那这个簇可能就是个伪簇考虑在聚类前就把这些基因高的细胞滤掉。3.2 打分法与阈值判定光看 top 基因列表有时候不够直观尤其是当一个簇的 marker 信号比较弱、top 基因里混着一堆核糖体基因的时候。这时候sc.tl.score_genes就派上用场了它能给每个细胞算一个基因集的综合得分然后你按簇看这个得分的分布。sc.tl.score_genes( adata, gene_list[CD3D, CD3E, CD2, IL7R], score_nameT_cell_score, ctrl_size50, n_bins25, use_rawTrue )几个参数的含义值得说清楚。ctrl_size是从表达量分箱里随机抽多少个对照基因默认 50调大一点结果更稳但计算慢n_bins是按平均表达量把基因分成多少个箱子默认 25主要目的是让对照基因的表达水平和你目标基因集匹配避免因为目标基因本身表达高就天然得分高。use_rawTrue这里同样重要用 scale 后的数据算出来的得分没法解释。得分算完怎么看我一般用sc.pl.violin按簇画出来看哪个簇的得分分布整体上移。判定标准不是“得分大于 0 就是阳性”而是看分布的形状如果某个簇的得分分布明显双峰说明这个簇内部可能混了两类细胞值得提高分辨率再分一次如果整个簇的分布都很平那说明这套 marker 对这个簇没有区分度。打分法还有一个经典用途是细胞周期评分这个放到后面第 5 节讲。3.3 自动注释与人工复核的配合现在基于参考数据集的自动注释工具不少思路基本一致拿一个已经注释好的参考数据集把你的数据映射过去用相似度或者分类器给出候选标签。这类方法的好处是快坏处是它一定会给你一个答案哪怕你的数据里存在参考集里没有的细胞类型。我的用法是自动注释出初稿人工复核定终稿。具体操作上先用自动方法给每个簇一个候选标签然后拿这个标签去反查差异表达结果如果自动注释说是 B 细胞那这个簇的 top 基因里应该出现 CD79A、MS4A1 这些对得上就通过对不上就手动改。这个交叉验证的步骤花不了多少时间但能挡掉大部分明显的错误。还有一种情况是自动注释给出的标签特别细比如把 T 细胞分成了七八个亚群但你的数据里每个亚群只有几十个细胞。这种细粒度标签在统计上站不住脚我一般会往上合并一层统一叫 T 细胞然后在讨论里注明“未进一步细分亚群”。别为了显得分析深入就硬保留细分标签细胞数不够的亚群差异表达结果基本是噪声。3.4 把注释写回obs并保持可追溯注释结果最终要落到adata.obs里通常是新增一列cell_type。这里我强烈建议同时保留原始的leiden列不要覆盖。后续做任何分组可视化、差异表达你都需要能在“聚类编号”和“细胞类型”之间来回切换。比如你想看某个细胞类型内部是否还有异质性就得回到 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而不用循环赋值好处是映射关系一目了然而且没映射到的簇会自动变成缺失值你能立刻发现漏了哪个。映射完记得看一眼adata.obs[cell_type].value_counts()确认没有大量 NaN。另外把这份字典单独存成一个文件比如 json 或者 csv跟分析脚本放一起下次别人接手或者你自己三个月后回看能立刻搞清楚当时的判断依据。4. 可视化细节从能看到能放进汇报图和图之间的差距在单细胞分析里体现得特别明显。同一份数据有人画出来的点图一眼就能看出细胞类型分布有人画出来糊成一团什么都看不清。这一节专讲参数。4.1 dotplot的关键参数点图是 marker 展示的主力横轴基因、纵轴簇点的大小表示表达比例颜色表示平均表达量。核心参数就那么几个但每个都能显著改变可读性。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是按基因做标准化把每个基因的表达量缩放到 0 到 1 之间。这个参数非常关键不同基因的绝对表达量差异巨大CD3D 可能均值 3PPBP 可能均值 8不标准化的话颜色分布完全被高表达基因主导低表达但特异的基因就看不出来了。按基因标准化之后每个基因的颜色都是相对自己的尺度特异性一目了然。dot_max和dot_min控制点大小的映射范围。默认dot_max是 1.0也就是表达比例最高的那个点会撑满如果你的数据里有个基因在所有簇里都接近 100% 表达那其他点就全被压得很小。把dot_max设成 0.8 或者 0.6能让中等表达的点的差异更明显。cmapRdBu_r是红蓝配色配合standard_scalevar的时候红色代表高表达、蓝色代表低表达中间是白色。我不建议用默认的 viridis因为它在表示“有正有负的相对高低”时不如发散型配色直观。当然如果你要投的期刊有配色要求另说。还有一个实用技巧基因分组。var_names可以直接传一个字典键是细胞类型值是对应的基因列表这样点图会自动按组排列基因并加上分组标签比自己手动排顺序省事得多。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 的取舍这三种图都能展示多基因多组的表达格局但侧重不同选错了会事倍功半。sc.pl.stacked_violin展示的是每个基因在每个组里的完整表达分布能看出双峰、零膨胀这些结构信息。缺点是基因一多图会变得极长。我一般只在需要展示两三个关键基因的分布细节时用它比如确认某个 marker 是不是真的在目标簇里特异性高表达还是所有簇都有一小撮阳性细胞。sc.pl.matrixplot是点图的简化版只保留颜色不保留点的大小适合临床或者跨学科读者信息密度低但易懂。如果汇报对象不熟悉单细胞用这个比点图友好。sc.pl.rank_genes_groups_heatmap是直接基于差异表达结果画热图每个簇取 top N 个基因。这张图的信息量最大能同时看到多个簇的 marker 格局但基因太多时行标签会糊成一片。我一般取每个簇前 5 到 8 个基因n_genes5配合swap_axes调整方向出来的图比较清爽。sc.pl.rank_genes_groups_heatmap( adata, n_genes5, groupbycell_type, show_gene_labelsTrue, use_rawTrue, cmapRdBu_r )提示热图默认会对表达量做标准化如果图上的颜色和你预期的表达高低对不上先检查standard_scale参数和use_raw设置。这两个是热图翻车的高频原因。4.3 出图保存与配色sc.pl系列函数都支持save参数传个文件名就会自动存到figures/目录下格式默认是 png可以通过sc.settings.figdir改路径、sc.settings.file_format_figs改格式。我一般设成 pdf矢量图放进汇报里放大不失真。配色上连续变量用发散型或者顺序型配色分类变量比如 UMAP 上的细胞类型多用定性配色。scanpy 默认的分类配色在类别超过 10 种之后会开始重复或者相近容易混淆。我习惯手动指定把主要的细胞类型分配差异大的颜色稀有的类型用灰色系这样 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把标签直接标在簇中心比侧边图例清晰。但如果你的簇有重叠on data会挤在一起那就退回默认的right margin。5. 往前一步细胞周期、轨迹与通路差异表达和注释做完一个完整的描述性分析就闭环了。但如果想再往前走一步看看细胞状态和分化关系还有几件事值得做。5.1 细胞周期回归增殖活跃的细胞在转录组上会和同类型的静息细胞拉开很大距离聚类的时候经常出现一个“细胞周期簇”把本该在一起的细胞拆开。处理方式是给每个细胞算 S 期和 G2M 期得分然后在降维之前把这两个得分的影响回归掉。# 使用公开的 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之前调用而且要在 log 归一化之后。regress_out这一步在细胞数上万的时候会很慢因为它是逐基因做线性回归。如果数据量大可以考虑用sc.pp.scale之前先做sc.pp.highly_variable_genes只对高变基因做回归能省不少时间。还有一个反直觉的点细胞周期回归不是必须的。如果你的研究对象本身就是增殖过程比如发育或者肿瘤那细胞周期差异就是生物学信号不该回归掉。只有在研究细胞类型构成、且确认细胞周期主导了聚类结果的时候才需要这一步。5.2 PAGA与DPT想看细胞类型之间的连接关系PAGA 是个好选择。它把每个簇压缩成一个节点用簇间的连接强度作为边的权重能直观看出哪些细胞类型在转录组上更接近。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会同时画出 PAGA 图和 UMAP两张图并排看能验证 UMAP 上的邻近关系是否和 PAGA 的连接强度一致。如果不一致往往说明 UMAP 的局部结构有失真这时候 PAGA 的结果更可信。DPT 是沿着 PAGA 的图结构算伪时间的需要一个明确的起点细胞。这个起点怎么选取决于你的生物学问题通常选一个处于分化起始状态的细胞类型里的代表性细胞。adata.uns[iroot] 0 # 替换成你选定的起始细胞的索引 sc.tl.dpt(adata)DPT 的结果对起点极其敏感换个起点伪时间轴就完全反过来。所以我一般会试几个候选起点看哪个得到的伪时间排序符合已知的生物学常识再定下来。5.3 基因集打分看通路除了细胞周期score_genes可以用来评估任意基因集的活性比如炎症相关、氧化磷酸化、细胞毒性等等。做法和前面一样关键是基因集来源要可靠别自己随手拼。打分结果可以按细胞类型画成热图看哪类细胞在哪个通路上更活跃。这一步的价值在于它能把差异表达得到的“一堆基因”提升到“功能层面”让结论更容易解释。比如你发现某个单核细胞亚群在炎症基因集上得分显著更高这比列出一百个上调基因有说服力得多。需要注意的是通路打分和差异表达是两套逻辑前者是基因集整体水平后者是单基因水平结论有可能不一致。不一致的时候优先相信差异表达的单基因结果因为它更直接通路得分受基因集质量和大小影响很大适合作为辅助证据。6. 踩坑实录与速查表写了一堆方法最后还是落到实际会遇到什么问题。下面这张表是我自己和我带过的同学最常撞的几类按出现频率排序。6.1 典型问题排查表现象可能原因排查方法处理方式所有簇的 top 基因都是 MT- 或 RPS 开头质控没做好或差异表达用了未过滤数据检查adata.var_names里线粒体基因比例回到质控步骤重新过滤logfoldchanges 出现大量负值或异常大X是 scale 后的 z-scoreprint(adata.X.max())看是否超过 10改用use_rawTrue某个簇没有任何显著差异基因簇太小或该簇是双细胞/伪簇看该簇细胞数和 marker 分布合并到邻近簇或滤除点图上所有点颜色接近没做standard_scale检查参数加standard_scalevarUMAP 上同类细胞被拆成几块批次效应未校正按样本着色看分布做批次校正后重聚类use_rawTrue报错adata.raw为 Noneprint(adata.raw)在 scale 前补上adata.raw adataleiden 报错找不到模块缺少 igraph / leidenalgpip list检查安装对应依赖或改用flavorigraph6.2 几个容易忽略的实操细节第一个是随机种子。Leiden 聚类、UMAP、t-SNE 都有随机性不设种子的话每次重跑图都不一样。在脚本开头加一句sc.settings.seed 0能保证结果可复现。我做过的项目里有一半以上的“图怎么变了”问题都是没设种子导致的。第二个是数据类型的坑。adata.obs里的分组列必须是 categorical 或者至少是字符串如果是整数类型rank_genes_groups在某些版本里会把它们当连续变量处理结果完全不对。用astype(category)显式转换一下能省掉很多莫名其妙的报错。第三个是基因名的大小写和别名。同一个基因在不同注释版本里可能叫CD3D也可能叫CD3dmarker 列表直接拷过来匹配不上是常事。判断方法是用adata.var_names.isin(gene_list).sum()看能匹配上几个如果明显少于预期就得检查命名规范必要时做大写转换或者别名映射。第四个是保存中间结果。单细胞流程跑一遍动辄十几分钟做完注释之后的AnnData一定要adata.write(annotated.h5ad)存下来。我见过太多人注释完直接关了 notebook第二天想画个图发现得从头跑。存的时候注意write默认会保留raw文件会比较大可以用adata.write(x.h5ad, compressiongzip)压缩一下。最后一个也是我觉得最重要的注释结论要留证据。不要只留一个cell_type列把当时的点图、差异表达表格、映射字典都存到一个文件夹里。半年后你或者别人回头看能立刻知道每个标签是怎么来的。这件事在单独分析的时候感觉多余但在需要复现或者答辩的时候能救命。我个人在实际操作中的体会是单细胞分析真正花时间的从来不是写代码而是反复在“这个簇到底是什么”上纠结。差异表达给你的是线索不是答案marker 知识库给你的是参照也不是答案。真正的答案往往来自多个证据的交叉top 基因对得上、打分分布支持、自动注释也这么认为、文献里这个组织确实有这类细胞。四条里对上三条就可以定下来了剩下的那一条对不上记下来别硬圆。

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑