复现这篇论文，最容易犯的错只有一行代码。把细胞类型等于 Myeloid 的行数出来，患者年龄与髓系细胞比例的相关系数是 0.139。把另一批数据中使用的 McDC 标签也计入，相关系数变成 0.261。患者还是那 168 位，改变的是我们怎样识别不同研究里的同一类细胞。下面从原始公开对象开始，逐步解释这两个数字，也交代我们和论文究竟对上了哪一步。Peng 与 Zhan，npj Precision Oncology，2025，原文图 1
研究问题和方法亮点
研究者想知道，较年轻的结直肠癌患者，肿瘤中的免疫细胞组成是否不同。他们整合六项公开单细胞研究，覆盖 22 到 91 岁的 168 位患者。论文的图 1 具体问的是，每位患者捕获的肿瘤组织细胞中，髓系细胞占多少。图 1C 数各年龄组有几位患者，图 1G 把年龄当作连续变量，每位患者画一个点。论文随后在 TCGA 结肠腺癌的 bulk RNA-seq 中用 BayesPrism 估计细胞比例，观察同方向的关系。后者提供了另一个队列的证据，测量方式也随之改变。论文图 1 与方法
这套设计值得借鉴的地方，是把单细胞的分辨率和患者层面的比较接了起来。假如直接把几万个细胞当成独立患者，建库时多捕获了细胞的样本就会获得更大权重。下面所有比例都以单个患者为单位。分母包括该患者通过质控的上皮、免疫、基质等所有细胞。如果先取免疫细胞子集再算，回答的就变成“髓系细胞在免疫细胞中占多少”。
第一步 选一个适合自己的起点
想先看方法，下载本文脚本包，解压后在目录里执行 Rscript --vanilla quickcheck.R age_moments.tsv 即可。这个小文件保存每个年龄下的人数、患者比例之和及平方和，足以重算患者级 Pearson 相关和双侧 P 值。若只拿年龄组均值做相关，样本量和标准误都会变。想检查质控、重新定义细胞类型，则按下面步骤使用作者对象。下载本文 R 脚本、轻量数据与核对结果
准备 R、readxl 包以及两个原始文件。Zenodo 记录里只有一个约 3.87 GB 的 crc_seu_raw.rds。补充材料压缩包中，文件名以 MOESM2_ESM.xlsx 结尾的工作簿存有研究来源和患者临床表。下例在 Linux 或 macOS 的终端里运行，已解压本文脚本包的目录作为工作目录。作者的 Zenodo 对象 · 论文补充材料下载
curl -fL --retry 3 -C - -o crc_seu_raw.rds \
'https://zenodo.org/records/15267716/files/crc_seu_raw.rds?download=1'
curl -fL --retry 3 -o supp.zip \
'https://www.ebi.ac.uk/europepmc/webservices/rest/PMC12627461/supplementaryFiles'
unzip -p supp.zip '*MOESM2_ESM.xlsx' > supp.xlsx
Rscript -e 'stopifnot(unname(tools::md5sum("crc_seu_raw.rds")) == "5730f7ff6b48be474a4846d278fbfaac")'
Rscript -e 'if (!requireNamespace("readxl", quietly=TRUE)) install.packages("readxl", repos="https://cloud.r-project.org")'
Rscript --vanilla reproduce.R crc_seu_raw.rds supp.xlsx out
本次用 R 4.3.3 从完整 RDS 实跑，进程峰值约 18.8 GiB。文件下载大小和载入内存差别很大，因为读入整个对象时还会展开表达矩阵。脚本一取出元数据就释放大对象，后面的年龄统计只处理一张表。内存紧张时，可以先走轻量路线学习计算过程。
第二步 先对齐患者，再复现图 1C
脚本使用 Supplementary Data 2 的 patientID 对齐临床表和细胞表，检查患者编号唯一、两侧编号集合相同，并核对每个细胞记录的年龄。按补表 Group 计数，四组分别有 17、57、79、15 位患者，与图 1C 的四根柱子一致。论文以小于 40 岁作为 G1，读者若研究通常定义的小于 50 岁早发病例，应重新定义分组、重新计数。作者图 1C 和分组说明
clinical <- as.data.frame(readxl::read_excel("supp.xlsx", sheet = "Supplementary Data 2"))
stopifnot(!anyDuplicated(clinical$patientID))
table(clinical$Group)
with(subset(clinical, Age %in% c(60, 80)), table(Age, Group))
最后一行专查边界。补表把全部 8 位 60 岁患者列入 G3，2 位 80 岁患者列入 G4。这是值得单独运行的检查。仅按正文中的“40–60”“60–80”“>80”重新写条件，可能把同一位患者重复分组、放到相邻组或漏掉。要复现论文，直接用补表标签。要改成自己的临床问题，把新区间的端点和人数一并记下来。
第三步 逐个核对质控指标，找出第一处分歧
原始对象有 560,238 个细胞。论文方法写明剔除检测基因数小于 400、UMI 小于 1,000 或线粒体比例大于 30% 的细胞。对象中的 percent_mt 已经是百分数。下面用原有元数据列执行一次交集筛选。论文单细胞方法和数据来源 · Zenodo 对象版本
object <- readRDS("crc_seu_raw.rds")
m <- attr(object, "meta.data", exact = TRUE)
rm(object)
invisible(gc())
keep <- with(m, nFeature_RNA >= 400 & nCount_RNA >= 1000 & percent_mt <= 30)
table(keep)
with(m[keep, ], table(dataset, celltype))
实跑中，26,214 个细胞的基因数低于门槛，50,842 个的 UMI 低于门槛，25,942 个的线粒体比例超过门槛。同一细胞可以触发多项，三组剔除数不能直接相加。三项同时通过的有 479,270 个，比论文图 1D 报告的 554,930 个少 75,660 个。仅凭现有元数据和文字方法，无法定位差异出在哪个处理环节。我们保留这两个数，并在后续重算中使用实际通过的 479,270 个。
第四步 审查标签，重算图 1G 的年龄相关
论文用 Seurat 4.1.1 归一化，选 2,000 个高变基因，以前 30 个主成分做 Harmony 样本整合，用分辨率 0.3 聚类，再根据标记基因等信息标出九个大类。作者举出的标记包括髓系的 CD14 和浆细胞的 JCHAIN。Zenodo 对象的 celltype 列实际有 20 个标签，部分来源使用 Myeloid，另一些使用 McDC，还单列 pDC、Neutrophils、Mast 等。我们没有获得每个细胞对应图 1D 最终九类的映射，本文代码属于根据公开资料编写的重算脚本。论文方法、数据可用性说明 · 本文脚本
运行脚本会得到 out/labels_by_dataset.tsv。打开这张来源与标签交叉表，你会看到部分来源没有一个细胞叫 Myeloid，却有数千个叫 McDC。按字面筛选就会把这些来源的髓系细胞数记成零。我们保留字面口径作对照，再试一个明确可复查的合并口径，只合并这两个标签。这里的合并是一项敏感性分析，pDC、中性粒细胞和肥大细胞等单列标签均未并入。要建立自己的髓系定义，还需要核查来源研究的注释说明和标记表达。
id <- factor(m$patientID[keep], levels = clinical$patientID)
labels <- m$celltype[keep]
total <- as.integer(table(id))
stopifnot(all(total > 0))
myeloid <- as.integer(table(id[labels %in% c("Myeloid", "McDC")]))
fraction <- myeloid / total
cor.test(clinical$Age, fraction, method = "pearson", alternative = "two.sided")
factor 的 levels 固定了 168 位患者的顺序。某位患者通过质控、但该口径下没有髓系细胞时，分子会正确记为 0。保留这些零才能避免只对“检出过髓系细胞”的人作相关。下表的组均值也是先算每位患者的比例，再在组内取平均。
| 数据与口径 | G1 | G2 | G3 | G4 | 年龄相关 r | 双侧 P |
|---|---|---|---|---|---|---|
| 本文，仅计 Myeloid | 6.21% | 8.29% | 8.59% | 15.97% | 0.139 | 0.0731 |
| 本文，合并 Myeloid 与 McDC | 8.68% | 11.22% | 14.86% | 21.25% | 0.261 | 0.000625 |
| 论文图 1G 的最终注释 | — | — | — | — | 0.25 | 0.0011 |
合并口径下的方向和量级与原图接近，95% 置信区间为 0.114 至 0.397。质控后的细胞集合和注释口径都与论文有尚未解决的差异，因此不能把相近的系数叫作精确复现。脚本另存 per_patient.tsv 和 age_myeloid.pdf，便于检查散点和离群患者。原图 1H 的 TCGA 反卷积属于另一项分析，这里的对照值取自图 1G。论文图 1G、H
把这个方法用到自己的队列
先做两张小表。第一张数“研究来源 × 年龄段”的患者数，找出哪些队列同时有年轻和年长患者。第二张数“来源 × 原始细胞类型”的细胞数，找出哪几种名称要核对。接着按来源分别画年龄和比例的散点，再逐个去掉一项研究重算。如果去掉某个来源后趋势消失，就要查看该研究的取材部位、样本制备和年龄构成。Harmony 可以改善表达空间的批次混合，不能让不同队列的捕获效率自动一致。
最后把结论写回它的分母。这里观察到的是捕获细胞中的相对比例。即使髓系细胞绝对数量不变，其他谱系捕获减少，也可能让这个比例升高。最年轻组只有 17 位患者，且各人年龄、部位和临床特征不能随意互换。这张图能提出“年龄与肿瘤免疫组成相关”的问题，单凭它无法判断某位患者的免疫治疗反应。
如果你的队列也出现“换一个注释口径，结论就变了”，可以先整理脱敏的分组人数、各来源标签和一张患者级散点图。带这三样材料与云生信讨论，才好判断接下来该统一注释、调整比较设计，还是补一次分队列验证。云生信课题分析与研究路线评估
