一、 基本介绍

Trim Galore是对FastQC和Cutadapt的包装。适用于所有高通量测序，包括RRBS (Reduced Representation Bisulfite-Seq )、Illumina、Nextera 和smallRNA测序平台的双端和单端数据。主要功能包括两步：第一步首先去除低质量碱基，然后去除3' 末端的adapter，如果没有指定具体的adapter，程序会自动检测前1 million的序列，然后对比前12-13bp的序列是否符合以下类型的adapter：Illumina: AGATCGGAAGAGC；Small RNA: TGGAATTCTCGG; Nextera: CTGTCTCTTATA。

Trimmomatic是一个主要针对Illumina测序的reads修剪和过滤的软件。它支持多线程，处理数据速度快，过滤模式有两种，分别对应paired end (PE) reads和single end (SE) reads 测序数据，也支持 gzip 和 bzip2 压缩文件。Trimmomatic的好处在于，它不但可以用来切除illumina测序平台的接头序列，还可以去除由我们自己指定的特定接头序列，而且同时也能够过滤read末尾的低质量序列。具体的原理就是通过滑动一定长度的窗口，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除，注意不是挖掉read中的这部分低质量序列，而是像切菜一样，直接从低质量区域开始把这条read后面的所有其它碱基全部剁掉！否则就是在人为改变实际的基因组序列情况。

Trimmomatic接头序列的选择：一般来说，目前的HiSeq系列和MiSeq系列用的都是TruSeq3，TruSeq2是以前GA2系列的测序仪所用的，已经很少见了。这些信息都可以在illumina的官网上找到，至于具体该用PE（Pair End）还是SE（Single End）就按照具体的测序类型进行选择即可。如果用的不是illumina测序平台，那么我们也可以按照adapters文件夹下的这些文件的格式做一个新的接头序列，然后再作为参数传入。不过在自定义接头序列的时候，命名时有一些小的细节需要注意，可以参考Trimmomatic的主页文档（The Adapter Fasta）。需要明确指明质量值体系是Phred33还是Phred64，默认是Phred64，这需要特别注意，因为我们现在的测序数据基本都是Phred33的了，所以一定要指定这个参数。

ls /home/analysis/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic/adapters/
# NexteraPE-PE.fa  TruSeq2-PE.fa  TruSeq2-SE.fa  TruSeq3-PE-2.fa  TruSeq3-PE.fa  TruSeq3-SE.fa

# TruSeq3-PE.fa的内容是：
>PrefixPE/1
TACACTCTTTCCCTACACGACGCTCTTCCGATCT
>PrefixPE/2
GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
# TruSeq3-SE.fa的内容是：
>TruSeq3_IndexedAdapter
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC
>TruSeq3_UniversalAdapter
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
# TruSeq3-PE-2.fa的内容是：
>PrefixPE/1
TACACTCTTTCCCTACACGACGCTCTTCCGATCT
>PrefixPE/2
GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
>PE1
TACACTCTTTCCCTACACGACGCTCTTCCGATCT
>PE1_rc
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
>PE2
GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
>PE2_rc
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

二、 trim_galore使用方法

(1) 用法和参数

trim_galore  [options]  <filename(s)>
trim_galore  -j 4  -q 25  -e 0.1  --stringency 4  --length 20  --max_n 2  -o outputdir  test.fq.gz  #单端
trim_galore  -j 4  -q 25  -e 0.1  --paired  --stringency 4  --length 40  --max_n 4  -o outputdir  test_1.fq.gz  test_2.fq.gz  #双端

-j/--cores：线程数，默认1（最好设置为4，假设已使用Python 3并安装了Pigz，那么内核设置为4，实际使用内核是15）。
-q/--quality <INT>：除了移除接头之外，还要修剪掉reads的低质量末端。Default Phred score: 20。
-e <ERROR RATE>：最大允许错误率（错误数除以匹配区域的长度）。Default: 0.1
--paired：双端测序，如果其中一个不合格则两个都被剔除。
--stringency <INT>：严格性，设定可以忍受的前后adapter重叠的碱基数，默认1（even a single bp of overlapping sequence will be trimmed off from the 3' end of any read），为非常严格的设置，通常可以适当放宽，因为后一个adapter几乎不可能被测序仪测到。
--length <INT>：丢弃由于质量或接头调整而变得短于长度n的read。值“0”禁用了该行为。Default: 20 bp。在PE150下，可以设置50或36。
--max_n：一条read在被完全删除之前可能包含的N的总数。
-o/--output_dir <DIR>：输入目录，需要提前建立目录，否则运行会报错
--phred33/--phred64：默认phred33。
-a/--adapter <STRING>：要修剪的接头序列。如果没有明确指定，Trim Galore将尝试自动检测是否使用了Illumina universal、Nextera transposase或Illumina small RNA接头序列。
--max_length <INT>：丢弃修剪后长于n的read。仅建议用于小RNA测序，以去除非小RNA序列。
--gzip：用GZIP压缩输出文件。如果输入文件是GZIP压缩的，输出文件也会自动被GZIP压缩。
--dont_gzip：输出文件不会用GZIP压缩。该选项覆盖--gzip。
--fastqc：剪切结束后用默认选项对结果文件进行fastqc分析

