NCBI 下载 10X 单细胞数据到 Cell Ranger 全流程

做单细胞的人绕不开两件事一是在 NCBI 上把别人发表过的 10X 单细胞测序原始数据扒下来二是把它喂进 Cell Ranger 跑出表达矩阵。论文里那句轻描淡写的raw data are available in GEO under accession GSEXXXXXX背后其实是一整套从 SRA 到 FASTQ、从文件命名到样本表对齐的体力活。我前前后后处理过三十多个公开数据集从人的 PBMC 到小鼠的脑组织踩过的坑包括 barcode 被当成普通 read 丢掉、文件名让 Cell Ranger 直接报 sample 找不到、下到一半磁盘爆掉、跑了两天的结果质控一塌糊涂。这篇就把 NCBI 下载 10X 单细胞测序原始数据到 Cell Ranger 分析这条完整链路拆开来聊从检索、下载、格式转换、输入准备一直到运行和排错尽量把每一步的为什么也讲清楚。不管你是刚进实验室的研一新生还是已经会跑 Cell Ranger 但每次下数据都手忙脚乱的老手这套流程应该都能直接用。1. 这条链路到底长什么样先拆解再动手很多人第一次从 NCBI 下 10X 数据是直接在 GEO 页面找Download按钮结果点进去发现全是 .sra、.bam、.tar 之类的文件根本不知道怎么变成 Cell Ranger 要的 fastq 目录。问题的根源在于GEO、SRA 和 Cell Ranger 三者对原始数据的理解完全不一样。1.1 为什么 NCBI 下载这件事比想象中麻烦GEO 是一个数据集展示层它关心的是样本分组、实验设计、表型信息SRA 是一个原始序列仓库它只负责把测序仪吐出来的碱基流保存成一种紧凑的二进制格式而 Cell Ranger 只认一种东西命名符合规范、read 结构正确的 gzip 压缩 fastq。这三个体系之间的转换没有任何一个按钮能一键完成。更麻烦的是 10X 数据本身的特殊性。常规的 WGS 或 RNA-seq 是双端测序R1 和 R2 各自就是一条真实的序列片段但 10X Chromium 的 read 结构是错位的——R1 装的是 16bp barcode 加 10bp UMIv2 版本 26bpv3 版本 28bpR2 才是真正的 cDNA 序列中间还可能夹着一条 I1 的 sample index。如果下载工具按常规的双端逻辑处理极容易把 barcode 那条 read 直接扔掉。我见过最典型的翻车场景是这样的用fastq-dump不带任何参数跑出来两个文件扔进 Cell Ranger 报错Read 1 is too short因为下载工具默认只取了 biological reads把 barcode 所在的 index read 当成了技术性 read 过滤掉了。搞清楚这一点后面所有操作就都有了依据。1.2 三种常见的数据交付形态在 GEO 上翻数据集的时候原始数据一般以三种形态存在处理难度差别很大形态典型文件处理难度推荐做法交付 processed 矩阵.h5、.mtx、.tsv最低直接用 Seurat/Scanpy 读不必跑 Cell Ranger交付 bam 文件possorted_genome_bam.bam中等用cellranger bamtofastq反推 fastq只挂 SRA 号SRR/SRP/PRJNA最高SRA Toolkit 或 ENA 下载第一种情况最省事很多 2019 年之后发表的文章会把 filtered_feature_bc_matrix.h5 直接放在 Supplementary Files 里。但注意这类矩阵往往已经被作者做过一轮过滤你想复现原文的 Cell Ranger 参数或者重新做 QC 就无能为力了。真正要做方法学复现、要自己调--expect-cells这类参数的还是得拿原始 fastq。1.3 完整链路与耗时预估把整条链路铺开大致是这么几步在 GEO 定位到研究数据集顺着链接找到对应的 SRA/BioProject accession用 SRA Toolkit 或者 ENA 下载成 fastq检查 read 长度确认 barcode 在不在按 Cell Ranger 规则重命名准备参考基因组最后cellranger count跑完拿矩阵。时间上一个常规的 10X 3 样本约 4 亿 reads下载取决于网络和源站ENA 走直连大概半小时到两小时fasterq-dump 转换比下载更慢磁盘 IO 是瓶颈20 到 60 分钟Cell Ranger count 在 16 核 64G 的机器上大约 2 到 5 小时。也就是说单个样本从零到矩阵留出 8 小时的窗口比较稳妥。批量处理 10 个样本的话建议头一天晚上挂上下载第二天白天跑转换和比对排期才不至于太紧张。2. 在 NCBI 上精准定位可用的 10X 数据集下载慢一点还能忍最亏的是辛辛苦苦下完才发现数据集根本不适合自己的问题。检索阶段多花二十分钟能省掉后面一整天的返工。2.1 GEO 检索从关键词到 accession 号打开 GEO DataSets 的检索框不要只输入研究主题要把平台特征一起加上。我常用的组合是(10x OR Chromium OR single cell RNA) AND 你的组织或疾病关键词。括号和 OR 的写法在 GEO 里是支持的能比空格分词精准不少。拿到候选 GSE 号之后进详情页重点看三处顶部的 Platform 字段、Samples 列表里的 Organism 和 Source Name、底部的 Series 关联的 BioProject。Platform 如果是 GPL24676 或者 GPL18573 这类带 10x 字样的基本可以确定是 Chromium 3 或 5 试剂盒如果 Platform 描述里出现10x Genomics Visium那是空间转录组跟单细胞完全是两套流程别下错了。从 GSE 页面左下角的Relations里点进 BioProject再到 SRA Run Selector 页面才能看到真正的测序数据。SRA Run Selector 是我最推荐的中转站它把每个 run 的样本名、reads 数、碱基量、文件大小都列成一张可下载的表格还能直接勾选导出 Accession List 和 Metadata 两份 txt。2.2 判断数据集是不是 Chromium 平台这里有个血泪教训GEO 上标的single cell不一定就是 10X早期用 Smart-seq2、CEL-Seq、Drop-seq 平台做的也是单细胞这些数据的 read 结构和 10X 完全不同喂给 Cell Ranger 会直接崩。判断方法有三条看 README 或者 Methods 部分有没有出现Chromium Controller、10x Genomics字样看 SRA 里 run 的 read 数量10X 一个 run 动辄上亿Smart-seq2 单细胞通常只有几百万看 read 长度10X 的 R1 通常是 26bp 或 28bp这个长度非常具有辨识度。如果三条里两条对不上先别急着下用fastq-dump -X 10 --split-files抓前 10 条序列看一眼长度几秒钟的事比下完 40G 再发现不对划算太多。2.3 样本数与元信息核对清单决定开下之前我会在表格里核对这几项run 总数和样本数是否一致有些数据集一个样本拆成了多个 lane需要合并每个 run 的 spots 数也就是细胞数有没有明显偏小低于 500 个细胞的样本后续分析价值有限物种和参考基因组版本人、鼠、还是混合物种是否用了 GRCh38 还是 hg19组织处理方式是新鲜样本还是冻存是细胞悬液还是细胞核后者在做 QC 时阈值要放宽是否有配套的 Cell Ranger 输出文件有的话可以用来做下游结果的交叉验证。把这些信息整理成一张表后面构建样本表和排查问题时随时能翻比反复回 GEO 页面找要快得多。3. 下载环节SRA Toolkit 与 ENA 两条路怎么选真正开始下载其实有两条成熟路线一是用 SRA Toolkit 从 NCBI 官方源拉二是绕道 ENA 直接下 fastq。两条路我都长期在用各有适用场景。3.1 SRA Toolkit 安装与初始配置首推 conda 安装避免自己编译的麻烦conda create -n sra python3.10 conda activate sra conda install -c bioconda sra-tools pigz装完之后立刻做一件事把默认的缓存目录改掉。SRA Toolkit 默认把数据放在~/ncbi/public/sra/这个目录直接吃 home 分区空间家里分区一满整个系统都可能出问题。在~/.ncbi/user-settings.mkfg里加两行echo /repository/user/main/public/root /data/sra_cache ~/.ncbi/user-settings.mkfg echo /repository/user/main/public/cache-enabled true ~/.ncbi/user-settings.mkfg另外建议把prefetch的最大文件限制调大否则遇到大样本会直接拒绝下载vdb-config --set /sratoolkit/refseq/download-max-size 100000000000这些配置看着琐碎但每一条都是我实际踩坑之后补上的——缓存目录设在 home 分区被撑爆、大文件被默认限制卡住都是新手最常见的两种翻车方式。3.2 prefetch 与 fasterq-dump 的标准组合正式下载分两步走先用prefetch把 .sra 二进制文件拉到本地再用fasterq-dump转换成 fastq。之所以不一步到位是因为 prefetch 支持断点续传fasterq-dump 不支持网络一抖动前功尽弃。# 第一步下载 sra 文件支持断点续传 prefetch SRR1234567 -O /data/sra_cache批量下载的话可以把 accession 写进 txt 每个一行然后cat srr_list.txt | while read id; do prefetch $id -O /data/sra_cache; done第二步转换关键是加对参数fasterq-dump /data/sra_cache/SRR1234567.sra \ --split-files \ --include-technical \ --threads 8 \ --temp /data/tmp \ -O /data/fastq/SRR1234567--split-files保证 paired reads 被拆成独立的文件这是 Cell Ranger 的硬性要求--include-technical是 10X 数据的救命参数加了它才会把 index read 一起输出。如果这个 run 确实没有保存 index加了也不会有副作用只是不会多出第三个文件而已。--temp指定临时目录也很重要fasterq-dump 生成过程中临时文件体积能到最终结果的近两倍默认走 /tmp 很容易把系统盘塞满。转换完的 fastq 默认不压缩逐个用 pigz 压一下pigz -p 8 /data/fastq/SRR1234567/*.fastqpigz 是多线程版的 gzip同样大小的文件压缩速度比单线程 gzip 快好几倍这个细节在大批量处理时省下的时间相当可观。3.3 ENA 直连下载更省事的一条捷径如果只是想拿到 fastq其实 ENA欧洲核酸档案库比 NCBI 贴心很多它把 SRA 里的数据预先转成了 fastq.gz还拆好了文件、起好了名字。拿到 BioProject 号之后先拉一份文件清单curl https://www.ebi.ac.uk/ena/portal/api/filereport?accessionPRJNA123456resultread_runfieldsrun_accession,fastq_ftp,fastq_bytes,read_countformattsv ena_manifest.tsv打开这个 tsvfastq_ftp一列就是下载地址。批量下载很简单tail -n 2 ena_manifest.tsv | cut -f2 | tr ; \n | sed s|^|https://| urls.txt aria2c -i urls.txt -j 4 -x 8 -d /data/fastq/enaENA 的优势非常明显免去了 fasterq-dump 那一步源站带宽也普遍比 NCBI 更友好实测同样的数据集速度能差三到五倍。但它也不是万能——有些提交者上传的 SRA 数据经过了额外处理ENA 上未必 100% 还原原始 read 结构所以下完仍然要做第 4 章的长度检查。我的习惯是ENA 有就直接用 ENAENA 缺数据再回 SRA Toolkit。4. 整理成 Cell Ranger 能吃的输入下载完成只是拿到了原材料真正决定成败的是把它们整理成 Cell Ranger 认识的样子。这一步不出错后面的比对基本就稳了。4.1 先看清 read 结构barcode 到底在哪个文件里拿到一堆 fastq 之后第一件事是抽前几条序列看长度head -n 8 /data/fastq/SRR1234567/SRR1234567_1.fastq head -n 8 /data/fastq/SRR1234567/SRR1234567_2.fastq对照下表判断平台版本R1 长度R2 长度是否有 I1Chromium 3 v22698可选Chromium 3 v32891 或 150可选Chromium 5 v1/v22698可选Chromium 3 v42890可选如果 _1 文件里长度是 26 或 28说明 barcode 结构完好如果 _1 的长度接近 100 甚至更长那大概率是提交者把 barcode 拼到了 cDNA 前面或者数据被重新处理过这时候需要跟原始发表文章核实必要时联系作者要原始数据。如果发现根本没出 _1 文件只有 _2那说明--include-technical没加对或者这个 SRA 条目本身没保留 index。还有一个隐形的坑有些数据集的 _1 长度正确但用错了 read 方向测序时把 cDNA 放在了 _1 里。这种情况很少但存在判断方法是看 _1 的碱基分布barcode 那部分因为要落入 whitelist四种碱基的分布会比较均匀而 cDNA 有明显偏向性。4.2 文件重命名规则与 sample 参数对齐Cell Ranger 对文件名的要求非常死板必须符合这个模式[SAMPLE_NAME]_S[NUM]_L[LANE]_R[1-2]_[NUM].fastq.gz从 SRA 出来的是SRR1234567_1.fastq.gz和SRR1234567_2.fastq.gz完全不匹配。重命名示例cd /data/fastq/SRR1234567 mv SRR1234567_1.fastq.gz Pbmc_Donor1_S1_L001_R1_001.fastq.gz mv SRR1234567_2.fastq.gz Pbmc_Donor1_S1_L001_R2_001.fastq.gz这里的Pbmc_Donor1就是后面--sample参数要填的值必须完全一致。我强烈建议用见名知意的样本名而不是 SRR 号因为后面出图的时候轴标签直接就是这个名字用 SRR 号看着非常难受。如果同一个样本被拆成多个 run比如 SRR111111、SRR222222 属于同一个样本重命名时把它们的 sample 部分写成一致S 和 L 编号不同即可Pbmc_Donor1_S1_L001_R1_001.fastq.gz Pbmc_Donor1_S1_L002_R1_001.fastq.gzCell Ranger 会自动识别属于同一个 sample 的多个 lane 并合并处理不用手动 cat。4.3 参考基因组下载与版本选择Cell Ranger 需要 10x 官方预构建的参考基因组不是随便一个 FASTA 加 GTF。下载地址在 10x Genomics 支持页面的Reference Genomes部分常见的两类构建refdata-gex-GRCh38-2020-A人类2020 版最通用refdata-gex-mm10-2020-A小鼠对应 GRCm38/mm10。下载的是一个约 11G 的 tar.gzwget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz tar -xzvf refdata-gex-GRCh38-2020-A.tar.gz -C /data/refs/版本选择上2020-A 是经过最多验证的版本除非原始文章明确说自己用了旧版否则优先选它。如果用 GRCh38-2020-A 而原文用的是 hg19 或 GRCh38-1.2.0基因注释和基因名会有细微差别做出来的下游结果和原文对不上号不能算 bug只能说是版本差异这一点写方法学的时候要说清楚。5. Cell Ranger count 实操与资源规划输入都备好之后就进入真正跑分析的环节。这个环节最贵的是时间参数设错一次可能白跑半天所以值得把每个参数都过一遍。5.1 命令逐参数拆解cellranger count \ --idPbmc_Donor1_run1 \ --transcriptome/data/refs/refdata-gex-GRCh38-2020-A \ --fastqs/data/fastq \ --samplePbmc_Donor1 \ --expect-cells5000 \ --localcores16 \ --localmem64 \ --create-bamtrue每个参数的作用和踩坑点--id输出目录名如果目录已存在会直接报错想重跑要加--force或者换个 id别傻乎乎删了又重跑。--fastqsfastq 所在目录Cell Ranger 会自动递归查找所以目录层级不必完全精确但样本名匹配必须精确。--sample这是最容易出错的地方它必须和文件名_S1_前面的部分完全一致。大小写、下划线数量都要对上差一个字符就会报no input FASTQs were found。--expect-cells告诉算法预期细胞数。10X 官方的自动估计在细胞数差异大的样本上不太靠谱如果已知样本上机时目标细胞数直接填比如说 5000 或 10000如果不知道可以先不填跑一次看估计值再补跑。--localcores和--localmem必须显式写否则 Cell Ranger 会按机器全部核数去要资源在共享服务器上很容易把别人挤下去也容易被系统 OOM killer 干掉。--create-bam默认 true输出 bam 文件占空间很大如果只想要矩阵可以设成 false 省空间。5.2 时间与内存的实际估算Cell Ranger count 是内存密集型内存需求的经验公式大致是每 1000 个细胞配 1G 内存但对于 reads 数非常高的样本这个公式会低估。实际经验值细胞数reads 数推荐内存16 核耗时30001 亿32G1.5-2 小时50002 亿48G2-3 小时100004 亿64G3-5 小时200008 亿96-128G6-10 小时内存不够的典型症状是任务在中途毫无征兆地被 kill日志里只有一行 Killed。这种情况下先看系统 OOM 日志然后按上表加内存再跑。共享服务器上跑之前最好先free -h看一眼可用内存以及用squeue或htop确认没人跟你抢。5.3 输出目录与质控指标怎么看跑完之后输出目录里重要的东西有这些outs/web_summary.html可视化报告第一个该看的东西outs/metrics_summary.csv关键质控指标的表格版outs/filtered_feature_bc_matrix/过滤后的表达矩阵下游分析入口outs/raw_feature_bc_matrix/原始矩阵包含空液滴做 SoupX 之类的背景校正会用到outs/possorted_genome_bam.bam排序后的比对结果占空间最多。web_summary 里最该盯的几个数字和大致阈值指标理想值说明Estimated Number of Cells与预期接近差一个数量级要警惕Median Genes per Cell 500低于 200 一般说明样本质量差Median UMI per Cell 1000低于 500 得上游找原因Reads Mapped Confidently to Transcriptome 60%低于 50% 检查参考版本Fraction Reads in Cells 70%偏低说明空液滴占比高Q30 Bases in Barcode 85%偏低是测序问题得换数据实际评估的时候我不会只看单项而是把几个指标一起看。比如同时出现细胞数正常但 median genes 只有 300和Fraction Reads in Cells 只有 50%那大概率是样本本身 RNA 降解严重得回原文确认他们用的是什么处理方法是不是冷冻组织或者细胞核样本后者的中位基因数本来就会偏低一些。6. 报错排查与实操心得跑流程的过程中报错是常态。下面把遇到的典型问题和处理方式整理出来遇到问题直接查表能省不少时间。6.1 下载与转换阶段的常见问题问题一prefetch 报maximum file size exceeded。原因是 SRA Toolkit 默认的文件大小限制太小解决方案是在配置里把download-max-size调大具体命令前面 3.1 节已经给了。问题二fasterq-dump 跑到一半提示磁盘空间不足。fasterq-dump 在转换过程中会产生几乎和最终结果等体积的临时文件然后把它们重新拼装。实际预留的空间要比最终结果大两倍以上。--temp参数一定要显式指定到空间充足的分区。问题三转换完发现只有一个 fastq 文件。说明这个 run 可能本身就是单端或者--split-files没加或者数据在 SRA 里的元信息把两条 read 标记成了同一条 spot。先加--split-spot试试实在不行用fastq-dump老命令加--split-files再跑一次两个工具的内部逻辑不完全一样。问题四ENA 下载下来的 fastq.gz 解压后发现文件截断。多半是下载中断了aria2c支持断点续传重跑一次同样的命令会接着下。如果重跑后还是坏的用gzip -t校验一遍确认完整性再往下一步走。6.2 运行阶段的常见问题问题一The chemistry of this run was not detected。Cell Ranger 无法从 read 长度自动推断试剂盒版本。这通常发生在 read 长度不太标准的情况下。解决办法是手动指定--chemistrySC3Pv3之类的参数具体是 v2 还是 v3 靠 R1 是 26 还是 28 判断。问题二no input FASTQs were found for sample xxx。回到第 4.2 节检查三件事文件后缀是不是 .fastq.gz不能是 .fq.gz、文件名里的_R1__R2_有没有漏、--sample参数和文件名前缀是否完全一致。我遇到过一次是文件名里用了中文下划线肉眼看不出来折腾了一小时才发现。问题三Cell Ranger 跑到比对阶段内存被 kill。如果是集群环境先看 SLURM 或 PBS 的资源申请是否设置对如果是本地机器看dmesg | grep -i killed确认是不是 OOM。加内存或者先用--subsample参数小规模跑一遍验证流程。问题四web_summary 里 Estimated Number of Cells 只有几百但实际材料明明有几千个细胞。两个常见原因一是--expect-cells没有设置或者默认值太低导致算法偏向保守二是原始数据质量问题RNA 含量低导致细胞识别不出来。可以设一个明确的--expect-cells重跑一次看看结果是否变化。6.3 几张私藏的经验卡片用了这么多年攒下这么几条不太写在文档里的心得都用得上。卡片一永远保留一套原始 SRA 文件。fasterq-dump 转换一次可能要几小时如果后面发现参数错了要重新转换没有原始 .sra 就得重新下载。我习惯在 /data/sra_cache 里留一份 .sra等所有分析都确认没问题了再统一清理。卡片二样本名里不要出现空格和特殊字符。Cell Ranger 对文件名里的特殊字符非常敏感用下划线连接是最安全的。样本名里的群体、处理批次信息用Ctrl_、Treat_这种前缀别用C-1、T#2这种符号。卡片三md5 校验一定要做。大文件下载完之后跑一次md5sum -c花不了几分钟但能避免后面在错误数据上浪费一整天。尤其是走 ENA 或者第三方镜像的时候文件缺失或者截断的概率并不低。卡片四下游分析之前先做一次聚类目视检查。拿到 filtered_feature_bc_matrix 之后别急着做差异分析先跑个 PCA 和 UMAP 看一眼有没有明显异常样本。有些时候 Cell Ranger 的指标看起来正常但实际细胞群的分布非常奇怪这时候就该回上游重新审视参数设置。卡片五方法学描述要写清楚版本号。NCBI、SRA Toolkit、Cell Ranger、参考基因组这四样东西的版本号都要记下来投稿的时候审稿人很可能会问。Cell Ranger 在 6.0 之后对 UMI 的处理有一处调整不同版本跑出来的矩阵在个别基因上会有差异写清楚版本号既是对读者负责也是避免日后自己都说不清楚。最后再补一个我用了很多年的小技巧把整个流程写成一个 shell 脚本把 accession 号做成参数传进去从 prefetch 到 cellranger count 一条龙跑完。第一个样本手工走通之后后面几十个样本照着脚本批量跑犯错的概率会大幅下降也方便日后复现。脚本里每跑完一步就写一行日志到文件出了问题回头翻日志定位比盯着终端输出回溯要高效得多。