Trimmomatic文献解读

《Trimmomatic:一款应用于Illumina测序数据的裁剪工具》2014年

摘要

1、目前的NGS处理软件不能满足灵活、准确处理 PE 测序数据的要求，因此开发了Trimmomatic 软件；
2、Trimmomatic的NGS预处理性能在有参、无参任务中进行了验证。在多数测试场景中，Trimmomatic 性能优于其他软件；
3、软件使用和说明：http://www.usadellab.org/cms/index.php?page=trimmomatic

1 引言

1、NGS数据需要进行适当的预处理。如果不进行处理，低质量序列和所含的接头会影响下游的数据分析；过量处理，又会增加生信分析的成本
2、预处理方法也不能干扰下游数据的分析。例如，NGS数据多采用PE测序，正反向reads存储在两个单独的FASTQ文件中，并且顺序相同。许多下游数据处理工具会使用这种对位关系，因此预处理序列数据时必须维护这种关系
3、目前存在的NGS处理工具（无论单独或者组合使用），没有一个能同时兼顾灵敏度和性能，而且大多数都不是为处理PE数据设计。Trimmomatic能解决上述问题

2 算法

Trimmomatic 包含一系列 reads 裁剪和过滤步骤，但主要的算法创新是确认接头序列和质量过滤，详细过程如下。其他一些处理过程，见补充材料。

2.1 去除人工序列

Trimmomatic 采用两种方式去检测 reads 中的人工序列。第一种方式，称为 ‘simple mode’ ，它主要通过在 reads中 匹配用户提供的序列来查询。它的优势是可以适用于任何人工序列，包括接头，PCR引物。这个种方式可以检测 reads 中任何位置和方向的人工序列，但是该方法需要一个最小的重叠区域来避免假阳性。但是，一些较短的接头序列一般存在于 reads 的末端，不可避免的达不到最短重合区域的要求，因此不能被检测到。

第二种方式，称为 ‘palindrome mode’，专门用于检测 ‘adapter read-through’ 场景，当测序DNA片段小于 reads 长度，那么reads的末端就会被接头污染。这种情况主要应用于 reads 更长的 Miseq 测序。虽然过短的插入片段应该在文库制备的过程中删除，但是实际过程中并不能完全删除，因此许多文库都存在这个样的问题。 ‘palindrome mode’ 只适合应用于 PE 测序，但是相比于 ‘simple mode’ 有其自身的优越性。

2.1.1 simple mode

在简单模型中，会对每条read进行5'至3'端遍历，去判断是否有用户提供的接头存在。会使用种子序列定位及延伸算法（seed and extend ）定位人工序列和reads的初始匹配。种子序列不需要完全匹配，用户可以定义mismatch的数量。进行局部比对，如果比对的分数超过一定的阈值，则删除比对后的区域。Fig1A-D 展示了人工序列从5'至3'端与 reads 局部比对的情况。A,B 情况会删除整条read，C,D情况会保留5'端未比对的reads。这种模式对于匹配的人工序列较短时，不能有效的检测。
Simple mode是根据与污染序列的比对情况判断切除序列：
A) 污染序列3'端的一部分与read5'端有overlap，去除整条read；
B) 污染序列从read的5'端开始能与read完全overlap, 去除整条read；
C) 污染序列在read中间有overlap，5'端一部分没有overlap的保留，overlap的部分及之后的进行切除；
D) 污染序列3'端的一部分与read有overlap，切除overlap部分

2.1.2 Palindrome mode

当 ‘read-through（read测穿）’时，配对reads中会含有相同长度的有效碱基，还有相反方向的污染接头，并且有碱基序列是互补配对的，当同时检测到这三种情况， ‘read-through’ 情况可以高灵敏度和特异性的确定。采用种子序列定位及延伸算法进行比对，匹配分数较高的部分是read的互补配对部分，剩余的部分是接头部分。回文模式的主要优点是较长的对齐长度，这确保了接头序列可以被可靠的检测到。由于回文模式只限于检测read测穿的情况，需要综合‘simple’和‘palindrome’模式共同进行分析。
Palindrome mode: 是根据a-R1, a-R2的overlap判断切除的序列：
A) overlap是在接头序列与reads的起始之间，即overlap只有接头序列，没有有用的信息，这种情形将这两条read都整条去除；
B) overlap有一部分是DNA插入片段，一部分是整个接头序列，该情形会将这对reads从overlap序列及剩余序列(3'端)切除；
C) overlap中有接头序列的一小部分，与B)相同方式切除；
D) overlap中不存在接头序列

2.1.3 比对检测和打分

match：匹配一个碱基加分，分值为(约等于0.602)；
mismatch: 错配一个碱基减分，分值为Q/10；
其中，Q为对应碱基的碱基质量值，由于一般碱基质量值区间为: [0-40]，即错配的罚分值区间[0-4]。
该分值公式说明碱基质量越高的碱基出现错配时罚分越多

例如：
12bp的接头序列能够完全比对到read上的分值为：12 * 0.602 ≈ 7；
25bp的接头序列完全比对到read上的分值为：25 * 0.602 ≈ 15；
a-R1和a-R2有50bp完全匹配的分值为50 * 0.602 ≈ 30

