用cutadapt软件来对双端测序数据去除接头

一般来讲,我们对测序数据进行QC,就三个大的方向:Quality trimming, Adapter removal, Contaminant filtering,当我们是双端测序数据的时候,去除接头时,也会丢掉太短的reads,就容易导致左右两端测序文件reads数量不平衡,有一个比较好的软件能解决这个问题,我比较喜欢的是cutadapt软件的PE模式来去除接头!尤其是做基因组或者转录组de novo 组装的时候,尤其要去掉接头,去的干干净净!
cutadapt是经典的python软件,但是因为我的linux服务器有点问题 ,可能是root权限问题,没有用pip install cutadapt 安装成功,我懒得搞这些了,其实可以自己去下载cutadapt的源码,然后进入源码文件夹里面 python setup.py install –user 到自己的 ~/.local/bin下面。
所以我用conda安装了cutadapt软件,http://www.bio-info-trainee.com/1906.html 所以我需要 python ~/miniconda2/pkgs/cutadapt-1.10-py27_0/bin/cutadapt –help 才能调用这个软件,不过,问题不大,我也就是试用一下。

首先用fastqc软件对测序数据进行简单检测,看看有什么接头:

既然fastqc能探测到你的接头,说明它里面内置了所有的接头序列,在github可以查到: https://github.com/csf-ngs/fastqc/blob/master/Contaminants/contaminant_list.txt 或者:Download common Illumina adapters from    https://github.com/vsbuffalo/scythe/blob/master/illumina_adapters.fa TruSeq Universal Adapter AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT Illumina Small RNA 3p Adapter 1 ATCTCGTATGCCGTCTTCTGCTTG 最严重的一般是TruSeq Universal Adapter , 而且它检测到你其它的接头可能就是这个 TruSeq Universal Adapter 的一部分而已~! 我们可以先用cutadapt去除试一下 cutadapt软件支持对PE 测序数据的处理,基本的用法是: cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq -a和-A是左右端测序数据的3端接头,-g和-G是左右端测序数据的5端接头。 支持fastq和fasta格式的gz压缩文件,必要时用-f参数指定测序文件数据格式即可。 我自己的实际例子如下:

python ~/miniconda2/pkgs/cutadapt-1.10-py27_0/bin/cutadapt -a AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT -a ATCTCGTATGCCGTCTTCTGCTTG -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -A CAAGCAGAAGACGGCATACGAGAT -e 0.1 -O 5 -m 50  -n  2 –pair-filter=both -o read_PE1.fq -p read_PE2.fq test1.fastq test2.fastq >& log.txt

