为什么要进行数据过滤

二代测序的下机数据经碱基识别得到的原始序列称为raw-data。raw-data中可能含有接头（adapter）、低质量序列及N序列（N 表示无法确定碱基信息）。过滤的目的就是通过设定一定的阈值，去除去除含N比例高，接头污染或者低质量占比高的reads。数据过滤的工具很多，这里我们简单介绍Trimmomatic。
Trimmomatic可以去除fastq数据中的接头，也可以用来切除reads中低质量的序列，软件的适用范围广：同时支持单端（SE）和双端（PE）测序数据；同时支持 gzip 和 bzip2 压缩文件；同时支持 phred-33 和 phred-64 格式。另外软件支持多线程，处理数据速度快。

资源

官网：USADELLAB.org - Trimmomatic: A flexible read trimming tool for Illumina NGS data
Github：GitHub - usadellab/Trimmomatic
说明文档：TrimmomaticManual_V0.32.pdf (usadellab.org)

使用示例

以双端测序（PE）为例，过滤逻辑如下：

image.png

使用示例如下：

java -Djava.io.tmpdir=./ -Xmx20g -jar /work/app/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 5 -phred33 S.raw.1.fastq.gz S.raw.2.fastq.gz S.1.clean.fq.gz S.1.unpaired.fq.gz S.2.clean.fq.gz S.2.unpaired.fq.gz ILLUMINACLIP:/work/app/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:51
#其中：
#PE ：双端测序（也可SE）
# -threads 5 ：线程数
# -phred33 ：转成phred33质量体系（也可是phred64）
#S.raw.1.fastq.gz S.raw.2.fastq.gz ：raw data，双端测序有两个文件，正向和反向测序序列。
#S.1.clean.fq.gz S.1.unpaired.fq.gz S.2.clean.fq.gz S.2.unpaired.fq.gz ：输出文件名，双端测序有4个文件（过滤之后双端序列都保留的就是 paired，反之如果其中一端序列过滤之后被丢弃了另一端序列保留下来了就是 unpaired ）
#ILLUMINACLIP:/work/app/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10:8:true ：过滤 reads 中的 Illumina 测序接头和引物序列，参数格式如下：ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads><keepBothReads>。“ILLUMINACLIP”表示去除illuminate的接头和引物：“fastaWithAdaptersEtc”接头和引物的序列：“seed mismatches”seed 搜索允许2个碱基错配：“palindrome clip threshold”palindrome 比对分值阈值 30：“simple clip threshold”simple clip 比对分值阈值 10：“minAdapterLength”palindrome 模式允许切除的最短接头序列为 8bp：“keepBothReads”palindrome 模式去除与 R1 完全反向互补的 R2。其中palindrome clip是 Trimmomatic 利用 PE 测序进行短接头序列去除的典范，如果文库插入片段比测序读长短，利用正反向测序 reads 中一段碱基可以完全反向互补的特点，将两个接头序列与 reads 进行比对，同时两条 reads 之间也互相比对，可以将 3' 末端哪怕只有 1bp 的接头序列都可以被准确去除，该模式去除接头污染更彻底；simple clip包含两种模式，A 模式：测序 reads 从起始位置开始就包含了完整的接头序列，那么根据 Illumina 测序原理，这整条 reads 都不可能包含有用序列了，整条 reads 被丢弃。B 模式：这种相对常见，由于文库插入片段比测序读长短，会在 reads 末端包含部分接头序列，若是这部分接头序列足够长是可以识别并去除的，但如果接头序列太短，比接头匹配参数设置的最短长度还短，那么就无法去除。但是，如果是 PE 测序，可以按照 palindrome clip 模式去除 reads 末端的很短的接头序列。
#LEADING:3 ：切除reads 5'端质量值低于阈值（3）的碱基。
#TRAILING:3 ：切除 reads 3'端质量值低于阈值（3）的碱基。
#SLIDINGWINDOW:4:20 ：从 reads 的 5' 端开始，以n个碱基为窗口进行滑窗（此处n=4），切如果窗口碱基质量平均值低于阈值（20），则切掉从当前窗口到3'端的序列。
#MINLEN:51：如果经过剪切后 reads 的长度低于阈值（51）则丢弃这条 reads。

注：命令的顺序决定了分析的顺序（如去除接头序列和去除低质量碱基），建议先进行接头切除