1.用法

Trimmomatic基本使用见https://my.oschina.net/u/4150954/blog/3065527

2.原理讲解

两种不同的模式

1. Simple Model

每个Read都会被从5‘端到3’端进行浏览，为了确认是否其中含有已知的接头序列（在程序中已经设定了，有相关的文件专门记录这些接头序列，接头序列信息可以通过https://github.com/csf-ngs/fastqc/blob/master/Contaminants/contaminant_list.txt进行查看）。软件中本身内置的接头文件为adapter目录中，可在使用时进行文件指定。

在进行这些人造序列和reads之间的比对时，用到的是seed and extend策略(参考文献https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2943993/），seed不用完全匹配，可以根据人为设定的阈值来限制允许错配的碱基数，之后会通过局部比对的方法进行人造序列（或者理解为污染序列）和read的匹配。

image

A. 如果人造序列的3‘端和read的5’端有重合，则整个read都会被扔掉，因为整个read的浏览过程都是从5‘端开始的。

B. 如果人造序列的与5’端开始都后面全部匹配的话，整个read被扔掉。

C. 如果人造序列在read原序列中是被包含关系，那么从read的5‘端到人造序列的开头这一段会被保留。

D. 人造序列还会沿着这个方向继续比对，只有当如图D中时，如果人造序列直到read的3’端结束之后才匹配上的话，那这个read才会被完整保留。

如果污染序列（任何人为造成的）匹配到read足够长时，则可以精确地检测出污染序列，但是如果污染序列跟read匹配到的较短时，则在A和D可能中，就可能被忽略了。

2. Palindrome Mode（回文模式）

回文模式专门对于“接头 read贯穿”（adapter read-through：可以理解为read测穿的情况）的情况做了优化，当这种read 测穿模式出现时，意味着一对中的两条read应该含有相同数量的有效碱基，且这些有效碱基应该是反向互补的，也就意味着，这两条read中可能伴随着由两个完全相反（反向互补）的接头造成的污染序列。所以如果在算法中同时根据这三个条件进行测定的话，能够提高特异度和敏感性。

在这个模式中，接头序列会先与read形成混合序列，然后与双端测序中的另一条反向链形成的混合链相互比对。如果比对结果打分较高的话，说明两条链的read的第一部分（应该可以理解为，read本身序列部分吧）都很好的反向互补了，而read的剩余部分匹配到另一个的adapter上去了。

同样，在这种模式中用到的也是seed and extend策略，不过这种情况下会进行全局比对。

image

图中红色的序列会被切除。

这种模式相较为上面的简单模式而言，有一点较大的优势，能够比对更长的序列（因为简单模式中是adapter和 整个read比，而这个模式下是整个），这样就能保证就算存在因为测序错误导致接头序列只出现了一部分的这种情况下也能够准确的测出接头。如果真的有需求的话，这种模式还可以在保证较低假阳率的情况下用切除单个接头碱基。

但是这种模式的最大弊端就是只用在测穿的情况中。

对比中的算法

• 在最开始的比对中，首先会取每个序列中的16个碱基，之后将这16个碱基转化为64bit，A = 0001, T = 0010, C = 0100 and G = 1000；然后进行位异或运算，比如A-T错配时：XOR（0001，0010）=0011，比对完成后会对1进行计数，得到的数将是比对结果中的错配数的2倍。

• 如果是匹配到的话，则会打分log10(4)~=0.602

• 如果是错配的话，则会根据碱基质量了进行罚分：Q/10=[1,4],因为碱基的质量值就在10-40之间，所以碱基质量值越高的碱基错配后的罚分越厉害。

质量过滤

Trimomatic中提供两种质量过滤方式，这两种情况中一般都是会根据碱基质量值保留5'端的数据，而剪切掉3’端的数据，这是符合Illumina的测序情况的，因为在Illumina中一般都会是3'端的数据质量较差，接下来看两种不同的方法：

• 1.Sliding Window quality filering

  这种方法中会从5'端向3'端读取read，如果当这一组碱基中的平均质量值低于域质的时候就会切除3'端的碱基。

• 2.Maximum Information quality filtering

  • 在这种方法中提到了保留额外碱基的概念，因为有一些端的read是提供不了任何有效筛选信息，因为他们可以在read中出现多次，所以提供的也是一些无用信息，为了我们能够更精确地对数据进行过滤，那应该在早期保留一些额外的碱基来保证能够让read的足够长，能够提供一些有效信息以供低质量碱基的过滤。

  • 但是有一些read本身就是足够长了。如果在早期又保留一些额外的base只会让这个read质量更差。

  所以，综合以上，我们在根据read的信息进行过滤的时候应该进行多方面考虑，而不是简单的给个保留base的阈值。

  • 在软件中通过三个因素的综合考虑来确定碱基保留长度：

  1. "长度阈值"，假设给定一个目标长度为t，最后需要剪切到长度l，那么最后会根据这个长度阈值来给出一个长度阈值打分：

  
  image.png

  2. “覆盖度”，会在假定保留的read长度的基础上给出一个覆盖度打分

  image.png

  3. “错误率”模型，根据read中base的质量分数，会将每个碱基的错误率相乘得到这个模型的打分!

  image.png

  将会根据这三个因素打分来综合确定最终需要剪切掉多少的read。用户可以设置此值的阈值s（0-1），因该均衡考虑覆盖度（当s=0时最大）和错误率因素（当s=1时最大）

  计算方法如图

image.png

示例：

image.png

图中第四个图的peak值就是需要剪切的read点。