2. RNA-seq数据获取与预处理：原始测序数据到表达矩阵

RNA-seq的下机数据是一堆FASTQ文件，看起来就像天书。但说白了，我们要的只是一个表格——基因表达矩阵。怎么从原始测序数据一步步走到这个矩阵？这里面有不少坑。

我在刚接触RNA-seq时也踩过不少雷。有一回跑完整个流程才发现，原来某个样本的测序质量根本不过关，白算了三天。所以这一章，我把整个流程拆开来讲，每一步都告诉你：为什么要做、怎么做、容易在哪儿翻车。

2.1 整体流程概览

从FASTQ到表达矩阵，核心就四步：质量控制、比对、定量、汇总。听起来简单，但每一步都有选择。你想想看，光是比对软件就有STAR、HISAT2、Bowtie2好几个，选哪个？默认参数够用吗？

下面这张图把这四步串起来了，建议你先看一遍，心里有个全局感。

2.2 数据获取：从哪找RNA-seq数据？

大部分公开RNA-seq数据都集中在两个库里：GEO（Gene Expression Omnibus）和TCGA（The Cancer Genome Atlas）。我个人习惯优先去TCGA拿肿瘤数据，因为它样本量大、临床信息全。

但要注意，TCGA的数据是经过处理的（比如HTSeq计数），而GEO上很多是原始FASTQ。你需要看清楚描述。

举个例子，用SRA Toolkit下载FASTQ的命令长这样：

# 安装SRA Toolkit后，下载单个SRA条目
prefetch SRR12345678
fasterq-dump --split-files SRR12345678

2.3 质量控制：别让垃圾数据污染结果

拿到FASTQ第一件事，不是直接比对，而是看质量。FastQC是行业标准。跑起来很快：

# 对每个FASTQ跑FastQC
fastqc sample1_R1.fastq.gz sample1_R2.fastq.gz -o qc_reports/

跑完后你会拿到一个HTML报告。重点关注这几项：

如果发现质量不好怎么办？Trimmomatic或cutadapt是常用的清洗工具。我一般这样处理：

# Trimmomatic 双端数据清洗示例
java -jar trimmomatic-0.39.jar PE \
  sample1_R1.fastq.gz sample1_R2.fastq.gz \
  sample1_R1_clean.fastq.gz sample1_R1_unpaired.fastq.gz \
  sample1_R2_clean.fastq.gz sample1_R2_unpaired.fastq.gz \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

multiqc qc_reports/ -o multiqc_report/

我曾经遇到过一个项目，32个样本里有3个的GC含量分布明显异常。后来查实验记录才发现，那3个样本的文库构建批次不同，引入了污染。如果不做QC直接往下走，那3个样本的差异表达基因会全是假阳性。

2.4 序列比对：把读段贴回基因组

清洗完FASTQ，下一步是把每一条read"贴"到参考基因组上。这一步叫比对（alignment）。

目前最主流的工具是STAR和HISAT2。我个人偏爱STAR，因为它速度极快，而且对剪接 Junction 的支持很成熟。但STAR比较吃内存——人类基因组索引大概需要30GB RAM。

STAR的用法分两步：建索引 + 比对。

# 第一步：建立基因组索引（只需做一次）
STAR --runMode genomeGenerate \
     --genomeDir star_index/ \
     --genomeFastaFiles hg38.fa \
     --sjdbGTFfile hg38.gtf \
     --sjdbOverhang 100 \
     --runThreadN 16

# 第二步：比对每个样本
STAR --genomeDir star_index/ \
     --readFilesIn sample1_R1_clean.fastq.gz sample1_R2_clean.fastq.gz \
     --readFilesCommand gunzip -c \
     --outFileNamePrefix sample1_ \
     --outSAMtype BAM SortedByCoordinate \
     --runThreadN 16

STAR会在输出日志里报告这些数字。我一般是拿到BAM后先看一眼Log.final.out，确认比对率没问题再继续。

2.5 基因定量：从BAM到Counts矩阵

有了比对好的BAM文件，下一步就是数一数每个基因上有多少条read。这一步叫定量（quantification）。

最常用的工具是featureCounts（来自subread包）和htseq-count。featureCounts速度更快、内存更省，我基本只用它。

用法很简单：

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

输出文件counts.txt长这样：

Geneid  Sample1.bam  Sample2.bam  Sample3.bam
ENSG000001.1  1234  1098  1456
ENSG000002.2  56    78    92
ENSG000003.3  0     2     1

这就是最原始的表达矩阵——行是基因，列是样本，单元格是整数counts。

2.6 表达矩阵的常见格式与选择

拿到原始counts后，你可能会看到三种表达量单位。它们之间的区别一定要搞清楚：

简单说：如果你要做差异表达，用counts；如果要做热图或聚类，用TPM。我一般会同时保留counts和TPM两份矩阵。

2.7 整个流程的自动化

当你手上有几十个样本，一个个跑STAR和featureCounts会疯掉。我建议用Snakemake或Nextflow来编排流程。下面是一个极简的Snakemake示例：

rule all:
    input:
        "results/counts.txt"

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

rule feature_counts:
    input:
        expand("results/bam/{sample}.bam", sample=SAMPLES)
    output:
        "results/counts.txt"
    shell:
        "featureCounts -T 8 -a hg38.gtf"
        " -o {output} {input}"

把这个写进Snakefile，一行命令snakemake -j 8就能跑完整个流程。中途断了还能断点续跑——这是我最喜欢Snakemake的地方。

到这里，你就拿到了最核心的输入——基因表达矩阵。下一章我们会在这个矩阵上做质量评估和批次校正。但现在，先确保你的矩阵是干净、可靠的。