转录组分析中Trimmomatic和Trim_galore详细使用说明（1） Trimmomatic

最新推荐文章于 2024-10-01 10:38:58 发布

lazymark2

最新推荐文章于 2024-10-01 10:38:58 发布

阅读量8.7k

点赞数 9

分类专栏：生物信息学文章标签：生物信息学 perl 生物学

本文链接：https://blog.csdn.net/lazymark2/article/details/119224606

版权

生物信息学专栏收录该内容

8 篇文章 8 订阅

订阅专栏

优劣势分析

Trimmomatic

优势：
1、可使用参数更多，如滑窗剪切，可以直接选择使用内置的接头序列等等
2、默认可生成paired和unpaired两种文件，更利于下游分析
劣势：
1、代码非常长，而且容易写错，最好写在一个脚本里；
2、参数比较难记，像ILLUMINACLIP中的几个数字分别代表什么必须要对照说明书才能看懂
3、运行时间较长，建议在测序质量较好并且接头残余少的时候使用trim_galore

Trim_galore

优势：
1、安装和使用都非常简单；
2、代码较短
3、参数更直观，不用去死记硬背
4、默认下不区分paired和unpaired，运行速度较快
劣势：
1、可调参数较少

Trimmomatic使用说明

Trimmomatic文件相关参数

Trimmomatic主要就是针对illumina平台的，其是基于java的，因此在非conda安装的情况下，运行必须加上java -jar
命令格式为：

java -jar PE [-threads <threads] [-phred33 | -phred64] [-trimlog ] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ...

其中-threads 代表线程数；

-phred33 代表采用phred33编码系统，默认为phred64，具体使用哪个要看你测序的机器和测序方法；

-trimlog可以生成一个log文件，包括以下信息
the read name read编号
the surviving sequence length 处理后剩余片段长度
the location of the first surviving base, aka. the amount trimmed from the start 开头切掉的碱基数
the location of the last surviving base in the original read read中最后一个剩余碱基的位置
the amount trimmed from the end 尾端切掉的碱基数

输入文件
input 1 和 input 2
这里有两点要注意
1、input文件前是不需要加连接符 “-”的
也就是说类似

java -jar trimmomatic PE -threads 12 -phred 33 -trimlog sample1_R1.fq.gz sample1_R2.fq.gz

2、可以利用-basein参数

java -jar trimmomatic PE -threads 12 -phred 33 -trimlog -basein sample1_R1.fq.gz

就只需要输入一端的文件，程序会自动推测另一端文件的名称

输出文件
一般来说，需要根据顺序读取/逆序读取/paired/unpaired输出四个文件，这样代码会非常长。为了节约时间，程序同样提供了-baseout参数

java -jar trimmomatic PE -threads 12 -phred 33 -trimlog -basein sample1_R1.fq.gz  -baseout mySampleFiltered.fq.qz

这样程序会自动生成以下四个文件:

mySampleFiltered_1P.fq.gz - `for paired forward reads`
mySampleFiltered_1U.fq.gz - `for unpaired forward reads`
mySampleFiltered_2P.fq.gz - `for paired reverse reads`
mySampleFiltered_2U.fq.gz - `for unpaired reverse reads`

Trimmomatic剪切相关主要参数

举例说明

java -jar trimmomatic-0.30.jar PE s_1_1_sequence.txt.gz s_1_2_sequence.txt.gz
lane1_forward_paired.fq.gz lane1_forward_unpaired.fq.gz lane1_reverse_paired.fq.gz
lane1_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:1:true LEADING:3
TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:1:true
ILLUMINACLIP代表使用的是ILLUMINACLIP剪切；
TruSeq3-PE.fa是对应的adapter序列的fasta文件，是内置在trimmomatic包下的adapters文件夹中。
根据说明书，先是看fastqc报告中的overrepresented sequences中是否有显示对应的adapter序列，然后根据显示的序列选择对应的fasta文件，比如：

Illumina Single End	Illumina Paired End
TruSeq2-SE.fa	TruSeq2-PE.fa

TruSeq Universal Adapter Single End	TruSeq Universal Adapter Paired End
TruSeq3-SE.fa	TruSeq3-PE.fa

