过滤软件的比较与选择:cutadapt/fastp/trimmomatic
注:还没有完全搞明白,先总结一下特点和使用,之后再慢慢体会、总结经验
本次只针对双端PE
算法都没好好读,因为看不懂==
首先,我们对数据进行过滤,是为了:
去掉接头
去掉低质量reads
去掉污染序列
在尽量去掉上述序列的同时,保留尽可能多的有用数据,减少损失
CutAdapt,2010
- 基于Python,作者是个德国人,长得还挺帅气(✧◡✧) 不过都9年过去了,嗯。。
- 不仅支持illumina,还支持SOLID,454等平台产出的数据
- 支持输入.gz
- 需要自己先检测接头类型(fastqc等),然后搜索接头序列是啥,手动输入到参数里。但是有一个参数 -n,若是两种接头,也可以指定然后去除:-n 2
- 一般命令:
cutadapt -a -A #a是read1的3'接头,A是read2的3'接头(5'接头的反向互补序列)
-e 0.1 -0.5 -m 50 #去除接头后read长度大于50才保留
-o -p #生成文件:过滤后的R1 R2
read1.fastq read2.fastq #输入文件
- 本次分析没用,所以详细参数可以阅读--help
fastp,2018
- 基于c++这种强大的语言所以算法比较高效,中科院深圳所发的。还没用过,不过身边做RNA-Seq的俩师兄强烈推荐,有空可以test一下。
- 主题就是ultra-fast,all-in-one,而且是只处理FASTQ也就是illuminate下机数据
- 特点:
能进行质控,生成比fastqc美观、全面的报告,但是我看了一遍,不如fastqc直观、fresh-friendly
号称去除低质量序列的方法类似于trimmomatic但是更快
自动识别序列并去除
支持illuminate short read,也一定程度支持Pacbio/Nanopore long reads,具体支持到什么程度,需要试验。
参数众多,但是挺有条理的,而且很多都是默认不是必需参数,不会“新手退散”
- 最简单的命令:
fastp -i r1.fq -o rr1.fq -i r2.fq -o rr2.fq
- 这篇介绍写的不错:知乎
但他说一般下机数据要经过fastqc+cutadapt+trimmomatic,有点不太理解,要这么麻烦吗?
Trimmomatic,2014
也是很好用的,引用量超高,good at去除低质量reads,只针对illuminate数据
最重要的特点:对数据的处理步骤与参数的顺序有关!
所以建议先去接头,否则接头被剪更无法有效去除。PE数据常用参数:
ILLMINACLIP: 注意以下参数以:隔开
<fastaWithAdaptersEtc>: 指定包含接头和引物序列(所有被视为污染的序列)的 fasta 文件
<seed mismatches>: 第一步seed搜索时允许的mismatch个数,一般2。
<palindrome clip threshold>: 指定针对 PE的palindrome clip模式下,需要R1和 R2之间至少多少比对分值,才会进行接头切除,例如30。
<simple clip threshold>: 指定切除接头序列的最低比对分值,一般7-15之间。
<minAdapterLength>: 只对 PE 测序的 palindrome clip 模式有效,指定 palindrome 模式下可以切除的接头序列最短长度,默认值是 8。但实际上 palindrome 模式可以切除短至 1bp 的接头污染,所以可以设置为 1。
<keepBothReads> 重要参数!第一次做的时候没加这个参数,结果20%+的数据Unpaired,扔掉不现实,比对处理太麻烦!正确用法:只对 PE 测序的 palindrome clip 模式有效,R1 和 R2 在去除了接头序列之后剩余的部分是完全反向互补的,默认参数 false,意味着整条去除与 R1 完全反向互补的 R2,当做重复去除掉,但在有些情况下,例如需要用到 paired reads 的 bowtie2 流程,就要将这个参数改为 true,否则会损失一部分 paired reads。本次所用命令:(也是公司报告中所用的)
java -jar trimmomatic-0.38.jar PE -threads 2 #双端模式,两个线程
ILLUMINACLIP: #顾名思义,去illumina接头
TruSeq3-PE.fa: #接头文件,需要指定全路径
2:30:10 # 默认格式为 2:30:10:8:false,可改做:2:30:10:8:true
LEADING:20 #从reads的起始端开始切除质量值低于设定的阈值的碱基,直到有一个碱基其质量值达到阈值。一般用LEADING:3???
TRAILING:20 #一般用3,因为Illumina 平台有些低质量的碱基质量值被标记为2,所以设置为 3 可以过滤掉这部分低质量碱基
SLIDINGWINDOW:4:20 #滑窗剪切,统计滑窗口中所有碱基的平均质量值,如果低于设定的阈值,则切掉窗口。此处设置4bp窗口,阈值20,一般阈值用15。
MINLEN:50 #可被保留的最短 read 长度
trimmomatic PE模式默认处理2个文件,也就是说,sh脚本中使用本办法只能一次列举R1 R2两个文件,不能 In_R1 In_R2 IP_R1 IP_R2这样四个文件都列出来,事实证明会报错,trimmomatic有点傻傻的不知道第三个开始的文件该干嘛。
所以要批量做,需要写循环,或者是认真阅读使用说明的参数。
- trimmomatic的更多解读可以参考这个,写得很详细。目前我理解的是以上。
最后附一个图:
出自:Chen et al. Source Code for Biology and Medicine 2014, 9:8. Software for pre-processing Illumina nextgeneration sequencing short read sequences.