Trimmomatic对raw reads的过滤

Trimmomatic是一个针对 Illumina高通量测序的reads trim的工具。既能够针对paired-end 也能弄single ended。它能够利用FASTQ文件(phred + 33 或者是phred + 64 碱基质量格式,取决于Illumina测序的机器)。对于single-ended,一个输入文件和一个输出文件,加上参数。对于paired-end数据,两个输入文件,4个输出文件,分别为2个是’paired',2个是’unpaired'(一个为forward的,一个为reverse的)

 一、下载地址

http://www.usadellab.org/cms/?page=trimmomatic (注意里面download那里,选择source和binary)
<strong>说明书地址:</strong>
http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.30.pdf

#下载
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
#解压
unzip Trimmomatic-0.39.zip

二、方法

Trimmomatic用两种策略来去除adapter:

simple trimming是利用每一个adapter序列去跟reads匹配,如果匹配上,就删除read的这部分序列。

Palindrome trimming 是在adapter序列中reading through一个短的片段。在这个方法中, 相应的adapter序列在硅片上的连接被连接到了reads的开始,形成了adapter+read序列,正链和反链比对,如果比对显示read-through,正链剪切掉adapter,反链去掉(由于它没有包含新的数据)

Trimmomatic它包括如下功能:

 ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the read.
 SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average quality within the window falls below a threshold.
 LEADING: Cut bases off the start of a read, if below a threshold quality
 TRAILING: Cut bases off the end of a read, if below a threshold quality
 CROP: Cut the read to a specified length
 HEADCROP: Cut the specified number of bases from the start of the read
 MINLEN: Drop the read if it is below a specified length
 TOPHRED33: Convert quality scores to Phred-33将碱基质量转换为pred33格式
 TOPHRED64: Convert quality scores to Phred-64 将碱基质量转换为pred64格式
 AVGQUAL:丢掉质量低于这个值的reads

功能详解:

 ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>
 fastaWithAdaptersEtc:包含adapter的文件
 seedMismatches: 许的最大mismatch 数
 palindromeClipThreshold :palindrome模式下匹配碱基数阈值
 simpleClipThreshold: simple模式下的匹配碱基数阈值
 序列的命名暗示了他们将如何被利用,对于'Palindrome'剪切来说,序列的名字应该以'Prefix'开始,以'/1'结尾是针对正链adapter,以'/2'结尾是针对反链adapter。
其他的序列用'simple' 模式,不是以'/1' or '/2'结尾的序列将会跟正链和反链都比对。如果你要针对某一个read序列,你需要对你的adapter命名上能对应。
 阈值用来界定相似度的范围,每一个匹配上的碱基得分是大于0.6,每一个不匹配减少比对的分数Q/10。
因此,一个12个碱基比对最好的得分应该是大于7,25个碱基需要15分。所以我们推荐7-15分这个参数。
对于palindromic匹配来说,更长的匹配是有可能的,因此这个阈值需要高一点,一般设为30分。'seed mismatch'用来形容最大容许的匹配不上的个数,一般是1或者2
 SLIDINGWINDOW::
 windowSize:设定几个碱基为一个单位
 requiredQuality: 该单位最低值
 LEADING:
 quality: Specifies the minimum quality required to keep a base.
 TRAILING:
 quality: Specifies the minimum quality required to keep a base.
 CROP:
 length: The number of bases to keep, from the start of the read.
 HEADCROP:
 length: The number of bases to remove from the start of the read.
 MINLEN:
 length: Specifies the minimum length of reads to be kept.

三、使用例子

java -jar /sam/trimmomatic/trimmomatic-0.32.jar PE -threads 12 -phred64 -trimlog logfile input_forward.fq input_reverse.fq output_forward_paired.fq output_forward_unpaired.fq output_reverse_paired.fq output_reverse_unpaired.fq ILLUMINACLIP:1.adapter.list:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 
同一个指令下的不同参数可以用':'来界定;
phred33 设置碱基的质量格式,如果不设置,默认的是-phred64 
trimlog file就是产生日志,包括如下部分内容:read的名字,留下来的序列的长度,第一个碱基的起始位置,从开始trimmed的长度,最后的一个碱基位于初始read的位置,最后trimmed的数量,其他步骤产生的日志
LEADING:3 切除首端碱基质量小于3的碱基或者N 
TRAILING:3 切除末端碱基质量小于3的碱基或者 
ILLUMINACLIP: 1.adapter.lis:2:30:10 1.adapter.list为adapter文件,允许的最大mismatch 数,palindrome模式下匹配碱基数阈值:simple模式下的匹配碱基数阈值
SLIDINGWINDOW:4:15 Windows的size是4个碱基,其平均碱基质量小于15,则切除 
MINLEN:36 最低reads长度为36 CROP: 保留的reads长度 
HEADCROP: 在reads的首端切除指定的长度

四、讨论

1.这么知道我的reads是-phred33,还是-phred64?

Fastq格式–phred33/64 我的是64。

2.如果我是PE测序,我有两分别针对的adapter,我是否需要将这两个adaper合并为一个?

需要。

3.是否可以批量处理?

ls * |grep "fasta"|perl -lne's/_R1.fasta//g;s/_R2.fasta//g;print' | sort |uniq | perl -lne 'print"java -jar /sam/trimmomatic/trimmomatic-0.32.jar PE -threads 12 -phred33-trimlog ${_}_logfile ${_}_R1.fasta ${_}_R2.fasta ${_}_output_forward_paired${_}_output_forward_unpaired ${_}_output_reverse_paired${_}_output_reverse_unpaired LEADING:25 TRAILING:25 SLIDINGWINDOW:4:25MINLEN:150 "'>work.sh;
 ​nohupsh work.sh & ;
 nohup sh work.sh & ; 在终端机后台工作,即使再次登录,命令还在执行
 nohup sh work.sh ; 在终端机前台中工作

五、参考文献

药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn