信息速览
recount3 是一件不太像「数据集」的数据集:它没有去采一次样、没有开一台测序仪,而是把别人早就交上去的原始读段,用自己的一套管线从头重跑了一遍。2021 年 12 月,Johns Hopkins 团队在 Genome Biology 上公布了这项工作,官网页面对它的定位写得很直白——人类与小鼠公共 RNA-seq 数据的统一处理摘要。它要卖的不是数据本身,是「同一把尺子」。
这个尺子值钱,是因为公共测序档案里横着上万份各自为政的研究。每一份的作者都用自己的比对工具、自己的注释版本、自己的过滤阈值,于是两份研究里的同一批基因,数值放在一起没法比。recount3 把比对、剪接位点统计、覆盖曲线生成、基因与外显子汇总全部固定下来,任何人拿到的两个项目,计数是在同一条流水线上产出的。论文口径的处理规模是七十余万条测序 run,输入侧压着近一千 TB 的压缩读段。
代价同样清楚:你放弃了对处理细节的控制权,接受了一套别人定死的注释与版本,还要面对一个事实——这批数据里有大约一半是全长单细胞协议,混用会出事。本篇按「它是什么、结构长什么样、怎么拿到、哪里会翻车」的顺序展开,所有规模数字以论文正文表 1 与官方文档为准,来源分级见附录 B。一屏速览见下表。
| 维度 | 内容 |
|---|---|
| 官方名 | recount3(全小写,R 包同名) |
| 一句话定位 | 把公共 RNA-seq 原始读段用同一管线重跑,产出跨研究可比的计数 |
| 发布方 | Johns Hopkins University 团队(Langmead、Hansen、Collado-Torres 等) |
| 发表 | Genome Biology,2021-12,doi:10.1186/s13059-021-02533-6 |
| 处理管线 | Monorail(网格化调度,核心为 STAR 2.7.3a) |
| 测序 run 合计 | 论文表 1:七十五万余;摘要口径:七十五万以上 |
| 人类 run | 316,443(SRA)+ GTEx v8 与 TCGA 合计约三万 |
| 小鼠 run | 416,803(全部来自 SRA) |
| 数据来源 | SRA、GTEx v8、TCGA(GEO 经 SRA 间接覆盖,无 TCP 源) |
| 参考基因组 | 人类 hg38/GRCh38,小鼠 mm10/GRCm38 |
| 人类注释 | GENCODE v26、v29,RefSeq,FANTOM-CAT v6,另含 ERCC 与 SIRV 对照 |
| 小鼠注释 | 仅 GENCODE v23 |
| 输出粒度 | 基因级计数、外显子级计数、外显子-外显子剪接位点计数、碱基级 bigWig |
| 体量 | 输入约 990 TB 压缩读段,摘要输出约一百五十 TB(论文 146.31 TB) |
| 访问方式 | HTTPS 明文矩阵直接下载 / recount3 与 snapcount 两个 R 包 |
| 文件许可 | 软件 MIT;数据随来源,论文声明 CC0 适用于文内数据 |
§1 它要解决的具体麻烦
1.1 公共测序档案的真正痛点不是数据少,是没法比
公共测序档案的体量早已不是瓶颈。SRA 每十八个月翻一倍,任何人想找「人类肝脏的 RNA-seq」都能找到成百上千条记录。麻烦出在下一步:这些记录来自不同实验室、不同年份、不同测序平台,每一条的数值都是在各自的处理流程里算出来的。有的作者用 HISAT2 加 featureCounts,有的用 STAR 加 RSEM,有的用 Kallisto 直接对转录本做伪比对;注释版本从 GENCODE v19 到 v40 都有;有的过滤低质量读段,有的不过滤。
这些差别叠在一起,使「表达量」这个词失去了可比性。你把 A 研究的肝组织 counts 和 B 研究的肝组织 counts 放在一张表里,看到的大部分差异来自流程,而不是生物学。要去掉这部分差异,传统做法是逐篇读方法学、逐篇复现、逐篇重算,一篇论文的时间成本就是一个人周。上千篇就是上千个人周,没有实验室扛得住。
recount3 的回答是把这件事做一次、做到底:所有原始读段进同一台机器,出来的是同一套数值。它不关心某一篇研究原本用了什么方法,因为那些方法被全部丢弃了。这个决断看起来粗暴,实际上是把「可比」放在了「忠于原作者意图」之前。
1.2 三条常见的错误期待
第一条期待是「它收集了新数据」。recount3 一条新测序都没有做。它的全部输入都是别人已经提交到公共档案的读段。这意味着任何在 SRA 里不存在的研究,在 recount3 里也不会存在;任何被原作者撤下或替换的提交,recount3 也无从知晓。把它当「转录组的全集」会失手,它只是「某个时间点之前、SRA 里能找到的那些」的快照。
第二条期待是「它已经做好归一化」。交付的基因级矩阵是原始整数计数,没有任何标准化。文件头部只有两行注释(注释代号与生成时间),没有 TPM、没有 FPKM、没有文库大小列。要做跨样本比较,使用者必须自己上 DESeq2、edgeR 或类似的工具,甚至还要额外处理项目间的批次效应。把它当「现成的表达量矩阵」直接喂给模型,是最常见的一次翻车。
第三条期待是「同一项目的所有样本都是同一种实验」。这一条最隐蔽。论文原文写明,大约一半的 run 来自全长单细胞协议,指的是 Smart-seq 与 Smart-seq2 这一类。这类协议一条 run 往往只对应一个细胞或极少数细胞,基因检出率低、3’ 端没有偏倚但覆盖不均,跟批量组织样本放在同一个矩阵里,等于把两种完全不同的测量混作一谈。项目内部如此,跨项目更是如此。
1.3 它能给到的与给不到的
| 你的需求 | 它给的 | 够不够 |
|---|---|---|
| 同一研究内部的样本比较 | 同管线、同注释、同版本的整数计数 | 够,配上项目自带的协变量即可 |
| 跨多个研究的合并分析 | 全部研究共用一套计数口径 | 够方向,但须自行做批次校正 |
| 剪接事件层面的分析 | 外显子-外显子剪接位点计数,含未注释位点 | 够,且这是它相对同类最厚的部分 |
| 想看某个位点的读段覆盖形状 | 每个样本一份全基因组 bigWig | 够,代价是文件多、体量大 |
| 想做转录本水平的定量 | 只到基因与外显子,没有转录本 | 不够,须回原始读段自行定量 |
| 想要已归一化、可直接建模的表达矩阵 | 原始整数计数 | 不够,须自行变换 |
| 想要样本的临床与表型细节 | 只有从来源抽取的少量字段 | 不够,须回原研究补 |
| 想要最新提交的样本 | 冻结在某个处理批次 | 不够,须回 SRA 或自己跑管线 |
1.4 你的处境与它能帮到的事
| 你的处境 | 它能帮你什么 | 下一步该做什么 |
|---|---|---|
| 手上有一组小样本转录组,找不到合适对照 | 从同一组织的大量公共样本中挑出可比的背景集 | 先用元数据筛组织与实验类型,再看 QC 字段剔除异常样本 |
| 想验证一个基因在多种组织里是否普遍表达 | 直接读跨项目汇总的计数,无需逐篇下载 | 用 R 包按组织拼接,先做归一化再看趋势 |
| 在做一个表达基础模型,缺大规模预训练语料 | 提供同一口径的数十万样本矩阵,格式统一 | 注意单细胞与批量混编问题,训练前按协议分层 |
| 自己刚做完一批 RNA-seq,想和公共数据比 | 用同一管线处理自己的数据即可直接并列 | 走 Monorail 本地或云端流程,输出可对齐 |
| 想找某个剪接位点在哪些样本里出现 | 用查询接口按区域直接取,不必下载整项目 | 用 snapcount 的查询接口,按样本数过滤 |
| 需要复现别人的差异表达结果 | 提供标准化口径下的近似输入 | 先读原研究的方法学,确认它是否也用了统一管线 |
1.5 一句话判断
recount3 的价值不在数据量,而在它把「可比」这件原本要逐篇手工完成的事,变成了一次性的基础设施。
1.6 不适合谁:四条边界
不适合需要转录本异构体定量的人。 全套输出止步于基因与外显子,剪接位点计数能告诉你某个连接处有没有读段,但不能直接折算成异构体比例。要做异构体水平的分析,得回到原始读段,用 Salmon、Kallisto 或 RSEM 重新跑。
不适合把样本量等同队列规模的人。 七十余万条 run 听起来像一个巨型队列,实际上它们来自上万份彼此独立的研究,同一个体的多次测序、同一细胞系的多个批次都在里面。把它当作「七十万人」来算统计功效,会得到一个完全错误的结论。
不适合需要最新样本的时效性任务。 处理批次冻结在某个时点,SRA 在那之后的新提交不在其中。做流行病学式的实时监测,它帮不上忙。
不适合要求原始读段或个体级数据的人。 它发放的是汇总层。GTEx 与 TCGA 的原始读段仍受 dbGaP 与 GDC 的访问条款约束,recount3 的开放性不覆盖那一层。
§2 数据构成与规格
2.1 三层目录与两个物种
recount3 的托管结构是一棵很浅的树。顶层按物种分 human 与 mouse,第二层按数据来源分 sra、gtex、tcga 三个目录(小鼠只有 sra),第三层按输出类型分 gene_sums、exon_sums、junctions、base_sums、metadata 五类。再往下是项目代号,项目代号再往下才是具体文件。
这个设计的用意是让路径本身携带语义。你看到一个路径里有 gene_sums,就知道拿到的是一张基因乘样本的矩阵;看到 junctions,就知道拿到的是稀疏的三元组。文件名则是三段式:来源、输出类型、项目代号,再接注释代号或矩阵片段后缀。这种命名法在批量脚本里可以直接用正则解析,不需要额外查表。
两个物种的成熟度并不对称。人类侧有四套注释、三个来源、额外提供 UNIQUE 类型的剪接位点矩阵;小鼠侧只有一套注释、一个来源、只有 ALL 类型的剪接位点。做跨物种比较时必须记住这一点,否则会以为小鼠数据「缺了一块」,实际上是本来就没提供。
2.2 文件与字段的完整字典
| 层级 | 名称 | 内容 | 规模或形态 |
|---|---|---|---|
| 目录 | human / mouse |
物种分区 | 两个 |
| 目录 | data_sources/sra |
来自 SRA 的项目 | 人类与小鼠都有 |
| 目录 | data_sources/gtex |
GTEx v8,按组织切分 | 人类专属 |
| 目录 | data_sources/tcga |
TCGA,按癌种切分 | 人类专属 |
| 目录 | gene_sums |
基因级计数矩阵 | 每项目每注释一份 |
| 目录 | exon_sums |
外显子级计数矩阵 | 每项目每注释一份 |
| 目录 | junctions |
外显子-外显子剪接位点计数 | 每项目三份(MM/ID/RR) |
| 目录 | base_sums |
碱基级覆盖 | 每样本一份 bigWig |
| 目录 | metadata |
五张项目级与样本级元数据表 | 每项目五份 |
| 矩阵文件 | *.gene_sums.*.G026.gz |
注释代号编码在文件名中 | 制表符分隔文本 |
| 矩阵头 | ##annotation= |
该矩阵使用的注释代号 | 单行 |
| 矩阵头 | ##date.generated= |
该矩阵的生成时间戳 | 单行,精确到微秒 |
| 矩阵体 | 首列 gene_id |
基因标识,带版本号后缀 | 每行一个基因 |
| 矩阵体 | 其余列 | 每列一个样本,列名即 run accession | 整数计数 |
| 剪接文件 | *.ALL.MM.gz |
Matrix Market 格式的稀疏计数 | 三元组 |
| 剪接文件 | *.ALL.ID.gz |
样本标识符顺序表 | 与 MM 的列一一对应 |
| 剪接文件 | *.ALL.RR.gz |
每个剪接位点的染色体、起止、链 | 与 MM 的行一一对应 |
| 元数据 | recount_project |
来源侧项目信息 | 八列 |
| 元数据 | recount_qc |
比对与汇总质量指标 | 数十列 |
| 元数据 | recount_seq_qc |
读段层面质量指标 | 数十列 |
| 元数据 | recount_pred |
预测与人工整理的组织、类型字段 | 十列 |
| 元数据 | project_meta |
定位文件所需的内部信息 | 内部使用 |
| 注释文件 | human.gene_sums.G026.gtf.gz |
基因区间定义 | GTF |
| 注释文件 | human.exon_sums.G026.gtf.gz |
外显子区间定义 | GTF |
2.3 一条记录的完整解剖
拿一个真实存在的人类项目做解剖。SRP009615 是数据的常用示例,属于 SRA 来源、人类、使用 GENCODE v26 注释。它的基因级矩阵文件头部有两行注释,第三行才是表头。
##annotation=G026
##date.generated=2020-05-12 01:01:49.411644
gene_id SRR389077 SRR387777 SRR387778 SRR389078 SRR387779 SRR389079
ENSG00000278704.1 0 0 0 0 0 0
ENSG00000277400.1 0 0 0 0 0 0
ENSG00000274847.1 0 0 0 0 0 0
三行足以说明几个结构事实。第一,矩阵是稠密的制表符分隔文本,行数约六万三千余行,对应注释里的全部基因,包括那些全零的基因,并不做稀疏化,也不剔除低表达基因。第二,列名是 run accession,与 SRA 的 SRR 编号一一对应,这是跨库连接的天然外键。第三,头部只有两行元信息,没有任何归一化系数,也没有任何样本级的质量标记。
剪接位点则是另一套结构。它被拆成三份文件:一份是 Matrix Market 格式的稀疏矩阵,里面每一行是一个剪接位点、每一列是一个样本、值是跨越该位点的读段数;一份是样本标识符列表,保证矩阵列的顺序可还原;一份是位点坐标表,给出每个位点的染色体、起点、终点与链方向。三份文件必须一起取,单独取任何一份都没有意义。这种拆分是为了让矩阵可以按需切块,而不是被当成一张巨大的稠密表。
外显子级矩阵的结构与基因级几乎相同,差别只在于行是外显子区间而不是基因,行数会多出数倍。碱基级覆盖则不落成表,而是每个样本压缩成一个 bigWig,读取需要专门的索引支持,这也是它体量最大的输出层。
2.4 各层的规模对照
| 输出层 | 组织形态 | 读取方式 | 体量感受 |
|---|---|---|---|
| 基因级计数 | 稠密文本矩阵,行是基因、列是样本 | 任何能读制表符分隔文件的工具 | 单个小项目一两兆,大项目上百兆 |
| 外显子级计数 | 稠密文本矩阵,行是外显子 | 同上,但内存需求高数倍 | 约为基因级的数倍 |
| 剪接位点计数 | 稀疏三元组拆分三文件 | 需要按 MM/ID/RR 拼装 | 小项目几兆,大项目数十兆 |
| 碱基级覆盖 | 每样本一份 bigWig | 需要 bigWig 读取库 | 样本数乘单个文件,总量级最大 |
| 元数据 | 五张制表符分隔表 | 直接读表 | 每张表几十到几百千字节 |
| 注释定义 | GTF | 区间工具 | 每套注释数十兆 |
论文给出的处理账本可以从另一个角度看这批数据的体量:输入侧压缩读段接近一千 TB,输出摘要约一百五十 TB,压缩比大致是六到七比一。计算侧用了约三万节点小时,折算约三年单节点时间。这三个数字合起来说明一件事——重跑一遍的成本不高到一个国家的超算中心才做得起,但也不低到可以随手复现。这正是它作为基础设施存在的理由。
2.5 与同类资源的字段差异
| 字段或能力 | recount3 | 只做基因级汇总的资源 | 只做原始档案的资源 |
|---|---|---|---|
| 基因级计数 | 有,多注释 | 有,通常单一注释 | 无 |
| 外显子级计数 | 有 | 多数无 | 无 |
| 剪接位点计数 | 有,含未注释位点 | 多数无 | 无 |
| 碱基级覆盖 | 有,每样本 bigWig | 多数无 | 需自行生成 |
| 跨研究统一管线 | 是 | 视资源而定 | 不适用 |
| 转录本级定量 | 无 | 部分有 | 无 |
| 原始读段 | 无,仅汇总 | 无 | 有 |
| 样本临床细节 | 极少 | 少 | 取决于提交者 |
§3 二次加工:它为什么不是「又一个 SRA 镜像」
3.1 「不产生新数据」是特性不是缺陷
把 recount3 说清楚,关键是把「二次加工」这四个字的位置摆正。它不是 SRA 的镜像,不是 GEO 的重新打包,也不是某个联盟新测的一批样本。它是同一批原始读段在另一套处理器下的产物。这个关系类似于同一批 RAW 文件被两套不同的显影参数冲洗——底片没变,成像变了。
这层定位带来一个容易被忽视的后果:recount3 与原研究之间不是替代关系,而是对照关系。原作者算出的 counts 反映了他当时的判断,recount3 的 counts 反映统一管线的判断。两者不一致的地方,恰恰是最有信息量的地方——可能是注释版本差异,可能是过滤策略差异,也可能是原研究真的处理错了。很多方法学论文正是利用这种对照来检验自己的结论是否稳健。
理解这一点,就不会对「它为什么不包含某个新研究」感到困惑。新研究只要进了 SRA,理论上就能进下一批;它没进,只是批次还没跑到,或者跑的时候那批数据还不在。
3.2 Monorail:为「重跑一遍」而生的调度系统
支撑这次重跑的系统叫 Monorail,是同一团队为这个项目新写的分析管线。它的设计思路是网格计算:一个中心调度器把上百万个独立的工作单元分派到各计算节点,每个单元对应一次测序 run 的比对与统计,单元之间不共享状态,失败只影响自己。
之所以要这样设计,是因为这类任务的失败率天然很高。公共档案里的读段质量参差,有的文件损坏、有的格式不规范、有的只有一端可读。如果整批任务是一个大作业,任何一条坏数据都会拖垮全局。拆成一run一单元之后,坏的那一条单独失败、记入日志,其余继续。这个工程选择直接决定了项目能不能在可接受的时间内跑完七十余万条 run。
管线内部的具体工具在论文方法部分列得很清楚:比对用 STAR 2.7.3a,基因与外显子计数用 featureCounts 1.6.3,剪接位点与覆盖用 MegaDepth 0.4.0,另有 samtools、seqtk、tabix 等辅助工具。整个工具链的版本被固定下来,这是「可比」承诺的技术基础。任何版本漂移都会破坏可比性,所以它宁可冻结在旧版本,也不追新。
3.3 注释无关的第一步与注释相关的第二步
管线在结构上分成两步,这个划分是理解 recount3 输出的钥匙。
第一步是注释无关的:STAR 比对到参考基因组,比对时不提供任何基因注释作为剪接提示。这一点在同类资源里并不普遍,多数流程会把 GTF 喂给比对器以提高剪接位点识别率。recount3 刻意不这么做,目的是让比对结果不带任何注释偏向。比对完成后,MegaDepth 直接从比对文件里数剪接位点,包括那些不在任何注释里的、由数据本身支撑的新连接。这一步的产物是剪接位点计数与覆盖曲线,与注释无关。
第二步才是注释相关的:拿第一步的比对结果,按选定的注释用 featureCounts 汇总基因与外显子计数。同一次比对可以配多套注释,所以人类样本会同时有 GENCODE v26、v29、RefSeq、FANTOM-CAT 四套基因矩阵。这一步的计算量远小于第一步,所以日后新增一套注释并不需要重跑比对。
这个两段式设计的实际好处是:使用者如果只信任某套注释,可以直接取那套的矩阵;如果不信任任何注释,可以绕过第二步,直接用剪接位点与覆盖做数据驱动的外显子发现。两条路都通。
3.4 参考序列里为什么混了 ERCC 与 SIRV
论文提到,人类与小鼠的 STAR 索引在构建时都额外加入了 ERCC 对照序列与 SIRV 合成转录本的拼接外显子。这不是随手添加的:ERCC 是外源掺入的已知浓度 RNA 标准品,SIRV 是已知结构的合成剪接异构体。把它们放进索引,可以让管线的输出里保留这些对照的计数,方便使用者判断某次实验的定量下限与剪接检测灵敏度。
对使用者而言,这意味着矩阵里会出现一些并非来自样本基因组的行。做常规分析时这些行通常被忽略,但在评估数据质量或做方法学基准时,它们是不可多得的内部参照。很多人拿到矩阵后直接按基因名过滤,顺手就把这部分信息丢掉了。
3.5 输出为什么是文本而不是二进制
recount3 的矩阵是 gzip 压缩的制表符分隔文本,不是 HDF5、不是 Parquet、不是任何列式格式。这个选择在 2021 年有其合理性:文本格式对任何语言都是零门槛,R、Python、awk、甚至 Excel 都能打开;文本也能被 HTTP 的 Range 请求部分读取,配合 gzip 的分块可解压特性,前端可以做流式解析。
代价是体积与解析速度。同一个矩阵用列式格式存,体积通常能再降数倍,读取速度快一个量级。对做全量训练的使用者,这个代价是实打实的:你下载的是文本,解析成内存数组还要再花时间。这也是为什么很多人宁愿用 R 包读,至少它把解析逻辑封装好了。真要上规模,通常的做法是先用脚本把文本流转成每一列一个二进制数组的紧凑格式,训练时按列读,磁盘占用与解析时间都能降下来。
§4 质量的那条线划在哪里
4.1 统一管线不等于统一质量
「统一处理」这个词容易被读成「质量统一」。事实并非如此。管线统一的是处理步骤,不是输入质量。一条深度只有几百万读段的 run 和一条有上亿读段的 run,走过同一套管线之后,产出的仍然是两种质量的数据。recount3 没有替使用者做剔除,只是把判断所需的指标一并交出来。
这些指标在元数据的 QC 表里。字段包括比对到常染色体、X 染色体、Y 染色体的读段比例,外显子与外显子-外显子连接的读段占比,以及按所有读段与唯一读段两种口径分别统计的覆盖指标。这些字段的设计意图很明确:让使用者可以按自己的标准筛。官方没有给出一个「合格线」,也没有标出哪些样本被剔除。
这个处理方式在工程上是诚实的,在使用上是危险的。没有默认过滤意味着默认全收,而全收的结果里混着大量不该进分析的低质量样本。任何跨项目分析的第一步,都应该是按 QC 字段做一遍粗筛,而不是直接把整张矩阵丢进去。
4.2 快一半的单细胞:最大的技术分层
论文里那句「大约一半的 run 来自全长单细胞协议」,是整份文档里最该被高亮的一句。它的含义是:如果在不分层的情况下把 recount3 的全部小鼠 run 当作一个整体,你实际得到的是一个由批量组织样本与单细胞样本混合而成的集合,两者的技术特征差异远大于大多数生物学差异。
这类协议的特点是一条 run 对应一个细胞或极少数细胞,读段总数通常远低于批量样本,基因检出以「有或没有」为主,表达矩阵极度稀疏。混进批量样本的矩阵里,它们会在归一化步骤中拉低整体的文库大小分布,在差异分析中制造大量假阳性。
好在元数据里有分层线索。预测表提供了类型字段与组织字段,其中类型字段会标出 rna-seq 这样的粗分类;质量表提供了文库大小的统计量。做分层时至少要按「是否单细胞协议」与「文库大小量级」两条线切开,再考虑组织与项目。
4.3 元数据是这套数据最软的一环
recount3 的样本注释不来自人工整理,而来自对提交记录的自动抽取与预测。论文在讨论部分专门提到,元数据可能不可靠,并把它列为后续工作的重点之一:需要更好的模型来预测缺失字段、纠正错误字段。
具体表现是:组织字段可能缺失,可能写成同义词的变体,也可能是从标题里猜出来的;细胞类型字段经常为空;样本类型(是否疾病、是否处理)多数没有结构化。这些字段在探索性分析里够用,在做严格的病例对照设计时不够。
实际做法是把 recount3 的元数据当作索引而不是真值:用它快速定位候选项目,再回到原研究的提交记录或论文里核对关键分组。多花这一步,能省掉后面大量的返工。
4.4 案例:SRP009615 的一条 run 如何被追溯
拿 SRP009615 里的 SRR389077 走一遍完整的追溯链路,可以把这套数据的连接关系具象化。这条 run 在基因级矩阵里是第二列,读数由矩阵直接给出;在预测表里,它的样本标识是 SRS282369,类型字段标为 rna-seq,组织字段标为 Blood;在质量表里,它的 rail_id 与基因矩阵的列顺序对应,比对到常染色体的读段比例、外显子读段占比等指标都能查到。
这条链路说明三件事。第一,矩阵与元数据通过 run accession 连接,不需要额外的映射表。第二,一个 run 在矩阵里的位置与在元数据表里的行并不必然相同,索引要靠 rail_id 或 accession 对齐,不能靠行号。第三,追溯到最上游时,最终还是要回到 SRA 的 run 页面看原始提交信息,recount3 本身不复制那些信息。
这个案例的第三点在 §7.3 的偏差讨论里会再次出现,因为绝大多数误用都出在「以为 recount3 的元数据就是全部元数据」这一步。§4.4 的案例在这里的作用是给出一个可复现的最小追溯路径,而不是一个特例。
4.5 可复现性的两个层次
recount3 提供的可复现性分两个层次。第一层是结果可复现:任何人下载同一个文件,得到的矩阵逐字节相同,因为文件是静态托管的。这一层很硬,不依赖任何环境。
第二层是过程可复现:要从原始读段重新生成同一个文件,需要完全一致的 STAR 版本、索引构建参数、featureCounts 参数与注释文件。论文把这些版本号都写进了方法部分,但完整复现仍然需要相当的计算资源与工程投入。多数使用者只需要第一层,做方法学研究的人才会触及第二层。
有第三层,是很多人没意识到的:把自己的数据用同一管线跑一遍,然后与 recount3 并列分析。官方为此提供了本地与云端两种运行方式,论文也专门讨论了用这套管线处理受控访问数据的路径。这层能力让 recount3 从「一个数据集」变成了「一个坐标系」。
§5 取数的四条路
5.1 先查有哪些项目
任何取数动作的第一步都是确认项目存在。人类与小鼠各有一份项目级清单文件,含 rail_id、外部标识、研究代号、项目路径、物种、来源与处理日期八列。这份清单是全部后续操作的索引。
#-- 取人类与小鼠的项目级清单,并统计项目数与来源分布
cd /tmp
curl -sL -o h_proj.gz "https://duffel.rail.bio/recount3/human/data_sources/sra/metadata/sra.recount_project.MD.gz"
curl -sL -o m_proj.gz "https://duffel.rail.bio/recount3/mouse/data_sources/sra/metadata/sra.recount_project.MD.gz"
#-- 项目代号在第三列,物种之外的来源在第六列
zcat h_proj.gz | awk -F'\t' 'NR>1{print $3}' | sort -u | wc -l
zcat m_proj.gz | awk -F'\t' 'NR>1{print $3}' | sort -u | wc -l
zcat h_proj.gz | awk -F'\t' 'NR>1{print $6}' | sort -u
清单之外还有 GTEx 与 TCGA 两份同类文件,路径把 sra 换成 gtex 或 tcga 即可。这两份的项目代号是组织名或癌种代号,不是 SRP 编号,这一点在写脚本时最容易搞错。
5.2 用 R 包按项目取
官方首推的路径是 R 包的函数式接口,它把路径拼装、矩阵下载、行名对齐、元数据合并全部封装,返回一个区间化的实验对象,可以直接进 Bioconductor 的下游工具链。
#-- 列出全部人类项目,挑出示例项目,取基因级矩阵
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("recount3")
library("recount3")
human_projects <- available_projects()
proj_info <- subset(human_projects, project == "SRP009615" & project_type == "data_sources")
rse_gene <- create_rse(proj_info)
dim(rse_gene)
#-- 取外显子级与剪接位点级,参数只需改一个词
rse_exon <- create_rse(proj_info, type = "exon")
rse_jxn <- create_rse(proj_info, type = "junctions")
这个接口的三个参数值得记住:type 控制输出粒度,annotation 控制注释版本,assay_transform 控制是否做变换。默认返回原始计数,不做任何变换。
5.3 不用 R 也能取:直接拼 URL
所有文件都是 HTTP 可取的静态对象,路径规则在官方文档里写得很明确。人类 SRA 项目的基因矩阵路径形如 human/data_sources/sra/gene_sums/<项目代号末两位>/<项目代号>/sra.gene_sums.<项目代号>.<注释代号>.gz。掌握这条规则,用任何语言都能取。
#-- 不用 R,直接按路径规则拼 URL 并读出基因级矩阵
import gzip, io, urllib.request
BASE = "https://duffel.rail.bio/recount3"
def gene_matrix(organism, project, annotation="G026", source="sra"):
url = (f"{BASE}/{organism}/data_sources/{source}/gene_sums/"
f"{project[-2:]}/{project}/{source}.gene_sums.{project}.{annotation}.gz")
raw = urllib.request.urlopen(url).read()
text = gzip.decompress(raw).decode("utf-8")
lines = [ln for ln in text.splitlines() if not ln.startswith("##")]
header = lines[0].split("\t")
sample_cols = header[1:]
body = [ln.split("\t") for ln in lines[1:]]
return sample_cols, body
samples, rows = gene_matrix("human", "SRP009615")
print(len(samples), "samples,", len(rows), "genes")
注意两点:压缩流要整体解压后再切行,不能直接按行读 gzip;头部两行以 ## 开头,必须先剔除再解析表头。
5.4 按区域查询:不必下载整项目
如果你的问题只是「某个基因或某个区间在哪些样本里有读段」,下载整项目是浪费。查询型接口可以直接按区域取,返回的是命中记录而不是整张矩阵。
#-- 用查询接口按基因名或区间取剪接位点,按覆盖样本数过滤
BiocManager::install("snapcount")
library(snapcount)
q <- query_jx(compilation = "sra", regions = "CD99")
q_filtered <- query_jx(compilation = "sra", regions = "CD99",
range_filters = samples_count == 10)
dim(q)
compilation 参数决定查哪一批汇总:人类 SRA、小鼠 SRA、GTEx、TCGA 各是一批。这个接口的优势在于它把这批数据的索引放在了服务端,你只需要传一个区间。
5.5 取碱基级覆盖
覆盖层是这套数据里体量最大的部分,每个样本一个 bigWig。取它时要克制:按需取特定区间,不要整文件拉。
#-- 取单个样本的 bigWig 并只看目标区间,避免整文件下载
cd /tmp
curl -sL -o s1.bw "https://duffel.rail.bio/recount3/human/data_sources/sra/base_sums/15/SRP009615/77/sra.base_sums.SRP009615_SRR387777.ALL.bw"
python3 - <<'PY'
import pyBigWig
bw = pyBigWig.open("/tmp/s1.bw")
print(bw.chroms().get("chr1"))
print(bw.values("chr1", 1000000, 1001000))
bw.close()
PY
bigWig 的随机访问特性让它即使在本地也能秒开,真正贵的是下载时间。做全基因组分析时,更经济的做法是先取基因矩阵定位候选区间,再只对候选区间取覆盖。另一个省事的做法是先用它自带的质量字段看出某个样本的比对率,比对率明显偏低的样本,与其花时间取覆盖,不如直接跳过。
5.6 把注释定义也取下来
计数矩阵给的是基因或外显子的标识,要知道它们在基因组上的位置,还得取注释定义文件。每套注释都有基因级与外显子级两份 GTF。
#-- 取 GENCODE v26 的基因与外显子定义
curl -sL -o g026_gene.gtf.gz "https://duffel.rail.bio/recount3/human/annotations/gene_sums/human.gene_sums.G026.gtf.gz"
curl -sL -o g026_exon.gtf.gz "https://duffel.rail.bio/recount3/human/annotations/exon_sums/human.exon_sums.G026.gtf.gz"
zcat g026_gene.gtf.gz | head -3
这一步常被跳过,但它对做区间合并至关重要:不同注释版本里同一个基因的区间可能不同,做跨注释对比时必须先把区间对齐。GTF 里的属性列还需单独解析,直接从文本里按字符串找基因名并不可靠,因为同名基因可能出现在不同的转录本行上。另一个常被忽略的细节是染色体命名:同一套坐标里,有的文件写 chr1,有的写 1,直接按字符串合并会让整条染色体的计数对不上,动手前应先统一命名前缀。
§6 与库内条目的关系
6.1 关系图谱
| 库内相关词条 | 与本词条的关系 | 一句话区隔 |
|---|---|---|
| GTEx(人类组织基因表达调控图谱) | 上游来源之一 | GTEx 是一个联盟自己测的一批样本、自己发布的图谱;recount3 只是把它的读段拿来回炉重算,并与其他来源并列 |
| NCBI SRA(测序档案) | 最大上游来源 | SRA 存原始读段、是管道源头;recount3 存处理后的汇总、是管道出水口 |
| TCGA(癌症基因组图谱) | 上游来源之一 | TCGA 以癌种为单位组织临床与分子数据;recount3 只取其读段,按癌种切分后重算计数,不携带临床终点 |
| ENCODE(功能元件百科全书) | 平行但不同取向 | ENCODE 追求调控元件的广度与实验设计;recount3 追求存量转录组的可比性,覆盖的是常规转录组而非功能实验 |
| ARCHS4(若已收录) | 功能最接近的同侪 | ARCHS4 同样做公共数据统一重处理,但输出集中在基因级大规模矩阵;recount3 多出剪接位点与覆盖层,注释选择更多 |
| Expression Atlas(若已收录) | 功能相近,选样逻辑不同 | Atlas 以人工挑选的少数高质量研究为主,样本量小但注释好;recount3 几乎全收,样本量大但注释弱 |
| Single Cell Expression Atlas(若已收录) | 在单细胞侧部分重叠 | 后者专做单细胞实验汇编;recount3 的单细胞部分来自 Smart-seq 类全长协议,不覆盖液滴法 3’ 端数据 |
关系可以用一句话概括:SRA 是原料仓,GTEx 与 TCGA 是两个高质量专区,recount3 是这三处原料进同一台机器加工后的成品仓。
6.2 检索路径
任务不同,入口不同。做「跨研究找组织」去找元数据清单;做「某个基因在什么样本里表达」去找查询接口;做「大规模预训练语料」去找 R 包或直接遍历项目矩阵;做「剪接事件的新发现」去找剪接位点文件。
| 你的问题 | 该去哪个词条或入口 | 理由 |
|---|---|---|
| 想找某个 GEO 系列对应的原始数据 | SRA 词条 | recount3 以 SRP 或 ERP 为项目单元,不直接用 GSE 编号 |
| 想看正常组织的表达参考 | GTEx 词条 | GTEx 的词条讲清了组织与供体结构,recount3 只是重算 |
| 想看肿瘤的表达与临床关联 | TCGA 词条 | recount3 不带临床字段 |
| 想要跨研究可比的基因计数 | 本词条 | 这是它独占的能力 |
| 想要单细胞液滴法的数据 | 单细胞表达图谱类词条 | recount3 的长尾单细胞是 Smart-seq 类全长协议 |
| 想要原始 FASTQ 自己做流程 | SRA 词条 | recount3 不提供读段 |
§7 质量评估与局限性
7.0 阅读本节前的一个前置校正
看评分之前要先校正一个误判:很多人会拿「样本数七十余万」去对标「队列研究几十万人」,得出「这是史上最大的转录组队列」的结论。这个类比是错的。队列研究的样本数是独立个体数,而这里的计数单位是测序 run,同一个体可以被测多次、同一细胞系可以有多个批次、同一个项目可能对同一批样本做了重复测序。
正确的读法是把它看成「上万份独立研究、七十余万次测量」的集合。这个集合的价值在于广度与口径统一,不在于统计功效。任何按个体数解释它的分析都会高估自己手上的独立信息量。§4.4 的案例已经演示了单条 run 的追溯路径,可以看到它最终的落点是「某个项目里的一次测量」,而不是「某个人的一次访视」。
7.1 五维评分
DAIMS 框架下的五个维度分别打分,满分五分,评分依据是官方文档与论文能支撑的事实,而不是使用体验。
| 维度 | 分数 | 依据 |
|---|---|---|
| 覆盖广度 | 5 | 人类与小鼠双物种、三个来源、上万项目,同类中极少有对等规模 |
| 口径一致性 | 5 | 全量共用一套冻结版本的工具链,这是它存在的理由本身 |
| 注释质量 | 3 | 元数据来自自动抽取与预测,组织与类型字段缺失错标普遍,需回源核对 |
| 粒度完整性 | 4 | 基因、外显子、剪接、覆盖四层齐备,缺转录本层 |
| 可获得性 | 5 | 明文矩阵 HTTP 直取,无申请、无密钥、无白名单 |
总分二十四分(满分二十五)。这个分数里最低的一项是注释质量,而它恰好是使用者最容易踩坑的地方。
7.2 三个可以自己验证的点
第一,验证项目清单与文件是否存在。清单里有项目代号,不等于对应文件一定可取。改写脚本批量探测时会发现少数项目只部分可下载,原因通常是处理失败或尚未完成,这一点官方文档没有逐条说明。
第二,验证矩阵列数与元数据行数是否一致。同一个项目的基因矩阵列数是样本数,元数据表的行数也应是样本数,但两者排序不一定相同。用 accession 对齐而不是行号,是最稳的做法。
第三,验证注释版本带来的差异。同一个项目取两套注释的基因矩阵,行数会不同、同名基因的数值也可能不同。这不是错误,而是注释定义差异的直接体现;做跨版本对比时必须在方法里说明用的是哪一版。
7.3 已知的、可枚举的偏差
| 偏差 | 表现 | 后果 |
|---|---|---|
| 提交自选择 | 只有被提交到公共档案的研究才在集合里 | 阴性结果与工业界数据系统性缺席 |
| 地域失衡 | 提交量以美国与欧洲为主 | 特定人群的表达多样性代表性不足 |
| 单细胞混编 | 约一半 run 来自全长单细胞协议 | 不分层直接合并会产生伪差异 |
| 注释自动抽取 | 组织与类型字段有缺失与错标 | 按组织检索会漏检与误检 |
| 无默认过滤 | 低质量样本全部保留 | 不做 QC 粗筛会显著拖累下游 |
| 批次冻结 | 停在某个处理时点 | 新提交样本缺席,时间上不可比 |
| 物种注释不对称 | 小鼠只有一套注释 | 跨物种对比时注释分辨率不对等 |
7.4 外部验证矩阵
| 验证角度 | 可用手段 | 预期结果 |
|---|---|---|
| 与原作者计数对比 | 取同一项目的原研究补充矩阵 | 数值高度相关但非完全一致,差异来自流程 |
| 与 GTEx 官方矩阵对比 | 取 GTEx v8 的同组织样本 | 相关系数应在高水平,但批次处理差异会留下残差 |
| 与独立重处理资源对比 | 找另一个统一重处理的公共资源 | 同源码的结论可复现,跨源则有系统性偏移 |
| 内部复现 | 同注释同项目重复下载 | 逐字节一致 |
| 跨注释验证 | 同项目两套注释的基因级结果 | 基因集合不同、共有基因趋势一致 |
| 对照序列验证 | 查 ERCC 与 SIRV 行的计数 | 应随掺入量呈梯度,可用于判断定量下限 |
7.5 它为何能持续:把一次性重跑变成可复用的坐标系
这套资源能长期存在,靠的不是不断更新,而是把自己变成了一个别人离不开的参照系。一旦足够多的方法学论文、基准评测与下游工具都拿它的计数当基准,新的工作就必须与它对齐;对齐的成本又低于自己重跑一遍,于是它被反复引用、反复依赖。这种「参照系锁定」是基础设施类资源的典型生存逻辑,也是它即使更新缓慢也不会被替代的原因。
§8 生态位与历史
8.1 从 Rail-RNA 到 Monorail 的三代演进
这条技术线不是从 recount3 开始的。第一代是 recount2,用 Rail-RNA 做全量重处理,覆盖人类与小鼠的公共 RNA-seq,重点在基因与剪接位点的联合汇总,并首次把覆盖曲线做成可复用的资源。第二代是 Snaptron,把剪接位点做成可查询的索引,让「按区间查」成为可能。第三代才是 recount3,把前两代的能力合并,用新写的 Monorail 重跑更大的集合。
三代之间的关键变化在于规模与工具链。Rail-RNA 时代样本量在数万量级,重跑一次的成本可控;到了 SRA 爆发之后,样本量上到数十万,旧流程的单位成本撑不住了。论文里有一处对比说明了这个压力:新流程的单条 run 成本比旧流程低了一个数量级以上。成本降下来,才有底气把「全部重跑」当成常规动作。
8.2 与同期统一重处理项目的分工
2021 年前后,公共转录组统一重处理是一个热点,至少有几家在做类似的事。它们的差别主要在选样策略与输出粒度。
一类是「精选型」:人工挑选高质量的少数研究,重点保证注释准确,样本量不大。另一类是「全收型」:把档案里的数据几乎全部纳入,重点保证广度,注释靠自动抽取。recount3 明确站在全收型这一侧,并在论文里把自己与精选型资源做了对照,指出后者的局限是「只有依赖注释的基因级摘要」。
这个分工对使用者是好消息:需要注释可靠时选精选型,需要覆盖面时选全收型,两者并不冲突,甚至可以互补——用精选型做验证集,用全收型做发现集。
8.3 生态里的下游工具
围绕这套数据长出了一批下游工具。取数侧有 R 包与查询包两个入口,前者面向项目级分析,后者面向区域级查询。分析侧有专门处理这套输出格式的差异分析与可视化包,把行名对齐、协变量提取这些琐事封装起来。处理侧则是同一套管线本身,允许把自己的数据跑成同样的格式。
这个生态的规模不大但结构完整:取数、分析、复现三个环节各有工具,且都指向同一套文件格式。这种「格式即接口」的设计,是它能在不同语言、不同团队之间流通的原因。
8.4 在中国语境下的可用性
从国内访问这套资源,主要障碍是传输而非权限。文件托管在境外学术服务器上,全部是公开可读,不需要任何账号。下载速度取决于链路,大项目的基因矩阵往往在几十到上百兆,剪接与覆盖层则动辄数吉字节。
实践上通常的做法是先取元数据与清单,在本地做筛选,再只下载筛出来的那几十个项目。全量镜像在多数场景下没有必要,也没有必要——这套数据的价值在「按需取用」,不在「攒一个本地副本」。若确实需要长期离线使用,建议只镜像基因级矩阵,覆盖层按需取。
§9 使用指引
9.1 三种角色的路径
做机制研究的人:从元数据出发,先按组织与实验类型筛出候选项目,再按 QC 字段剔除异常样本,最后取基因级与外显子级矩阵做差异分析。核心风险是批次效应,务必把项目作为协变量。一个稳妥的顺序是先在小项目上验证假设、再把同样的假设放到大集合上复核,两次结果方向一致才继续推进,否则多半是批次而非生物学在起作用。
做工具与基准的人:从剪接位点与覆盖层出发,用未注释位点构造数据驱动的外显子集合,用 ERCC 与 SIRV 行评估定量下限。核心优势是不必自己重跑公共数据,直接拿现成的统一结果做对照。写评测时要把「同源码」与「跨源码」分开报告,前者检验实现,后者才是真实场景。
做模型预训练的人:从清单出发做规模抽样,先按协议类型分层,再按量级分层,避免单细胞与批量混编。核心代价是文本格式的解析开销,建议先在本地转成列式格式再训练。抽样比例应当记录下来,因为不同分层比例会系统性地改变模型学到的表达先验。
9.2 一条最短的可用流程
从零到一张可分析的表达矩阵,最短路径是四步:查项目、取矩阵、取注释、对齐后归一化。
#-- 从项目清单到可分析矩阵的最短流程
import gzip, io, urllib.request
import numpy as np
BASE = "https://duffel.rail.bio/recount3"
def fetch(path):
return gzip.decompress(urllib.request.urlopen(f"{BASE}/{path}").read()).decode("utf-8")
#-- 第一步:列出人类 SRA 的项目代号与来源
listing = fetch("human/data_sources/sra/metadata/sra.recount_project.MD.gz")
rows = [ln.split("\t") for ln in listing.splitlines()[1:]]
projects = sorted({r[3] for r in rows if len(r) > 5})
target = "SRP009615"
assert target in projects
#-- 第二步:取该项目的基因级矩阵
mat_txt = fetch(f"human/data_sources/sra/gene_sums/{target[-2:]}/{target}/"
f"sra.gene_sums.{target}.G026.gz")
lines = [ln for ln in mat_txt.splitlines() if not ln.startswith("##")]
header = lines[0].split("\t")
samples = header[1:]
counts = np.array([[int(v) for v in ln.split("\t")[1:]] for ln in lines[1:]], dtype=np.int64)
print(counts.shape)
#-- 第三步:按文库大小做最简单的每百万缩放,跨样本可比的最低要求
lib = counts.sum(axis=0)
cpm = counts / lib * 1e6
print(np.median(lib), cpm.shape)
这段代码不做任何过滤与批次校正,只保证跑得通。真实分析里,第四步应当是把 cpm 换回整数计数、交给差异分析工具处理,归一化的具体形态由工具决定。
9.3 引用方式
学术引用应引论文,不该引站点。论文的作者列表较长,共十五人,第一作者 Wilks,通讯作者之一为 Langmead。正式引用格式在官方文档的引用章节给出了 BibTeX 条目,可直接取用。
软件引用应引两个 R 包各自的 DOI,因为包的接口与版本是独立维护的。数据引用则应同时给出论文 DOI 与所取文件的具体路径与生成时间戳——文件头部的 ##date.generated 那一行,是唯一能证明你取了哪个版本的信息。
§10 AI 使用声明
| 项目 | 说明 |
|---|---|
| 数据性质 | 公共测序档案的二次加工汇总,不含新采集的受试者数据 |
| 训练可用性 | 矩阵为整数计数,可直接作为回归或分类目标,但需按协议与项目分层 |
| 标签可靠性 | 组织与类型标签来自自动预测,存在缺失与错标,不应当作金标准 |
| 隐私风险 | 汇总层不含个体标识,但 GTEx 与 TCGA 的原始数据仍受原条款约束 |
| 偏倚风险 | 提交自选择与地域失衡会传导到任何基于它的模型 |
| 建议用途 | 预训练、跨研究对照、剪接事件发现、方法学基准 |
| 不建议用途 | 临床决策、个体诊断、把样本数当个体数的统计推断 |
| 版本固定 | 引用时须记录文件生成时间戳,否则结果不可复现 |
常见误读(避坑清单)
| 常见说法 | 实际 | 证据 |
|---|---|---|
| 它有七十多万个病人 | 计数单位是测序 run,不是个体 | 论文表 1 的列名是 Sequencing Runs Processed |
| 它收集了新数据 | 全部输入来自 SRA、GTEx、TCGA 的既有读段 | 论文开头即说明是让归档数据可查询 |
| 它的矩阵已经归一化 | 交付的是原始整数计数,无 TPM 无文库标准化 | 矩阵头部只有注释代号与生成时间两行 |
| 同一项目里的样本可以直接比较 | 约一半 run 来自全长单细胞协议,须分层 | 论文正文明确写出「大约一半」 |
| 小鼠数据也支持多套注释 | 小鼠只有 GENCODE v23 一套 | 官方文档的小鼠注释章节只列一项 |
| 它是 GEO 的重新打包 | 项目单元是 SRA 的 SRP/ERP 编号,不按 GSE 组织 | 项目清单的 study 列取值形态 |
| 元数据的组织字段可以直接分组 | 字段来自文本抽取与预测,有缺失与错标 | 论文讨论部分承认元数据不可靠 |
| 它包含转录本水平的定量 | 输出止步于基因、外显子、剪接位点、覆盖 | 官方文档的数据章节只列这四类 |
| 它是实时更新的 | 冻结在某个处理批次 | 矩阵头部的生成时间戳固定不变 |
| 它只能通过 R 使用 | 全部文件是公开 HTTP 静态对象 | 官方文档的原始文件章节给出完整路径 |
FAQ
Q1:recount3 和直接下载 SRA 数据有什么区别?
SRA 存的是原始读段,比对、注释、过滤全要自己做,而且不同人的做法互不相同。recount3 已经用同一套管线把这些步骤固定下来,同一研究、不同研究之间的计数因此落在同一把尺子上。代价是你放弃了处理细节的控制权。
Q2:样本数到底是七十五万还是七十六万?
论文摘要写的是 over 750,000,正文表 1 的合计是七十五万八千二百,两个口径都出自同一篇论文。差异来自摘要取整,处理过程中也可能重新核对过部分 run。引用时写摘要口径更稳妥。
Q3:能不能直接拿 recount3 的基因计数做差异表达?
可以,但必须先做归一化和批次校正。recount3 交付的是原始整数计数,没有 TPM、没有文库大小标准化,也没有去除项目效应。同一项目内部的比较相对安全,跨项目比较需要自己引入协变量。
Q4:recount3 支持小鼠吗?
支持,而且小鼠的 run 数比人类更多。小鼠只有一条注释通路(GENCODE v23,基于 mm10),人类有 GENCODE v26、v29、RefSeq、FANTOM-CAT 等多套,跨物种比较时要记住这个不对称。
Q5:为什么说它一半的数据是单细胞?
论文原文写明大约一半 run 来自全长单细胞协议,主要指 Smart-seq 与 Smart-seq2。这类数据每一条 run 对应一个单细胞或极少量细胞,基因检出率与批量组织样本差距很大,混用会引入严重的技术批次。
Q6:它包含 GEO 的数据吗?
包含,但不是以 GEO 的系列编号组织。大量 GEO 研究在提交时同时上传到 SRA,recount3 以 SRA 的 SRP 或 ERP 编号为项目单元,所以同一份研究在两边都有,只是键不同。想按 GSE 查,需要先做一次编号映射。
Q7:能用它做单细胞分析吗?
能取到数据,但不能替代专门的单细胞资源。它只覆盖全长单细胞协议,没有液滴法的细胞条形码与 UMI 结构,单个细胞的表达矩阵在这里表现为一条独立的 run,而不是细胞乘基因矩阵。真正的单细胞分析应去专门的单细胞资源。
Q8:文件这么大,有没有更省的办法?
有。先用元数据把项目筛到几十个,再只取基因级矩阵;剪接位点与覆盖层按需取;如果只需要某个基因或区间,直接用查询接口而不是下载整项目。全量镜像在多数场景下没有必要。
Q9:为什么同一个基因在不同注释下数值不同?
注释版本决定了基因的边界定义。GENCODE v26 与 v29 对同一个基因的外显子边界可能不同,落在边界上的读段归属就会变化。这不是错误,是定义差异;跨版本对比时必须在方法里说明用的是哪一版。
Q10:它之后还会更新吗?
站点与 R 包持续可用,但与 SRA 的增量相比更新速度有限,处理批次的时间戳决定了它冻结在某个时点。需要最新样本时应回到 SRA,或自行用同一套 Monorail 管线处理,这样输出仍可对齐。
Q11:怎么确认自己取到的是哪个版本?
看矩阵头部的两行注释。第一行是注释代号,第二行是生成时间戳,精确到微秒。把这两行记进方法部分,别人才能复现你取到的确切文件。
术语表
| 术语 | 含义 |
|---|---|
| recount3 | 本词条描述的资源,公共 RNA-seq 的统一重处理汇总 |
| Monorail | 支撑该资源的分析管线,网格化调度,单 run 一单元 |
| STAR | 剪接感知的比对工具,本资源使用 2.7.3a 版本 |
| featureCounts | 按注释汇总读段到基因与外显子的工具,本资源使用 1.6.3 |
| MegaDepth | 生成剪接位点计数与覆盖曲线的工具,本资源使用 0.4.0 |
| run | 测序档案的最小记录单位,一次测序产出的读段集合 |
| rail_id | recount3 内部给每条 run 分配的整数标识 |
| accession | 外部编号,如 SRR 开头的 run 号、SRP 开头的项目号 |
| SRP / ERP | SRA 与 ENA 的项目级编号前缀,本资源的项目单元 |
| GTEx | 基因型-组织表达计划,本资源的人类来源之一 |
| TCGA | 癌症基因组图谱,本资源的人类来源之一 |
| GENCODE | 常用的基因注释集合,人类侧有多个版本 |
| RefSeq | NCBI 维护的参考序列与注释集合 |
| FANTOM-CAT | FANTOM 联盟的注释集合,强调长链非编码 |
| ERCC | 外源掺入的已知浓度 RNA 对照 |
| SIRV | 已知结构的合成剪接异构体对照 |
| bigWig | 压缩的基因组区间数值格式,用于存覆盖曲线 |
| 外显子-外显子剪接位点 | 跨越两个外显子接头的读段所代表的连接位置 |
| Matrix Market | 稀疏矩阵的文本交换格式,用三元组表示非零元素 |
| RangedSummarizedExperiment | Bioconductor 的区间化矩阵容器,R 包返回的对象类型 |
| snapcount | 面向区域查询的 R 包,可查本资源与上一代资源 |
| Snaptron | 支撑区域查询的服务端索引系统 |
| 注释无关 | 不依赖外部注释、直接由数据推出结构的处理方式 |
| CPM | 每百万计数,最简单的跨样本缩放方式 |
| 批次效应 | 由处理时间、平台、项目等非生物学因素引入的系统差异 |
| 全长单细胞协议 | 覆盖整条转录本的单细胞测序方案,如 Smart-seq 系列 |
| 受控访问 | 需申请与密钥才能取用原始数据的访问等级,如 dbGaP 与 GDC |
| 数据源与集合 | 数据源指原始托管位置,集合指按主题拼起来的一组研究 |
| 表达区域 | 由覆盖曲线数据驱动划出的转录活跃区段,不依赖注释 |
| 解码表 | 稀疏矩阵的行列说明文件,用于还原三元组对应的样本与坐标 |
附录 A:核心数字核对表
| 项 | 值 | 来源 | 来源层级 |
|---|---|---|---|
| 测序 run 合计 | 七十五万余(表 1 合计值) | 论文表 1 | 论文正文 |
| 摘要口径规模 | 七十五万以上 | 论文摘要 | 论文摘要 |
| 人类 SRA run | 316,443 | 论文正文 | 论文正文 |
| 小鼠 SRA run | 416,803 | 论文正文 | 论文正文 |
| SRA 来源 run | 732,246 | 论文正文 | 论文正文 |
| GTEx v8 run | 19,214 | 论文表 1 | 论文正文 |
| TCGA run | 11,348 | 论文表 1 | 论文正文 |
| 输入压缩读段 | 990.133 TB | 论文表 1 | 论文正文 |
| 输出摘要体量 | 146.31 TB | 论文表 1 | 论文正文 |
| 计算消耗 | 28,965 节点小时 | 论文表 1 | 论文正文 |
| STAR 版本 | 2.7.3a | 论文方法部分 | 论文正文 |
| featureCounts 版本 | 1.6.3 | 论文方法部分 | 论文正文 |
| MegaDepth 版本 | 0.4.0 | 论文方法部分 | 论文正文 |
| 人类注释套数 | 4 套主体加 2 套对照 | 官方文档 | 官方页面 |
| 小鼠注释套数 | 1 套 | 官方文档 | 官方页面 |
| TCGA 样本数 | 11,348(10,396 人,33 癌种) | 论文正文 | 论文正文 |
| 人类项目文件数 | 18,830 | 官方文档代码输出 | 官方页面 |
| GTEx 组织切片数 | 32 | 官方元数据清单 | 官方页面 |
| TCGA 癌种切片数 | 33 | 官方元数据清单 | 官方页面 |
附录 B:数字的来源分层
第一层,官方页面。 包括 rna.recount.bio 的文档章节、托管服务器上的元数据清单、矩阵文件头部。这一层提供的是结构信息与可核对的计数:注释套数、项目代号形态、文件命名规则、矩阵列数。凡是能从文件里直接数出来的,都归这一层。
第二层,论文摘要与正文。 包括 Genome Biology 那篇的摘要、表 1、方法部分。这一层提供的是规模信息与版本信息:run 总数、输入输出体量、计算消耗、各工具版本号。摘要与正文口径不一致时,以正文表 1 的精确值为准,摘要的取整值作为常用引用口径。
第三层,未核实。 凡是需要下载全量数据才能得出的统计,例如各组织的样本分布、各项目的样本量中位数、基因检出率的分布,均归入这一层,本篇不写。同样,任何依赖站点当前状态而非文档记载的数字,也归入这一层。
附录 C:完整参考文献
- Wilks C, Zheng SC, Chen FY, Charles R, Solomon B, Ling JP, Imada EL, Zhang D, Joseph L, Leek JT, Jaffe AE, Nellore A, Collado-Torres L, Hansen KD, Langmead B. recount3: summaries and queries for large-scale RNA-seq expression and splicing. Genome Biology. 2021;22:323. doi:10.1186/s13059-021-02533-6
- recount3 官方文档,数据章节与原始文件章节,rna.recount.bio。
- recount3 托管服务器目录清单,duffel.rail.bio/recount3。
- recount3 R/Bioconductor 包文档与快速上手教程。
- snapcount R/Bioconductor 包文档。
- Snaptron 服务端索引说明。
- GTEx 计划官方门户,本资源的人类来源之一。
- TCGA 项目说明,本资源的人类来源之一。
- NCBI SRA 文档,本资源的原始读段来源。
- GENCODE 注释发布说明,本资源的注释来源。
- STAR 比对工具论文与版本说明。
- featureCounts 汇总工具论文。
- MegaDepth 覆盖与剪接统计工具论文。
附录 D:常见任务的代码骨架
#-- 任务一:只取元数据,在本地筛选后再批量下载
cd /tmp
curl -sL -o proj.gz "https://duffel.rail.bio/recount3/human/data_sources/sra/metadata/sra.recount_project.MD.gz"
#-- 挑出样本数超过某个阈值的项目,仅打印代号
zcat proj.gz | awk -F'\t' 'NR>1{c[$3]++} END{for(p in c) if(c[p]>=100) print p}' | sort | head -20
#-- 任务二:批量探测项目矩阵是否存在,记录缺失清单
import urllib.request
BASE = "https://duffel.rail.bio/recount3"
def exists(organism, project, ann="G026"):
url = (f"{BASE}/{organism}/data_sources/sra/gene_sums/"
f"{project[-2:]}/{project}/sra.gene_sums.{project}.{ann}.gz")
req = urllib.request.Request(url, method="HEAD")
try:
with urllib.request.urlopen(req, timeout=30) as r:
return r.status == 200
except Exception:
return False
candidates = ["SRP009615", "SRP107565", "SRP066781"]
missing = [p for p in candidates if not exists("human", p)]
print("missing:", missing)
#-- 任务三:同一项目取两套注释,比较共有基因的相关性
library("recount3")
proj <- subset(available_projects(), project == "SRP009615" & project_type == "data_sources")
rse26 <- create_rse(proj, annotation = "gencode_v26")
rse29 <- create_rse(proj, annotation = "gencode_v29")
g26 <- assay(rse26, "counts")
g29 <- assay(rse29, "counts")
common <- intersect(rownames(g26), rownames(g29))
cor(log1p(g26[common, 1]), log1p(g29[common, 1]))
-- 任务四:把元数据导入 SQLite 后按组织统计样本数
CREATE TABLE proj (
rail_id INTEGER,
external_id TEXT,
study TEXT,
project TEXT,
organism TEXT,
file_source TEXT,
metadata_source TEXT,
date_processed TEXT
);
-- 导入时用 .import --skip 1 proj.tsv proj
SELECT organism, COUNT(*) AS n FROM proj GROUP BY organism ORDER BY n DESC;
SELECT study, COUNT(*) AS n FROM proj GROUP BY study ORDER BY n DESC LIMIT 10;
附录 E:如果只记五件事
第一,它不产生新数据。全部输入是 SRA、GTEx v8、TCGA 里的既有读段,用同一套管线重跑一遍,产出跨研究可比的计数。这一条决定了它的一切边界。
第二,规模口径有两个。摘要写的是七十五万以上,正文表 1 是七十五万八千二百。引用写摘要值,精确讨论写表 1 值。
第三,输出四层:基因、外显子、剪接位点、覆盖。没有转录本层。要转录本就得回原始读段。
第四,注释是软的,管线是硬的。元数据的组织与类型字段靠自动抽取,别当金标准;工具链版本被冻结,这才是它值钱的地方。
第五,约一半是全长单细胞。不分层直接合并,是这类数据上最常见也最致命的一次误用。
编辑说明
本篇的结构事实来自官方文档的数据章与原始文件章,规模数字来自 Genome Biology 论文的摘要与表 1,项目单元的具体形态来自托管服务器上的元数据清单。凡是需要下载全量数据才能得出的分布型统计,一律未纳入。
核验时间:2026 年 10 月。核验方式:读取官方文档页面、拉取项目级元数据清单、抽查单个项目的矩阵文件头部、比对论文摘要与正文的数字口径。未覆盖的内容包括:各组织的样本分布、各项目的质量分布、基因检出率的统计特征——这些都留给使用者在拿到数据后自行判断。
对数字口径的处理原则是:官方内部不一致时不做裁决,只并列陈述两个口径与各自的出处。摘要与正文的七十余万之差即按此处理。
相关数据集导航
以下为站内 AI-Ready 数据集百科中与本词条共享多个主题标签的相关数据集,按相关度降序排列:
- geo — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- fantom5 — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- gencode — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- astra — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- archs4 — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- ddbj-gea — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- gtex-legacy — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- gtex — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- tcga-luad — 共享标签:基因组学与多组学 / 转录组 / 测序数据
- cellxgene-census — 共享标签:基因组学与多组学 / 转录组 / 测序数据
导航说明:本章节由全站统一标签体系自动计算生成(标签重合度算法),双向可达;点击链接可跳转至对应数据集词条。

