📚 课程目录第3章序列质量控制基础

序列质量控制基础：FASTA/Q格式解析、Phred质量评分标准、低质量序列过滤方法

做序列分析的人，有一半的时间都在跟数据质量较劲。

我见过太多人，上来就直接跑BLAST或InterProScan，结果出不来就怀疑软件有问题。其实吧，大多数情况下，是数据本身质量不行。

这一章，我把序列质量控制的核心内容拆成三块讲：格式怎么读、质量值怎么算、低质量的序列怎么处理。

3.1 FASTA与FASTQ格式解析

先说格式。你手里的序列文件，要么是.fasta，要么是.fastq。它们长得很像，但用途完全不同。

3.1.1 FASTA格式

FASTA是最简单的序列存储格式。就两样东西：标题行和序列行。

>seq_id description
ATCGATCGATCGATCGATCG
ATCGATCGATCGATCGATCG

标题行以>开头，后面跟着序列ID，空格后可以写注释。序列行就是碱基或氨基酸字母，通常每行60-80个字符。

实操要点：

序列ID不能有空格（用下划线代替）
碱基用大写字母更规范
避免使用>之外的非法字符

3.1.2 FASTQ格式

FASTQ比FASTA多了质量值。每四条记录代表一条序列：

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

第1行：以@开头，跟FASTA的>类似
第2行：序列本身
第3行：只有一个+（有时也带ID）
第4行：质量值字符串，长度必须跟序列一致

我这里要提醒一句：序列行里不能有换行。FASTA可以分行写，但FASTQ不行。我在项目中遇到过有人把FASTQ里的序列分行写了，结果解析器全乱了。这个坑，踩过就不会忘。

特征	FASTA	FASTQ
包含质量值	否	是
序列换行	允许	不允许
每序列行数	2+	4（固定）
标题行前缀	`>`	`@`
适用场景	参考序列、结果输出	原始测序数据

3.2 Phred质量评分标准

质量值是FASTQ的灵魂。它告诉你，某个碱基测错的概率有多大。

3.2.1 评分公式

Phred质量分数Q跟错误概率P的关系很简单：

Q = -10 * log10(P)

比如：

Q20：错误概率为1%，即每100个碱基里可能错1个
Q30：错误概率为0.1%，每1000个碱基里错1个
Q40：错误概率为0.01%，相当高了

我个人习惯把Q30作为质量合格线。低于这个值，就得留意了。

Q值	错误概率	碱基准确率	常见含义
Q10	10%	90%	极差
Q20	1%	99%	可接受
Q30	0.1%	99.9%	良好
Q40	0.01%	99.99%	极好

3.2.2 ASCII码与质量值的映射

质量值不能直接写数字，那样文件太大了。所以用ASCII字符来编码。

最常见的编码是Sanger格式（Phred+33）：

质量值 = ASCII码 - 33

举个例子：

ASCII字符 ! 的十进制是33，对应质量值0
ASCII字符 5 的十进制是53，对应质量值20
ASCII字符 ? 的十进制是63，对应质量值30

还有一种Solexa/Illumina 1.3+格式（Phred+64），不过现在基本不用了。你看到质量值全是B、C这类大写字母，大概率就是格式没搞对。

注意：千万不要把Phred+33的文件当成Phred+64去解析。质量值会全乱。我曾经踩过这个坑，花了半天才发现是编码格式选错了。

3.2.3 用Python快速检查质量编码

想确认文件是哪种编码？写几行代码跑一下就知道了：

def guess_phred_encoding(qual_string):
    ascii_values = [ord(c) for c in qual_string]
    min_val = min(ascii_values)
    max_val = max(ascii_values)
    
    if min_val >= 64:
        return "Phred+64"
    elif min_val >= 33:
        return "Phred+33"
    else:
        return "unknown"

# 示例
print(guess_phred_encoding("!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65"))
# 输出: Phred+33

小技巧：看质量值的第一个字符。如果第一个字符是!、"、#、$这类符号，就是Phred+33。如果第一个字符是@或更大的字母，可能是Phred+64，但也可能是质量值很高。

3.3 低质量序列过滤方法

拿到原始数据，第一件事就是过滤。说白了，就是把那些质量不行的序列扔掉或修剪掉。

