一般来讲,我们对测序数据进行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 <marcel.martin@scilifelab.se>
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.
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.
去接头对基因组或者转录组组装是非常重要的: