Trimmomatic基本介绍
FastQC做质控时有的raw reads质量不好,会存在overrepresented reads,也有人反馈给小编说FastQC只是一个可视化的软件,并不包含QC (quality control ) 的过程,不会改变reads序列本身,这话没毛病。这时候就轮到trimmomatic登场了。Trimmomatic是一个主要针对Illumina测序的reads修剪和过滤的软件,包括paired end (PE) reads和single end (SE) reads。它支持多线程,处理数据速度快,过滤模式也有两种,分别对应 SE 和 PE 测序数据,也支持 gzip 和 bzip2 压缩文件。
Trimmomatic下载安装
mkdir Trimmomatic && cd Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip
unzip Trimmomatic-0.36.zip
java -jar ~/biosoft/Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h
Trimmomatic过滤步骤
Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关,通常的过滤步骤如下:
ILLUMINACLIP: 过滤 reads 中的 Illumina 测序接头和引物序列,并决定是否去除反向互补的 R1/R2 中的 R2。
SLIDINGWINDOW: 从 reads 的 5’ 端开始,进行滑窗质量过滤,切掉碱基质量平均值低于阈值的滑窗。
MAXINFO: 一个自动调整的过滤选项,在保证 reads 长度的情况下尽量降低测序错误率,最大化 reads 的使用价值。
LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
CROP: 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
HEADCROP: 从 reads 的开头切掉指定数量的碱基。
MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。
Trimmomatic简单用法
由于 Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关,因此,如果需要去接头,建议第一步就去接头,否则接头序列被其他的过滤参数剪切掉部分之后就更难匹配更难去除干净。
single end 测序
在 SE 模式下,只有一个输入文件和一个过滤之后的输出文件:
java -jar <path to trimmomatic jar> SE -threads <threads> [-trimlog<logFile>] <input> <output> <step 1> <step 2> ...
-trimlog 参数指定了过滤日志文件名,日志中包含以下四列内容:
read ID
过滤之后剩余序列长度
过滤之后的序列起始碱基位置(序列开头处被切掉了多少个碱基)
过滤之后的序列末端碱基位置
序列末端处被剪切掉的碱基数
由于生成的 trimlog 文件中包含了每一条 reads 的处理记录,因此文件体积巨大(GB 级别),如果后面不会用到 trim 日志,建议不要使用这个参数。
paired end 测序
在 PE 模式下,有两个输入文件,正向测序序列和反向测序序列,但是过滤之后输出文件有四个,过滤之后双端序列都保留的就是 paired,反之如果其中一端序列过滤之后被丢弃了另一端序列保留下来了就是 unpaired 。
${name}1.fq.gz ${name}2.trim.fq.gz 为输入文件;${name}R1.clean.fq.gz ${name}R1.unpaired.fq.gz ${name}R2.clean.fq.gz ${name}R2.unpaired.fq.gz 为对应的四个输出文件;
ILLUMINACLIP:"$adapter"/Exome.fa:2:30:9:1:TRUE,这部分指定 2 种去接头模式的参数:"$adapter"/Exome.fa 指明需要匹配的接头文件,2 代表 16 个碱基长度的种子序列中可以有 2 个错配,30 代表采用回文模式时匹配得分至少为30 (约50个碱基),10 代表采用简单模式时匹配得分至少为10 (约17 个碱基);
LEADING: 20,从序列的开头开始去掉质量值小于 20 的碱基;
TRAILING:20,从序列的末尾开始去掉质量值小于 20 的碱基;
SLIDINGWINDOW:4:15,从 5' 端开始以 4 bp 的窗口计算碱基平均质量,如果此平均值低于 15,则从这个位置截断 read;
MINLEN:36, 如果 reads 长度小于 36 bp 则扔掉整条 read。
Trimmomatic过滤原理
去除接头以及引物序列看似简单,但需要权衡灵敏度(保证接头和引物去除干净)和特异性(保证不是接头和引物的序列不被误切除),由于测序中可能存在的随机错误让去接头这样一个简单的操作变的复杂。虽然理论上接头序列和引物序列可能出现在 reads 中的任何位置,但实际上序列中出现接头和引物大部分情况下都是由于文库插入片段比测序读长短导致的,这种情况在 reads 的开头部分是有一段可用序列的,末端包含了接头的全长或部分序列,如果末端只有接头的一部分序列,那么去除这残缺的接头序列也不是容易的事。然而,在 PE 测序模式下如果文库的插入片段比测序读长短,那么 read1 和 read2 中非接头序列的那部分会完全反向互补,Trimmomatic 有一个 ‘palindrome’ 模式会利用这个特点进行接头序列的去除。
下图中 A、B、C、D 四种情况就是 Trimmomatic 去除接头和引物的四种模式:
红色条形:被切除的序列
绿色条形:保留下来的有效读长
深蓝色条形:接头序列
浅蓝色条形:引物序列
A 模式:测序 reads 从起始位置开始就包含了完整的接头序列,那么根据 Illumina 测序原理,这整条 reads 都不可能包含有用序列了,整条 reads 被丢弃。
B 模式:这种相对常见,由于文库插入片段比测序读长短,会在 reads 末端包含部分接头序列,若是这部分接头序列足够长是可以识别并去除的,但如果接头序列太短,比接头匹配参数设置的最短长度还短,那么就无法去除。但是,如果是 PE 测序,可以按照 D 模式去除 reads 末端的很短的接头序列。
C 模式:PE 测序可能出现这种情况,正向测序和反向测序有部分完全反向互补,但是空载的文库,两个接头直接互连,这样的 reads 不包含任何有用序列,正反向测序 reads 都被丢弃。
D 模式:是 Trimmomatic 利用 PE 测序进行短接头序列去除的典范,如果文库插入片段比测序读长短,利用正反向测序 reads 中一段碱基可以完全反向互补的特点,将两个接头序列与 reads 进行比对,同时两条 reads 之间也互相比对,可以将 3’ 末端哪怕只有 1bp 的接头序列都可以被准确去除,相对 B 模式去除接头污染更彻底。
References
http://www.biotrainee.com/thread-1484-1-1.html
http://www.biotrainee.com/thread-856-1-1.html
还有更多文章,请移步公众号阅读
如果你生信基本技能已经入门,需要提高自己,请关注上面的生信技能树,看我们是如何完善生信技能,成为一个生信全栈工程师。
如果你是初学者,请关注下面的生信菜鸟团,了解生信基础名词,概念,扎实的打好基础,争取早日入门。