4. 基本统计指标：序列数量、长度、GC含量

做测序质控，你拿到手的第一个东西就是 fastqc_report.html。打开一看，绿色对勾、黄色感叹号、红色叉……嗯，先别急着下结论。我习惯从最基础的三个指标开始扫一遍——序列数量、序列长度、GC 含量。这三个数字要是出了毛病，后面的模块都不用看，直接打回重跑都可能。

说白了，这几个指标就像体检报告里的身高、体重、体温。不正常往往不是孤立问题，而是整批数据从建库、测序到数据产出哪个环节出了岔子。我带团队时，新人经常被那十几个模块花眼，我总说：先把前三项盘清楚，再碰后面的。

4.1 序列数量：多少条才算够？

FastQC 会直接告诉你总序列数（Total Sequences）。这个数字跟你测序下机的原始 Reads 数一致。但麻烦在于——它不告诉你哪种「不够」。

我记得有一次，客户给的 RNA-Seq 数据每个样本只有 800 万条 read。Q30 再高，基因检测覆盖度也是笑话。从那以后我养成了习惯：先核对测序承诺的 reads 数，再对比 FastQC 报出来的数字。差太多？跑前就出问题了。

阈值参考：
对于单端 50bp 的 ChIP-Seq，我一般要求 ≥ 2000 万条；双端 150bp WGS ≥ 5 亿条。当然，不同的数据类型差别很大。关键是——FastQC 不设硬阈值，黄色警告只根据 read 数奇偶性判断（偶数时 normal），所以你需要自己定一个业务阈值。

实操中我会这样校验：

  • 查看原始测序报告中的 total reads / total bases
  • 用 seqkit stats 或 fastq-stats 交叉验证
  • 对比 FastQC 的 Total Sequences 是否一致（差 > 5% 基本就是处理过程中丢数据了）
# 命令行交叉验证 FastQ 统计
seqkit stats sample_R1.fastq.gz sample_R2.fastq.gz

输出示例：

file              format  type   num_seqs   sum_len  min_len  avg_len  max_len
sample_R1.fastq.gz  FASTQ   DNA   12,345,678  1.9 G    150      150      150
sample_R2.fastq.gz  FASTQ   DNA   12,345,678  1.9 G    150      150      150

两个文件 num_seqs 完全一致而且和 FastQC 报告的 Total Sequences（12345678）相同，那基本放心。

注意：如果 FastQC 跑的是合并后文件（例如 cat 过的），Total Sequences 会是所有样本总数。审计时一定要拿到单样本的原始 FastQ，避免被合并数据误导。

4.2 序列长度：建库碎片还是酶切事故？

FastQC 的 Sequence Length 模块会显示所有 read 的长度分布。正常情况下，如果你的测序是固定读长模式（比如 150bp PE），那么所有 read 的长度应该完全相同。一点波动都没有。

一旦出现长度不一致，警报就响了。我自己踩过最大的坑是做 microRNA 测序时，小 RNA 文库 insert 本身就很短，加上 adapter 也没用完，结果 FastQC 报「长度红色警告」。我当时以为测序故障，后来才明白——miniRNA 文库就那样，阈值不能套用 WGS。

所以对于长度指标，我建议分两种情况看：

  • 常规 固定读长 测序（WGS/WES/RNA-Seq）： 所有 read 长度必须 100% 等于标称长度。出现 > 1% 的短 read → 大概率是接头污染、试剂故障或机器光学识别错误。
  • 可变长度文库（如 miRNA、cfDNA 片段组 或 靶向捕获后切胶回收）： FastQC 红色警告是正常的。你要关注的是长度分布峰形是否符合预期。比如 miRNA 应该在 21-23nt 出现主峰，如果变成宽峰，那才需要排查。
我个人的检查清单（审计时）：
1. 非可变长度文库 → 长度必须 100% 等于设置值，否则标记为 FAIL
2. 可变长度文库 → 允许 WARN，但需查看 per base sequence content 是否异常（短的片段更容易出现前几个碱基的偏好性）
3. 双端文件两个 R1/R2 长度必须一致（除非特意做不对称测序）

4.3 GC 含量：比想象中敏感得多

GC 含量（%GC）是 FastQC 中非常敏感的一个模块。它直接画一条曲线，对比理论正态分布（根据你的基因组 GC 平均含量计算）。偏离太远 → 红色警告。但你得知道·GC 含量·并不是「越接近 50% 越好」。

举个例子，人类基因组平均 GC 含量约 41%（取决于参考版本），但 AT 丰富区（比如 chr13、chr21）能低到 35%。如果测的是这种区域富集的数据，GC 曲线左偏完全正常。

