NGS数据过滤的详细描述
NGS原始数据的过滤对于后续分析非常重要，去除一些无用序列也可以提高后续分析的精度和效率。Trimmomatic是一款功能强大的数据过滤软件。
Trimmomatic发表的文章至今已被引用2810次，是Illumina平台上比较受欢迎的数据过滤工具。来自其他平台的数据，如Iron torrent和PGM测序数据，可以通过fastx_toolkit和NGSQC toolkit进行过滤。
Trimmomatic支持多线程，数据处理速度快。主要用于去除Illumina平台Fastq序列中的连接子，根据基质量值修剪Fastq。软件有两种过滤模式，分别对应SE和PE测序数据，同时支持gzip和bzip2压缩文件。
此外，它还支持phred-33和phred-64格式之间的相互转换。现在phred-33和phred-64格式的混淆是Illumina的锅造成的(Daminyou，Illumina！)，但现在Illumina平台的大部分输出数据也转换为phred-33格式。
Trimmomatic过滤数据的步骤与命令行中过滤参数的顺序有关。通常的过滤步骤如下:
因为Trimmomatic过滤数据的步骤与命令行中过滤参数的顺序有关，所以如果有必要的话，建议第一步就去连接器，否则连接器序列被其他过滤参数切断后，匹配和移除连接器序列会更加困难。
在SE模式下，只有一个输入文件和一个过滤后的输出文件:
-trimlog参数指定过滤器日志的文件名，它包含以下四列:
由于生成的trimlog文件包含了每次读取的处理记录，文件大小巨大(GB级)。如果以后不再使用Trimlog，建议不要使用trim参数。
在PE模式下，有两个输入文件，正向排序序列和反向排序序列，但过滤后有四个输出文件。过滤后，两个末端序列保持成对，而如果过滤后丢弃一个末端序列，另一个末端序列保持不成对。
其中，-phred33和-phred64参数指定了fastq的质量编码格式。如果未设置该参数，软件将自动确定输入文件的格式(v0.32以后的版本支持)。虽然软件默认的参数是phred64，但是如果不确定序列是什么质量编码格式，就不需要设置这个参数。
PE模式的两个输入文件:sample _ r 1 . fastq sample _ R2 . fastq和四个输出文件:sample _ paired _ r 1 . clean . fastq sample _ unpaired _ r 1 . clean . fastq sample _ pair。ed _ r 1 . clean . fastq sample _ unpaired _ r 1 . clean . fastq
通常PE测序的两个文件R1和R2的文件名差不多，可以用-basein参数指定R1的文件名，软件会推断出R2的文件名，但是这个函数不好用，因为软件只能自动识别推断出-basein:
建议指定两个文件名(R1和R2)作为输入，而不使用-basein参数。
有四个输出文件。当然可以如上指定四个文件名，但是参数太长，有点麻烦。有一个省心的方法。使用-baseout参数指定输出文件的基本名称，软件将自动命名四个输出文件。比如，-base out mySampleFiltered.fq.gz，在文件名后面加一个. gz后缀，软件会自动用gzip压缩输出结果。四个输出文件将自动命名为:
此外，如果直接指定输入和输出文件名，在文件名后添加一个. gz后缀会告诉软件输入文件是一个. gz压缩文件，输出文件需要用gzip压缩。
如果过滤的每一步都需要多个参数，那么通常使用冒号来分隔参数。当然，参数的顺序是需要的。
顾名思义，这一步是拆下illumina连接器。这个软件其实是为illumina平台数据设计的。
为了更好的理解测序读数中为什么会有引物和接头序列，我画了一个加入接头后文库的结构示意图，也大致标出了引物结合位点:
理解了这个文库的示意图结构，就很容易理解测序过程了。
去除接头和引物序列看似简单，但需要平衡灵敏度(保证接头和引物去除干净)和特异性(保证接头和引物以外的序列不被误删除)。由于测序中可能的随机错误，去除接头的简单操作变得复杂。
虽然理论上，接头序列和引物序列可能出现在阅读中的任何地方，但事实上，在大多数情况下，接头和引物出现在序列中，因为文库中插入的片段比测序阅读的片段短。在这种情况下，在阅读的开始处有一个可用的序列，并且结尾包含接头序列的全部或部分。如果末端只含有部分接头序列，就不容易去除这个不完整的接头序列。
然而，在PE测序模式中，如果文库的插入片段比测序读数短，那么read1和read2中的非接头部分将完全反向互补，并且Trimmomatic具有“回文”模式，该模式将利用该特征去除接头序列。
下图中的四种情况A、B、C和D是接头和引物的四种修剪去除模式:
模式A:测序读数从起始位置开始包含一个完整的接头序列，所以根据Illumina测序原理，整个读数不可能包含有用的序列，整个读数被丢弃。
模式B:这个比较常见。因为文库中插入的片段比测序读数短，它将在读数的末端包含一些接头序列。如果这个连接子序列足够长，它可以被识别并移除，但是如果连接子序列太短，短于连接子匹配参数设置的最短长度，它就不能被移除。然而，在PE测序的情况下，可以根据D模式去除阅读结束时的短接头序列。
c模式:这可能发生在PE测序中。一些正向测序和反向测序是完全互补的，但是卸载的文库和两个接头是直接互连的。这样的读数不包含任何有用的序列，正向和反向测序读数都被丢弃。
d模式:这是Trimmomatic使用PE测序去除短接头序列的模式。如果文库中插入的片段比测序读数短，可以利用正向和反向测序读数中的一个碱基可以完全互补的特点，将两个接头序列与读数进行比较，同时也将两个读数进行相互比较，这样即使是1bp的3’端的接头序列也可以准确地去除，这比B模式下更彻底。
Trimmomatic使用类似于序列比对软件(如超快速比对软件Isaac aligner)的两步策略来搜索潜在的联合序列。首先，将接头序列中的一段种子序列(种子长度不超过16bp)与测序读数进行比较。如果种子序列在测序读数中具有良好的比较结果(具体由种子错配参数决定)，则开始与读数进行接头长度比较的第二步。种子搜索的第一步非常快，可以过滤掉没有联合污染的读数。这种两步搜索方法使得联合序列的搜索非常高效。
在第二步中，当比较接头序列的全长和测序读数时，惩罚策略考虑测序碱基的质量值Q，比较中每个碱基为0.6，每个错配的碱基为Q/10。考虑碱基的质量值可以减少错配的低质量碱基(测序错误率高)对整体比较分数的影响。根据这一规则，12bp接头序列的片段与读数完全比对，得分为7.2，25bp接头序列的片段与读数完全比对，得分为15。因此，建议ILLUMINACLIP参数中简单剪辑阈值的取值在7-15之间(即上图中A/B比较模式的比较分数阈值)。
对于回文模式比对(上图中的D模式)，可以比对的序列长度会更长，比对分数的阈值会更高，以保证识别接头序列的准确性。比如R1和R2中有50bp的序列可以反向匹配，分值为30。在这种模式下，Trimmomatic可以识别并删除读数中非常短的接头序列。
ILLUMINACLIP参数说明:按照指定的顺序列出了ILLUMINACLIP的参数如下(每个参数之间用冒号隔开)，PE排序时要注意最后一个参数。SE序列的最后两个参数可能无法设置。
FastaWithAdaptersEtc:指定包含链接器和引物序列(所有序列都被视为污染)的fasta文件路径。Trimmomatic自带fasta文件，包含Illumina平台的连接子和引物序列，可以直接使用。
SeedMismatches:指定种子搜索第一步中允许的错配碱基数，例如2。
回文剪辑阈值:指定回文剪辑模式下PE需要多少R1和R2(上图中模式D)的比较分数，例如30。
简单剪辑阈值:指定切除关节序列的最低比较分数(上图中的A/B模式)，通常在7-15之间。
MinAdapterLength:只对PE测序的回文clip模式有效，指定回文模式下可以切除的接头序列的最短长度。由于历史原因，默认值为8，但实际上回文模式可以切除短至1bp的连接子污染，所以可以设置为1。
KeepBothReads:只在PE测序的回文剪辑模式下有效。这个参数非常重要。在上图的D模式中，去除接头序列后，R1和R2的剩余部分完全互补。默认参数为false，这意味着与R1完全互补的R2被完全去除。将其作为重复项删除，但在某些情况下，如bowtie2过程的成对读取，您应该将此参数更改为true，否则您将丢失一些成对读取。
看一个PE150数据的测试，就知道keepBothReads参数的重要性了:
滑动窗口切割，统计滑动窗口内所有碱基的平均质量值，如果低于设定的阈值，则切断窗口。
滑动窗口参数如下:
WidowSize:设置窗口大小。
RequiredQuality:在窗口中设置平均基础质量的阈值。
它包含一个可以自动调整的过滤条件，并在尽可能长时间地保留序列和尽可能低地保留碱基测序错误率之间取得平衡，以便在修剪后最大化保留序列的价值。
对于不同的应用场景，读取序列的值由以下三个因素决定:
MAXINFO有两个参数。第一目标读取长度控制上述第一个因素，即最短的允许读取长度。第二个参数，严格性，是控制因子2和因子3的平衡，即在满足最短阅读长度的条件下，是保留尽可能多的碱基，还是保证最低的测序错误率。
MAXINFO filtering从读数的3’端开始，在考虑以上三个因素的情况下，统计所有可能的trim模式的信息分数(如此干净的读数值),这三个因素以不同的方式影响最终的读数信息分数:
INFO分数是为任何可能的读数切割模式计算的，读数的最终长度和要切割的碱基由INFO的最大值决定。事实上，这三个影响因素以不同的方式发挥作用:
参数描述:
TargetLength:能使读数映射到参考序列上唯一位置的最短长度。
严格度:一个介于0-1之间的小数，决定了如何平衡读取长度的最大化或读取错误概率的最小化。当参数设置小于0.2时，倾向于最大化读取长度，当参数设置大于0.8时，倾向于最小化读取中的测序错误概率。
从读取开始，质量值低于设定阈值的碱基被切除，直到质量值达到阈值。
质量:设置基本质量值的阈值，低于该阈值将被切断。
从读数结束开始，质量值低于设定阈值的碱基被切除，直到一个碱基达到阈值。Illumina平台上的一些低质碱基标记为2，那么设置为3就可以过滤掉这些低质碱基。官方推荐使用滑动窗口或MaxInfo，而不是LEADING和TAILING。
质量:设置基本质量值的阈值，低于该阈值将被切断。
不考虑碱基质量，从阅读开始保留设定的碱基长度，所有其他的都被切除。将所有读数切割成相同的长度。
Length:从末端开始分割后，reads保留的序列长度。
不考虑碱基质量，一些碱基从阅读开始就被直接去除。
Length:从读取开始时移除的碱基数。
设置最小读取长度。在前面过滤读取时，如果保留长度低于此阈值，则整个读取将被丢弃。丢弃的读取数将计入Trimmomatic日志中的丢弃读取数。
Length:可以保留的最短读取长度。
该选项可以将过滤后的Fastq文件中的质量值行转换为phred-33格式。
此选项可以将过滤后的Fastq文件中的质量线值转换为phred-64格式。
Trimmomatic也可以制作自己的fasta文件，包含连接子和引物序列，格式可以参考软件自带的adapters文件夹中的格式。
adapters文件夹包含用于se和PE的illumina SEquencing TruSeq2和TruSeq3通用接头和引物序列。