2:30:10:1:true分别代表的含义：

<seed mismatches>:<palindrome clipthreshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>

seed mimatches 种子错配数：palindrome clipthreshold 双端剪切阈值：simple clip threshold单端剪切阈值：minAdapterLength 接头序列最小配对碱基数：keepBothReads 是否总是保留互补序列

种子错配数：
在trimmomatic针对双端测序的回文算法中，会进行2步检测；
1、首先是用所谓的种子序列（16bp）在reads中进行滑窗查找，找相似序列
2、找到相似序列后才进行比对评分，评分标准为每个配对碱基0.6分，每个错配减去Q/10，因此质量分数（Q）越低的碱基影响也越小。
因此，一般单端序列阈值定为7-15分，也就是12-25个碱基配对成功，判断该段为接头序列；而双端序列根据回文算法，该数值可以提高到30左右，也就是50个碱基左右。
那么这里的种子错配数就是在进行第1步检测的滑窗查找时最多可以允许16个种子碱基序列中出现2个错配。
双端剪切阈值：
就是上文中提到的一般为30，也就是配对达到50个碱基左右就判断为接头序列并进行剪切；
单端剪切阈值：
一般定义为10，也就是17个碱基左右；
接头序列最小配对碱基数：
简单来说就是最多能允许末端残留多少个接头序列的碱基，默认值为8，极端点可以设为3甚至1；该参数与trim_galore中的-stringency含义相同
是否总是保留互补序列：
如下图所示，红色区域为接头序列，在剪切掉接头序列后，两条互补的剩余序列实际上能提供的碱基信息是一样的，默认会丢弃其中一条序列，而选择总是保留互补序列则会同时保留两条序列，这样的好处就是1、提高了保留序列的比例；2、保持每条序列都是成对的，以便于下游的分析（某些下游软件可能无法识别unpaired 的序列）

实际上互补序列的序列信息是一样的
SLIDINGWINDOW:4:15
SLIDINGWINDOW 滑窗剪切:<windowSize 窗宽>:<requiredQuality 质量分数阈值>

窗宽: 每次滑动所检测的碱基数
质量分数阈值：最小质量分数，出现低于该分数的碱基则该窗宽段的碱基被剪切

LEADING:<quality质量分数>
TRAILING:<quality质量分数>
分别是剪切开头和末尾的低于质量分数阈值的碱基

MINLEN:<length最小长度>
最小保留长度为36，低于该长度的片段舍弃

其他参数：

MAXINFO:最大信息<targetLength 目标长度> <strictness 严格度>
简单来说该参数就是综合衡量了最小读取长度、额外读取长度以及错误敏感性三点的一个剪切方法。
它会根据规定的读取长度和实际读取长度进行比较（逻辑回归），短于规定长度的reads会趋向于被舍弃；根据所谓的宽容度（1-严格度），对reads进行加权；剩余的碱基的质量分数会结合进行运算，计算出该reads不出错的概率（正确性），然后再根据严格度进行加权。

目标长度：规定的读取长度
严格度：该值应该在0-1之间，两个极端分别代表尽量保留reads长度/去除所有错误的碱基。一般低于0.2的值认为是倾向于保留长度，高于0.8的值倾向于正确性。

Trimmomatic 批量处理命令：

假设在./raw下有4个fastq文件，将结果保存到…/cleandata中；fastqc显示adapter content主要是Nextera Transposase Sequence：

sample1_R1.fastq.gz
sample1_R2.fastq.gz
sample2_R1.fastq.gz
sample2_R2.fastq.gz

批量处理

ls *.gz | cut -d '_' -f 1 | sort -u | while read id; do \ 
##样本量较少可以nohup到后台；或者加一个echo sample${id} is running，确保程序完全运行
trimmomatic PE -threads 8 -phred33 -trimlog trim.log ${id}_R1.fastq.gz ${id}_R2.fastq.gz -baseout ../cleandata/${id}.fastq.gz \
ILLUMINACLIP:/Your/path/to/Trimmomatic/adapters/NexteraPE-PE.fa:2:30:10:1:true SLIDINGWINDOW:4:20 LEADING:3 TRAILING:3 MINLEN:36 TOPHRED33 & done