📚 课程目录第3章数据预处理

数据预处理：质量评估、比对与定量

数据预处理这一步，说白了就是把原始测序数据“洗”干净，再把它转换成我们能用的基因表达矩阵。我个人习惯分成三个环节：质量评估、比对、定量。三件事串起来，才算完成从FASTQ到counts矩阵的蜕变。

做这个步骤时，我曾遇到过一个坑——有批样品比对率奇低，查了半天才发现是接头污染没去掉。从那以后，我再也不敢跳过FASTQC了。

1. 质量评估：FASTQC

拿到FASTQ文件，第一件事就是看看数据质量怎么样。FASTQC是最常用的工具，它能给你一份HTML报告，告诉你碱基质量、GC含量、接头污染、重复序列等信息。

运行很简单，一行命令搞定：

fastqc sample1_R1.fastq.gz sample1_R2.fastq.gz -o ./qc_results/

输出的HTML报告里，有几个模块是你必须仔细看的：

Per base sequence quality：每个位置的碱基质量分数。如果尾巴上掉得厉害(Q<20)，就要考虑截断。
Adapter Content：接头残留的百分比。一旦>5%，必须用cutadapt或trimmomatic清理。
Overrepresented sequences：高丰度的未知序列，可能是接头二聚体或rRNA残留。
GC content distribution：正常的正态分布，如果出现双峰可能意味着污染。

我记得有一次，一个项目所有样本的GC含量曲线都呈双峰，后来发现是文库制备时混入了另一批引物。处理过后，曲线才恢复正常。

警告：FASTQC报告不能只看绿色✓。就算所有项目都是绿色，也要检查一下"Per tile sequence quality"——如果某些tile出现偏差，说明测序仪可能有局部问题。

如果你的数据来自公共数据库，比如SRA，通常已经过初步过滤。但我建议还是自己跑一遍FASTQC——你永远不知道别人上传的数据到底干不干净。

2. 比对：STAR vs HISAT2

质量合格后，就该做比对了。把reads贴回参考基因组，我们才能知道每个read到底来自哪个基因。

选哪个比对工具？我个人倾向：

工具	优点	缺点	适用场景
STAR	速度极快，几乎是最快的；支持双端；能处理长read；内存大但可调	内存占用高（一般要求30GB+）；对服务器要求高
HISAT2	内存占用低（~5GB）；支持剪接比对；适合个人笔记本	速度慢于STAR；对低质量reads处理不够好

我通常在服务器上做项目就用STAR。你知道的，RNA-seq数据量动不动就几亿条reads，STAR能省下大把时间。但如果是个人电脑或者内存受限的环境，HISAT2是更好的选择。

STAR的用法分为两步：建立索引和比对。索引需要基因组FASTA和GTF文件：

# 建立索引（一次）
STAR --runThreadN 8 \
     --runMode genomeGenerate \
     --genomeDir ./star_index/ \
     --genomeFastaFiles Homo_sapiens.GRCh38.fa \
     --sjdbGTFfile Homo_sapiens.GRCh38.gtf \
     --sjdbOverhang 100

# 比对
STAR --genomeDir ./star_index/ \
     --readFilesIn sample1_R1.fastq.gz sample1_R2.fastq.gz \
     --readFilesCommand zcat \
     --runThreadN 8 \
     --outFileNamePrefix ./sample1/ \
     --outSAMtype BAM SortedByCoordinate

HISAT2也是类似，不过索引需要额外准备：

# 建立索引（一个命令即可）
hisat2-build Homo_sapiens.GRCh38.fa hisat2_index

# 比对
hisat2 -q -x hisat2_index \
       -1 sample1_R1.fastq.gz \
       -2 sample1_R2.fastq.gz \
       -S sample1.sam

提示：STAR用--outSAMtype BAM SortedByCoordinate会直接输出排序后的BAM，省去后续samtools排序的步骤。我强烈推荐这个参数。

3. 定量：featureCounts

比对之后，我们有了BAM文件。下一步就是定量——统计每个基因被多少条read覆盖。featureCounts是Subread包里的主力工具，速度快，内存占用低，而且支持多种基因模型。

基本命令：

featureCounts -T 8 -p -t exon -g gene_id \
              -a Homo_sapiens.GRCh38.gtf \
              -o counts.txt \
              sample1.bam sample2.bam sample3.bam

参数说明：

-T：线程数，根据CPU核数设置。
-p：双端测序标记，让featureCounts正确处理pair。
-t：特征类型，默认是exon。如果你想看基因水平的定量，保留就好。
-g：分组ID，通常用gene_id。这个要和GTF文件里的属性对应上。

输出的counts.txt是制表符分隔的，直接可以用R或Python读入。列顺序、基因名、每个样本的counts一清二楚。

我曾遇到过一个情况：同一个样本用STAR比对后，再用featureCounts计数，结果发现很多基因counts为0。排查了半天，原来是GTF版本和基因组FASTA版本不一致——基因名对不上号。从那以后，我再也不混用不同版本的注释文件了。

关键要点：GTF文件和参考基因组必须来自同一个版本（比如GRCh38.110）。如果你用了Ensembl的GTF，就别用UCSC的FASTA，反之亦然。

整体流程

从原始FASTQ到表达矩阵，整个数据预处理流程可以浓缩成这张图：

这个流程图中，我强调了一个分支路径：如果FASTQC里发现接头或低质量区域太多，就要先做过滤再比对。不要跳过这步——你想想看，带着一堆垃圾reads去做STAR，不仅浪费计算资源，还可能产生假的跨外显子连接。

小技巧：如果你想自动化整个流程，可以用Snakemake或Nextflow写流水线。我个人偏好Snakemake，Python语法写起来顺手。一个简单的规则就能把FASTQC、比对、定量串起来。

常见问题与避坑

比对率低：检查物种是否匹配，参考基因组是否选对。我有个同事比对小鼠RNA-seq用了人基因组，结果比对率只有5%。
内存不够：STAR可以用--genomeSAindexNbases 12减小索引大小，但会影响速度。
featureCounts输出不全：检查你是否用了-p参数处理双端reads。如果忘了，counts会翻倍。
多线程冲突：不要在同一个目录下并行跑多个STAR，会写冲突。每个样本单独目录。

记得我之前接手一个公开数据集，用HISAT2比对完后featureCounts报了一堆错——原因是GTF文件里基因ID格式不统一（有些带点，有些不带）。最后用sed做了批处理才解决。嗯，这些细节才是真正耗时间的部分。

做完这三步，你应该已经拿到一个整洁的counts矩阵了：行是基因名，列是样本名，每个格子是个整数。这就是下游差异表达分析的“原料”。