高通量测序数据质控神器Trimmomatic
2019-03-19 本文已影响28人
kongxx
简介
高通量测序下机的原始数据中存在一些低质量数据、接头以及barcode序列等，为消除其对后续分析准确性产生的影响，在数据下机以后对原始数据进行质控处理就成了至关重要的环节。Trimmomatic就是一个高通量测序数据质控神器，可以对测序数据进行过滤。
Trimmomatic 支持多线程，处理数据速度快，主要用来去除 Illumina 平台的 Fastq 序列中的接头，并根据碱基质量值对 Fastq 进行修剪。软件有两种过滤模式，分别对应 SE（单末端测序模式） 和 PE（双末端测序模式） 测序数据，同时支持 gzip 和 bzip2 压缩文件。另外，也支持 phred-33 和 phred-64 格式互相转化。
准备
下载测试数据
$ curl -O ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR571/ERR571271/ERR571271.sra
将 sra 文件转换成 fastq 文件
为了将sra文件转换成fastq格式，我们需要使用 fastq-dump 工具，这个工具被打包在 sratoolkit 工具包中。
sratoolkit 工具包下载地址
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software
Linux环境可以执行选择下面地址直接下载并解压安装
# CentOS
$ wget -c https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-centos_linux64.tar.gz
# Ubuntu
$ wget -c https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-ubuntu64.tar.gz
$ tar zxvf sratoolkit.2.9.2-centos_linux64.tar.gz
将 sra 文件转换成 fastq 文件
./sratoolkit.2.9.2-centos_linux64/bin/fastq-dump --split-files ERR571271.sra
解压后产生 ERR571271_1.fastq 和 ERR571271_2.fastq 文件。
下载 Docker 镜像
为了测试方便，这里没有从源码进行编译，而是直接使用的 Docker 容器来做测试。
为了从 quay.io 获取 spades 镜像，需要先注册一个账号，注册比较简单，就是填个表格就行了。
$ sudo docker login quay.io
$ sudo docker pull quay.io/biocontainers/trimmomatic:0.38--1
运行 Docker 容器
在上面下载的数据文件同级目录下运行下面命令来启动一个容器，且后续的测试命令都需要在此容器里运行。
# 启动一个容器
$ sudo docker run -it --rm -v `pwd`:/trimmomatic quay.io/biocontainers/trimmomatic:0.38--1 bash
# 运行一下trimmomatic命令验证容器可用
$ trimmomatic
运行
$ cd /trimmomatic
$ trimmomatic PE \
-phred33 \
/trimmomatic/ERR571271_1.fastq /trimmomatic/ERR571271_2.fastq \
/trimmomatic/ERR571271_1_paired.fq /trimmomatic/ERR571271_1_unpaired.fq \
/trimmomatic/ERR571271_2_paired.fq /trimmomatic/ERR571271_2_unpaired.fq \
ILLUMINACLIP:/usr/local/share/trimmomatic-0.38-1/adapters/TruSeq3-PE.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 HEADCROP:8 MINLEN:36
TrimmomaticPE: Started with arguments:
-phred33 /trimmomatic/ERR571271_1.fastq /trimmomatic/ERR571271_2.fastq /trimmomatic/ERR571271_1_paired.fq /trimmomatic/ERR571271_1_unpaired.fq /trimmomatic/ERR571271_2_paired.fq /trimmomatic/ERR571271_2_unpaired.fq ILLUMINACLIP:/usr/local/share/trimmomatic-0.38-1/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 HEADCROP:8 MINLEN:36
Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'
ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 1812467 Both Surviving: 1557974 (85.96%) Forward Only Surviving: 248988 (13.74%) Reverse Only Surviving: 3135 (0.17%) Dropped: 2370 (0.13%)
TrimmomaticPE: Completed successfully
参数说明：
- PE/SE - 设定对Paired-End或Single-End的reads进行处理，其输入和输出参数稍有不一样。
- threads - 设置多线程运行数
- phred33 - 设置碱基的质量格式，可选pred64
- ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 - 切除adapter序列。参数后面分别接adapter序列的fasta文件：允许的最大mismatch数：palindrome模式下匹配碱基数阈值：simple模式下的匹配碱基数阈值。
- LEADING:3 - 切除首端碱基质量小于3的碱基
- TRAILING:3 - 切除尾端碱基质量小于3的碱基
- SLIDINGWINDOW:4:15 - 从5'端开始进行滑动，当滑动位点周围一段序列(window)的平均碱基低于阈值，则从该处进行切除。Windows的size是4个碱基，其平均碱基质量小于15，则切除。
- MINLEN:50 - 最小的reads长度
- CROP:<length> - 保留reads到指定的长度
- HEADCROP:<length> - 在reads的首端切除指定的长度
- TOPHRED33 - 将碱基质量转换为pred33格式
- TOPHRED64 - 将碱基质量转换为pred64格式