(2) 举个例子

trim_galore  -j 4  -q 20  -e 0.1  --stringency 4  --length 40  --max_n 4  -o ./02_cleandata  00_rawdata/WTF5.Input_1.fastq.gz  #单端
trim_galore  -j 4  -q 20  -e 0.1  --paired  --stringency 4  --length 40  --max_n 4  -o ./02_cleandata  00_rawdata/WTF5.Input_1.fastq.gz 00_rawdata/WTF5.Input_2.fastq.gz  #双端

(3) 批量处理

#以双端为例：trim_galore.sh
for i in `ls 00_rawdata/*_1.fastq.gz`
do

    i=`basename $i`
    i=${i%_1.fastq.gz}
    echo $i >>log/trim_galore_log.txt

    trim_galore  -j 4  -q 20  -e 0.1  --paired  --stringency 4  --length 40  --max_n 4  -o ./02_cleandata  00_rawdata/$i\_1.fastq.gz  00_rawdata/$i\_2.fastq.gz  2>> log/trim_galore_log.txt

done

# sh trim_galore.sh &

三、 Trimmomatic使用方法

(1) 用法和参数

trimmomatic PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...  #双端
trimmomatic SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...  #单端

trimmomatic  PE  -phred33  test_1.fq test_2.fq \ test_1_paired.fq test_1_unpaired.fq \ test_2_paired.fq test_2_unpaired.fq \ ILLUMINACLIP:/usr/local/src/Trimmomatic/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 \ LEADING:5  TRAILING:5  SLIDINGWINDOW:3:20  MINLEN:50  #双端
trimmomatic  SE  -phred33  test.fq test.out.fq  ILLUMINACLIP:/home/analysis/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic/adapters/TruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:5:20  LEADING:5  TRAILING:5  MINLEN:50  #单端

<inputFile1> <inputFile2>：为输入文件，fastq或fq.gz。
<outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>：为对应的四个输出文件，双端序列都保留的paired序列文件和只保留一段序列的unpaired序列文件，unpaired是不需要的。单端模式不需要指定文件来存放被过滤掉的read信息。
-threads：线程数。
-phred33：设置碱基的质量格式（默认是64，但现在常用33）。
LEADING:5：切除首端碱基质量小于5的碱基或者N。
TRAILING:5：切除末端碱基质量小于5的碱基或者N。
SLIDINGWINDOW:3:20：从5'端开始进行滑动，滑动窗口Windows的长度是3个碱基，其平均碱基质量小于20则切除。
MINLEN:50：最低reads长度为50。
ILLUMINACLIP:adapters/TruSeq3-PE.fa:2:30:10：切除adapter序列。参数ILLUMINACLIP后面分别接: adapter序列的fasta文件: 比对时接头序列时所允许的最大错误mismatch 数: palindrome模式下匹配碱基数阈值（PE的两条read同时和PE的adapter序列比对，匹配度加起来超30%，就认为这对PE的read含有adapter）: simple模式下的匹配碱基数阈值（只要这条read的某部分和adapter序列有超过10%的匹配率，就代表含有adapter了）。推荐30%和10%。

(2) 举个例子

trimmomatic  PE  -phred33  -threads 4  00_rawdata/WTF4.Input_1.fastq.gz  00_rawdata/WTF4.Input_1.fastq.gz  02_cleandata/WTF4.Input_1.clean.fastq.gz  02_cleandata/WTF4.Input_1.unpaired.fastq.gz  02_cleandata/WTF4.Input_2.clean.fastq.gz  02_cleandata/WTF4.Input_2.unpaired.fastq.gz  ILLUMINACLIP:/home/analysis/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic/adapters/TruSeq3-PE.fa:2:30:10  SLIDINGWINDOW:3:20  LEADING:5  TRAILING:5  MINLEN:50  #双端

