这两年眼科和衰老领域的单细胞研究越来越密集,但真正把中国人群视网膜衰老单细胞图谱做出来、还把数据和代码全部公开的,确实不多。这篇文章我要聊的就是2026年刚发布的这张图谱——从数据下载、环境配置到核心分析流程的完整复现记录,包括我在实际操作中踩过的坑和排查思路。如果你是刚接触单细胞数据分析的生信入门者,或者想用公共数据做视网膜衰老方向的研究,这篇指南可以帮你少走不少弯路。
我拿到这个项目的第一反应是:图谱本身的价值不用多说,但更难得的是作者把全流程代码和数据一起放出来了,这意味着我们可以不依赖论文补充材料里那些零散的描述,直接把整个分析链路在本地跑通。这篇文章会围绕三个核心问题展开:这张图谱解决了什么问题、数据长什么样、如何从原始数据一步步复现出论文里的细胞分群和衰老特征。
1. 这张中国人群视网膜衰老单细胞图谱的特别之处
1.1 为什么盯着视网膜衰老做单细胞图谱
视网膜作为中枢神经系统的延伸部分,一直是衰老研究的焦点。它的神经细胞不可再生,一旦损伤就很难修复,而年龄相关的视网膜病变——比如老年黄斑变性、糖尿病视网膜病变、青光眼——本质上都和细胞衰老脱不开干系。过去几年里,小鼠视网膜的单细胞图谱已经有好几套了,但人和小鼠在免疫细胞组成、Müller胶质细胞反应模式、感光细胞基因表达上差异非常大,直接用小鼠结果推算人类病理机制,经常对不上号。
这套图谱填补的正是这个空白:它把中国人群不同年龄阶段的视网膜组织做了单细胞转录组测序,覆盖了从胎儿期到老年期的多个时间点,构建了一套完整的视网膜细胞类型老化轨迹。中国人群的数据尤其重要,因为不同族裔之间的遗传背景、环境暴露、饮食结构差异很大,视网膜衰老的表型特征可能存在群体特异性。这不仅是数据层面的补充,更是对以往以欧美人群为主的视网膜单细胞图谱的一个重要对照。
1.2 图谱提供的数据形态与代码结构
作者公开的内容分为两大部分。第一部分是原始数据,通常以filtered_feature_bc_matrix格式存放,也就是10X Genomics平台标准的三个文件:barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz。每个样本对应一个目录,按年龄段和供体编号命名。第二部分是分析代码,按步骤编号组织,从质控过滤到拟时序分析一个不少。
我建议你拿到数据后先不要急着跑代码,而是把目录结构完整看一遍。很多复现失败的问题,根源不在代码本身,而是数据路径和文件命名和你预期的不一致。下载清单里一般还会附带一个metadata.csv,记录了每个样本的年龄、性别、组织来源、测序深度等关键信息,后续做整合分析、批次效应校正时都要反复用到这张表。
1.3 复现这套流程需要什么基础
如果你想完整复现这套分析,需要具备的硬件和软件条件如下:一台内存不低于32GB的工作站(后面我会解释为什么这个要求不过分)、R 4.3以上版本、Python 3.9以上环境。核心R包包括Seurat、SingleR、monocle3、CellChat、SCENIC,Python端主要用到scanpy和scVelo做RNA速率分析。
如果你是第一次接触单细胞数据的生信新手,我的建议是先不要急着全流程复现,而是先跑通我最前面给出来的核心流程——质控、降维聚类、细胞注释。这三个步骤跑通了,整个分析框架基本就掌握了。
2. 数据下载与项目准备:别在第一步就翻车
2.1 公共数据库的检索策略
作者把数据提交到了公共数据库,下载方式通常以GEO编号或GSA编号的形式在论文的Data Availability部分标注。检索的时候我建议不要在数据库首页直接搜关键词,而是先去论文里找到完整的accession编号。单细胞数据经常按照GEO系列编号(GSE开头)访问,但实际测序文件可能存放在同一个编号下的多个子平台里,点进去之后要仔细核对样本数是否和论文一致。
下载前还要确认数据版本。有些作者会在修订期间更新数据,你下载到手的可能是v2版本,但论文正文和代码注释里写的是v1的统计口径,这会造成后续分析的数字对不上。下载后第一时间比对文件的总大小和GEO页面标注的大小,不一致就重新下载,不要抱着侥幸心理继续往下走。
2.2 本地目录结构的组织方式
下载完成后,强烈建议按照下面的结构整理数据:
retina_aging/ ├── data/ │ ├── raw/ │ │ ├── sample01/ │ │ │ ├── barcodes.tsv.gz │ │ │ ├── features.tsv.gz │ │ │ ├── matrix.mtx.gz │ │ └── ... │ └── metadata.csv ├── scripts/ │ ├── 01_qc_filtering.R │ ├── 02_normalization_pca.R │ ├── 03_clustering_annotation.R │ └── 04_de_aging_signature.R └── results/这个结构的好处是路径清晰、脚本可复用。特别是scripts目录,后续你要改参数、换样本集重新跑,只需要修改脚本里对应的文件路径,不需要在RStudio里反复手敲路径。
2.3 环境配置与依赖版本锁定
单细胞分析最头疼的问题就是包版本冲突。Seurat从v4升级到v5,很多函数的默认参数都变了,FindAllMarkers的返回结果格式也变了。我的做法是用renv锁定版本:
install.packages("renv") renv::init() renv::install("Seurat", version = "5.0.1") renv::install("SingleR", version = "2.4.1")renv会把当前项目的包版本记录到一个lock文件中,别人拿到之后直接renv::restore()就能还原出完全一致的环境。这个问题在代码复现中非常关键——大量"代码跑不通"的反馈,最后发现都是包版本不一致导致的,而不是代码本身有问题。
Python环境我用conda单独建了一个虚拟环境,避免和系统Python冲突:
conda create -n retina_aging python=3.9 conda activate retina_aging pip install scanpy==1.9.3 scvelo==0.2.5 cellchat==0.2.0这一步看起来是准备工作,实际是整个复现过程中最值得投入时间的环节。环境稳定了,后面所有分析都会顺畅很多。
3. 核心分析流程:从count矩阵到细胞注释的完整实操
3.1 数据读入与Suerat对象构建
数据下载好、环境配好之后,第一步就是把10X标准输出读入Seurat。这里有个细节要注意:如果下载的是sample01目录下的三个文件,直接用Read10X函数指定目录即可,不需要手动拼接:
library(Seurat) # 读取单个样本 data_dir <- "data/raw/sample01" counts <- Read10X(data.dir = data_dir) # 创建Seurat对象,project命名尽量带上样本ID,方便后续merge时追踪来源 obj <- CreateSeuratObject( counts = counts, project = "Retina_sample01", min.cells = 3, # 至少在3个细胞中表达的基因才保留 min.features = 200 # 至少表达200个基因的细胞才保留 )min.cells = 3和min.features = 200这两个参数对单细胞来说是经验值。min.cells过滤掉的是那些只在极少数细胞中出现、大概率是环境RNA污染或测序错误导致的低质量基因;min.features过滤掉的则是测序深度太低、根本无法支撑下游分析的细胞。不同组织类型的过滤阈值其实不一样,上皮组织细胞基因数通常比较高,而血液细胞相对低。视网膜组织里感光细胞的线粒体基因比例普遍偏高,因为感光细胞的外段富含线粒体,这一点和普通组织的质控标准要区分开。
3.2 质控参数的选择逻辑
质控是单细胞分析中最主观、也最容易影响最终细胞注释结果的一步。常见的标准是每个细胞检测到的基因数在200到6000之间,线粒体基因比例低于20%,但视网膜组织有自己的特殊性。
我在复现时发现,感光细胞和视网膜色素上皮细胞的线粒体基因比例天然比其他细胞高。如果你一刀切用10%的线粒体阈值,会把大量真实的感光细胞当作低质量细胞过滤掉,导致后续分析里感光细胞亚群缺失,这会直接毁掉整个图谱的完整性。我的建议是先用可视化的方式比较不同阈值下的细胞分布,看看你的数据到底应该卡在哪个位置:
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-") VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)如果看到感光细胞特异性marker基因(比如RHO、ARR3)在某个线粒体阈值区间内仍然大量表达,说明这部分细胞是有生物学意义的,不应该被过滤。这一步需要一定的生物学判断力,不能完全依赖默认参数自动跑。
质控完成后,把所有样本合并成一个Seurat对象。合并前要注意,不同样本的细胞编号会冲突,Seurat会自动加上前缀。你需要在合并后手动检查一下每个样本的细胞数是否和原始数据一致,防止漏样本:
merged_obj <- merge(obj1, y = c(obj2, obj3, obj4), project = "Retina_Aging_Full") table(merged_obj$orig.ident)3.3 标准化与降维聚类:分辨率参数怎么选
合并之后,先做标准化和方差基因筛选,然后进入PCA降维。这里我不展开原理,只强调几个实操参数:
merged_obj <- NormalizeData(merged_obj, normalization.method = "LogNormalize", scale.factor = 10000) merged_obj <- FindVariableFeatures(merged_obj, selection.method = "vst", nfeatures = 2000) merged_obj <- ScaleData(merged_obj) merged_obj <- RunPCA(merged_obj, npcs = 30)npcs的选择可以通过ElbowPlot判断,一般选到拐点附近的主成分数量。但注意,视网膜组织的细胞类型异质性很大,感光细胞和胶质细胞之间的转录差异本来就非常明显,所以主成分的拐点可能不如其他组织那么陡峭。我这次用了20个PC,跑出来的UMAP分群效果很清楚。
聚类阶段最重要的是resolution参数。resolution越低,分群数越少;越高则分得越细。视网膜组织里Müller胶质细胞、星形胶质细胞、小胶质细胞之间边界比较清晰,但感光细胞里的视锥和视杆亚型之间差异很微妙,需要逐步提高分辨率来区分。我建议从resolution = 0.5开始跑,先看全局结构,再针对具体感兴趣的亚群单独做一次子聚类,而不是一开始就用高分辨率。resolution的调整会直接影响后续注释的颗粒度。
merged_obj <- FindNeighbors(merged_obj, dims = 1:20) merged_obj <- FindClusters(merged_obj, resolution = 0.5) merged_obj <- RunUMAP(merged_obj, dims = 1:20)这里有个非常实用的经验:每调整一次resolution,不要只看UMAP图上的分群数量,要结合marker基因的分布去判断新分出来的群是否真的有生物学意义。如果某个亚群只是把同一类细胞按表达量高低切开了,而没有任何独特的marker基因,就说明分辨率拉过头了。
3.4 细胞注释:marker基因与SingleR联合判定
细胞注释是整个流程中最依赖经验的一步,也是决定复现质量的关键。视网膜组织的细胞类型相对固定,下表是复现这套图谱时最常用到的marker基因:
| 细胞类型 | 核心marker基因 |
|---|---|
| 视杆细胞 | RHO, NRL |
| 视锥细胞 | ARR3, OPN1SW, OPN1MW |
| 双极细胞 | VSX2, GRM6, OTX2 |
| 无长突细胞 | GAD1, SLC6A9, TFAP2A |
| 水平细胞 | ONECUT1, LHX1, CALB1 |
| Müller胶质细胞 | RLBP1, GLUL, VIM |
| 星形胶质细胞 | GFAP, S100B, AQP4 |
| 小胶质细胞 | C1QA, C1QB, P2RY12 |
| 血管内皮细胞 | PECAM1, CLDN5, FLT1 |
实际操作中,我用三种方法交叉验证。先用SingleR自动注释——它是基于参考转录组数据集,通过计算每个细胞与参考集已知细胞类型的相关性来打标签,非常快,但精度有限。然后手动检查关键marker基因在各个cluster中的表达情况。最后再把两组结果做对比,冲突明显的cluster单独回到UMAP图上看位置。
冲突最常出现在Müller胶质细胞和小胶质细胞之间。两者在应激状态下都会高表达一些炎症相关基因,容易被SingleR混淆。这时候就要靠形态学知识了——小胶质细胞的经典marker是P2RY12、TMEM119,而Müller胶质细胞有RLBP1。用FeaturePlot把两组marker同时画出来,基本一眼就能分辨。
3.5 衰老特征分析:差异表达与基因调控网络
细胞注释完成后,接下来就是整个图谱的核心卖点——衰老特征分析。作者并不只是简单地画了不同年龄细胞的UMAP分布图,而是把每个年龄段差异表达的基因做了系统比较,定位出随年龄增长在特定细胞类群中表达水平改变的老化关键基因。
我在复现时关注的是Müller胶质细胞亚群。这个细胞类型在衰老过程中扮演的角色很特别:它既负责维持视网膜结构,又在损伤后发生胶质活化反应。作者在各个年龄段之间做差异分析时,Müller胶质细胞差异基因数量明显增多,这个现象我在自己的聚类结果里也看到了。
# 差异分析示例:比较老年组与年轻组的Muller胶质细胞表达差异 markers_muller <- FindMarkers(merged_obj, ident.1 = "Elderly_Muller", ident.2 = "Young_Muller")完整的衰老特征分析还需要做基因集富集分析、转录因子活性分析,甚至细胞通讯分析。这些分析虽然代码量不大,但每一步的输出都需要结合生物学背景解读。我强烈建议你复现时不要只盯着最终的p值和log2FC,把每一层输出的结果都保存下来,后续写论文的时候这些都是重要的补充材料。
4. 复现过程中最容易踩的坑与排查思路
4.1 内存溢出:32GB内存为什么仍可能不够
这是所有做单细胞分析的人都会遇到的问题。视网膜样本虽然不像肿瘤样本那样动辄几十万个细胞,但全套样本合并之后,如果同时加载大量基因和细胞,R的内存占用很容易飙到30GB以上。
如果代码跑到ScaleData时内存溢出,不要急着加内存条,有几个更高效的处理方式:一是用Seurat的Plan设置多层磁盘存储;二是分批次处理,每次只对部分细胞做标准化和scale;三是降低nfeatures的筛选数量,比如从2000降到1500。实际测试下来,分层存储的效果最明显,运行时间没有显著增加,内存占用却下降了接近一半。
另一个容易忽略的问题是R版本和BLAS库的搭配。R 4.2以上版本在某些Linux发行版上默认使用多线程BLAS,看起来是好事,但内存占用会随线程数翻倍。你可以通过options(mc.cores = 1)限制并行线程数,或者安装RhpcBLASctl控制BLAS线程数。
4.2 批次效应:不同样本放在一起,UMAP图出现"样本聚类"而非"细胞聚类"
这是单细胞整合分析里最常见的问题。如果你在UMAP图上看见相同来源的样本各自聚成一团,而不是相同细胞类型的细胞跨样本聚在一起,说明批次效应没有处理好,而不是生物学差异真的那么大。
Seurat经典的整合流程是SCTransform配合FindIntegrationAnchors:
obj_list <- SplitObject(merged_obj, split.by = "orig.ident") obj_list <- lapply(obj_list, SCTransform) anchors <- FindIntegrationAnchors(object.list = obj_list, dims = 1:20) integrated <- IntegrateData(anchorset = anchors, dims = 1:20)整合之后,UMAP图明显正常多了,同类型细胞跨样本聚到一起。这个环节有个很关键的经验:FindIntegrationAnchors会在每个样本内部先做一次锚点筛选,如果某些样本的细胞组成差异过大(比如某个样本刚好只测到了感光细胞),锚点的数量就会不足,整合效果不理想。解决方法是适当调整k.anchor参数,或者考虑用harmony这种基于PCA嵌入的快速整合方法,跑出来的效果在视网膜组织数据上差异不大,但速度要快很多。
4.3 注释偏差:SingleR的参考数据集不适合视网膜组织
SingleR内置的参考集里,免疫细胞的参考数据非常丰富,但视网膜特异性细胞类型的数据覆盖不足。我跑下来发现,Müller胶质细胞和小胶质细胞这两个类型的注释稳定性最差。这是因为两者都高表达一些共同的应激反应基因和炎症相关基因,在转录组特征上的区分度不够高。
遇到这种情况,我的应对策略是建立一个本地参考集:从已发表的质量较高的视网膜单细胞图谱里提取各细胞类型的marker,作为自定义参考集喂给SingleR。这个操作并不复杂,只需要把你的marker基因列表整理成SingleR能接受的List格式,效果立竿见影。
4.4 代码版本差异:函数重名导致的连环报错
Seurat升级到v5之后,FindAllMarkers的返回对象结构和v4不再兼容,ECT里调用的SplitObject行为也有了变化。如果你的分析流程是网上搜来的代码片段,很可能混杂了几个不同版本时代的写法,跑起来之后报错的一个接一个,排查非常耗时。
最稳妥的办法是打开代码后先看头部导入的包版本号,确定整段代码是基于哪个版本写的。如果是v4写的流程,就在v4环境里重跑;如果只有v5环境,就要逐行检查API变化。我在复现过程中遇到过RunPCA之后fetch.data格式变化的问题,排查了一整个下午才发现是版本差异导致的。
4.5 数据下载不完整:矩阵维度对不上
matrix.mtx.gz下载不完全时,Seurat读取不会直接报错,而是会在后续分析中莫名出现无法转换稀疏矩阵的错误。排查方法是在读取数据之后手动检查矩阵维度:
dim(counts) # 预期输出应该是 基因数 x 细胞数,和你从GEO页面上看到的数字一致如果维度不对,优先重新下载,而不要去修改数据处理逻辑。数据源本身就是坏的,后续再补救都是徒劳。
5. 这张图谱还能怎么用
5.1 跨人群、跨物种比对
图谱发布后,最直接的应用方向就是和其他物种的视网膜单细胞数据做整合。比较不同物种的衰老特征,可以区分哪些衰老机制是哺乳动物共有的,哪些是人类特有的。这种跨物种分析对药物靶点的筛选尤其有价值——如果某个衰老相关基因在小鼠和人类视网膜里的变化方向一致,那它在临床转化中的可信度会更高。
跨人群比对方面,这套图谱给出了中国人群的基线数据,未来如果有其他族裔的视网膜衰老数据集发布,可以做跨族裔的差异分析。这是个典型的信号:这个方向的公共数据会越来越多,谁先跑通分析流程,后续出成果的概率就更大。
5.2 衰老相关疾病的风险基因定位
把图谱里的衰老特征基因和已有的GWAS数据结合起来,是一种很高效的疾病机制研究策略。比如已知某些位点和老年黄斑变性高度相关,你可以定位这些位点附近的基因,看它们在视网膜图谱的哪些细胞类型中高表达。如果某个风险基因恰好在衰老过程中表达变化显著的细胞类型里富集,那么这个基因的功能研究就有了明确的细胞学基础。
这种分析不需要重新做实验,完全是数据挖掘的工作,但前提是你对单细胞图谱的细胞注释质量有把握。所以前面反复强调注释环节要仔细,就是因为这个环节的质量会直接影响下游所有分析的可信度。
5.3 药物靶点与干预策略
从应用层面讲,图谱能帮我们找到针对特定细胞类型起作用的药物靶点。传统药物筛选往往只关注靶点在大组织中的表达水平,但单细胞数据告诉我们,视网膜里不同细胞类型的基因表达差异非常大。一个在Müller胶质细胞里高表达的靶点,如果在感光细胞里不表达,那么靶向这个位点的药物对感光细胞的影响就会比较小,副作用也可能更小。
一个简单的操作路径:直接在图谱数据里筛选在目标细胞类型中高表达、且在衰老过程中表达量显著改变的基因,再比对已知的药物靶点数据库。这样得到候选基因的临床转化依据比盲筛要扎实得多,这也是这类公共图谱数据未来大量被引用的原因。
6. 我的实操体会与扩展建议
整个复现流程走下来,我最大的体会是:这套图谱的价值不只在数据本身,更在于代码和分析逻辑的示范意义。作者把每个步骤的参数选择逻辑都解释得很清楚,特别是质控阈值和细胞注释这两个最容易主观化的环节,给出了可复现的筛选依据。这对整个视网膜衰老领域的基础建设都有推动作用。
最后再分享一个小技巧。我在复现时习惯把关键分析节点的可视化结果按步骤编号保存下来,这样不仅方便自查,而且给组里的学生做培训时特别有用——他们可以直接对照标准的中间结果判断自己跑到哪一步出了问题。如果你打算在这条数据上做深入分析,建议额外保存每个cluster的平均表达矩阵,后续做基因集评分、细胞状态转换分析时都会用得到。
对于想在这个方向做研究的读者,我的建议很直接:先把这套数据完整跑通一遍,然后挑一个你真正感兴趣的细胞类型,专门做它的亚群分析和衰老轨迹分析。单细胞图谱的数据只是起点,真正的科学发现往往藏在二次分析和深度挖掘里。