2. RNA-seq技术原理与数据特性

RNA-seq，也就是转录组测序，说白了就是给细胞里所有的RNA分子拍一张「全家福」。你想想看，一个细胞在不同状态下，哪些基因在活跃表达，哪些基因被关闭了，RNA-seq都能告诉你。我刚开始做这个的时候，总以为测序仪出来的数据直接就能用，结果被各种批次效应和低质量数据坑过好几次。嗯，今天咱们就把这底层逻辑彻底掰扯清楚。

2.1 RNA-seq的实验流程

整个RNA-seq流程，可以从拿到样本开始，一直到我们手里有个表达矩阵为止。我习惯把它拆成三大块：湿实验、测序、数据分析。这里咱们重点看前两块，因为数据特性很大程度上是由它们决定的。

2.2 文库制备策略

文库制备这一步，决定了你测到的到底是哪些RNA种类。三种最常见的策略，差异还挺大的。

策略	富集目标	适用场景
poly(A) 捕获	mRNA（带polyA尾巴）	常规蛋白编码基因差异表达分析
rRNA 去除	去除rRNA，保留lncRNA、mRNA等	非编码RNA研究、降解样本
总RNA-seq	包含rRNA	特殊需求（非主流）

2.3 测序平台与数据特性

目前绝大多数RNA-seq数据都来自Illumina平台。它会产出一个FASTQ文件，里面记录了每条读段的碱基序列和对应的质量分数。质量分数的格式是这样的：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTTTTCTTGATACTTACCAGTTCTTAATTTCCAAGATCCTCCGCAATATCTT
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII9

第三行是碱基序列，第四行是质量打分。感叹号「I」代表高质量，数字「9」嘛，就比较差了。我个人习惯，拿到数据第一件事就是跑fastqc看看整体质量分布。

参数	典型值	影响
读长	SE 50 / PE 150	PE150更有利于比对到重复区域
测序深度	20-50M 读段/样本	深度不够，低表达基因的差异检测不出来
碱基质量	Q30 > 80%	Q30即每1000个碱基错1个

2.4 从FASTQ到表达矩阵

拿到原始序列之后，我们需要把它「翻译」成基因表达矩阵。这一步常用的工具包括STAR、Salmon、kallisto。我这里用STAR举个简单的例子：

# 构建参考基因组索引（只需做一次）
STAR --runThreadN 8 \
     --runMode genomeGenerate \
     --genomeDir ./star_index \
     --genomeFastaFiles Homo_sapiens.GRCh38.fa \
     --sjdbGTFfile Homo_sapiens.GRCh38.gtf

# 比对样本
STAR --genomeDir ./star_index \
     --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outFileNamePrefix sample_

比对完成后，会得到一个BAM文件。然后用featureCounts或者htseq-count对每个基因上比对到的读段进行计数。最终生成一个像这样的矩阵文件：

GeneID  sample1  sample2  sample3  sample4
ENSG000001  1234     998     1456    2033
ENSG000002  56       23       78      92
ENSG000003  0         1        0       0

每一行是一个基因，每一列是一个样本。数字就是原始count。注意这个数字是离散的、非负的整数。正是因为这个特点，DESeq2这类基于负二项分布的方法才会派上用场。你想想看，要是用普通的t检验去分析这个计数数据，那不翻车才怪。

2.5 为什么原始count不能用做直接比较？

这里有个新手最容易踩的坑。表里的counts为什么不能直接用来比较基因A在样本1和样本2中的表达量？原因有三个：

• 测序深度不同：样本1测了30M reads，样本2只测了15M，那样本2里所有基因的count天然就低一半左右。
• 基因长度不同：长基因比短基因更容易被随机打断测到，count必然偏高。
• RNA组成偏倚：少数高表达基因（比如核糖体蛋白）如果占了大部分reads，会「挤压」其他基因的计数空间。

为了解决前两个问题，常见的做法是计算FPKM或TPM。但DESeq2不走寻常路——它通过内部的一个缩放因子（size factor）来校正测序深度，完全不关心基因长度。为什么？因为差异分析比较的是同一个基因在不同样本间的变化，基因长度是个常数，会约掉。这个逻辑其实很巧妙，我第一次理解到这点时，真的有种豁然开朗的感觉。

2.6 数据的离散特性——为什么需要DESeq2？

RNA-seq的数据天生就有个毛病：方差大于均值。你做个技术重复，测同一个样本两次，每个基因的count都不会完全一样。这种过离散（over-dispersion）特性，用泊松分布是描述不了的。记得我最早用edgeR做分析时，总调参调得头大，后来深入理解后才知道，负二项分布才是这类数据的正确打开方式。