2. RNA-seq数据分析基础：从原始数据到表达矩阵的标准流程

RNA-seq 分析，说白了就是把测序仪吐出来的原始数据，变成一张你能拿去画火山图、做富集分析的表达矩阵。这套流程在生信领域已经相当成熟，但每一步都有不少坑。我自己当年第一次跑完整的 RNA-seq 流程，就因为跳过某个参数检查，导致后面差异分析全是垃圾——嗯，后来再也不敢偷懒了。

通常我们把 RNA-seq 分析分成“上游”和“下游”。上游是从 FASTQ 到基因/转录本表达量定量，下游是差异表达、功能注释这些。本章只讲上游，也就是从原始数据到表达矩阵的标准流程。

2.1 数据从哪里来：FASTQ 格式

测序仪下机数据一般是 FASTQ 文件。每条记录有 4 行：标识符、序列碱基、分隔符（+）、质量值。我见过不少新人直接去看序列文件，但其实第一件事应该是检查文件完整性——曾经有一次我发现某个样本的 FASTQ 文件少了两百条记录，原因是传输过程中断了。所以拿到数据后，先跑 md5sum 校验一下。

我习惯用 zcat 直接看压缩包里的头几行：

zcat sample_1.fastq.gz | head -n 8

这样能快速确认读长、碱基组成、质量编码版本（一般是 Phred+33）。

2.2 质量评估 FastQC

质控是绕不开的第一道关。FastQC 会生成一个 HTML 报告，告诉你每个碱基位置的质量分布、GC含量、接头污染等。我个人习惯把所有样本的 FastQC 报告用 MultiQC 汇总成一个页面，这样一眼就能看出异常样本。

比如，我曾在一次项目中看到某个样本的 Per base sequence quality 在末尾急剧下降，结果发现是测序仪跑着跑着试剂用完了……这个样本后期修剪也救不回来，只能重测。所以别跳过这一步。

2.3 接头修剪与低质量过滤

做完 FastQC 发现接头残留或低质量尾巴，就得修剪。我用得比较多的是 cutadapt 和 Trimmomatic。这里要注意一点：双端测序的修剪参数要成对处理。我建议先用 cutadapt 去接头，再用滑动窗口截除低质量碱基。

cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC \
         -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz \
         raw_R1.fastq.gz raw_R2.fastq.gz

实际项目中，我曾经因为参数 -m 36（最小长度）设得太高，把很多有效短片段都扔了，导致比对率低。后来我改成 -m 20，数据利用率上升了不少。具体阈值取决于你的读长和实验设计。

2.4 比对：STAR 与 HISAT2 的选择

比对器里 STAR 应该是 RNA-seq 的标配了。它速度快、敏感度高，还直接输出拼接后的 junction 信息。我第一次用 STAR 时被它的索引命令吓到了——需要几十 G 内存，但其实只需要跑一次：

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

我的个人经验：如果你想做可变剪接分析，一定要提供 GTF 文件生成索引，这样 STAR 才知道已知的剪切位点。另外，双端数据跑比对时加上 --outSAMtype BAM SortedByCoordinate，直接输出排序后的 BAM，省后续步骤。

STAR --genomeDir /ref/star_index \
     --readFilesIn trimmed_R1.fastq.gz trimmed_R2.fastq.gz \
     --readFilesCommand zcat \
     --outFileNamePrefix sample_ \
     --outSAMtype BAM SortedByCoordinate \
     --runThreadN 8

2.5 基因表达定量：featureCounts vs htseq-count

比对后的 BAM 文件要转换成基因表达量。主流的工具有 featureCounts（subread 包）和 htseq-count。我偏爱 featureCounts，因为它快，而且能同时统计基因、外显子、转录本级别。命令长这样：

featureCounts -T 8 -p -t exon -g gene_id \
              -a /ref/GRCh38.gtf \
              -o counts.txt \
              sample_Aligned.sortedByCoord.out.bam

参数 -p 是在双端模式下计数片段而不是读段；-t exon 指定只统计外显子区域。我踩过一个坑：如果你用 htseq-count，它默认是 union 模式，但如果有重叠基因（比如酵母基因组），得改成 intersection-nonempty，否则会丢掉大量计数。

2.6 构建表达矩阵

featureCounts 输出是一个文本文件，包含每行一个基因，每列一个样本。咱可以用 R 或 Python 去读，但我还是建议在管道里直接把它整理成标准的表达矩阵——行名是基因 ID，列名是样本名，数值是 raw count。

举个例子，在 bash 里合并多个样本的 counts：

# 假设每个样本输出单独文件 sample1.count.txt
# 可以用类似下面的脚本做合并（这是示意）
echo "gene_id" > merged_counts.txt
for sample in sample*; do
    cut -f7 $sample/${sample}.count.txt | tail -n +2 | paste -d'\t' merged_counts.txt - > tmp
    mv tmp merged_counts.txt
done

当然在实际 Nextflow 管道里，我们直接用流程自带的 collect 和 join 组合操作，不会这么手动处理。但理解这个逻辑很重要——表达矩阵是一切下游分析的起点。

2.7 流程中的常见坑与避坑指南

• 命名不一致：样本名里带空格、特殊字符，会导致很多工具报错。我习惯全部用下划线+数字，如 CRC_001_R1.fastq.gz。
• 比对参考版本：同一个物种可能有多个参考基因组版本（比如 hg19 vs GRCh38），一定要统一。我曾经混用导致 30% 的 reads 没比对上去。
• 计数时忘记双端参数：featureCounts 不加 -p 会把一对 reads 算成两个计数，直接翻倍。
• 染色体命名：有些数据库用 “chr1”，有些用 “1”。不匹配时，比对或计数会直接跳过。我有个同事就因为这个问题排查了整整一天。