3.3.1 常见过滤策略

我常用的策略有四种，可以组合使用：

基于平均质量分数：计算整条序列的Q值平均值，低于阈值（如Q20）就整条丢弃
基于窗口滑动：从5'端开始，每次取一个窗口（比如10bp），计算平均质量。如果低于阈值，从该位置开始往后的序列全部切除
基于最大N含量：序列里含有太多模糊碱基（N），说明测序质量不好。一般要求N的比例不超过5%-10%
基于序列长度：太短的序列（比如<50bp）没法可靠比对，直接过滤掉

为什么会有这样的文本？举个极端的例子。一条序列开头质量很好，但末尾都是低质量，你直接用平均分数过滤，可能会把整条好序列都扔了。用滑动窗口就很灵活——只切掉尾部不好的部分。

我个人习惯的参数组合：

滑动窗口大小：10bp
窗口平均质量阈值：Q20
最小序列长度：50bp
最大N比例：5%

3.3.2 实战：用Fastp做过滤

Fastp是我最常用的工具之一。它管着快，而且自动识别质量编码格式。

基本用法：

fastp -i raw_R1.fastq -I raw_R2.fastq \
      -o clean_R1.fastq -O clean_R2.fastq \
      --qualified_quality_phred 20 \
      --unqualified_percent_limit 40 \
      --length_required 50 \
      --cut_front \
      --cut_front_window_size 1 \
      --cut_front_mean_quality 20 \
      --cut_tail \
      --cut_tail_window_size 1 \
      --cut_tail_mean_quality 20

参数说明：

--qualified_quality_phred 20：合格线设为Q20
--unqualified_percent_limit 40：允许最多40%的低质量碱基，超过就整条丢弃
--length_required 50：最短序列50bp
--cut_front/--cut_tail：从5'端和3'端各切一个碱基做滑动窗口

3.3.3 手动过滤：自己写脚本

如果你需要更灵活的过滤逻辑，或者做二次开发，可以自己写Python脚本：

from Bio import SeqIO
import gzip

def filter_fastq(input_fastq, output_fastq,
                min_length=50, min_qual=20, max_n=0.05):
    
    count = 0
    passed = 0
    
    with gzip.open(input_fastq, "rt") if input_fastq.endswith(".gz") \
         else open(input_fastq) as in_handle:
        with gzip.open(output_fastq, "wt") if output_fastq.endswith(".gz") \
             else open(output_fastq, "w") as out_handle:
             
            for record in SeqIO.parse(in_handle, "fastq"):
                count += 1
                
                # 长度过滤
                if len(record.seq) < min_length:
                    continue
                
                # N比例过滤
                n_count = record.seq.count("N") + record.seq.count("n")
                if n_count / len(record.seq) > max_n:
                    continue
                
                # 平均质量过滤
                qualities = record.letter_annotations["phred_quality"]
                if sum(qualities) / len(qualities) < min_qual:
                    continue
                
                SeqIO.write(record, out_handle, "fastq")
                passed += 1
    
    print(f"总计: {count}, 通过: {passed}, 过滤率: {((count-passed)/count)*100:.1f}%")

filter_fastq("raw.fastq.gz", "clean.fastq.gz")

注意：如果你用的是Biopython，要注意版本。旧版本SeqIO.parse在处理大文件时可能内存占用偏高。加个iterator选项或者用Bio.SeqIO.QualityIO里的低层级API可以解决。

3.4 本章知识结构图

下面这张图，把序列质量控制的核心逻辑串起来了：

3.5 课后练习：动手试试

光看不练假把式。我从实际项目里抽了个小任务：

找一份公共的FASTQ测试数据（比如NCBI的SRA样本，SRR号里随便挑一个小的）
用BioPython读取前1000条序列，输出它们的平均质量值和序列长度
写一个滑动窗口过滤函数，将尾部质量低于Q20的碱基切掉
统计过滤前后序列长度的变化

写完了，你把过滤前后的统计结果跟原始文件对比一下，就能直观感受到质量控制带来的变化。数据质量好了，后面跑BLAST和InterProScan才能出靠谱的结果。

嗯，这一章的核心就这么多。你先把格式和质量值搞明白，过滤方法熟练用起来，后面几章的内容学起来就轻松多了。