Cutadapt对测序数据质控(去接头)

简单的说,cutadapt是一款去接头并且也能去低质量碱基的软件。至于接头的产生,一般来说是由于在随机打断时,有些片段的长度过于短,导致插入片段不够,使得测序时会读到测序引物。因此理论上来说接头是不应该被测序仪读到的,但是个别情况读到引物,我们需要通过cutadapt进行去除。

网上多为早些时候的cutadapt教程讲解,对于双端测序的实践较少。因此做个笔记记录下cutadapt对双端测序数据的处理。

  1. 软件下载及安装

    现在这个软件已经更新到1.14版本了,如果下载的是老版本,会导致部分功能无法使用(特别对于双端测序去接头来说),由于cutadapt是python写的,因此下载地址可以是:https://pypi.python.org/pypi/cutadapt#downloads,我选择下载的是cutadapt-1.14.tar.gz

    tar zxvf cutadapt-1.14.tar.gz
    cd cutadapt-1.14
    python setup.py install --user
    

    然后开始看cutadapt的说明书吧,蛮详细的https://cutadapt.readthedocs.io/en/stable/index.html

  2. 对双端数据进行质控

    一般现在测序以双端测序为止,那么双端测序的接头是如何产生的呢?可以参考下面这篇博文http://blog.sina.com.cn/s/blog_14d1975e90102x5o1.html

    先看下cutadapt是如何去接头的

    cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq
    

    这里我们主要看到-a 接头,-A 接头的反向互补序列,并且-a和-A均是去除3’端的接头序列,那么以一个接头序列(NEBNext® UltraTM RNA Library Prep Kit RNA)为例:

    5′ Adapter(5’端接头)
    5’AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT

    3′ Adapter(3’端接头)
    5’GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCT
    GCTTG

    按照我的理解,-a去除的是read1的3’的接头序列,因此-a参数后接3′ Adapter(3’端接头);-A 去除的是read2的3’的接头序列,那么-A参数后应该接5′ Adapter(5’端接头) 的反向互补序列,至于为什么是这样,可以简单的理解为由于插入片段过短,测序仪不仅测了插入片段,还测到了接头序列,侧穿了(不知道对不对,查了一些别人的解释,测序原理还不是非常清楚)

    所以命令如下:

    cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG \
    -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
    -o out.1.fq -p out.2.fq read.1.fq read.2.fq
    

    下面几个参数跟去接头比较相关

    -e 最大错配比例,比如cutadapt在某条序列上检测的接头有15bp长,那么允许这个匹配上的15bp接头中有15*0.1约为1个碱基的错配

    -m –minimum-length 切除接头后的序列长度的最小值

    -O –overlap 默认必须至少有3个碱基匹配时才会认为是adapter序列,但有时可以适当的调大

    –discard-trimmed 去除掉有检测到接头的序列(默认cutadapt只是截掉接头序列以及接头序列以后的序列)

    –untrimmed-output 将没有接头的序列输出到目标文件中(但是必须要跟-o 一起用)

    –untrimmed-paired-output 将没有接头的paired序列输出到目标文件中(也要跟-p 一起用)

    –pair-filter=(any|both) 这个参数很好用,对于双端测序而言,read1和read2都有可能检测到接头。如果选择any,则只要两个中其中一个检测到接头,read1和read2均舍弃;如果选择both,则必须两个都检测到接头,read1和read2才舍弃

    因此如果我们只想要输出没有接头的序列并且paired也同样舍弃,总结上述参数,那么可以如下:

    cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG \
    -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
    -e 0.1 -O 5 -m 50 --discard-trimmed --paired-filter=any
    -o out.1.fq -p out.2.fq read.1.fq read.2.fq
    

    cutadapt除了去接头序列,也能去除低质量碱基。其具体原理可看https://cutadapt.readthedocs.io/en/stable/guide.html#quality-trimming-algorithm,使用-p 20参数即可,默认是对3’端进行去除低质量的碱基;如果想对5’端进行处理的话,使用-p 20,20 也就是用逗号将3’端和5’端的阈值分开表示

  3. 总体上来说,Cutadapt也是一款不错的去接头软件,也有人将其与Trimmomatic的去接头效果进行比较,差别并不大;但是个人觉得Trimmomatic去低质量的算法更加给力,去的更加彻底。单从去接头的来说,确实也差不多,都可以用。Trimmomatic只能去除Illumina平台的接头,而cutadapt则可以去除任意指定的接头

参考文章:
http://www.bio-info-trainee.com/1920.html
http://blog.sina.com.cn/s/blog_14d1975e90102x5o1.html
https://cutadapt.readthedocs.io/en/stable/guide.html#