表达矩阵构建：从fastq到count matrix的完整流程

从下机数据到基因表达矩阵，这条路我走了不下百遍。说白了就是把几亿条短读序列，一步步浓缩成一个数字矩阵——每个基因在哪个样本里表达了多少。流程不复杂，但坑不少。今天带你捋一遍核心环节。

整体流程概览

我把完整流程拆成五个阶段：质量控制 → 比对 → 定量 → 格式转换 → 整合。每个阶段都有成熟工具，选对组合能省一半时间。

fastq（原始数据） 
   ↓ 质量控制（FastQC + Trimmomatic / fastp）
clean_fastq（清洁数据） 
   ↓ 比对（STAR / Hisat2 / Salmon）
bam / sam / quant.sf（比对结果）
   ↓ 定量（featureCounts / HTSeq / tximport）
count matrix（基因×样本矩阵）
   ↓ 整合 + 初步过滤
最终表达矩阵（用于下游分析）

第一步：质量控制与清洗

拿到fastq第一件事不是直接比对，而是先看看数据底子怎么样。我个人习惯用FastQC做前检查，跑个批量看所有样本的碱基质量、GC含量、接头污染。

# 批量跑 FastQC
fastqc -o qc_reports/ *.fastq.gz

看完报告，再用Trimmomatic或fastp做修剪。记得以前有个项目，数据里混了一堆Illumina接头，跑完STAR比对率只有50%。后来加了接头切除，比对率直接到90%。嗯，这一步不能省。

# Trimmomatic 示例（PE reads）
trimmomatic PE sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_trimmed_R1.fastq.gz sample_unpaired_R1.fastq.gz \
  sample_trimmed_R2.fastq.gz sample_unpaired_R2.fastq.gz \
  ILLUMINACLIP:adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

第二步：比对与伪比对

这是整个流程的算力大头。选择哪种比对器，取决于你的数据特点和下游需求。下表是我的常用组合：

我个人最常用STAR + Salmon组合：先用STAR生成BAM，再用Salmon做定量。STAR的索引建立比较耗时，但比对速度很快。

# 建立STAR索引（只需一次）
STAR --runThreadN 8 \
  --runMode genomeGenerate \
  --genomeDir star_index \
  --genomeFastaFiles genome.fa \
  --sjdbGTFfile annotation.gtf

# 比对
STAR --genomeDir star_index \
  --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
  --readFilesCommand zcat \
  --runThreadN 8 \
  --outFileNamePrefix sample/ \
  --outSAMtype BAM SortedByCoordinate

第三步：定量——从BAM到counts

获得BAM后，你需要把他转化成基因水平的read counts。主流工具是featureCounts和HTSeq-count。我倾向于featureCounts，因为快，而且能处理多线程。

# featureCounts 批量处理
featureCounts -T 8 -t exon -g gene_id \
  -a annotation.gtf \
  -o counts.txt \
  sample1.bam sample2.bam sample3.bam

输出的 counts.txt 长这样：

Geneid  sample1.bam  sample2.bam  sample3.bam
ENSG000001.1    120    98    145
ENSG000002.2    5    2    11
...

你想想看，一个几GB的BAM文件，最后变成了一个只有几MB的计数表。这就是数据压缩的魅力。

第四步：整合成count matrix

如果用了Salmon/kallisto这类伪比对工具，输出的不是直接counts，而是转录本丰度（TPM、估计counts）。这时候需要用tximport把结果汇总到基因水平。

# R代码：tximport 整合
library(tximport)
files <- file.path("quant", samples, "quant.sf")
names(files) <- samples
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "no")
count_matrix <- txi$counts

第五步：初筛与过滤

拿到count matrix别急着跑下游。先做两件事：低表达基因过滤和样本均一化。低表达基因是指那些在大部分样本里count数为0或很小的基因。我通常保留至少在50%样本里count≥10的基因。

# R代码：过滤低表达基因
keep <- rowSums(count_matrix >= 10) >= ncol(count_matrix) * 0.5
count_matrix_filtered <- count_matrix[keep, ]

均一化常用DESeq2的median of ratios或edgeR的TMM。这一步是为了消除文库大小和RNA组成差异带来的偏差。

本章知识结构图

下面的SVG概括了从fastq到count matrix的完整路径，以及每阶段的关键工具：

常见对比：比对 vs 伪比对

新手经常纠结该用STAR还是Salmon。我直接列个对比表：

我个人观点：如果你只关心基因表达差异，且实验设计里有大量重复，Salmon足够好。但如果你还要做剪切分析或寻找融合基因，那STAR是必须的。

几个关键避坑点

• 参考资料一致性：从构建索引到定量，必须使用同一个GTF/GFF文件。换版本会导致基因ID对不上。我曾在不同版本之间来回切换，结果下游分析跑出一堆NA——浪费了三天时间。
• 链特异性：建库时如果是链特异性文库（例如dUTP法），定量时一定要指定链方向参数。featureCounts加 -s 1 或 -s 2，否则计数会翻倍或丢失。我曾经漏了这个参数，最后发现基因表达量整体偏高两倍——调成链特异性后恢复正常。
• 多比对read处理：featureCounts默认忽略多比对的read（mapQ次要比对）。如果想保留多比对，可以加 -M 选项，但要小心：多比对read往往来自重复区域或同源基因，会引入噪声。
• 并行与资源：STAR的线程数建议不超过物理核心数。我在32核服务器上跑STAR，设了32线程，结果内存耗尽直接崩了——后来减到16线程，稳稳跑完。

自动化管线的搭建思路

手动一个个跑命令太累。我通常用Snakemake或Nextflow构建流程。下面是个Snakemake的最小示例骨架：

# Snakefile 片段
rule all:
    input: "results/counts.txt"

rule fastqc:
    input: "data/{sample}_R1.fastq.gz"
    output: "qc/{sample}_fastqc.zip"
    shell: "fastqc {input} -o qc/"

rule star_align:
    input: "data/{sample}_R1.fastq.gz", "data/{sample}_R2.fastq.gz"
    output: "aligned/{sample}.bam"
    params: index = "star_index"
    shell: "STAR --genomeDir {params.index} --readFilesIn {input} --outFileNamePrefix aligned/{sample}_ --outSAMtype BAM SortedByCoordinate"

rule quantify:
    input: "aligned/{sample}.bam"
    output: "counts/{sample}.txt"
    shell: "featureCounts -a annotation.gtf -o {output} {input}"

有了这个骨架，再复杂的流程都能自动化跑起来。我强烈建议你用这类工具管理流程——不仅可重复，还能自动跳过已完成的步骤。