2.2 质量值过滤

Trimmomatic 提供了两种质量过滤策略。两种方法利用每个碱基的Illumina质量值来确定read的哪些区域应该被切除，保留5'端部分，切除3'端部分。这与典型的Illumina数据非常相符，它通常在3'末端的质量较差。下面部分将对这两种方法进行详细描述。

2.2.1 滑动窗口碱基质量过滤

滑动窗口利用了一个相对标准来进行质量过滤。它会遍历read 的5'端至3'端，当 group 中碱基的平均质量值低于设定的阈值时，会切割掉3'端的read。这种方法可以防止单个低质量碱基引起的后续高质量碱基的删除，同时也能确保连续的低质量碱基能被删除。

2.2.2 最大信息碱基质量过滤

在许多应用中，保留 read 中额外碱基的原因是与read 的长度有关。因为，短的 read 几乎毫无价值，它在目标序列中多次出现，只能提供模棱两可的信息。即使冒着引入错误的风险，在 read 过滤早期保留额外的低质量碱基也是值得的，这样修剪之后的 read 就足够长，可以提供尽可能多的信息。当然，超过一定的 read 长度，保留额外的碱基可能是有害的。因此必须权衡好保留 read 长度和去除低质量碱基之间的关系。因此，在遍历 read 过程中，修剪过程越来越严格是值得的，而不是应用固定的质量值过滤。

最大信息碱基质量值过滤 ，它结合三个因素来决定每次 read 应该保留哪些内容。
1、 第一个因素模拟了 “长度阈值” 的概念，即 read 必须达到最小长度才能应用于下游程序。 非常短的 read 可能在基因组上有多个位置，会造成信息模糊混乱。另一方面，一些长 read 也可以比对到基因组的多个位置，因为它们来自于一个重复的区域，那么少量额外碱基就不可能解决这个问题。对于这两种极端情况，少量额外碱基的边际收益是相当可观的，因为这些额外的碱基会造成 信息模糊 和 有效read 之间的差异。
选择逻辑曲线来模拟这种评分行为，它在极端值时较为平坦，在用户指定的阈值附近较为陡峭。如下公式，t 代表给定目标长度，l 代表修剪的长度，score 代表最终评分。

2、第二个因素模拟了“coverage”的概念，根据保留的序列长度提供线性评分。这反映了较长的 read 包含对大多数应用程序有用的信息。

3、第三个因素模拟了 “错误率” 的概念，根据 read 中质量分数对应的错误概率来确定 read 中累积错误的可能性。为了计算这个分数，这里简单地计算每个碱基正确概率的乘积。每个碱基的正确概率值为 Pcorr ，错误分值通常在 read 的开始位置较高，后续达到某个点时迅速下降。

最大信息碱基质量过滤算法综合考虑了三个因素，最终得分决定了 read 的修剪方式。用户可以设定严格参数 s（在0和1之间） , 它可控制第二因素和第三因素的平衡。最终公式如下：

Fig3 展示了3个因素合并成最终结果的过程。

3 执行

Trimmatic 使用pipeline架构，各个步骤（去除接头，质量值过滤）可以独立或者组合分析。

4 结果

为了展示预处理步骤的作用，使用两个不同的场景进行测试：基于参考基因组的比对Escherichia coli K-12/MG1655（Bowite2 和 BWA 软件） ，从头组装E.coli K-12/MG1655 datasets （Velvet 软件）。

4.1 基于参考基因组比对

4.2 从头组装

Trimmatic处理之后，组装的N50和最大contig值有提升，adapter进行了去除，不影响组装。

4.3 与现存工具的对比

5 讨论

5.1 read进行预处理的必要性

文章的上面的数据展示了NGS 预处理步骤在参考基因组比对和从头组装中的优势。对于高质量的数据，在基于基因组比对的应用中，预处理的作用好像有限。唯一比对的reads只有小部分的提升（~1.5%）。然而，在实践中，这样的高质量数据集，对下游应用程序（变异检测）的提升可能是比较小的。
第二个数据集的 read 质量要低得多，它表明即使是基于参考基因组比对的任务也可以从 read 预处理中受益匪浅。未经预处理的 BWA 只能对不到25%的读取进行比对。通过预处理步骤能提高到近80%，即使在严格的设置下也有近78%的一致性。
在从头组装的场景下，需要修剪以确保接头序列不会被合并到新组装的基因组中。在两个数据集上，N50分别有58%/77%的显著改善，最大contig 大小分别有28%/55%的显著改善。统计数据的实质性改进进一步证明了对从头组装进行 read 预处理的合理性。

5.2 与其他工具的比较

在测试中，Trimmomatic优于所有其他工具。当测试数据集质量较高时，各种工具差异不大，当数据集质量较低时，Trimmomatic 有明显的优势，特别是“最大信息碱基质量过滤”模式。

6 参考文献

[1] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 Aug 1;30(15):2114-20. doi: 10.1093/bioinformatics/btu170. Epub 2014 Apr 1. PMID: 24695404; PMCID: PMC4103590.