第3节 数据质控
一、正式流程的搭建，整个完整的流程分为以下6部分：
① 原始测序数据 fastq 的质控 QC；
② read比对，排序和去除重复序列；
③ Indel区域重（“重新”的“重”）比对；
④ 碱基质量值重校正BQSR；
⑤ 变异检测；
⑥ 变异结果质控和过滤VQSR；
注意：
① 当前以 Illumina 为代表的新一代测序 (NGS) 平台，主要采用边合成边测序（SBS, Sequencing by Synthesis）的技术。该方法通过化学反应完成碱基的合成，使 DNA 链从 5' 端逐步延伸至 3' 端。然而，随着合成链长度的增加，DNA 聚合酶的效率会逐渐降低，其特异性也会变差，这就导致了碱基合成的错误率在链延伸后期逐步升高。这种现象是目前 NGS 测序读长普遍较短的主要原因之一。此外，在测序开始阶段，由于化学反应尚未稳定，可能出现质量值波动的情况。但这种波动通常局限于高质量值区域，影响相对较小。
② 而测序数据的质量直接影响下游分析的准确性，而不同测序平台的错误率分布和偏差特性存在显著差异。因此，在分析测序数据前，必须清晰了解以下几点：
（1）测序平台与错误率分布
- 确定数据来源的测序平台及其错误率分布特征。
- 判断测序数据是否存在偏向性或局限性，例如是否显著受 GC 含量的影响。
GC含量是指DNA或RNA序列中鸟嘌呤（G）和胞嘧啶（C）这两种核苷酸的比例。GC含量通常用以下公式计算：
A代表腺嘌呤（Adenine），T代表胸腺嘧啶（Thymine）
- 了解平台特性的同时，结合官方文档与自身数据分析进行验证。
（2） 错误率对下游分析的潜在影响
- 变异检测：测序深度是否充足？变异位点是否过度集中在 read 的末端？是否存在正反链偏向性？
- 基因组组装：组装过程中对错误碱基的敏感性较高，需要更严格的剪切策略，以避免计算资源和时间的额外消耗。
（3）实际策略与建议
- 全面了解数据：分析数据质量图谱，评估错误率分布和 GC 偏向性，必要时进行质量剪切和纠错。
- 视具体需求调整策略：如变异分析中评估深度与偏向性，或基因组组装时优化数据预处理流程。
二、如何全面认识原始测序数据（FASTQ 数据）
原始测序数据的分析是了解实验质量和评估下游分析潜在问题的第一步。一般来说，我们可以从以下几个方面入手：
（1）Read 各个位置的碱基质量值分布
通过评估每个位置上的质量值，可以直观反映测序错误率的趋势。
（2）碱基的总体质量值分布
分析所有 read 的质量值概览，了解数据的整体可靠性。
（3）Read 各位置的碱基分布比例
目的是分析各碱基 (A/T/G/C) 的分离程度是否正常，例如是否存在某些碱基过多或过少的现象。
（4）GC 含量分布
检查测序数据是否存在 GC 偏向性或异常波动，这可能会影响下游分析。
（5）Read 各位置的 N 含量
N 表示不确定碱基，其分布能反映测序错误的随机性或局部错误的集中性。
（6）Read 是否包含测序接头序列
接头序列的存在会干扰比对和组装，需检测并剪切。
（7）Read 重复率
评估因 PCR 扩增引入的重复情况，重复率过高可能表明样本复杂度不足或实验问题。
通过全面分析以上几个方面，可以对原始数据的基本质量和潜在问题有一个清晰的认识。
三、Read 各位置碱基质量值分布的分析
在前文中已经展示了如何使用 Python 代码计算单条 read 的质量值和测序错误率。然而，当面对成千上万条 read 时，逐条分析显然不现实，此时更直观的表达方式是数据可视化。
FastQC是目前广泛使用的工具之一，它能够高效生成测序数据的质量控制 (QC) 报告。其报告涵盖了上述多个方面，并提示数据可能存在的问题。例如，FastQC 会通过综合分析 FastQ 文件中所有 read 的数据，为 read 各位置的碱基质量值分布绘制图表，如下图所示：
测序质量分布图
注：在这个分布图中，横轴表示 read 的位置，纵轴表示质量值。对于一个测序质量较好的数据集，大多数碱基的质量值应在高分区域（通常为绿色）。FastQC 的优势在于快速生成报告并帮助识别潜在问题，但其缺点是图表美观度较差啦~
除了上面read各位置的碱基质量值分布之外，FastQC还会为我们计算其他几个非常有价值的统计结果，包括：
1）碱基总体质量值分布，只要大部分都高于20，那么就比较正常。
一般来说，对于二代测序，最好是达到Q20的碱基要在95%以上（最差不低于90%），Q30要求大于85%（最差也不要低于80%）。
2）read各个位置上碱基比例分布
碱基分离是指在测序过程中，A与T、C与G的比例应大致相当。如果测序过程是随机的（即理想状态），那么在每个位置上，A和T的比例、C和G的比例应该相近。理想情况下，两者之间的偏差应控制在1%以内。如果偏差过大，除非有合理解释（如特定的捕获测序），否则需要关注测序过程是否存在偏差。
3）GC含量分布图
GC含量是指G和C两种碱基在总碱基中的比例。二代测序平台通常存在一定的测序偏向性，通过分析GC含量可以判断测序过程的随机性。人类基因组的GC含量一般约为40%。如果GC含量的图谱明显偏离这个值，说明测序过程中存在较高的序列偏向性，可能导致某些特定区域被反复测序的几率高于平均水平。这种偏差不仅会影响覆盖度，还可能对下游的变异检测和CNV分析产生负面影响。
4）N含量分布图
N在测序数据中一般是不应该出现的，如果出现则意味着，测序的光学信号无法被清晰分辨，如果这种情况多的话，往往意味着测序系统或者测序试剂的错误。
5）接头序列
（1）测序接头与其在数据中的影响
在第 1 节中提到，测序文库构建是测序流程中的关键一步，而**测序接头（adapter）**正是在这个过程中被添加的。接头的作用有以下两点：
① 固定与区分样本
接头用于将 DNA 片段固定到测序平台的流动槽（flowcell）上。
在多样本混合测序时，接头上还包含特定的条形码序列（barcode），以便在测序后区分不同样本。
② 连接与扩增
接头为插入片段的末端提供连接点，便于后续扩增和测序。
（2）接头序列的测到与其影响
当测序 read 的长度大于目标 DNA 片段（插入片段）的长度时，read 的末尾可能会延伸到接头序列。
① WGS（全基因组测序）中的情况
在 WGS 测序中，这种现象较少发生。构建 WGS 文库时，插入片段通常较长（几百 bp），而测序 read 的长度一般为 100–150 bp。因此，read 通常不会测到接头序列。
② RNA 测序中的情况
在一些 RNA 测序（如 miRNA 测序）中，由于目标 RNA 片段本身较短（通常几十 bp），read 很容易测完整个片段并延伸到接头序列。这种现象在短片段 RNA 测序中尤为常见。
测到的接头序列会干扰下游分析，其影响类似于低质量碱基。未切除的接头序列可能导致：
比对错误，降低比对率和精确性；
错误的变异检测或基因表达量估计。
（3）接头序列的清理
在正式数据分析前，接头序列与低质量碱基一样，需要通过预处理软件清理。常用的清理工具包括：
① Trimmomatic：高效去除接头序列与低质量片段；
② Cutadapt：能够快速识别并剪切指定的接头序列；
③ FASTP：支持自动检测接头序列并进行质量控制。
清理后的数据可以显著提高下游分析的准确性，避免因接头干扰而消耗计算资源。
通过预处理确保测序数据的高质量是后续研究的基础，也是数据可靠性的保障。
四、FastQC该怎么用呢？
（1）通过网页搜索或者直接到它的主页上下载最新的版本。
也可在终端通过wget命令下载： $ wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip ./ 解压之后，修改文件夹中fastqc的权限，就可以直接运行了： $ unzip fastqc_v0.11.5.zip $ cd FastQC $ chmod 755 fastqc
（2）FastQC的运行非常简单，直接在终端通过命令行是最有效直接的，例如：
$ /path_to_fastqc/FastQC/fastqc untreated.fq -o fastqc_out_dir/
注意：-o 参数用于指定FastQC报告的输出目录，这个目录需要事先创建好，如果不指定特定的目录，那么FastQC的结果会默认输出到文件untreated.fq的同一个目录下。它输出结果只有两个，一个html和一个.zip压缩包。
$ tree fastqc_out_dir/
（3）我们可以直接通过浏览器打开html，就可以看到FastQC给出的所有结果，zip压缩包解压后，从中我们也可以在对应的目录下找到所有的QC图表和Summary数据。
除了上述用法之外，FastQC支持同时输入多个fq文件（或者以通配符的形式输入fq），当我们的fq文件比较多时，这种用法会比较方便，如：
$ /path_to_fastqc/FastQC/fastqc /path_to_fq/*.fq -o fastqc_out_dir/
这个链接是FastQC的一个结果模板：online report
五、切除测序接头序列和read的低质量序列
（1）接头序列与低质量碱基的切除工具及原理
目前市面上已有多种工具用于切除测序接头序列和低质量碱基，其中主流工具包括SOAPnuke、cutadapt、untrimmed等十余种，但在便捷性和功能性上，以下几种工具尤为突出：
① Trimmomatic：功能全面，支持接头序列切除与低质量碱基过滤。
② Sickle：专注于低质量碱基过滤，支持双端（PE）和单端（SE）数据。
③ Seqtk：快速处理低质量碱基，但仅支持单端（SE）数据。
（2）Trimmomatic 的优势
Trimmomatic不仅可以切除 Illumina 平台默认的接头序列，还支持切除用户自定义的接头序列，同时过滤 read 末尾的低质量碱基。相比之下，Sickle和Seqtk仅用于去除低质量片段，不支持接头切除。
① 工具工作原理
以Trimmomatic为例，其工作原理基于滑动窗口法：
- 在指定窗口长度内，计算窗口内碱基的平均质量值。
- 如果窗口的平均质量值低于阈值，从该窗口开始，将 read 的尾部全部切除，而不是挖掉局部低质量片段。
- 这种方式保证了 read 的连续性，避免人为改变实际的基因组序列。
特别注意：切除后的序列可能会显著缩短，因此在实际操作中需要设定最低 read 长度以避免无效数据进入下游分析。
② 工具选择建议
- 如果下机数据中包含接头序列，推荐使用Trimmomatic，它功能全面，且对多种场景适用。
- 如果不包含接头序列，且对计算资源要求较高，推荐使用Sickle（支持 PE 和 SE）或Seqtk（仅支持 SE），两者处理速度更快且资源占用较低。
(3) 使用 Trimmomatic 过滤数据
安装与准备
1. 获取 Trimmomatic
- 最新版本可从其官网下载。
- 如需查看源代码，可在GitHub上找到项目页面。
2. 安装与运行
- 下载打包的二进制文件（如 Trimmomatic-0.36.zip），解压后可看到trimmomatic-*.jar文件。
- 使用 Java 运行程序：java -jar trimmomatic-0.36.jar
3. 构建过滤流程
使用 Trimmomatic 时，通常通过命令行指定输入文件、接头序列文件和过滤参数。以下是一个基本流程示例：
ILLUMINACLIP：指定接头序列文件及参数（最大错配数、精确性评分等）。
SLIDINGWINDOW：设置滑动窗口长度和质量阈值。
MINLEN：过滤掉低于指定长度的序列。
4.输出文件解释
output_R1_paired.fastq.gz和output_R2_paired.fastq.gz：过滤后仍成对的 reads。
output_R1_unpaired.fastq.gz和output_R2_unpaired.fastq.gz：未成对的 reads，通常在后续分析中被丢弃。
接头序列的选择与设置
（1）接头序列的选择
① 接头序列
- HiSeq和MiSeq系列通常使用TruSeq3接头序列。
- TruSeq2接头序列主要用于较老的 GA2 测序仪，目前已较少使用。
- 这些默认接头序列信息可以在 Illumina 官网查询。
② 自定义接头序列
- 如果所用测序平台非 Illumina，或者需要切除自定义接头序列，可以自行创建符合格式的接头序列文件，并作为参数传入 Trimmomatic。
- 具体命名和格式要求可参考 Trimmomatic 文档（The Adapter Fasta）。
（2）PE 与 SE 模式选择
① 接头序列的使用需根据具体测序类型选择：
- **PE 模式（Pair End）：**双端测序，使用TruSeq3-PE.fa文件。
- **SE 模式（Single End）：**单端测序，使用TruSeq3-SE.fa文件。
Trimmomatic 的运行模式（两种）
① PE 模式：适用于双端测序数据。
PE 模式命令：java -jar trimmomatic-0.36.jar PE [参数] 输入文件 输出文件 过滤参数
示例命令
PE 模式（HiSeq PE 测序）
java -jar /path/Trimmomatic/trimmomatic-0.36.jar PE \ -phred33 -trimlog logfile \ reads_1.fq.gz reads_2.fq.gz \ out.read_1.fq.gz out.trim.read_1.fq.gz \ out.read_2.fq.gz out.trim.read_2.fq.gz \ ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq3-PE.fa:2:30:10 \ SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
② SE 模式：适用于单端测序数据。
SE 模式命令：java -jar trimmomatic-0.36.jar SE [参数] 输入文件 输出文件 过滤参数
SE 模式（HiSeq SE 测序）
java -jar /path/Trimmomatic/trimmomatic-0.36.jar SE \ -phred33 -trimlog se.logfile \ raw_data/untreated.fq out.untreated.fq.gz \ ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq3-SE.fa:2:30:10 \ SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
参数说明
常用参数
- phred33 或 -phred64 指定碱基质量值体系。
- Phred33是当前测序数据的主流格式，必须显式指定此参数。
- 若未指定，默认使用Phred64，可能导致错误处理。
- trimlog 指定修剪日志文件，记录每个 read 的修剪信息。
- ILLUMINACLIP
指定接头序列文件和参数，格式为：ILLUMINACLIP:<文件路径>:<最大错配数>:<精确性阈值>:<剪切长度>
- SLIDINGWINDOW
滑动窗口参数，格式为：SLIDINGWINDOW:<窗口大小>:<质量阈值>
- LEADING 和 TRAILING 分别指定去除 read 起始或末端低质量碱基的质量阈值。
- MINLEN 设置修剪后 read 的最低长度，低于此值的 read 将被丢弃。
PE 与 SE 模式的差异
PE 模式
- 需要输出过滤后和未过滤的双端 read 数据：
- out.read_1.fq.gz和out.read_2.fq.gz：过滤后的配对 reads。
- out.trim.read_1.fq.gz和out.trim.read_2.fq.gz：未通过过滤的 reads。
SE 模式
- 不需要指定未通过过滤的 reads 文件。过滤后直接生成一个输出文件。
- 例如，在 SE 模式下：raw_data/untreated.fq --> out.untreated.fq.gz
- 通过灵活设置接头序列和运行模式，Trimmomatic可以高效完成测序数据的预处理，为下游分析提供高质量的输入数据。
以下是对Trimmomatic参数的整理：
Trimmomatic 参数的整理
ILLUMINACLIP 参数解释
- 接头序列文件：指定用于匹配的接头序列文件（如TruSeq3-PE.fa或自定义文件）。
- 最大错配数：允许接头序列匹配时的最大错配数，值越小匹配越严格。
- 双端 reads 匹配度：两条 reads 同时与接头序列匹配时，总体匹配度的最低阈值（推荐值为 30%）。
- 单端匹配度：单条 read 与接头序列部分匹配时的最低匹配度（推荐值为 10%）。
**注：**实际测序数据可能只覆盖接头序列的一部分，因此不要求完全匹配。上述推荐值（30% 和 10%）是常用参数。
以上就是数据质控的内容结束。
生物信息学领域非常广泛，难以一次说尽。我们下次继续更新，一起深入学习生物信息学的内容！
喜欢的宝子们点个赞吧~码字不易，且行且珍惜~