参数解释如下: 两个-a 参数后面接的是两种接头,两个-A参加后面接的是同样的两个接头的反向互补序列! 加上 -n 2 是因为有两个接头,我的测序数据里面有可能一个read里面含有两种接头,所以需要检测两次。 -e 0.1 -O 5 -m 50 是标准参数,自己看readme就好了,其中-m 设置为50 是表示去除接头后如果read长度小于50我就不要了,因为我是PE150测序的,这也就是为什么要用PE模式来去除接头,保证过滤后的reads还是数量继续平衡的。 去除了接头之后的文件再用fastqc软件跑一下,很明显可以看到 接头去除成功啦! 当然,这个软件还有很多其它的功能,我就不一一介绍了: 参考: http://cutadapt.readthedocs.io/en/stable/guide.html 帮助文档还是蛮长的: cutadapt version 1.9.1 Copyright (C) 2010-2015 Marcel Martin cutadapt removes adapter sequences from high-throughput sequencing reads. Usage:     cutadapt -a ADAPTER [options] [-o output.fastq] input.fastq For paired-end reads:     cutadapt -a ADAPT1 -A ADAPT2 [options] -o out1.fastq -p out2.fastq in1.fastq in2.fastq Replace “ADAPTER” with the actual sequence of your 3’ adapter. IUPAC wildcard characters are supported. The reverse complement is *not* automatically searched. All reads from input.fastq will be written to output.fastq with the adapter sequence removed. Adapter matching is error-tolerant. Multiple adapter sequences can be given (use further -a options), but only the best-matching adapter will be removed. Input may also be in FASTA format. Compressed input and output is supported and auto-detected from the file name (.gz, .xz, .bz2). Use the file name ‘-’ for standard input/output. Without the -o option, output is sent to standard output. Citation: Marcel Martin. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.Journal, 17(1):10-12, May 2011. http://dx.doi.org/10.14806/ej.17.1.200 Use “cutadapt –help” to see all command-line options. See http://cutadapt.readthedocs.org/ for full documentation. Options:   –version             show program’s version number and exit   -h, –help            show this help message and exit   –debug               Print debugging information.   -f FORMAT, –format=FORMAT                         Input file format; can be either ‘fasta’, ‘fastq’ or                         ‘sra-fastq’. Ignored when reading csfasta/qual files                         (default: auto-detect from file name extension).   Options that influence how the adapters are found:     Each of the three parameters -a, -b, -g can be used multiple times and     in any combination to search for an entire set of adapters of possibly     different types. Only the best matching adapter is trimmed from each     read (but see the –times option). Instead of giving an adapter     directly, you can also write file:FILE and the adapter sequences will     be read from the given FASTA FILE.     -a ADAPTER, –adapter=ADAPTER                         Sequence of an adapter that was ligated to the 3’ end.                         The adapter itself and anything that follows is                         trimmed. If the adapter sequence ends with the ‘$’                         character, the adapter is anchored to the end of the                         read and only found if it is a suffix of the read.     -g ADAPTER, –front=ADAPTER                         Sequence of an adapter that was ligated to the 5’ end.                         If the adapter sequence starts with the character ‘^’,                         the adapter is ‘anchored’. An anchored adapter must                         appear in its entirety at the 5’ end of the read (it                         is a prefix of the read). A non-anchored adapter may                         appear partially at the 5’ end, or it may occur within                         the read. If it is found within a read, the sequence                         preceding the adapter is also trimmed. In all cases,                         the adapter itself is trimmed.     -b ADAPTER, –anywhere=ADAPTER                         Sequence of an adapter that was ligated to the 5’ or                         3’ end. If the adapter is found within the read or                         overlapping the 3’ end of the read, the behavior is                         the same as for the -a option. If the adapter overlaps                         the 5’ end (beginning of the read), the initial                         portion of the read matching the adapter is trimmed,                         but anything that follows is kept.     -e ERROR_RATE, –error-rate=ERROR_RATE                         Maximum allowed error rate (no. of errors divided by                         the length of the matching region) (default: 0.1)     –no-indels         Do not allow indels in the alignments (allow only                         mismatches). (default: allow both mismatches and                         indels)     -n COUNT, –times=COUNT                         Remove up to COUNT adapters from each read (default:                         1)     -O LENGTH, –overlap=LENGTH                         Minimum overlap length. If the overlap between the                         read and the adapter is shorter than LENGTH, the read                         is not modified. This reduces the no. of bases trimmed                         purely due to short random adapter matches (default:                         3).     –match-read-wildcards                         Allow IUPAC wildcards in reads (default: False).     -N, –no-match-adapter-wildcards                         Do not interpret IUPAC wildcards in adapters.     –no-trim           Match and redirect reads to output/untrimmed-output as                         usual, but do not remove adapters.     –mask-adapter      Mask adapters with ‘N’ characters instead of trimming                         them.   Additional read modifications:     -u LENGTH, –cut=LENGTH                         Remove LENGTH bases from the beginning or end of each                         read. If LENGTH is positive, bases are removed from                         the beginning of each read. If LENGTH is negative,                         bases are removed from the end of each read. This                         option can be specified twice if the LENGTHs have                         different signs.     -q [5’CUTOFF,]3’CUTOFF, –quality-cutoff=[5’CUTOFF,]3’CUTOFF                         Trim low-quality bases from 5’ and/or 3’ ends of reads                         before adapter removal. If one value is given, only                         the 3’ end is trimmed. If two comma-separated cutoffs                         are given, the 5’ end is trimmed with the first                         cutoff, the 3’ end with the second. See documentation                         for the algorithm. (default: no trimming)     –quality-base=QUALITY_BASE                         Assume that quality values in FASTQ are encoded as                         ascii(quality + QUALITY_BASE). This needs to be set to                         64 for some old Illumina FASTQ files. Default: 33     –trim-n            Trim N’s on ends of reads.     -x PREFIX, –prefix=PREFIX                         Add this prefix to read names. Use {name} to insert                         the name of the matching adapter.     -y SUFFIX, –suffix=SUFFIX                         Add this suffix to read names; can also include {name}     –strip-suffix=STRIP_SUFFIX                         Remove this suffix from read names if present. Can be                         given multiple times.     –length-tag=TAG    Search for TAG followed by a decimal number in the                         description field of the read. Replace the decimal                         number with the correct length of the trimmed read.                         For example, use –length-tag ‘length=’ to correct                         fields like ‘length=123’.   Options for filtering of processed reads:     –discard-trimmed, –discard                         Discard reads that contain an adapter. Also use -O to                         avoid discarding too many randomly matching reads!     –discard-untrimmed, –trimmed-only                         Discard reads that do not contain the adapter.     -m LENGTH, –minimum-length=LENGTH                         Discard trimmed reads that are shorter than LENGTH.                         Reads that are too short even before adapter removal                         are also discarded. In colorspace, an initial primer                         is not counted (default: 0).     -M LENGTH, –maximum-length=LENGTH                         Discard trimmed reads that are longer than LENGTH.                         Reads that are too long even before adapter removal                         are also discarded. In colorspace, an initial primer                         is not counted (default: no limit).     –max-n=COUNT       Discard reads with too many N bases. If COUNT is an                         integer, it is treated as the absolute number of N                         bases. If it is between 0 and 1, it is treated as the                         proportion of N’s allowed in a read.   Options that influence what gets output to where:     –quiet             Do not print a report at the end.     -o FILE, –output=FILE                         Write modified reads to FILE. FASTQ or FASTA format is                         chosen depending on input. The summary report is sent                         to standard output. Use ‘{name}’ in FILE to                         demultiplex reads into multiple files. (default:                         trimmed reads are written to standard output)     –info-file=FILE    Write information about each read and its adapter                         matches into FILE. See the documentation for the file                         format.     -r FILE, –rest-file=FILE                         When the adapter matches in the middle of a read,                         write the rest (after the adapter) into FILE.     –wildcard-file=FILE                         When the adapter has N bases (wildcards), write                         adapter bases matching wildcard positions to FILE.                         When there are indels in the alignment, this will                         often not be accurate.     –too-short-output=FILE                         Write reads that are too short (according to length                         specified by -m) to FILE. (default: discard reads)     –too-long-output=FILE                         Write reads that are too long (according to length                         specified by -M) to FILE. (default: discard reads)     –untrimmed-output=FILE                         Write reads that do not contain the adapter to FILE.                         (default: output to same file as trimmed reads)   Colorspace options:     -c, –colorspace    Enable colorspace mode: Also trim the color that is                         adjacent to the found adapter.     -d, –double-encode                         Double-encode colors (map 0,1,2,3,4 to A,C,G,T,N).     -t, –trim-primer   Trim primer base and the first color (which is the                         transition to the first nucleotide)     –strip-f3          Strip the _F3 suffix of read names     –maq, –bwa        MAQ- and BWA-compatible colorspace output. This                         enables -c, -d, -t, –strip-f3 and -y ‘/1’.     –no-zero-cap       Do not change negative quality values to zero in                         colorspace data. By default, they are changed to zero                         since many tools have problems with negative                         qualities.     -z, –zero-cap      Change negative quality values to zero. This is                         enabled by default when -c/–colorspace is also                         enabled. Use the above option to disable it.   Paired-end options:     The -A/-G/-B/-U options work like their -a/-b/-g/-u counterparts.     -A ADAPTER          3’ adapter to be removed from second read in a pair.     -G ADAPTER          5’ adapter to be removed from second read in a pair.     -B ADAPTER          5’/3 adapter to be removed from second read in a pair.     -U LENGTH           Remove LENGTH bases from the beginning or end of each                         second read (see –cut).     -p FILE, –paired-output=FILE                         Write second read in a pair to FILE.     –pair-filter=(any|both)                         Which of the reads in a paired-end read have to match                         the filtering criterion in order for it to be                         filtered. Default: any.     –interleaved       Read and write interleaved paired-end reads.     –untrimmed-paired-output=FILE                         Write second read in a pair to this FILE when no                         adapter was found in the first read. Use this option                         together with –untrimmed-output when trimming paired-                         end reads. (Default: output to same file as trimmed                         reads.)     –too-short-paired-output=FILE                         Write second read in a pair to this file if pair is                         too short. Use together with –too-short-output.     –too-long-paired-output=FILE                         Write second read in a pair to this file if pair is                         too long. Use together with –too-long-output.

 

I think assembly is one of the things were you definitely want to remove adapters, if not you’ll have spurious overlaps creating erroneous contigs/transcripts, just because they share an adapter. 去接头对基因组或者转录组组装是非常重要的: http://www.opiniomics.org/we-need-to-stop-making-this-simple-fcking-mistake/ http://grahametherington.blogspot.jp/2014/09/why-you-should-qc-your-reads-and-your.html

声明:本站部分文章及图片源自用户投稿,如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢!

上一篇 2018年7月26日
下一篇 2018年7月26日

相关推荐