Bulk RNA-seq 分析入门（一）：从实验设计、FASTQ 到可信的计数表

Perry2026-10-062026-10-06

RNA-seq 分析的起点不是某条比对命令，而是一句话：我想比较什么样的生物学对象，比较能够排除哪些其他解释？先把这句话说清楚，后面的文件、模型与图才会围绕同一个问题组织起来。

这是三篇连续教程的第一篇，讨论常规 bulk RNA-seq：每个文库代表一个组织、细胞群体或其他整体样本。它不是单细胞教程，不能直接套到带有细胞条形码和 UMI 的单细胞矩阵上。本文中的六样本设计是原创的假设练习；命令只是准备好真实数据之后的执行模板，写作过程中没有在本机下载 FASTQ、安装软件或运行测序流程。

读完这一篇，应该能交付样本表、参考文件清单、处理参数、质控记录，以及来源明确的计数或定量结果。下一篇再用官方公开的小数据集练习差异表达，而不是让入门电脑先承担完整人类基因组比对。

1. 先画出实验单位，而不是先数文件

假设我们研究某种培养细胞在处理前后是否出现转录变化。三个独立培养来源进入对照组，另三个独立培养来源进入处理组，这里共有六个生物学重复。若同一个培养来源被分装后重复建库，重复的文库提供的是技术层面的信息，不能自动变成多个独立生物学重复。测了两个 lane、得到 R1 和 R2、把一个样本重复测序，也都不会把生物学样本数量翻倍。

我会先在纸上给每个样本写四项：它来自谁、接受什么处理、在哪个批次建库、与其他样本共享什么来源。这样可以提早发现两个危险情况：两组其实各只有一个独立来源；或者所有对照都在第一批处理、所有处理都在第二批处理。第一种缺少组内变异信息，第二种把处理与批次完全绑在一起。再复杂的统计软件也不能凭空补出缺失的实验信息。

“每组至少三个”可以作为入门设计练习的起点，却不是通用功效保证。正式实验应根据预期效应、变异、研究目标和资源做功效评估，并尽量在处理、提取、建库与测序环节平衡或随机化分组。技术重复的处理方式以及配对设计也要提前记录，不应看到结果以后才决定。关于实验单位与伪重复，可参阅 Lazic 等人的原始方法文章。

2. 一张表给流水线，一张表给科学问题

流水线样本表负责告诉程序文件在哪里；生物学元数据负责告诉分析者这些文件代表什么。两张表通过稳定的 sample_id 连接，不能靠“应该是这个顺序”连接。

下面是六个假设样本的元数据。subject_id 在这里表示独立来源，不是患者姓名；所有字段都只是教学标签。B1 与 B2 都有两个条件，但分配并不完全均衡，正式设计还可以进一步改善。

sample_id,subject_id,condition,batch
C1,demo_01,control,B1
C2,demo_02,control,B2
C3,demo_03,control,B1
T1,demo_04,treated,B1
T2,demo_05,treated,B2
T3,demo_06,treated,B1

本系列固定参考 nf-core/rnaseq 3.27.0 文档，不用“最新版”代替版本记录。该版本 FASTQ 样本表必需的四列是 sample,fastq_1,fastq_2,strandedness。双端填写 R1/R2，单端的 fastq_2 留空。同一样本多个 lane 可以重复使用同一个 sample 标识，但不同生物学样本不能因此被意外合并。未知链特异性可设 auto，并在输出后复查推断与协议是否一致。3.27.0 官方样本表说明

sample,fastq_1,fastq_2,strandedness
C1,/data/fastq/C1_R1.fastq.gz,/data/fastq/C1_R2.fastq.gz,auto
C2,/data/fastq/C2_R1.fastq.gz,/data/fastq/C2_R2.fastq.gz,auto
C3,/data/fastq/C3_R1.fastq.gz,/data/fastq/C3_R2.fastq.gz,auto
T1,/data/fastq/T1_R1.fastq.gz,/data/fastq/T1_R2.fastq.gz,auto
T2,/data/fastq/T2_R1.fastq.gz,/data/fastq/T2_R2.fastq.gz,auto
T3,/data/fastq/T3_R1.fastq.gz,/data/fastq/T3_R2.fastq.gz,auto

/data/fastq/ 是需要自行替换的路径，博客并未提供这些 FASTQ。真实元数据还应携带组织、采样时间、建库方式、poly(A) 富集或 rRNA 去除、读长、性别等与研究有关的字段。只填写实际已知的信息，不要用猜测补空格。涉及人类数据时，先确认数据使用范围；未经允许不要把 FASTQ、真实身份映射或临床字段上传到云服务及公开仓库。

