最近生信圈在转的一件事，应该就是这份“中国人群视网膜衰老单细胞图谱”的公开。做眼科研究的人在等它，做衰老研究的人在等它，不少想系统练一遍单细胞流程的初学者也在等它。等它发布后大家发现，作者不仅把数据放了上去，还把完整的注释结果、marker基因列表、核心绘图代码一起公开，等于把一份带参考答案的习题集摆在了桌上。这篇文章我就围绕这份图谱展开，把数据背后做了什么、为什么值得复现、获取数据后怎么从零开始跑通全流程、以及实际复现中容易掉进去的坑全部梳理一遍。哪怕你之前没跑过单细胞项目，对照这份图谱也能把Seurat那套核心分析吃透。

1. 这份图谱凭什么引人关注：项目定位与核心价值拆解

1.1 视网膜为什么是衰老研究的理想模型

视网膜是中枢神经系统里少数可以直接被观察到的组织。它结构分层清晰，细胞类型明确，主要由感光细胞（视杆和视锥）、双极细胞、无长突细胞、水平细胞、神经节细胞（RGC）、Müller胶质细胞、小胶质细胞、血管内皮细胞等组成。不同细胞各司其职，又彼此通过突触和旁分泌信号紧密耦合，任何一个细胞亚群的衰老变化都可能牵动整条视觉传导链。

更重要的是，视网膜的衰老和多种高发疾病直接相关。年龄相关性黄斑变性（AMD）是老年人视力丧失的主要原因，糖尿病视网膜病变、青光眼的发病率也随年龄显著上升。研究视网膜衰老，本质上是在为这些退行性疾病寻找早期标记和干预窗口。而之前公开的单细胞视网膜数据大多来自欧美人群，东亚人群的参考图谱长期缺乏。这份中国人群视网膜衰老图谱，恰好补上了这个缺口。

1.2 图谱数据集的整体构成

从发布资料来看，这份图谱覆盖了不同年龄阶段的中国人群视网膜样本，经过严格的单细胞转录组测序之后，系统性地完成了细胞注释和衰老相关分析。数据公开部分主要包括以下几个层级：

• 原始测序数据或比对后的表达矩阵，用于后续独立分析；
• 经过质控、降维聚类和注释后的细胞亚群信息，每一类细胞都有对应的marker基因列表；
• 衰老相关的差异表达结果、功能富集结果；
• 用于生成论文核心图表的完整分析代码。

对于复现者来说，表达矩阵加上注释信息就是最值钱的部分。你不用重复跑一遍cell ranger比对，可以直接从质控后的矩阵起步，节省大量时间和计算资源。这一点对没有服务器集群的个人研究者特别友好。

1.3 这份图谱适合谁去复现

我大致把适合复现的人群分成三类。

第一类是做眼科和衰老基础研究的人。他们关心的是：中国人群视网膜里到底有哪些细胞亚群在衰老中变化最大，哪些基因可以作为候选靶点。复现一遍注释和差异分析后，可以得到自己的分析结果，也可以在这些结果基础上直接做后续功能实验。

第二类是刚入门的单细胞生信学习者。这份数据注释完整、代码公开、图谱结构清晰，比自建课题数据练手要舒服得多。你可以拿着它完整走一遍“质控-聚类-注释-差异-可视化”的标准流程，遇到问题还有论文结果可以对照。

第三类是做生信工具开发或流程搭建的人。有了标准数据集，就能在本地验证新方法、做多数据集整合，或者用来跑自动注释工具的性能评测。公开图谱某种程度上承担了benchmark数据集的角色。

2. 图谱背后的技术逻辑：单细胞分析流程是怎么设计的

2.1 建库方式与表达矩阵的格式判断

视网膜组织比较特殊，细胞种类多、组织致密度高，而且感光细胞等大细胞比例高。从这类公开图谱的通例来看，大多采用10x Genomics Chromium平台进行单细胞转录组建库。你拿到手的矩阵通常是一个高度稀疏的计数矩阵，行是基因，列是细胞，数值表示每个细胞中检测到的转录本count数。

判断建库方式最直接的方法看文件格式。如果下载目录里包含filtered_feature_bc_matrix这样带barcode和feature的文件夹，基本就是10x平台的产物。如果拿到的是.csv.gz或者matrix.mtx.gz，也大概率是10x格式转换过来的。复现前先确认矩阵格式，避免后面读取时反复报错。

2.2 质控参数的选择逻辑

单细胞质控的核心是筛掉低质量细胞和空液滴。常规做法看三个指标：

• 每个细胞检测到的基因数（nFeature_RNA），过低说明细胞破裂或转录本捕获失败，过高可能是一个液滴包含了两个细胞；
• UMI总数（nCount_RNA），是测序深度的直接体现；
• 线粒体基因比例（percent.mt），过高说明细胞状态差，胞质RNA流失而线粒体RNA相对富集。

很多人直接在Seurat里写subset(nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 20)，但这份图谱发布在不同组织上，视网膜感光细胞转录本数量本身偏高，Müller胶质细胞代谢活跃，线粒体比例分布也和外周血完全不同。我的建议是先画violin图和QC散点图，观察分布边界，再结合实际确定的阈值。数据发布方如果给出了QC阈值，优先以原始阈值为主，复现时再去套通用参数很容易得到和原文不一致的细胞数。

2.3 降维聚类与注释方法的组合思路

拿到质控后的细胞后，标准路径是归一化、高变基因筛选、PCA降维，然后用UMAP或t-SNE展示，最后做聚类分群。Seurat里的FindClusters使用的Louvain算法需要指定分辨率，分辨率越高分群越多。不要一个分辨率跑到底，可以按0.5、0.8、1.0分别试，对比marker表达后选择生物学上最合理的分群数。

细胞注释是整条流程中最依赖经验的环节。这个项目里我推测采用的是“marker基因打分+参考图谱映射+人工核对”的组合方式。视网膜的经典marker比较明确：视杆细胞看RHO、NRL，视锥细胞看ARR3、OPN1MW，双极细胞看VSX2，无长突细胞看TFAP2A、GAD1，RGC看RBPMS、SNCG，Müller胶质看RLBP1、GLUL，小胶质看P2RY12、C1QA，血管内皮看CLDN5、PECAM1。先拿这些经典marker做DoHeatmap和FeaturePlot，按表达模式把大类分出来，再往下分亚群。

2.4 衰老相关分析的常用方法

图谱既然叫“衰老图谱”，分析重点自然落在年龄相关的细胞状态变化上。常见做法有几类：

• 分年龄组做差异表达分析，看每个细胞亚群中随年龄上调或下调的基因；
• 用基因集打分（AddModuleScore）评估衰老相关通路活性，比如炎症、氧化应激、DNA损伤修复、线粒体功能障碍相关基因集；
• 做细胞通讯分析，比较不同年龄组之间配体-受体相互作用的强度变化，细胞通讯常用的工具是CellChat；
• 做拟时序分析，用Monocle3等工具推断细胞分化或状态转变轨迹。

这些分析在公开代码里大概率都有对应脚本。复现的时候不需要每个脚本都跑，先挑和你研究问题最相关的部分，跑通了再补其他分析。

3. 复现之前的关键准备：数据获取、环境搭建与资源评估

3.1 数据去哪里找

先回答一个很多人私信问的问题：数据从哪下。

这类公开图谱数据的存放通常分两个地方。其一，国际通用数据库GEO（Gene Expression Omnibus），检索关键词用“retina aging single cell Chinese”之类，找到对应编号的GSE数据即可下载。其二，国内的国家基因组科学数据中心（NGDC）下的GSA子库，检索“视网膜 衰老 单细胞”或项目名称也可以找到。作者通常还会在论文的Data availability段落写明两个数据库的访问编号，优先按论文给的编号去检索，最准确。

下载时注意区分两个概念：raw data是原始测序下机数据，体积很大，按T计算；processed data是处理后的表达矩阵，通常几十到几百MB，复现分析用它就够了。非必要不下raw data，除非你要自己重新比对。

