2. RNA-seq 数据基础：从 FASTQ 到 count 矩阵

拿到 RNA-seq 项目，第一件事不是跑 DESeq2，是先把原始数据变成 count 矩阵。

没有 count 矩阵，后面所有差异分析和标记筛选都是空谈。

这一章，我就带你走一遍这个标准流程。

说白了，RNA-seq 数据分析有三个核心步骤：

  • 原始数据：测序仪下机的 FASTQ 文件
  • 比对/定量：把 reads 定位到参考基因组或转录本
  • count 矩阵：汇总每个基因被多少条 reads 覆盖
RNA-seq 数据流程：从 FASTQ 到 count 矩阵 RNA-seq 数据流程：从 FASTQ 到 count 矩阵 原始 FASTQ 文件 质量检查 + 过滤 (FastQC/Trimmomatic) 比对到基因组 STAR / HISAT2 直接定量到转录本 Salmon / kallisto count 矩阵 核心思路：将测序片段 → 基因表达定量 → 形成行是基因、列是样本的计数表

2.1 FASTQ 文件，长什么样？

FASTQ 是测序仪最原始的产出格式。

每条 read 有 4 行：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
GTGCTTGCATCTGATCTGCTGTGCCTGACCGTCAGTGTGTG
+
BBBB@B@BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
  • 行 1：序列标识符，@ 开头
  • 行 2：碱基序列
  • 行 3：分隔符，通常就一个 +
  • 行 4：质量评分，每个字符对应一个 Phred 分数

我的一点习惯：拿到新项目，我会先跑一遍 FastQC 看看质量分布。有一次项目数据质量很差，3' 端全是低质量碱基，如果不提前修剪，后面比对率会非常低。

注意：FASTQ 的质量编码有 Phred+33 和 Phred+64 两种。现在绝大多数数据是 Phred+33，但老数据可能是 +64。FastQC 能自动识别，不过如果你写脚本自己解析，一定要确认。

2.2 比对，到底比到哪？

有两种主流策略，我个人都接触过，各有适用场景。

策略代表工具输出典型场景
有参比对到基因组STAR、HISAT2BAM 文件有高质量参考基因组
直接定量到转录本Salmon、kallisto量化的 count/TPM速度快，无需完整基因组

STAR 是我用得最多的比对工具。它的核心是「种子搜索 + 拼接」两步，速度非常快。

一个典型命令长这样：

STAR --genomeDir /path/to/genome_index \
    --readFilesIn sample_R1.fastq sample_R2.fastq \
    --outFileNamePrefix sample_ \
    --outSAMtype BAM SortedByCoordinate

如果你用 featureCounts 做定量，它会直接读取 BAM 文件，统计落到每个基因外显子区域的 reads 数量。

命令非常简单：

featureCounts -T 8 -t exon -g gene_id \
    -a annotation.gtf \
    -o counts.txt \
    sample_sorted.bam

-T 8 用 8 个线程，-t exon 指定统计 exon 特征，-g gene_id 按基因 ID 汇总。

一个坑：featureCounts 默认用 --primary 只统计 primary 比对。我曾经遇到过一批 BAM 文件没经过 properly paired 过滤，结果计数偏高。后来加上了 -s 0（不区分链）重新跑才正常。

2.3 count 矩阵的三种生成方式

我把它分成三类，你根据项目情况选一个就行。

方式一：HTSeq-count

Python 写的，基于 GTF 注释，逐条 reads 判断落在哪个基因上。

htseq-count -f bam -r pos -s no \
    sample_sorted.bam \
    annotation.gtf > sample_counts.txt
  • -s no：不区分链特异性
  • -r pos：BAM 已按位置排序
  • 优点是结果稳定，缺点是速度慢

方式二：featureCounts

subread 包里的 C 程序，速度比 HTSeq 快很多。

featureCounts -T 4 -t exon -g gene_id \
    -a annotation.gtf \
    -o all_counts.txt \
    sample1.bam sample2.bam sample3.bam

一次跑多个样本，输出一个合并的 count 矩阵。

我的建议：项目样本数超过 20 个，我肯定选 featureCounts。HTSeq 在样本少、需要仔细排查比对质量时用。

方式三：STAR --quantMode GeneCounts

STAR 在比对的同时可以直接输出 count 数量，不需要额外工具。

STAR --genomeDir index \
    --readFilesIn R1.fastq R2.fastq \
    --quantMode GeneCounts \
    --outFileNamePrefix sample_

它会生成一个 ReadsPerGene.out.tab 文件，第四列是 unstranded 的 count。

三种方式对比：

工具速度内存多线程多样本合并
HTSeq-count慢低否需手动
featureCounts快中等是一次性输出
STAR GeneCounts极快高是每样本独立文件

注意：STAR 的内置定量模式没有多映射 reads 的校正。如果你的项目是重复序列富集的物种（比如植物），我还是建议用 featureCounts 配合 -M --fraction 参数。

2.4 从多个样本到一张 count 表

无论用哪个工具，最终要得到一个矩阵。行是基因 ID，列是样本名，单元格是整数。

一个标准的 count 表长这样：

gene_id control_1 control_2 treatment_1 treatment_2
ENSG000001 120  145  200  180
ENSG000002 50   60   22   30
ENSG000003 900  850  1100 1050

注意：count 必须是整数，不能是归一化值。DESeq2 自己会做内部归一化，所以给原始整数就行。

一个实用操作：featureCounts 输出文件第一行以 # 开头，第二行是列名。导入 R 时用 read.table(..., skip=1) 跳过注释行。

# R 代码示例
counts <- read.table("all_counts.txt", header=TRUE, row.names=1, skip=1)
counts <- counts[, -c(1:5)]  # 去掉注释列，只保留 count 列

2.5 本章小结

好了，你应该清楚了：

  • FASTQ 是原始数据，要检查质量
  • 比对或定量是关键步骤，STAR + featureCounts 是经典组合
  • count 矩阵是 DESeq2 的输入基础
← 上一章 📖 返回目录 下一章 →