trimmomatic  SE  -phred33  -threads 4  00_rawdata/WTF4.Input_1.fastq.gz  02_cleandata/WTF4.Input_1.clean.fastq.gz  ILLUMINACLIP:/home/analysis/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic/adapters/TruSeq3-SE.fa:2:30:10  SLIDINGWINDOW:3:20  LEADING:5  TRAILING:5  MINLEN:50  #单端

(3) 批量处理

#以双端为例：trimmomatic.sh
for i in `ls 00_rawdata/*_1.fastq.gz`
do

    i=`basename $i`
    i=${i%_1.fastq.gz}
    echo $i >>log/trimmomatic_log.txt

    trimmomatic  PE  -phred33  -threads 6  00_rawdata/$i\_1.fastq.gz  00_rawdata/$i\_2.fastq.gz  02_cleandata/$i\_1.clean.fq.gz  02_cleandata/$i\_1.unpaired.fq.gz  02_cleandata/$i\_2.clean.fq.gz 02_cleandata/$i\_2.unpaired.fq.gz  ILLUMINACLIP:/home/analysis/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic/adapters/TruSeq3-PE.fa:2:30:10  SLIDINGWINDOW:3:20  LEADING:5  TRAILING:5  MINLEN:50  2>>log/trimmomatic_log.txt

done

# sh trimmomatic.sh &

四、 Cutadapt

(1) color space数据

例如，测序仪AB 5500xl Genetic Analyzer得到的reads：

@SRR1449908.1 ugc_450_solid0105_20110426_FRAG_BC_ugc_450_3_48_545/1
T20103321021332323020223030003323033200023300033032
+
!05')<:-2(*4-(**(&(%%.)%,%(&((.4%(.%&%(*%2'3)%7,)%&

color space数据fastq里面的第二行不是序列，而是color的编码。Colors may be encoded either as numbers (0=blue, 1=green, 2=orange, 3=red) or as characters A/C/G/T (A=blue, C=green, G=orange, T=red).

Here, T is the primer base.（“T”是引物碱基）bowtie detects and handles primer bases properly (i.e., the primer base and the adjacent color are both trimmed away prior to alignment) as long as the rest of the read is encoded as numbers.

(2) cutadapt使用方法

SOLiD colorspace format不支持使用trim_galore，可以用cutadapt进行数据修剪。

cutadapt  -a ADAPTER  [options]  [-o output.fastq]  input.fastq  #单端
cutadapt  -a ADAPT1  -A ADAPT2  [options]  -o out1.fastq  -p out2.fastq  in1.fastq  in2.fastq  #双端

cutadapt  -j 4  -a AACCGGTT  -o output.fastq  input.fastq  #去除3'接头
cutadapt  -a TGAGACACGCA  -a AGGCACACAGGG  -n 2  -o output.fastq  input.fastq  #多个3'接头
cutadapt  -q 10  -o output.fastq  input.fastq  #过滤3'端的低质量序列
cutadapt  -q 15,10  -o output.fastq  input.fastq  #过滤3'和5'端的低质量序列

输入或输出文件格式：fasta文件指定后缀：.fasta，.fa，.fna；fastq文件指定后缀：.fastq，.fq；或以上任意后缀+压缩格式（.gz，.bz2，.xz），cutadapt会自动解压或压缩。
-j：选择几个核
-a：3’接头，如果机器读长大于fragment 的长度，在测序得到的序列的3'末端会出现部分adapter序列；3'末端有完整的adapter序列，在adapter序列之后还有其他随机测到的序列。对上述两种情况，cutadapt 会把adapter 序列都清除干净。
-n：检测几次，有两种接头就检测两次
-o：输出reads.1.fastq去掉接头的结果
-p：是--paired-output 的缩写，输出reads.2.fastq去掉接头的结果
-q/--trim-qualities：过滤低质量序列，默认使用phred quality+33 的方式识别序列质量。默认只过滤3'端的低质量序列，如果想要过滤5'端低质量序列，需要用逗号隔开。
-m/--minimum-length：去除序列长度低于N的read
--max-n：单条reads中N碱基的最大数量
-A：3’端接头的反向互补序列。The -A/-G/-B/-U options work like their -a/-b/-g/-u counterparts, but are applied to the second read in each pair.
-g：5'接头
-G：5’端接头的反向互补序列
-u/--cut：在序列的3端和5端切除固定长度的序列，当值为正数时，从序列开头切除， 当值为负数时，从序列末端开始切除。
-e：最大允许错误率（错误数除以匹配区域的长度，默认值0.1）
-O：最小重叠长度。如果read和adapter之间的重叠小于该长度，则不修改read。这减少了纯粹由于短随机接头匹配而修剪的碱基数量。（默认值3）
--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.

