📚 课程目录第4章 Conda软件包管理

4. RNA-seq分析流程概述：从原始测序数据到差异表达基因的典型步骤

说实话，很多刚入行的朋友问我，RNA-seq分析到底要做哪些事。

我通常回答：一场数据清洗与信号提取的接力赛。

原始数据里充满了噪音、接头、低质量碱基。只有一步步把干扰信息去掉，才能找到真正的生物学信号。

我个人习惯把整个流程拆成5个核心阶段。每个阶段都有它的坑，也有它的巧劲。下面这张图能帮你快速建立全局观。

4.1 质量控制——数据的第一次体检

拿到FASTQ文件，我第一件事就是跑FastQC。

这不是例行公事。你想想看，一张模糊的照片，再怎么分析也看不出细节。数据质量不过关，后面全是白费功夫。

质量指标	正常范围	报警信号
每个碱基的测序质量	Q30 > 80%	Q30 < 70%，需切除低质量端
GC含量分布	均匀或接近理论分布	异常偏移提示污染或偏好性扩增
接头污染	接头含量 < 5%	> 10% 需要去除接头
重复序列比例	与实验类型匹配	RNA-seq重复率不宜过高，警惕PCR重复

如果发现质量不行，我会用Trimmomatic或fastp做一次清洗。命令大概是这样：

trimmomatic PE sample_R1.fastq sample_R2.fastq \
  sample_R1_clean.fastq sample_R1_unpaired.fastq \
  sample_R2_clean.fastq sample_R2_unpaired.fastq \
  ILLUMINACLIP:adapters.fa:2:30:10 \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

我的经验：不要只依赖FastQC的默认阈值。不同测序平台、不同文库类型，正常值范围有差异。多看几个样本，建立自己的质量基线。

清洗完记得再跑一遍FastQC。我习惯把多个样本的结果用MultiQC汇总成一个报告，一目了然。

multiqc . --title "RNA-seq QC Report"

4.2 序列比对——把reads送回它们的家

清洗干净的数据，接下来要跟参考基因组做比对。这一步说白了，就是把一条条短的reads，映射回它原本该在的基因组位置。

我个人偏好STAR。为什么？因为它快，而且对剪接点处理得很好。RNA-seq的特点是存在内含子，STAR的剪接感知比对能力是刚需。

STAR --runThreadN 8 \
     --genomeDir /path/to/STAR_index \
     --readFilesIn sample_R1_clean.fastq sample_R2_clean.fastq \
     --outFileNamePrefix sample_ \
     --outSAMtype BAM SortedByCoordinate \
     --outBAMsortingThreadN 4

注意：STAR索引很吃内存。人类基因组索引大约需要30GB内存。我曾经在16GB的笔记本上跑，直接崩了。建议提前规划计算资源。

比对完成后，会得到一个BAM文件。这里有一个我踩过的坑——比对率过低。比如RNA-seq的比对率通常应在85%以上。如果低于70%，先别急着往下走。看看是不是参考基因组版本不对，或者样本有污染。

另一个选择是HISAT2。它比STAR更省内存，但速度慢一些。看你的场景选。

4.3 表达定量——给每个基因算账

有了BAM文件，接下来就是统计每个基因有多少条reads覆盖。这就是表达定量。

常用的方法是featureCounts。它简单、快速、稳定。我90%的项目都用它。

featureCounts -T 4 -t exon -g gene_id \
  -a /path/to/annotation.gtf \
  -o counts.txt \
  sample1.bam sample2.bam sample3.bam

输出是一个矩阵文件，行是基因，列是样本，每个值就是count数。这就是后续差异分析的原材料。

  为什么不推荐直接看FPKM/RPKM做差异分析？ 因为在统计建模中，raw counts更符合负二项分布的假设。标准化后的值反而会丢失一些统计信息。记住，DESeq2需要raw counts作为输入。

如果想走无比对定量路线，可以用Salmon或kallisto。它们直接对fastq做准比对，速度很快，适合大规模筛选。

4.4 差异表达分析——找出关键基因

有了表达矩阵，正菜来了：比较两组样本，找出表达水平有显著差异的基因。

常用的工具有DESeq2、edgeR、limma。我个人偏爱DESeq2。它的稳健性很好，特别适合小样本量实验。

在R里跑DESeq2的标准步骤：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = meta,
                              design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treatment", "control"))
res <- res[order(res$padj), ]

结果里主要看两列：log2FoldChange（差异倍数）和padj（校正后的p值）。

通常筛选阈值：|log2FoldChange| > 1（即2倍差异），padj < 0.05。

避坑指南：我曾经发现一批差异基因在重复验证时全不成立。后来排查发现，样本分组有混淆变量。所以在跑差异分析前，一定要先做PCA或聚类，看样本是否按实验分组自然分开。如果对照组和治疗组混在一起，那结果不可信。

4.5 功能富集分析——从基因列表到生物学故事

拿到差异基因列表，我们不能只看清单。一个基因一个基因解读效率太低，也看不出整体趋势。功能富集分析就是帮我们把基因列表"翻译"成生物学通路和功能模块。

常用的工具是clusterProfiler，它支持GO富集、KEGG通路富集等。

library(clusterProfiler)
ego <- enrichGO(gene = deg_genes, 
                OrgDb = org.Hs.eg.db,
                keyType = "ENSEMBL",
                ont = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05)
dotplot(ego, showCategory = 20)

除了GO和KEGG，GSEA（基因集富集分析）也值得一试。它的好处是不需要设定阈值，所有基因按表达变化排序后，看整个基因集是否协同变化。

常见问题：有些朋友直接用DAVID在线工具富集，结果通路全是"癌症通路"或"代谢通路"，什么解释都做不了。这就是输入基因太多太杂导致的。建议先用|log2FC|和padj严格筛选，让富集更聚焦。

OK，到这里，从原始测序数据到差异表达基因的典型步骤，已经完整走了一遍。质量控制、比对、定量、差异分析、富集解释，一共5个阶段。每个阶段都有它的工具和陷阱。

我建议新手不要着急跑流程。先花一天理解每个步骤在做什么，再用小数据集跑一遍。等有了手感，再上真实数据。嗯，这样不容易翻车。