我曾经审计过一个微生物扩增子项目，样本是 GC 含量差异很大的几种菌混在一起，结果 FastQC 的 GC 曲线变成了一个双峰。客户吓坏了，觉得数据污染了。我让他们跑一遍 kraken2 确认物种组成，结果就是不同菌的天然 GC 差异导致的。所以，红色警告不一定是坏数据，但要搭配实验背景理解。

4.3.1 如何设定 GC 阈值？

  • WGS（全基因组）： 期望 GC 均值 ± 5% 以内。人类样本正常范围 36% – 46%。超出直接报警。
  • RNA-Seq： 因为 rRNA 去除和 polyA 捕获，GC 曲线往往右倾（高 GC 偏高）。我一般放宽容忍度到 ± 7%。
  • 靶向测序（Panel）： 必须根据目标区域计算理论 GC 分布，不能跟全基因组参考曲线硬比。

实际在 FastQC 里，你可以复现 GC 含量的计算：

def gc_fraction(seq):
    """计算单条序列的 GC 比例"""
    if not seq:
        return 0.0
    gc_count = seq.count('G') + seq.count('C')
    return gc_count / len(seq)

然后对全部 read 做统计，就能画出那条 %GC per read 分布的柱状图。FastQC 内部正是这么做的。

4.4 三个指标联合审计流程

讲到这里，我把它们串在一张图里给你看。这是我自己总结的判断流程，从序列数量开始，走不下去就亮红灯：

基本统计指标联合审计流程 读取 FastQC 报告 Total Sequences 达标？ FAIL：回测序中心确认 OK → 检查长度 序列长度一致吗？ FAIL：检查接头污染 GC 含量曲线偏离？ FAIL：需确认实验背景 PASS → 继续下游 说明：GC 红色警告 不一定错误，需结合文库类型
一个重要提醒：我见过不少团队在自动化 pipeline 里直接把 FastQC 的绿色当作 pass、红色当作 fail。这很危险。比如 GC 含量的红色警告，在 AT 丰富基因组、单细胞低深度数据、富GC文库中都非常正常。一定要根据 project 定制阈值，而不是看颜色。

4.5 实操：用 Python 批量提取三项指标

FastQC 的纯文本输出是 fastqc_data.txt，藏在每个样本的 .zip 包内部。我习惯写一个小脚本批量抓取：

import zipfile
import re

def extract_fastqc_basic(zip_path):
    """
    从 FastQC zip 中提取三个基本统计量
    """
    info = {}
    with zipfile.ZipFile(zip_path, 'r') as z:
        # 假设 zip 内结构: sample_fastqc/fastqc_data.txt
        for name in z.namelist():
            if name.endswith('fastqc_data.txt'):
                text = z.read(name).decode('utf-8')
                # 匹配 Total Sequences
                seq_match = re.search(r'Total Sequences\s+(\d+)', text)
                # 匹配 Sequence length
                len_match = re.search(r'Sequence length\s+(\d+)', text)
                # 匹配 %GC
                gc_match = re.search(r'%GC\s+([\d.]+)', text)
                if seq_match:
                    info['total_sequences'] = int(seq_match.group(1))
                if len_match:
                    info['sequence_length'] = int(len_match.group(1))
                if gc_match:
                    info['gc_percent'] = float(gc_match.group(1))
                break
    return info

# 用法示例
sample_data = extract_fastqc_basic('SRR123456_fastqc.zip')
print(sample_data)

输出：

{'total_sequences': 12345678, 'sequence_length': 150, 'gc_percent': 41.2}
审计时可以这样用：把整个项目数百个样本的 zip 跑一遍，汇总成 CSV，然后对着表格看异常点——哪些样本 Total Sequences 减半了？哪些 GC 低于 36%？一眼就能锁定。

4.6 开篇遇到的情况，我这样处理

有一次，一个 RNA-Seq 项目的 48 个样本里，有 6 个样本的序列数量只有平均值的 30%。按理说直接打回。但客户坚持说取样批次不一样。我让技术把 6 个样本的插入片段长度和 GC 曲线调出来看——结果发现这些样本的 GC 曲线都出现了诡异的“双峰”，而且 read 长度从 150bp 变成 50~150bp 不等。最后查明是 adapter dimer 污染 + 测序机器某通道故障。如果不是把三个基本指标联合起来看，光看数量就被误导了。

所以，记住我的话：序列数量、长度、GC 含量，这三个点必须同时看。任意两个出问题，基本就是实验或测序有问题。只有一个出问题？先查是不是算法或阈值选错了。

← 上一章 📖 返回目录 下一章 →