3.2 复现环境与依赖清单

跑这套流程，主流方案是R + Seurat。除了Seurat本身，还需要一批配套包。我建议按功能分组安装：

• 数据读取与矩阵格式转换：Seurat、SeuratData、Matrix；
• 质控与可视化：ggplot2、patchwork、dplyr、RColorBrewer；
• 批次整合与多样本合并：harmony（或者Seurat的IntegrateData）；
• 差异分析：Seurat内置的FindMarkers就够了，但有时需要MAST、DESeq2做补充；
• 功能富集：clusterProfiler、org.Hs.eg.db；
• 细胞通讯和拟时序：CellChat、monocle3，这两个包依赖比较多，建议单独花时间装。

R版本建议用4.2以上，Seurat用5.x版本。注意Seurat 5和Seurat 4的对象结构有差异，如果作者代码是基于旧版写的，读入后可能要用UpdateSeuratObject做一次对象升级，否则某些函数会报错。

3.3 计算资源与目录规划

很多人低估了单细胞分析对内存的需求。一个典型的视网膜图谱，细胞数量通常在几万到十几万之间，基因数两万以上。用Seurat做完整流程，16GB内存勉强能跑，32GB会比较舒服。如果你打算在本地笔记本上玩，建议先把读取矩阵后的对象用saveRDS存下来，后续聚类、差异分析都从RDS对象读取，不要反复从原始矩阵重来。

内存不够还有一个折中方案：用Python的Scanpy做完整流程，Scanpy对内存的管理比R好一些，不过可视化风格和Seurat差异比较大，和原文代码对照起来不如R方便。

目录规划上我建议按这个结构放文件：

retina_atlas/ ├── data/ │ ├── raw_matrix/ │ └── metadata/ ├── scripts/ │ ├── 01_qc.R │ ├── 02_cluster.R │ ├── 03_annotation.R │ └── 04_diff_exp.R ├── output/ │ ├── figures/ │ └── tables/ └── rds/ └── retina_seurat.rds

这样每一步的输入输出都很清晰，回头排查问题也方便。

4. 全流程复现实操：从只读矩阵到核心图表

4.1 第一步：读入矩阵与基础质控

确认你拿到的是filtered矩阵（只包含细胞barcode），就可以用Read10x读入。如果拿到的是其他格式，用readRDS或read.csv读入后转换成Seurat对象。

library(Seurat) library(dplyr) library(Matrix) # 读取10x格式矩阵 data_dir <- "data/raw_matrix/filtered_feature_bc_matrix" counts <- Read10X(data.dir = data_dir) # 创建Seurat对象 so <- CreateSeuratObject(counts = counts, project = "Retina_Aging_CN", min.cells = 3, min.features = 200) # 计算线粒体基因比例 so[["percent.mt"]] <- PercentageFeatureSet(so, pattern = "^MT-")

这里有两个参数值得说明。min.cells = 3要求基因至少在3个细胞中表达，用来过滤在组织中几乎不表达的基因。min.features = 200用来过滤基因检出数过少的细胞。这两个值属于通用设定，之后还会做更严格的质控过滤。

质控可视化是决定阈值的第一步：

VlnPlot(so, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.01) # 根据分布进行过滤，阈值需要结合数据分布调整 so <- subset(so, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 15)

我做完这一步通常会再把过滤前后的细胞数对比一下。如果过滤掉了超过30%的细胞，说明原始数据质量不好或阈值设得太严，此时不要强行继续，先回到分布图上重新判断。

4.2 第二步：归一化、高变基因与批次整合

视网膜图谱如果包含多个样本，样本之间必然存在测序深度、批次效应。处理批次最常用的方式就是harmony。它是Python的harmony-py和R的harmony包，本质是迭代聚类降维，把混杂信息从PCA嵌入中移除。

so <- NormalizeData(so) so <- FindVariableFeatures(so, selection.method = "vst", nfeatures = 3000) so <- ScaleData(so) so <- RunPCA(so, npcs = 30) # 如果数据来自多个个体，建议按样本ID做harmony整合 so <- RunHarmony(so, group.by.vars = "sample_id")