colorspace选项：https://cutadapt.readthedocs.io/en/v1.18/colorspace.html
-c, --colorspace：Enable colorspace mode: Also trim the color that is adjacent to the found adapter
--maq/--bwa：MAQ- and BWA-compatible colorspace output. This enables -c, -d, -t, --strip-f3 and -y '/1'. Enables colorspace mode (-c), double-encoding (-d), primer trimming (-t)
--no-zero-cap：不要将色彩空间数据中的负质量值更改为零。色彩空间读数的质量值有时是负的，BWA和Bowtie无法处理，默认情况下Cutadapt将色彩空间数据中的负质量值转换为0。
--format=sra-fastq：当使用来自sra-toolkit包的FASTQ-dump程序将这些.sra文件转换为FASTQ格式时，colorspace读取将在每次读取开始时获得额外的质量值，使质量值数目比序列数多1。为了让cutadapt忽略额外的质量基础，在命令行中添加--format=sra-fastq。

fastqc检测到的接头：
https://github.com/csf-ngs/fastqc/blob/master/Contaminants/contaminant_list.txt

(3) 举个例子

cutadapt  --bwa  -a 330201030313112312  -q 20  -m 20  --max-n 3  --format=sra-fastq  SRR1449908.fastq.gz  -o SRR1449908_clean.fastq.gz  >& log/cutadapt_log.txt  &
cutadapt  -c  -t  -a 330201030313112312  -q 20  -m 20  --max-n 3  --format=sra-fastq  SRR1449909.fastq.gz  -o SRR1449909_new.fastq.gz  >& log/cutadapt_log.txt  &

旧版本的bowtie，支持colorspace，参数稍有不同（需要加-C/--color，索引也是需要单独构建的）。

./shellscript/bowtie-0.12.3/bowtie-build  -C  -f  genome.fa  index/bowtie_c  &
./shellscript/bowtie-0.12.3/bowtie  -p 4  -C  index/bowtie_c  SRR1449908_clean.fastq  -S  SRR1449908.sam  2>>log/bowtie_log.txt  &

五、 fastx_toolkit

(1) 去接头

fastx_clipper  [-h] [-a ADAPTER] [-D] [-l N] [-n] [-d N] [-c] [-C] [-o] [-v] [-z] [-i INFILE] [-o OUTFILE] 
fastx_clipper  -i SRR620204.1.fastq  -o SRR620204.clipper.fastq  -Q 33 

[-a ADAPTER]：接头序列（默认为CCTTAAGG）
[-l N]：忽略那些碱基数目少于N的reads，默认为5
[-M N]：要求最小能匹配到接头的长度N，如果和接头匹配的长度小于N不修剪
[-i INFILE]：输入文件，不支持压缩格式的输入文件，不允许序列中存在N碱基，这样的序列会自动去除
[-o OUTFILE]：输出文件
[-Q]：illumina测序需要加上-Q 33，默认情况下认为fastq文件的碱基编码格式为phred64
[-z]：压缩输出
[-d N]：保留接头序列后的N个碱基，默认-d 0
[-c]：放弃那些没有接头的序列
[-C]：只保留没有接头的序列
[-k]：报告只有接头的序列
[-n]：保留有N多序列，默认不保留
[-v]：详细报告序列编号
[-D]：输出调试结果

(2) 去除低质量reads

fastq_quality_filter  [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE] 
fastq_quality_filter  -q 20  -p 80  -Q 33  -i SRR306394.fastq  -o SRR306394_filtered.fastq 
#要求一条reads中最少有80%以上的碱基质量值达到Q20才作保留

[-q N]：最小的需要留下的质量值
[-p N]：每个reads中最少有百分之多少的碱基需要有-q的质量值
[-Q]：illumina测序需要加上-Q 33
[-z]：压缩输出
[-v]：详细报告序列编号，如果使用了-o则报告会直接在STDOUT，如果没有则输入到STDERR

(3) FASTQ-to-FASTA

fastq_to_fasta  [-h] [-r] [-n] [-v] [-z] [-i INFILE] [-o OUTFILE] 

[-h]：This helpful help screen.
[-r]：Rename sequence identifiers to numbers.
[-n]：Keep sequences with unknown (N) nucleotides. Default is to discard such sequences.
[-v]：Verbose - report number of sequences.
If [-o] is specified, report will be printed to STDOUT.
If [-o] is not specified (and output goes to STDOUT), report will be printed to STDERR.
[-z]：Compress output with GZIP.
[-i INFILE]：FASTA/Q input file. default is STDIN.
[-o OUTFILE]：FASTA output file. default is STDOUT.