3. 完整小练习：仅检查设计，不处理测序数据

下面的 Python 标准库代码可以单独保存为 check_rnaseq_design.py，使用 python check_rnaseq_design.py 执行。它只读取代码里的教学 CSV，不访问磁盘中的实验数据、不安装包、不生成差异基因。

import csv
import io
from collections import Counter, defaultdict

text = """sample_id,subject_id,condition,batch
C1,demo_01,control,B1
C2,demo_02,control,B2
C3,demo_03,control,B1
T1,demo_04,treated,B1
T2,demo_05,treated,B2
T3,demo_06,treated,B1
"""
rows = list(csv.DictReader(io.StringIO(text)))
required = {"sample_id", "subject_id", "condition", "batch"}
if not rows or not required.issubset(rows[0]):
    raise ValueError("Missing metadata fields")
if any(not row[field].strip() for row in rows for field in required):
    raise ValueError("Blank metadata field")
ids = [row["sample_id"] for row in rows]
if len(ids) != len(set(ids)):
    raise ValueError("Duplicate sample_id")
groups = defaultdict(set)
batches = defaultdict(set)
for row in rows:
    if row["condition"] not in {"control", "treated"}:
        raise ValueError("Unexpected condition")
    groups[row["condition"]].add(row["subject_id"])
    batches[row["batch"]].add(row["condition"])
if set(groups) != {"control", "treated"}:
    raise ValueError("Both conditions are required")
if any(len(subjects) < 3 for subjects in groups.values()):
    raise ValueError("This teaching design requires 3 independent sources per group")
if any(n > 1 for n in Counter(row["subject_id"] for row in rows).values()):
    raise ValueError("Repeated source: review pairing or technical replication")
if all(len(conditions) == 1 for conditions in batches.values()):
    raise ValueError("Batch and condition may be fully confounded; review design")
print("Teaching design check passed")
print("Independent sources:", {k: len(v) for k, v in groups.items()})
print("Conditions represented in each batch:", dict(batches))

这段检查故意只接受“每个来源仅出现一次”的教学非配对设计。真实配对数据会被它要求人工复查，而不是被直接判定为坏数据。批次检查也是保守提示，不代替正式设计矩阵的秩检查。第二篇会展示配对建模与完整秩检查。

可以做三个自检：把 T3 的来源改成 T2 的来源；把全部 C 样本放到 B1、全部 T 样本放到 B2；把 C2 的 ID 改成 C1。观察程序在哪一步停止，并用自己的话解释停止原因。这个练习的价值是让“样本表错误”在昂贵的运算开始前被发现。

4. 参考基因组与注释必须是一套可追溯的东西

FASTA 告诉比对器参考序列是什么，GTF 告诉计数或转录本构建过程哪些区间属于哪些基因和转录本。二者不能只因为都写着“human”就视为匹配。需要核对物种、组装版本、注释发布版本、染色体命名、是否包含替代 contig，以及转录本 ID 是否与定量索引一致。

一个经常被忽略的错误是：比对使用带 chr 的 FASTA，注释却使用不带 chr 的 contig；另一个是用某个发布版本的转录组建索引，却拿另一个版本的转录本到基因对应表做汇总。程序有时可以继续运行，但有些读段或转录本已经悄悄失去了正确对应关系。不要在不理解来源的情况下统一删掉 ID 后缀，或者将缺失映射全部填成同名基因。

我建议保存一份参考清单，列出来源网址、下载日期、组装版本、注释版本、文件校验值、索引构建命令及工具版本。后面换参考文件时，把它当成一次新的分析决定，而不只是“更新了文件”。相关参数可在 nf-core/rnaseq 3.27.0 参数页逐项核对。

5. 从 FASTQ 到计数，各步骤究竟改变了什么

可以把流程记成一条有审计点的链：原始 FASTQ → 原始读段检查 → 必要的接头处理 → 比对或转录本定量 → 样本汇总 QC → 基因级输入检查。每个箭头都应该能回答“输入是什么、参数是什么、输出在哪里、为什么允许进入下一步”。

STAR 基因组比对保留读段的基因组位置，方便检查剪接、比对分布和其他 BAM 层面的指标；Salmon 估计转录本丰度，再依据转录本到基因的映射汇总。它们不是“哪个名字更高级”的关系，选择要对应目标、资源和实验协议。普通初学者可以先使用固定的标准工作流，掌握每个产物，再决定是否需要改变流程。nf-core 的比对和定量选项