多单样本整合后再聚类，得到的分群通常比直接用原始PCA聚类更干净。有人会直接跳过harmony，但我觉得对于这种跨年龄、多个体样本的图谱，整合这一步不要省，否则后续注释看到的可能不是生物学差异而是批次差异。

4.3 第三步：聚类分群与细胞注释

聚类分辨率的选择我前面提过，实际操作是从低到高扫一遍。

so <- FindNeighbors(so, reduction = "harmony", dims = 1:30) so <- FindClusters(so, resolution = c(0.5, 0.8, 1.0)) # 先看不同分辨率的聚类数 sapply(seq_along(levels(so@meta.data$RNA_snn_res.0.8)), function(i) i)

选定合适分辨率后，先用经典marker确认大类。比如我想确认感光细胞群，就画这几个基因的FeaturePlot和DotPlot：

FeaturePlot(so, features = c("RHO", "ARR3", "RBPMS", "RLBP1", "P2RY12"), cols = c("lightgrey", "red"), ncol = 3) DotPlot(so, features = c("RHO", "ARR3", "NRL", "VSX2", "GAD1", "RBPMS", "SNCG", "RLBP1", "GLUL", "P2RY12", "C1QA", "CLDN5"))

看到特征表达模式后，手动给每个cluster赋细胞类型标签。需要注意，亚群注释这一步主观性很强，不同人会得出不同粒度。比如Müller胶质细胞可能被分成一个群，也可能在不同状态下被分成应激态和静息态；小胶质细胞也可能存在homeostatic和activated两种状态。这些细节和论文报道不一定完全一致，我的处理原则是：大类一定要和原文一致，亚群按自己的聚类结果重新判断，并在方法部分说明差异。

4.4 第四步：衰老差异分析与可视化

注释完成后，就可以做年龄相关的差异分析。把metadata里加上年龄分组，比如young和old，然后按细胞类型分别跑FindMarkers。

so$age_group <- ifelse(so$age < 60, "young", "old") so$age_group <- factor(so$age_group, levels = c("young", "old")) Idents(so) <- "celltype" diff_list <- list() for (ct in levels(so$celltype)) { diff <- FindMarkers(so, ident.1 = "old", ident.2 = "young", subset.ident = ct, min.pct = 0.1, logfc.threshold = 0.25) diff$gene <- rownames(diff) diff_list[[ct]] <- diff }

这里subset.ident参数直接指定在特定细胞类型内部做比较，逻辑是“同一个细胞类型中，老年组相对年轻组有哪些基因变化”。如果样本来自多个个体，建议用test.use = "MAST"或者混合模型，可以更好处理个体间的随机效应。

拿到差异基因后，接着做功能富集：

library(clusterProfiler) library(org.Hs.eg.db) # 以Müller胶质细胞为例 genes_up <- diff_list$Muller$gene[diff_list$Muller$avg_log2FC > 0.5] ego <- enrichGO(gene = genes_up, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", pAdjustMethod = "BH")

最后把UMAP、marker气泡图、差异火山图、富集条形图一起输出，复现基本就完成了。论文里的核心结论一般都能在你自己跑出来的图上看到对应趋势，如果趋势完全相反，那就要回到前面的步骤找问题。

5. 复现过程中的常见问题与排查心得

5.1 数据读取阶段报错

这个阶段最常见的报错是Read10X找不到barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三个文件，或者目录名不对。10x新版输出目录是filtered_feature_bc_matrix，里面又分了一级目录。直接指向包含三个文件的文件夹即可，不要多指一层。

另一个高频问题是矩阵中基因名格式。人类视网膜数据如果上游用Ensembl ID做基因名，后面所有marker检测都会失效。看到ENSG00000169174这种格式，先做ID转换再继续，不要硬着头皮往下跑。

5.2 内存不足和会话崩溃

R跑单细胞流程，长期运行后内存碎片化严重，最后直接session terminated。我的经验是分步执行，每跑完一个大步骤就saveRDS一次，下次从头文件恢复，而不是把全部代码放在一个脚本里一口气跑完。

