单细胞测序这两年最不缺的就是公开数据,GEO 上随便检索一个关键词,动辄就是几十个 GSE 数据集。但真正让刚入门的人卡住的,往往不是分析本身,而是第一步:把 NCBI 上的 10X 原始数据拿下来,整理成 Cell Ranger 能认的格式。我自己第一次做这件事的时候,光是"这个 SRR 号下下来怎么有三个 fastq 文件""为什么 --sample 写的名字和文件名对不上就报 No fastq files found"就折腾了整整一个下午。这篇文章就把从 NCBI 找数据、判断数据形态、用 SRA Toolkit 转 FASTQ、必要时用 bcl2fastq 拆 BCL,一直到 cellranger count 跑通、看 web_summary 判断数据质量这一整条链路讲清楚。适合完全没碰过 SRA 数据库的新手,也适合跑过几次但总在格式问题上翻车的同学。
1. 先搞清楚数据在 NCBI 的哪个角落
1.1 GEO 页面上三个入口的区别
很多人拿到一篇单细胞文章,第一反应是打开 GEO 页面,然后看到一大堆文件就懵了。其实 GEO 页面上跟"能不能跑 Cell Ranger"相关的入口就那么几个,先分清楚它们的性质,后面能省掉大量无用下载。
打开一个 GSE 页面往下拉,你会看到Supplementary file这一栏。这里面的东西分两种:一种是GSE123456_RAW.tar这种打包的原始文件,里面可能是每个样本的 BCL 压缩包,也可能是已经切好的 FASTQ;另一种是形如GSMxxxxxx_sample_matrix.mtx.gz、barcodes.tsv.gz、features.tsv.gz的表达矩阵文件。前者是原始数据,后者是已经处理过的结果。
再往下或者往页面顶部的链接区,会有一个SRA Run Selector或者BioProject的入口。点进 SRA Run Selector,你会看到一个表格,每一行是一个 SRR 号,对应一次测序运行。这才是最"干净"的原始数据来源。
区分方法很简单:矩阵文件只能做下游分析,跑不了 Cell Ranger。Cell Ranger 需要的是原始 reads,也就是 FASTQ 或 BCL。你要是拿 mtx 去喂 Cell Ranger,它会直接告诉你找不到 fastq 文件。所以第一步判断:这个数据集有没有提供 raw data。很多 2020 年以后的文章会同时提供两者,也有不少老数据集只给了矩阵。
1.2 从 GEO 跳到 SRA Run Selector 的正确姿势
GEO 和 SRA 的关系是"样本元数据"和"测序原始记录"的关系。同一个项目在 GEO 叫 GSE 号,在 SRA 叫 SRP 号,样本层面 GEO 是 GSM,SRA 是 SRS,一次测序运行是 SRR。
从 GSE 页面跳转的路径通常是这样的:在页面右侧的Relations区域找到SRA链接,点进去就是该项目的 SRA 页面。或者直接改 URL:把https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123456换成https://www.ncbi.nlm.nih.gov/Traces/study/?acc=GSE123456,SRA Run Selector 会自动把该项目下的所有 run 列出来。
这里有个经验:优先用 SRP 号而不是 GSE 号去 Run Selector 检索。因为一个 GSE 有时候会关联多个 SRP,用 GSE 检索可能漏掉一部分 run。在 Run Selector 页面顶部有个Study下拉框,能看到该项目关联的所有 SRA study。
Run Selector 页面上最有用的是那个可以筛选和导出的表格,列包括 Run、BioSample、Experiment、LibraryLayout、Bytes、Bases、spots 等。你可以把这些行勾选后点Metadata导出成 CSV,拿到本地做批量下载脚本。这一步比一个个点网页高效太多,尤其是面对几十个 run 的项目。
1.3 先说结论:哪些项目根本跑不了 Cell Ranger
踩过几次坑之后,我养成了一个习惯:在下载之前先花五分钟判断这个数据集到底值不值得下,因为一个 10X 样本的 FASTQ 动辄 20 到 50 GB,下错了纯属浪费时间和磁盘。
判断标准我总结成这么几条:
| 情况 | 能否跑 Cell Ranger | 说明 |
|---|---|---|
| SRA 里有 3' 或 5' 转录组 reads | 可以 | 最常见,直接 fasterq-dump |
| GEO 提供 RAW.tar 且内含 BCL | 可以 | 需要 bcl2fastq 或 mkfastq |
| 只提供 filtered_feature_bc_matrix | 不可以 | 已有细胞过滤,无法回溯 |
| 只提供 raw_feature_bc_matrix | 不可以 | 只有计数矩阵,没有 reads |
| 提供了 FASTQ 但被 trim 过 | 视情况 | 若 barcode 端被截短则无法用 |
| 空间转录组 Visium | 可以但流程不同 | 需要 spaceranger 不是 cellranger |
| 单细胞 ATAC | 可以但流程不同 | 需要 cellranger-atac |
还有一个细节:有些项目做的是单细胞核测序(snRNA-seq),这个用 Cell Ranger 跑完全没问题,参考基因组和参数都不变,只是在解读 web_summary 的时候要注意细胞数会偏低、线粒体基因比例会明显更低。
2. 拿到 SRR 号之后,先别急着下:判断数据形态
2.1 RAW.tar、SRA Run、现成 FASTQ 三种形态
NCBI 上的 10X 数据,落到你手里通常有三种形态,每种对应的后续操作完全不同。
第一种是GEO 的 RAW.tar。解压之后一般是每个样本一个GSMxxxxxx.tar,再解压里面是bcl/目录或者fastq/目录。如果是 BCL,你需要 bcl2fastq 或者 cellranger mkfastq;如果是 FASTQ,直接跳到命名整理那一步。
第二种是SRA Run。这是最普遍的形态,一个 SRR 号对应一次 sequencing run,用 SRA Toolkit 转成 FASTQ。要注意的是,SRA 里的 run 粒度不一定等于样本粒度——有的项目一个样本拆成多个 run(多 lane),有的项目多个样本混在一个 run 里(用 index 区分)。
第三种是作者直接提供的 FASTQ 附件。这种情况在近两年的数据集里越来越常见,因为测序成本的下降让作者更愿意直接上传大文件。这类文件通常已经在压缩包里按样本分好目录,但文件名往往被改得乱七八糟,需要重新整理。
三种形态的判断方法:下载前先读 GEO 页面上Supplementary file的文件名和大小。如果看到*_RAW.tar且几 GB 以上,多半是 BCL 或 FASTQ;如果看到*_fastq.tar.gz之类的,那就是现成的。
2.2 读 Run Browser 里的读长和 Layout
在 SRA Run Selector 里点任意一个 SRR 号,会进到 SRA Run Browser 页面。这个页面里有两个字段非常关键:Layout和Read specification。
Layout一般显示PAIRED,说明是双端测序。10X 的数据基本都是双端,Read1 是 cell barcode 加 UMI,Read2 是 cDNA 片段。
Read specification会给出每一端的长度,比如28,91,8,91或者151,151。这个数字序列的解读方式是:如果只有两段,比如28,91,那 Read1 是 28bp 的 barcode+UMI,Read2 是 91bp 的 cDNA。如果是四段,前两段是 Read1 和 Read2,后两段是 Index Read 的两端。
这个信息决定了你后面拿到 FASTQ 之后哪个文件对应 R1、哪个对应 R2。我遇到过不止一次,fasterq-dump 出来的_1.fastq长度是 91,_2.fastq是 28,跟 Cell Ranger 的预期正好反了。判断方法很简单,用zcat file.fastq.gz | head -2 | tail -1 | wc -c看一眼第一条 read 的长度就知道了。
顺便说一句,如果是 10X v3 化学,R1 通常是 28bp(16bp barcode + 12bp UMI),v2 是 26bp(16 + 10)。Read2 长度跟测序配置有关,常见 91 或 98。这些数字记住之后,看到 read 长度就能大致判断化学版本。
2.3 一个容易被忽略的坑:SRA Lite 与质量值
这是近两年新出现的坑,值得单独拎出来说。NCBI 从 2023 年前后开始把一部分旧的 SRA 数据转换成了 "SRA Lite" 格式。这种格式的特点是不存储原始的碱基质量值,下载出来转成 FASTQ 之后,质量值那一列是固定的占位字符。
对 Cell Ranger 来说,如果你的数据本身质量很好,这个影响可能不大,但严格来说它会丢失 Q30 之类的统计信息,web_summary 里的质量相关指标会失真。
判断方法是在 SRA Run Browser 页面看有没有SRA Lite的标记,或者在 Run Selector 表格的列里找SRA-Lite字段。如果中招了,可以尝试用prefetch --type all拉取完整数据,或者去其他镜像数据库找同一个 run 的原始版本。我个人的做法是:如果这是个关键样本,宁可换个数据源,也不要拿一份没有质量值的 FASTQ 硬跑。
3. SRA Toolkit 把 SRA 变成 Cell Ranger 认的 FASTQ
3.1 prefetch 与 fasterq-dump 的分工
SRA Toolkit 装好之后,你会看到两个名字像双胞胎的命令:prefetch和fastq-dump(以及新版的fasterq-dump)。很多人搞不清该用哪个,其实它们的分工非常清楚。
prefetch负责下载。它把 SRA 文件从远端拉到本地缓存目录,默认在~/ncbi/public/sra/下面。这个命令最大的好处是支持断点续传,网络抖一下断了再跑一次,它会从断掉的地方继续,不会从头再来。
fasterq-dump负责转换。它把本地的.sra文件解压、拆分成 FASTQ。之所以叫 "faster",是因为它比老的fastq-dump快好几倍,支持多线程,是现在的主流选择。
最小可用的命令组合长这样:
# 下载 prefetch --max-size 100G SRR1234567 # 转换,-e 指定线程数 fasterq-dump --split-files --threads 8 --outdir ./fastq_raw SRR1234567--split-files是必须加的。不加的话,双端数据会被塞进一个文件里,后面完全没法用。
3.2 fasterq-dump --split-files 出来的文件到底哪个是 barcode
跑完上面的命令,./fastq_raw目录里会出现SRR1234567_1.fastq和SRR1234567_2.fastq,有时候还有_3.fastq。
这里就是最容易翻车的地方。默认情况下,_1对应 Read1,_2对应 Read2,_3一般是 index read。但并非所有提交者都按照这个顺序提交数据,尤其是早期项目或者作者自己用 bcl2fastq 转换后上传的情况。
判断方法就是前面提到的看长度:
# 看第一条序列的长度 zcat SRR1234567_1.fastq.gz | sed -n '2p' | wc -c zcat SRR1234567_2.fastq.gz | sed -n '2p' | wc -c如果_1是 28 左右,_2是 90 以上,那就是标准情况,不用动。如果反过来了,你就需要在后续命名的时候手动交换,把短的那个命名为 R1,长的命名为 R2。这个交换操作一定要在命名阶段做掉,不要去改文件内容。
另一个情况是出现_3.fastq。这个是 index read(I1),Cell Ranger 不需要它,因为它只看 FASTQ 内容不读 index 信息做拆分。你可以直接删掉以节省空间,但如果这个 run 里混了多个样本,那_3就是拆样本的唯一线索,得留着配合--lanes或者其他工具处理。
3.3 重命名成 _S1_L001_R1_001 规范
Cell Ranger 对输入 FASTQ 的文件名有严格的格式要求,这一点官网上写得很清楚但很多人第一次看会忽略。规范的命名格式是:
[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz举个例子,SampleA_S1_L001_R1_001.fastq.gz和SampleA_S1_L001_R2_001.fastq.gz。其中Sample Name必须和你命令行里--sample参数的值完全一致,L001是 lane 号,单 lane 数据固定写 L001 就行,R1/R2是读端。
所以转换完之后的关键一步就是重命名:
cd ./fastq_raw mv SRR1234567_1.fastq SampleA_S1_L001_R1_001.fastq mv SRR1234567_2.fastq SampleA_S1_L001_R2_001.fastq pigz -p 8 SampleA_S1_L001_R1_001.fastq pigz -p 8 SampleA_S1_L001_R2_001.fastq注意pigz这个工具,它是 gzip 的多线程版本,压缩 20G 的文件用它能快好几倍。原始 FASTQ 不压缩的话,一个样本可能占 60 到 80 GB,压完通常能到 15 到 25 GB,磁盘压力的差别相当大。
提示:如果你的样本有多个 lane,每个 lane 都要重命名成
_L002_、_L003_这样的形式,--sample只写样本名,Cell Ranger 会自动把所有 lane 合并。命名不一致是后面 "No fastq files found" 报错的第一大原因。
3.4 磁盘、内存与并发参数
这条链路里最容易在硬件上翻车。几个实测数字供参考:
--threads参数控制 fasterq-dump 的并发。这个值不是越大越好,因为每个线程都要占内存和磁盘 IO。在 16 核 64G 内存的机器上,我给--threads 8比较稳妥;32 核的机器可以给到 16。给太大反而会因为磁盘瓶颈导致整体变慢。
磁盘空间上,fasterq-dump在转换过程中会同时存在.sra文件、中间的临时文件和输出的 FASTQ,峰值占用可能是最终 FASTQ 体积的两倍多。所以转换一个样本前,最好预留 100 GB 以上的空间。我有个朋友就是在磁盘只剩 40G 的时候跑转换,跑到一半把系统盘写满了,最后连 log 文件都没写出来。
prefetch还有一个--max-size参数,默认是 20G,超过这个大小的 run 会被跳过。10X 的 run 很容易超过这个值,所以一定要显式给大一点,比如--max-size 200G,否则你会看到它下载了几秒钟就"完成"了,其实什么都没下。
3.5 断点续传与失败重试
大规模下载最怕的就是跑到 90% 断掉。prefetch在这点上做得很好,它的缓存机制让你可以直接重跑命令:
# 第一次 prefetch --max-size 200G SRR1234567 # 断了之后直接再来一次,它会续传 prefetch --max-size 200G SRR1234567下载完成后建议做一次完整性校验:
vdb-validate ~/ncbi/public/sra/SRR1234567.sra这个命令会检查本地 SRA 文件的完整性,输出ok才算真正下好了。别跳过这一步,否则你可能在 fasterq-dump 跑到一半才发现数据是坏的。
如果项目里有几十个 run,写个循环是必然的:
while read srr; do prefetch --max-size 200G "$srr" fasterq-dump --split-files --threads 8 --outdir ./fastq_raw "$srr" done < srr_list.txt这个循环可以放后台跑,用nohup或者screen挂着,第二天来看结果。
4. 如果 GEO 给的是 BCL 原始数据
4.1 判断手里的 RAW.tar 是不是 BCL
解压GSE123456_RAW.tar之后,你会看到一堆GSMxxxxxx.tar。随便挑一个解开,如果目录结构里有Data/Intensities/BaseCalls/这样的路径,并且有.bcl或者.bcl.gz文件,还有RunInfo.xml、runParameters.xml,那这就是标准的 Illumina BCL 输出目录。
如果是fastq/目录下面直接是.fastq.gz,那就跳回上一节的重命名流程。
BCL 格式的好处是保留了全部原始信息,理论上可以重新做 base calling;坏处是必须用 bcl2fastq 转换,多一道工序,而且 bcl2fastq 的安装在某些系统上有点折腾人。
4.2 bcl2fastq 与 mkfastq 的版本对应关系
Cell Ranger 里的mkfastq实际上是 bcl2fastq 的一个封装,它自己不带 bcl2fastq 的二进制,需要你在系统里单独装好。版本对应关系大致是:
| Cell Ranger 版本 | 需要的 bcl2fastq 版本 |
|---|---|
| 3.x | bcl2fastq2 v2.20 |
| 4.x | bcl2fastq2 v2.20 |
| 5.x | bcl2fastq2 v2.20 |
| 6.x | bcl2fastq2 v2.20 |
| 7.x | bcl2fastq2 v2.20 |
可以看到其实都是 v2.20,这个版本从 2017 年发布之后就没怎么变过。安装方式有两种:一是用官方提供的 rpm/deb 包,二是从源码编译。第二种比较麻烦,依赖一堆 boost 和 zlib,个人建议能用包管理器就用包管理器。
装完之后用bcl2fastq --version确认一下,然后在 Cell Ranger 里用--bcl2fastq或者--bcl2fastq2参数指定路径。
4.3 SampleSheet.csv 怎么写
BCL 转换的核心是 SampleSheet.csv,它告诉 bcl2fastq 怎么根据 index 把数据拆成不同样本。10X 的 SampleSheet 通常长这样:
[Header] IEMFileVersion,4 Investigator Name,xxx Experiment Name,xxx [Data] Lane,Sample_ID,Sample_Name,Index,Sample_Project 1,SampleA,SampleA,SI-GA-A1,Project1 1,SampleB,SampleB,SI-GA-A2,Project1关键点在Index这一列。10X 的 v2 和 v3 用的是成套的 index(比如 SI-GA-A1 到 SI-GA-H12 这一组),如果你的 SampleSheet 里 index 写错了,转换出来的样本会全部为空。最保险的做法是直接用 GEO 附件里自带的 SampleSheet.csv,不要自己重写。
然后跑:
cellranger mkfastq --id=run1 \ --run=/path/to/bcl_dir \ --csv=/path/to/SampleSheet.csv \ --localcores=16 --localmem=64出来的结果在run1/fastq/下面,目录结构自动就是样本名分子目录,文件名也已经是规范格式,可以直接喂给 count。
4.4 参考基因组的版本陷阱
这一步跟 BCL 无关但同样致命。Cell Ranger 需要一个参考基因组,10X 官方提供了预构建的包,比如refdata-gex-GRCh38-2020-A。这个包有版本,而且和 Cell Ranger 的版本有对应关系:
- 2020-A 系列可以用在 Cell Ranger 3.0 到 7.0
- 2024-A 需要用 Cell Ranger 8.0 及以上
如果你用的是老版本参考基因组配合新版本 Cell Ranger,多数情况能跑,但会提示警告;反过来,新参考配老软件则可能直接报错。
还有一个常见误解:GRCh38-2020-A 和 GRCh38-3.0.0 不是同一个东西。前者更新了基因注释,包含了更多的 lncRNA 和更新过的基因名。用不同版本跑出来的基因数会有差异,做横向比较的时候必须保证全项目用同一个版本。我的习惯是在项目开始就把参考基因组的 md5 记在 README 里,避免半年后回来看不知道自己用的哪个版本。
5. cellranger count 跑通与参数取舍
5.1 目录结构与 --sample 的匹配逻辑
假设你的 FASTQ 已经整理好了,目录结构长这样:
fastqs/ ├── SampleA_S1_L001_R1_001.fastq.gz ├── SampleA_S1_L001_R2_001.fastq.gz ├── SampleB_S1_L001_R1_001.fastq.gz └── SampleB_S1_L001_R2_001.fastq.gz那么跑 SampleA 的命令是:
cellranger count --id=SampleA_run \ --transcriptome=/ref/refdata-gex-GRCh38-2020-A \ --fastqs=./fastqs \ --sample=SampleA \ --expect-cells=5000 \ --localcores=16 \ --localmem=64--sample的值是SampleA,Cell Ranger 会自动去--fastqs目录里匹配所有以SampleA_开头的文件,收集成 R1/R2 配对。这就是为什么命名必须严格——它匹配靠的是前缀字符串,不是智能识别。
一个非常容易犯的错:--sample写成了文件名的完整前缀但带了_S1,比如写--sample=SampleA_S1。这样它匹配不到任何文件,报 "No fastq files found"。记住--sample只需要样本名前缀,_S1_L001_R1_001这部分是格式约定,不用写进去。
5.2 expect-cells 和资源参数怎么给
--expect-cells这个参数的作用是给 Cell Ranger 一个初始的细胞数量预期,它会根据这个值来调整 barcode 的过滤策略。给得太少,真实的细胞会被当成背景过滤掉;给得太多,背景 barcode 会被误判成细胞。
经验值是这样:先大致估一下样本上机时目标细胞数,然后按这个数来给。如果完全不知道,可以从 3000 到 5000 起步,跑完看 web_summary 里的Estimated Number of Cells,如果这个值跟你设的 expect-cells 差得很远,就调整参数重跑一次。常见的判断标准是估计值落在 expect-cells 的 0.5 到 2 倍区间内比较合理。
--localcores和--localmem控制资源。这两个参数填错会导致两个典型问题:cores 给太多而系统实际核数不够,任务会因为抢不到资源而卡住;mem 给超过物理内存,会在比对阶段被系统 OOM killer 干掉,日志里只会留一行被杀的记录,非常难排查。
我的一般做法是localcores给到系统核数的 80%,localmem给到物理内存的 80%。比如 32 核 128G 的机器,就给 24 核 100G。
5.3 跑完之后先看 web_summary 的哪几项指标
cellranger count跑完大约需要 2 到 8 小时,取决于数据量和机器性能。出来之后第一件事是打开outs/web_summary.html,重点看这几项:
| 指标 | 参考范围 | 说明 |
|---|---|---|
| Estimated Number of Cells | 与预期相符 | 和 expect-cells 差太多要警惕 |
| Mean Reads per Cell | > 20000 | 低于 10000 说明测序深度不足 |
| Median Genes per Cell | > 1000 | 太低可能是细胞活力差或 RNA 降解 |
| Fraction Reads in Cells | > 70% | 低于 50% 说明背景噪音大 |
| Sequencing Saturation | 60% - 90% | 过高说明测序过饱和,加测无意义 |
| Q30 Bases in Barcode | > 85% | 低于 80% 要检查数据质量 |
| Valid Barcodes | > 75% | 太低可能是 barcode 端有污染 |
这几项里我最看重的是Fraction Reads in Cells和Median Genes per Cell。前者反映样本的干净程度,后者反映 RNA 的完整性。如果这两个指标都正常,基本可以放心往下做。
5.4 输出目录里哪些文件后面还会用到
outs/目录里文件不少,但真正高频使用的就那几个:
filtered_feature_bc_matrix/是过滤后的矩阵,只包含判定为真实细胞的 barcode,是下游分析的主力输入。raw_feature_bc_matrix/包含全部 barcode,做背景估计或者空液滴分析时用。possorted_genome_bam.bam是排好序的 BAM 文件,做 CNV 分析、可变剪接、或者想自己重新定量的时候需要它,但体积很大,一个样本可能 30 到 50 GB。
molecule_info.hdf5是做aggr合并多样本时必须的文件,很多人跑完 count 就把它删了,等到要合并样本时又得重新跑,非常亏。cloupe.cloupe是给 Loupe Browser 用的可视化文件,做汇报的时候挺方便。
6. 报错排查:从日志第一行开始
6.1 No fastq files found
这是出现频率最高的报错,原因就那么几种:--fastqs路径写错了(相对路径和绝对路径混用最容易出问题);文件名不符合_S1_L001_R1_001规范;--sample和文件名前缀不一致;R1和R2标签写反。
排查方法是从内向外一层层确认。先ls一下--fastqs指的那个目录,确认文件在;再用ls fastqs/SampleA_*确认前缀匹配得上;最后确认每个样本都有成对的 R1 和 R2。
注意:Cell Ranger 对大小写敏感。
SampleA和samplea是两个不同的东西,命名的时候统一风格,别一会儿大写一会儿小写。
6.2 读长与chemistry不匹配
报错信息里可能出现Read 1 is too short或者The read length does not match the expected chemistry。这是因为 Cell Ranger 会自动检测化学版本,如果 read 长度偏离了它预期的范围就会报错。
v2 化学的 R1 是 26bp,v3 是 28bp,两者差 2bp。大部分情况下 Cell Ranger 能自动识别,但如果你之前在 fastq 处理中做过 trim,把 R1 截短了,它就会认不出来。
解决办法是用--chemistry参数显式指定,可选值包括SC3Pv2、SC3Pv3、SC3Pv3LT、SC5P-PE等。指定对了之后报错就消失了。但更好的做法是根本不要对 R1 做 trim,因为 barcode 和 UMI 就在这 28bp 里,截短任何一个碱基都会导致比对失败。
6.3 一个 SRR 里混了多个样本
有时候你会遇到一个 SRR 跑完,Estimated Number of Cells 高得离谱,比如预期 5000 结果出来 20000。这很可能是因为这个 run 里用 index 混了多个样本。
判断方法是看 SRA Run Browser 里的Library Strategy和样本描述。如果多个 BioSample 共享一个 Run,那就是混样了。
处理这种数据比较麻烦,需要根据 index 序列来拆分。可选方案是:把_3.fastq(index read)拿出来,用demultiplex类工具按 index 拆分,或者干脆用cellranger multi配合 feature barcode 的配置。这种情况我一般会优先去找作者是不是另