表达量化方法：从FASTQ到基因表达矩阵

基因表达矩阵是后续一切分析的起点。说白了，就是把每个基因在每个样本里有多少条 reads 给统计出来。我入行那会儿，有人拿 Excel 手算——嗯，那都是历史了。现在主流就两条路：一条是先比对再计数（STAR + featureCounts），另一条是伪比对直接定量（Salmon）。两条路我都跑过，今天把关键细节拆开讲。

核心思路：不管是哪个工具，最终都要得到一张表——行是基因，列是样本，格子是 count。流程上分四步：质量控制 → 比对/伪比对 → 定量 → 合并表达矩阵。
从FASTQ到表达矩阵：两条主流路线 FASTQ 原始数据 STAR 基因组比对 BAM 对齐文件 featureCounts 基因×样本 count 表 Salmon 伪比对 quant.sf 定量文件 tximport 汇总 基因×样本 count 表 两条路最终都汇入同一张表达矩阵，区别在于速度、内存和能否处理多比对等情况。

方法一：STAR + featureCounts

这是最经典的“先比对后计数”路线。STAR 负责把 reads 比对到参考基因组，速度极快但吃内存。featureCounts 则从 BAM 里把每个基因上的 read 覆盖数捞出来。

我个人习惯这么跑 STAR：

# 构建索引（只需要做一次）
STAR --runMode genomeGenerate \
     --genomeFastaFiles Homo_sapiens.GRCh38.fa \
     --sjdbGTFfile Homo_sapiens.GRCh38.88.gtf \
     --runThreadN 16 \
     --genomeDir star_index

# 比对
STAR --genomeDir star_index \
     --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix sample_ \
     --outSAMtype BAM SortedByCoordinate \
     --outFilterMultimapNmax 10 \
     --runThreadN 16
注意：--outFilterMultimapNmax 默认是 10。如果你的数据来自多拷贝基因家族，建议调低到 2-3，不然很多多比对 reads 会被错误分配。我曾经在一个重复序列高表达的样本上栽过跟头，调参前后差异能到 30%。

得到 BAM 后，用 featureCounts 计数：

featureCounts -T 16 -t exon -g gene_id \
  -a Homo_sapiens.GRCh38.88.gtf \
  -o counts.txt \
  sample_Aligned.sortedByCoord.out.bam

-t exon 指按 exon 为计数单元，-g gene_id 则按基因 ID 汇总。注意 GTF 版本必须和索引一致，否则计数会错乱。

方法二：Salmon + tximport

Salmon 走的是伪比对路线——不把 reads 完整拼到基因组上，而是快速判断每条 read 可能来自哪个转录本。好处是内存占用小得多，而且直接处理 fastq，省掉生成 BAM 的磁盘开销。

跑 Salmon 前需要先准备转录本索引：

salmon index -t transcripts.fasta -i salmon_index --gencode -k 31

--gencode 是给 GENCODE 格式的 fasta 用的，如果是 RefSeq 就别加。我之前遇到过一个坑：不加 --gencode 会导致转录本 ID 解析错误，定量出来的 gene-level 汇总全乱。

定量命令：

salmon quant -i salmon_index -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  -p 16 --validateMappings -o sample_quant

-l A 让 Salmon 自动检测文库类型；--validateMappings 可以过滤掉低质量的比对，建议始终加上。

每个样本会生成一个 quant.sf 文件，里面是转录本水平的 TPM 和有效计数。要汇总到基因水平，需要用 R 的 tximport 包：

library(tximport)
files <- file.path("sample_quant", "quant.sf")
names(files) <- c("sample1","sample2")
txi <- tximport(files, type="salmon", tx2gene=tx2gene, countsFromAbundance="no")
counts <- txi$counts

countsFromAbundance="no" 表示保留原始的 count 值，不转换成 TPM。如果是做 DESeq2 这类基于负二项分布的差异分析，原始 count 才是正确的输入。

STAR+featureCounts vs Salmon 核心差异
维度STAR + featureCountsSalmon + tximport
内存需求高（30G～50G）低（8G～16G）
运行速度较快（但索引大）快（尤其多次运行）
输出文件BAM（大）quant.sf（小）
多比对处理需谨慎调参内置概率模型更鲁棒
能否检测新剪接可以（借助 STAR）不能（依赖已知转录组）
基因水平汇总featureCounts 直接出需 tximport
我的建议：如果实验室服务器内存充裕，且后续要做剪接分析，选 STAR。如果只想做基因表达差异分析，Salmon 更轻量、更省时间。我自己两边都折腾过，现在个人项目几乎只用 Salmon，除非老板非要 BAM 文件。

合并表达矩阵

不管走哪条路，最终每个样本只有一个 count 文件。把它们拼成一张大表就是最后一步。我习惯用 Python/pandas 搞定：

import pandas as pd
import glob

files = glob.glob("*/quant.genes.sf")  # Salmon 产物
cols = []
for f in files:
    df = pd.read_csv(f, sep="\t", index_col=0, usecols=[0,5], header=0)
    sample = f.split("/")[0]
    df.columns = [sample]
    cols.append(df)
matrix = pd.concat(cols, axis=1)
matrix.to_csv("expression_matrix.csv")

如果是 featureCounts，输出自带 Geneid 和样本列，直接合并即可。

避坑指南：所有样本的 GTF/索引必须一致。我见过有人拿不同版本 GTF 分别运行，最后合并时基因名对不上——花了三天时间重跑，血的教训。

常见调参问题

  1. 多比对 reads 要不要留：STAR 默认滤掉 >10 个比对，featureCounts 可以设 -M 允许多比对。如果你研究的是基因家族或重复区，建议保留多比对并用 --fracOverlap 做更精细控制。
  2. 垃圾 reads 的影响：Salmon 的 --validateMappings 能过滤掉低质量映射，但无法替代原始读长修剪。我习惯先把 FASTQ 用 fastp 切一下 adapter 和低质量尾。
  3. 基因水平 vs 转录本水平：表达矩阵通常用基因水平做差异分析。如果后续做异构体分析，Salmon 的 transcript-level 输出可以直接用 tximport(..., type="salmon") 拿到转录本 count。

本章就到这里。记住一条原则：工具是手段，表达矩阵的质量才是命根子。多花十分钟理解参数，省下后面三天 debug 的时间。嗯，这样才划算。

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