如果16GB内存机器上还是OOM，可以限制工作线程数并开启gc：

options(future.globals.maxSize = 8 * 1024^3) options(future.rng.onMisuse = "ignore") # 减少并行线程 plan(strategy = "sequential") # 手动触发垃圾回收 gc()

另外，FindMarkers跑几十个细胞类型时非常占用资源，循环里每跑完一个类型就gc()一次，能明显降低峰值内存。

5.3 注释结果和论文不一致

拿到自己聚类结果后很可能发现cluster数量和论文不完全一致。这是正常现象。可能原因包括：过滤阈值略有差异、harmony参数设置不同、分辨率选择不同、甚至Seurat和Scanpy的聚类算法本身就有差异。我会把“大类注释一致、亚群合理解释了marker表达”作为复现成功的标准，而不是追求cluster编号和论文完全一样。

5.4 复现质量怎么验证

我常用的验证方法有三种。

第一，看关键marker是否只在对应细胞群表达，比如RBPMS只在RGC群表达，如果它在感光细胞群里也高表达，说明注释出了问题。

第二，看细胞比例趋势。论文报道中老年组某种细胞比例降低或升高，在你的复现结果里方向应该一致。

第三，看差异基因的特征。比如老年视网膜中炎症相关基因在小胶质和Müller胶质中上调，这类生物学趋势是稳定的。方向不对就回去检查分组变量是否设置反了。

我把常见问题整理成一个速查表，方便你对照排查：

6. 关于这份公开图谱的几点个人体会与扩展思路

6.1 复现不能只跑代码，更要读代码

我见过不少人下载公开数据后，把作者代码从头到尾跑一遍，图出来了就觉得自己会了。但单细胞分析的可复制性本来就有限，服务器环境、R包版本、数据路径都影响结果。复现的真正价值在于读懂每一步为什么要这么做。比如为什么要用harmony而不是Seurat的IntegrateData，为什么质控阈值在视网膜上和外周血不一样，为什么注释要分两个层次而不是一步到位。这些思考才是在你以后自己处理课题数据时真正能带走的。

我自己的习惯是拿到公开代码后先看README或代码注释，理清每个脚本的输入输出关系，再在代码里加自己的注释，改成适合本地路径的版本。这样一套流程下来，既复现了论文结果，也顺手搭好了一套以后可以直接套用的分析模板。

6.2 这份图谱还能做什么扩展

复现基础分析只是起点。这份数据真正值钱的地方在于后续可以衍生出很多新分析，我简要列几个方向：

一是跨数据集整合。把这份中国人群图谱和之前发布的欧美人群视网膜衰老数据整合，比较不同人种间的衰老共性变化和特异变化。这类分析可以直接写成一篇方法学或比较研究型论文。

二是衰老模型构建。用不同细胞类型的衰老特征基因做打分模型，在独立数据集上验证评分是否能区分年龄组，再和疾病状态做关联分析。

三是细胞通讯网络的时间变化。用CellChat比较年轻和老年组的配体受体网络，重点找“指向某个关键靶点”的通路，比如炎症通路、补体通路。这些通路往往就是后续湿实验验证的候选机制。

四是把图谱作为参考做反卷积。对大量bulk转录组数据做去卷积，估计样本中细胞类型比例的变化，看哪些组织性疾病样本的细胞组成向“衰老模式”偏移。这个方向对临床转化的价值很高。

6.3 最后给复现者的一句话

说实话，我踩过最深的坑，是不看数据规模就盲目开跑，结果内存爆掉、代码中断、注释返工，来回折腾一两周。先看数据说明、先确认矩阵格式、先把QC阈值画出来再过滤、每跑完一步就存RDS，这四件事做好，复现过程会顺畅太多。公开图谱的价值不仅在于结论本身，更在于它给了大家一套可以对照练习的标准流程。认认真真从头到尾复现一遍，你收获的不仅是一套图，而是整个单细胞转录组分析框架的完整认知。这套认知会在你以后处理任何组织来源的单细胞数据时反复用到。