下面模板选择 star_salmon。这不是对本机资源的可运行承诺。执行前需要可用的 Nextflow、与该版本要求匹配的 Java、Docker 或经管理员批准的容器运行环境、足够的 CPU/内存与磁盘，并确认镜像和参考下载权限。人类全基因组 STAR 工作不能按“普通小脚本”估计内存；具体资源应按参考、索引与集群政策规划。Nextflow 安装与运行前提以官方说明为准。

# Linux/HPC 执行模板：先替换路径并完成环境、权限和资源检查。
# 会读取真实FASTQ并产生大量文件；不是本篇已执行的练习。
nextflow run nf-core/rnaseq \
  -r 3.27.0 \
  -profile docker \
  --input /data/project/metadata/samplesheet.csv \
  --fasta /data/reference/genome.fa \
  --gtf /data/reference/annotation.gtf \
  --aligner star_salmon \
  --outdir /data/project/results/rnaseq_v1 \
  -work-dir /data/project/work/rnaseq_v1 \
  -with-report /data/project/nextflow_report.html \
  -with-trace /data/project/nextflow_trace.tsv \
  -with-timeline /data/project/nextflow_timeline.html

首次运行应使用新的结果目录，记录完整命令与退出状态。在 HPC 上不要在登录节点启动大任务；采用站点批准的执行器和资源配置。-resume 是恢复符合缓存条件的任务，不是自动确认旧参考、旧参数与新实验相容。正式分析还应记录实际 Nextflow、容器版本和镜像标识；只有 -r 固定了流程版本，其他外部配置也需要追踪。

6. QC 不是追求所有按钮都变绿

第一层看原始读段：数量、长度、质量下降位置、接头、过度代表序列。第二层看处理后信息：比对率、唯一/多重比对、rRNA 或其他污染、链特异性、插入片段、基因体覆盖与定量分配。第三层看样本之间：同一实验的样本是否有异常深度、异常分布，是否与建库记录一致。

RNA 读段来自不均匀表达的转录本，重复序列的解释不能机械照搬均匀覆盖的 DNA 测序。接头、低复杂度、过度扩增等也需要结合协议与其他证据区分，不要因一个警告就随意删掉全部重复读段。FastQC 提供的是待解释的线索，而不是替分析者完成决定。FastQC 官方概览、官方过度代表序列说明

MultiQC 把多个工具的报告集中展示，方便发现某个样本与同批其他样本不一致；它并不会替代实验协议、阈值依据和异常原因调查。MultiQC 官方说明

不要在这里给自己设置一个跨物种、跨协议的万能比对率阈值。应预先约定可接受范围和调查流程：发现问题后先查文件配对、参考、链特异性与污染，再结合原始记录判断是否重做。确有失败文库可以排除，但理由应与预先约定的 QC 依据相符，不能只因为排除后分组更好看。

7. 什么可以交给 DESeq2，什么不可以

基因整数计数、Salmon 的估计计数、TPM、VST 是不同对象。整数计数可以通过计数矩阵入口进入 DESeq2；Salmon 转录本估计通常用 tximport 或 tximeta 的专门入口，保留长度相关信息。TPM、FPKM、对数表达或 VST 值不能直接冒充原始计数。把 TPM 四舍五入也不会把它变回测序计数。DESeq2 输入说明

查看 Salmon 结果时，至少保留原始每样本 quant.sf、定量索引来源与 tx2gene 映射。计数矩阵必须附带生成方法；不能从几个名为 counts 的文件里凭名字选一个。nf-core 的输出包含原始估计计数、TPM、长度和不同缩放产物，选择下游入口时需要查看准确含义，而不是把所有列统称为“表达量”。3.27.0 输出文件说明

最终的交接包应含样本元数据、参考清单、原始计数或完整定量文件、流程与软件版本、QC 报告、异常处理记录，并明确哪些样本允许进入差异分析。若生物重复不足、关键标签存疑或处理与批次完全混杂，正确结果可能是“暂不进行差异检验”，不是强行产出火山图。

8. 本篇结束时的自查清单

我能说清实验单位，并区分生物学重复、技术重复、lane 和双端读段。
两张样本表使用相同稳定 ID，参考、索引和注释有精确来源。
链特异性有协议依据或经过推断复核，而不是从文件名猜测。
我知道每个输出是计数、估计丰度还是转换值，没有将 TPM 交给计数模型。
每次排除与重跑都有原因、版本和日志，可以追溯到原始输入。

下一篇：DESeq2 设计公式、对比与可复查的差异表达练习。