📚 课程目录第3章实验设计

实验设计：数据获取与参考基因组准备

做RNA-seq分析，第一步不是写代码，而是搞清楚你的数据从哪来，以及用什么当“地图”去比对。这一步要是搞错了，后面所有结果都是垃圾。我见过不少项目，就是卡在数据格式不匹配，或者参考基因组版本对不上，白白浪费几周时间。

数据来源与设计原则

先说数据获取。RNA-seq的原始数据，最主流的来源就是高通量测序平台产出的FASTQ文件。你拿到的数据，无论是自己的测序结果，还是从公共数据库下载的，都需要先确认几个关键点：

测序平台：Illumina还是华大？这决定了read长度和错误率分布。
文库类型：单端（Single-end）还是双端（Paired-end）？我个人习惯优先选双端，虽然贵一点，但比对准确率高很多。
测序深度：一般基因表达定量，10-20M reads就够用了。但你要是做可变剪接或者低频突变检测，那得往30-50M走。

公共数据库的话，最常用的就是GEO（Gene Expression Omnibus）和SRA（Sequence Read Archive）。下载时注意看有没有对应的

SRA Run Selector，可以批量获取下载链接。我以前遇到过从GEO下载的数据，跑完流程才发现是单端数据，但参考基因组索引是按双端建的…嗯，那叫一个头疼。

小提示： 下载FASTQ时，直接用fastq-dump或者prefetch工具。记得加--split-files参数，否则双端数据会合并成一个文件，后面分析直接崩。

实验设计上，还有个容易忽略的点：生物学重复。你想想看，如果只有一个样本，你拿什么评估基因表达的可变异性？我建议至少3个重复，找差异表达基因时才更有底气。

参考基因组的准备

接下来，就是准备“地图”——参考基因组。这一步看起来简单，但坑不少。你得先决定用哪个版本的基因组文件。

主流来源有几个：

Ensembl：注释比较全，更新快，我平时用得最多。
NCBI RefSeq：经过审查的，比较可靠，但基因编号可能和Ensembl不同。
GENCODE：联合项目，注重精准注释，适合人类和小鼠。

选版本的时候，我建议统一用同一套来源。比如你比对用Ensembl的基因组，那么后面定量时用的GTF注释文件也一定要从Ensembl下载，千万不能混搭。我记得有一次图省事，用了NCBI的基因组配Ensembl的GTF，结果比对率直接掉了15%，排查了好久才找到原因。

下载命令大概是这样：

# 下载人类参考基因组（GRCh38）
wget ftp://ftp.ensembl.org/pub/release-112/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 解压
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

注意： 不要下载dna_sm或dna_rm版本的基因组文件，那里面是小写的masked序列，会影响后续比对。要下dna.primary_assembly或者dna.toplevel。

GTF注释文件也得一起准备：

# 下载人类GTF注释文件
wget ftp://ftp.ensembl.org/pub/release-112/gtf/homo_sapiens/Homo_sapiens.GRCh38.112.gtf.gz

# 解压
gunzip Homo_sapiens.GRCh38.112.gtf.gz

这里我多说一句，有些人会直接下gtf文件，但我个人习惯先用gzip -d解压后检查一下文件有没有损坏。尤其是网络不稳定的情况下，下载的压缩包可能不完整，直接跑流程会报奇奇怪怪的错误。

索引构建：把地图变成导航系统

有了基因组FASTA和GTF，下一步就是构建索引。为什么需要索引？说白了，就是让比对工具能快速在几亿个碱基里找到每个read的位置。不用索引的话，比对一个样本可能要跑几天。

不同的比对工具，索引构建方式不一样。我这里以STAR和Salmon为例：

工具	索引命令	特点
STAR	`STAR --runMode genomeGenerate`	基于剪接比对，准确度高，适合基因表达和可变剪接分析
Salmon	`salmon index`	基于伪比对，速度极快，适合大规模差异表达分析

STAR的索引构建命令：

STAR --runMode genomeGenerate \
     --genomeDir /path/to/star_index \
     --genomeFastaFiles /path/to/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
     --sjdbGTFfile /path/to/Homo_sapiens.GRCh38.112.gtf \
     --runThreadN 8

Salmon的索引构建命令：

salmon index -t /path/to/transcriptome.fa \
             -i /path/to/salmon_index \
             -k 31

这里要注意，Salmon需要的是转录组序列，不是整个基因组。你可以用gffread从GTF和基因组FASTA中提取转录组：

gffread /path/to/Homo_sapiens.GRCh38.112.gtf \
        -g /path/to/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
        -w /path/to/transcriptome.fa

构建索引时，内存是关键。STAR很吃内存，人基因组索引大概需要30GB左右。我服务器如果内存不够，就加--limitGenomeGenerateRAM 31000000000参数。Salmon相对轻量一些，但也不能在内存紧张的时候跑。

  核心要点：
  数据获取：确认测序平台、文库类型、测序深度，建议至少3个生物学重复。
参考基因组：统一使用Ensembl或GENCODE的来源，避免混搭。
索引构建：STAR适合精细分析，Salmon适合高通量定量。

完整的实践流程

好，纸上谈兵完了，我们直接上手。假设你要分析一个人类双端RNA-seq数据集，下面是我惯用的完整流程：

下载公开数据：先用prefetch获取SRA文件，再用fastq-dump --split-files转成FASTQ。
质量检查：跑FastQC看数据质量，注意看GC含量、碱基质量分布、接头污染。
下载参考基因组：从Ensembl下载GRCh38版本。
构建STAR索引：用上面给的命令，记得调大内存限制。
比对：用STAR把FASTQ比对到索引，输出BAM文件。
定量：可以用featureCounts或者htseq-count从BAM文件统计每个基因的read count。

为什么我不直接用Salmon？看情况。如果你只是做差异表达分析，Salmon更快，结果和STAR差不多。但如果你想研究可变剪接或者新转录本，那就得上STAR。我觉得，两者结合使用才是最佳实践。

避坑点： 我之前在构建索引时，忘了加--sjdbGTFfile参数，结果STAR用的内置剪接位点，导致后续的junction reads比对率低了一大截。后来重新建索引，问题才解决。

最后，别忘了把下载的数据整理好。我习惯建一个data/目录，下面分raw、ref、index三个子目录。Nextflow的流程里，这些路径会作为参数传入，提前规划好能省去后面很多麻烦。

嗯，到这里，数据获取和参考基因组准备这部分就讲完了。你想想看，其实核心就三件事：搞清数据来源、选对参考版本、构建正确索引。把这三点做好，后面的定量分析就顺了。下一步，我会带你把这些环节用Nextflow串起来，让整个流程自动化跑起来。