RNA-seq 数据预处理：从 fastq 到 count 矩阵

说真的，这步要是搞砸了，后面用 DESeq2 还是 edgeR 都白搭。我自己刚入行时就吃过这个亏——拿着漂亮的 p-value 列表兴奋半天，结果发现是上游比对参数配错了。嗯，咱们从头捋一遍，把原始测序数据变成能用的 count 矩阵。

整体流程速览

把 fastq 文件变成表达矩阵，说白了就是四步：

  1. 质量控制 – 看看原始数据干不干净
  2. 修剪/过滤 – 把接头、低质量碱基砍掉
  3. 比对 – 把 reads 贴到参考基因组上
  4. 定量 – 统计每个基因有多少 reads 落上去

每一步都有对应的主流工具。我习惯先画张图把逻辑理清楚。

RNA-seq 数据预处理流程 原始 fastq 质量控制：FASTQC + MultiQC 修剪/过滤：fastp / Trimmomatic 比对：STAR / HISAT2 定量：featureCounts / htseq-count count 矩阵 常见问题： • 链特异性？STAR用--outSAMstrandField • 多比对reads怎么处理？ • 是否需要去重？ • 低质量碱基修剪标准

第一步：质量控制与修剪

我拿到 fastq 文件后，第一件事就是跑 FASTQC。别偷懒，直接对所有样本跑一遍，再用 MultiQC 汇总看一眼整体趋势。

fastqc *.fastq.gz -o ./qc_raw
multiqc ./qc_raw -o ./multiqc_report

看过报告后，你通常会决定要不要修剪。我个人偏好 fastp，因为它快，而且能自动检测接头。这是典型命令：

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o sample_clean_R1.fastq.gz -O sample_clean_R2.fastq.gz \
      --detect_adapter_for_pe \
      --trim_poly_g \
      --cut_front --cut_tail --cut_mean_quality 20 \
      --length_required 36
我的经验：遇到过用 Trimmomatic 切得太狠，导致很多 reads 变太短，反而浪费了数据。可以先用 fastp 的默认参数试跑一个样本，看看看通过率再批量跑。

第二步：比对到参考基因组

比对这一步，选对工具很关键。我个人一直用 STAR，因为它快，而且对剪接 junction 支持好（对于真核生物是必须的）。Hisat2 也可以，但速度差一些。建索引是前提：

STAR --runMode genomeGenerate \
     --genomeDir ./star_index \
     --genomeFastaFiles hg38.fa \
     --sjdbGTFfile hg38.gtf \
     --runThreadN 16

比对命令：

STAR --genomeDir ./star_index \
     --readFilesIn sample_clean_R1.fastq.gz sample_clean_R2.fastq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outSAMunmapped Within \
     --outSAMattributes Standard \
     --runThreadN 16
注意：如果你做的是链特异性文库（比如 dUTP 法、SMARTer），一定要在 STAR 里加 --outSAMstrandField intronMotif，否则后续定量会出错。我曾经在这上面栽过，链搞反了，结果基因表达方向全反。

第三步：定量得到 count 矩阵

拿到排序后的 BAM 文件，就该数数了。我最常用 featureCounts（来自 subread 包），因为它快而且灵活。当然 htseq-count 也很经典，但速度慢一些，尤其对于大样本。

featureCounts -T 8 -p --countReadPairs \
              -t exon -g gene_id \
              -a hg38.gtf \
              -o counts.txt \
              sample1.bam sample2.bam sample3.bam

输出的 counts.txt 长这样：

Geneid      sample1.bam   sample2.bam   sample3.bam
ENSG000001  1234          5678          9012
ENSG000002  234           345           456
ENSG000003  6789          1011          1213

这不就是我们要的 count 矩阵吗？直接拿去做 DESeq2 或 edgeR 的输入。

关键决策：如何处理多比对的 reads？
featureCounts 默认只统计唯一比对的 reads（-M 参数可以改为包括多比对）。我一般保持默认——唯一比对 reads 最可靠。如果实验设计中有大量重复区域或多基因家族，那可以考虑用 -M --fraction 分配分数。

我踩过的坑（避坑指南）

  • 忘记处理 adapter 残留 – 我曾经用 STAR 直接比对未修剪的 fastq，结果 junction 和 mapping rate 惨不忍睹。后来养成习惯：先跑 fastp 再进 STAR。
  • 基因注释版本不一致 – 建 STAR 索引用的 GTF 和 featureCounts 用的 GTF 必须完全一致。我吃过一次闷亏：建索引用 gencode v41，定量时用了 ensembl v108，导致一堆基因 ID 对不上。
  • 盲目去重 – 对于 RNA-seq，PCR 重复是生物学需要还是实验噪音？我强烈建议不要去重，除非你做的是 UMI 标记的文库。否则会严重低估低表达基因。
  • 未检查链特异性 – 先跑一个样本，用 infer_experiment.py（来自 RSeQC）看看实际链特异性。我遇到过文库标记为 "unstranded" 但实际上有 95% 的 reads 落在正义链，featureCounts 默认无链模式导致基因计数翻倍。

工具选择速查表

环节主流工具我的推荐一句话理由
质量控制FASTQC, MultiQCFASTQC + MultiQC一个看细节，一个汇总，绝配
修剪fastp, Trimmomatic, cutadaptfastp快，自动检测接头，输出报告清晰
比对STAR, HISAT2, Salmon (pseudoalignment)STAR速度与准确性平衡，junction 支持好
定量featureCounts, htseq-count, RSEMfeatureCounts快，内存占用少，支持多线程

当然，也有人用 Salmon 或 kallisto 直接做 pseudoalignment 跳过比对，但那就不是本文讨论的传统流程了。如果你做的是非常规样本（比如单细胞），那些工具更合适。

小技巧：写完 pipeline 后，用 samtools flagstat 检查每个样本的比对统计，比如 mapping rate、paired reads 比例。我曾经一个样本 mapping rate 只有 60%，发现是样本污染，及时剔除了，避免了后续差异分析的偏差。

到了这一步，你已经有了一个干净的 count 矩阵。下一件事就是导入 R 开始做差异分析了。不过别急——拿到矩阵后最好再检查一下基因名的格式（ Ensembl ID 还是 gene symbol？），必要时做个 基因ID-基因名 转换。我自己一般在 featureCounts 命令里就用 GTF 里的 gene_id，后续用 biomaRt 或者自带注释表转换。

← 上一章 📖 返回目录 下一章 →