📚 课程目录第4章基因组学质控（一）

4. 基因组学质控（一）：碱基质量分数阈值与过滤策略

做生信的人，每天跟测序数据打交道，最常问的一个问题就是：这批数据质量行不行？

质量行不行，看的第一个指标就是碱基质量分数。说白了，就是每个碱基的靠谱程度。你想想，一个碱基叫错了，后面变异检测、组装、注释全跟着跑偏。我在项目中遇到过好几次，数据看起来量很大，一跑FastQC发现尾巴质量稀烂，最后只好狠狠切掉一截——这就是阈值没设好的后果。

4.1 碱基质量分数：你读到的每个碱基有多可信？

测序仪在读取每个碱基的同时，会给出一个Phred质量分数（Q值）。公式很简单：

Q = -10 × log10(P_error)

P_error 是这个碱基被读错的概率。Q30 对应错误率千分之一，Q20 对应百分之一。我习惯记住几个常用档位：

Q值	错误概率	碱基正确率	通俗理解
Q10	1 / 10	90%	每10个碱基就有1个可能错
Q20	1 / 100	99%	基本可用，但保守区域要小心
Q30	1 / 1000	99.9%	多数分析的标准门槛
Q40	1 / 10000	99.99%	高精度测序的目标

注意一点：Q值不是越高越好。Q40 确实漂亮，但为了追求它而把数据切得太短，反而损失了覆盖度。这个平衡点，就是咱们要设的阈值。

4.2 先看分布再动手：别一上来就定阈值

我个人习惯，拿到FASTQ第一件事不是跑过滤，而是看质量分布图。FastQC 的 Per base sequence quality 图能告诉你每个位置的Q值分布。

举个例子，Illumina 的早期数据经常出现前几个碱基质量偏低，然后中间平稳，尾巴逐渐下滑。这是正常的。但如果你看到整段都在Q20以下晃荡，那就要警惕了——可能是测序跑崩了，或者文库有问题。

怎么判断？我一般看三个点：

整体趋势：是平稳下降还是断崖式下跌？
离散程度：箱线图的范围宽不宽？宽说明数据一致性差。
尾部表现：最后几十个碱基的Q值中位数有没有跌破阈值？

我的小习惯：在跑过滤前，先对每个样品单独看质量分布。不同批次、不同文库的样品，质量特征经常不一样。统一用一套参数过滤，有时候会误伤。

4.3 阈值怎么设？Q20 还是 Q30？

这个问题没有标准答案。但我可以给你几个参考场景：

应用场景	推荐最低阈值	理由
全基因组测序（WGS）	Q20	覆盖度深，少量错误可被覆盖修正
全外显子测序（WES）	Q30	目标区域有限，要求高准确度
RNA-seq	Q20 ~ Q25	表达量分析对单碱基错误容忍度较高
肿瘤体细胞突变检测	Q30	低频突变需要极低假阳性
单细胞测序	Q25	数据稀疏，但每个read都珍贵

我在处理肿瘤配对样本时，曾经把阈值从Q30降到Q25，结果多保留了15%的reads，而且后续突变检测的假阳性率并没有明显上升。你想想看，有时候稍微放一点点，数据利用率就能提升一大截。

4.4 过滤策略：不是一刀切

设好阈值只是第一步。怎么切，讲究挺多。我常用的过滤策略有这三种：

滑动窗口过滤——Trimmomatic 的经典做法。从5'端开始，每次滑4~5个碱基，计算窗口平均质量。低于阈值就从这个位置切掉尾巴。好处是能保留读长中质量好的部分。
直接截断——从3'端往5'端找第一个质量低于阈值的位置，一刀切。简单粗暴，适合质量分布稳定的数据。
全长低质量丢弃——如果整条read的平均质量或最低质量低于阈值，整条扔掉。适用于短读长测序（比如100bp以下）。

我自己的经验：滑动窗口最稳妥。直接截断有时候会切掉一些质量还凑合的碱基，造成不必要的浪费。而全长丢弃太浪费了，除非数据质量实在太差。

核心参数速查

# Trimmomatic 滑动窗口过滤示例参数
SLIDINGWINDOW:4:20
# 含义：窗口大小4bp，窗口平均质量不低于Q20

4.5 动手实践：用Python实现碱基质量过滤

有时候你需要自己写过滤代码，特别是处理非标准数据或者做定制化分析。下面是一段我常用的过滤逻辑：

import gzip
from Bio import SeqIO

def sliding_window_quality_filter(record, window_size=5, min_qual=20):
    """
    滑动窗口过滤：找到第一个质量低于阈值的窗口位置，截断read
    """
    quals = record.letter_annotations["phred_quality"]
    length = len(quals)
    if length < window_size:
        return None  # 太短的read直接丢弃
    
    for i in range(length - window_size + 1):
        window_avg = sum(quals[i:i+window_size]) / window_size
        if window_avg < min_qual:
            # 截断到窗口起始位置
            trimmed = record[:i]
            if len(trimmed) >= 30:  # 过滤后长度不低于30bp
                return trimmed
            else:
                return None
    return record  # 全部通过

input_file = "sample_R1.fastq.gz"
output_file = "sample_R1_filtered.fastq.gz"

with gzip.open(input_file, "rt") as f_in, \
     gzip.open(output_file, "wt") as f_out:
    for record in SeqIO.parse(f_in, "fastq"):
        filtered = sliding_window_quality_filter(record, 5, 20)
        if filtered:
            SeqIO.write(filtered, f_out, "fastq")

代码很简单，但有个细节要注意：letter_annotations["phred_quality"] 这个属性存的是每个碱基的Q值列表。如果你的数据是Sanger格式（大部分Illumina数据都是），Q值范围是0~40，直接拿来算就行。

我曾经踩过的坑：不同版本FastQ编码格式不一样。老旧的Solexa格式（Q值范围-5~40）直接用的话，算出来的平均质量会偏高，导致过滤不够严格。读取时一定用 SeqIO.parse(..., "fastq")，它会自动处理编码格式。如果你自己解析字符串，记得检查第一个读段的ASCII码——Phred+33 还是 Phred+64。

4.6 避坑指南：阈值设定的几个常见错误

阈值设得太高（比如Q35）：数据保留率可能不到50%，很多有效信息被扔掉了。尤其对于低起始量的样本，这是灾难。
阈值设得太低（比如Q10）：错误碱基太多，后续变异检测的假阳性会飙升。我曾经在一个项目中用Q10过滤，结果SNP calling出来一堆假阳性，花了两周才排查出来。
不分位置一刀切：read前10bp的质量通常偏低，但那些碱基其实很重要（比如用于比对锚定）。如果因为前10bp质量低就整条丢弃，太可惜了。建议用滑动窗口，或者对前几个碱基单独设宽松阈值。
忽略双端数据的配对关系：过滤完R1不管R2，导致配对信息丢失。过滤R1和R2时要保持同步——要么都保留，要么都丢弃。

推荐的快速检查方法：过滤完成后，用FastQC再跑一遍过滤后的数据，看看质量分布图有没有明显改善。如果改善了但数据量损失不大，说明阈值设得合适。如果改善不明显但数据量掉了一大截，说明阈值过严了。

好了，关于碱基质量分数阈值和过滤策略，咱们就聊到这儿。下面这张图总结了整个决策流程，你可以照着这个思路走，基本不会跑偏。