生物信息学/trim galore

使用trim_galore进行质量控制 编辑

官网:http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/

参考文档:https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md

trim_galore的安装 编辑

使用bioconda安装 编辑

conda install -y trim-galore

官网下载安装 编辑

cd ~/software
wget https://github.com/FelixKrueger/TrimGalore/archive/0.6.2.tar.gz
tar zxvf 0.6.2.tar.gz -C ~/software/
echo 'PATH=$PATH:~/software/TrimGalore-0.6.2/' >> ~/.bashrc
source ~/.bashrc

trim-galore的用法 编辑

USAGE:

trim_galore [options] <filename(s)>

trim-galore的选项和参数 编辑

-h/--help               打印帮助信息并退出。
-v/--version            打印版本信息并退出。
-q/--quality <INT>      去除低质量的碱基和接头。对RRBS样本,第一轮先进行去除低质量碱基,第二轮去除接头。其他类型样本,这两步是一起进行。算法与BWA一样(提取所有质量值转换为数值;计算从序列所有索引到末端的部分和;切除和最小的索引代表的序列)。默认为20。
--phred33               Cutadapt使用 ASCII+33 质量值分数体系(Sanger/Illumina 1.9+ encoding)作为Phred scores来进行序列修剪。默认:ON
--phred64               Cutadapt使用 ASCII+64 质量值分数体系(Illumina 1.5 encoding)作为Phred scores来进行序列修剪。可以从Fastqc的报告得到确定使用那种质量值体系。
--fastqc                质控完成后运行FastQC。
--fastqc_args "<ARGS>"  传递给FastQC的参数,如果需要传递的参数超过一个,使用 "arg1 arg2 etc."形式传递。一个示例:--fastqc_args "--nogroup --outdir /home/"。传递参数时,会自动调用 FastQC, --fastqc不需要再设置
-a/--adapter <STRING>   输入需要裁剪的adapter序列。如果不指定,Trim Galore会自动寻找可能性最高的平台对应的adapter,默认会从Illumina universal, Nextera transposase或Illuminasmall RNA adapter的三个平台中搜索最适合的,也可以直接使用参数指定这三种平台,即--illumina、--nextera和--small_rna。
如果在指定的第一个文件的前一百万个序列中未检测到adapter,或者在多个adapter序列之间没有关联,则Trim Galore默认为'--illumina'(默认是Illumina adapter, 否则“ --nextera”为默认值)。
单碱基也可以这样赋值:-a A{10}, 即-a AAAAAAAAAA
-a2/--adapter2 <STRING> 针对paired-end文件read 2的接头序列,该选项需要和 '--paired' 一起使用。该选项需要设置好,如果被裁剪的是smallRNA,a2会自动设置成Illumina small RNA 5' adapter(GATCGTCGGACT)。
单碱基也可以这样赋值:-a2 A{10}, 即-a2 AAAAAAAAAAA
--illumina              程序会裁剪read起始13bp是Illumian通用接头'AGATCGGAAGAGC'的序列,代替自动检测接头序列。

--nextera               程序会裁剪read起始12bp是Nextera adapter 'CTGTCTCTTATA'的序列,代替自动检测接头序列。

--small_rna             程序会裁剪read起始12bp是Illumina Small RNA 3' Adapter 'TGGAATTCTCGG'的序列,代替自动检测接头序列。设置--small_RNA接头也会将--length值降低到18bp。 如果smallRNA文库是配对末端,除非已明确定义-a2,否则a2将自动设置为Illumina samallRNA 5' adapter(GATCGTCGGACT)。
--consider_already_trimmed <INT>     在接头自动检测过程中,允许用户设置一个<INT> 阈值,来判断接头是否已被修剪。如果没有adapter序列超过该阈值, 不再执行额外的adapter修剪(从技术角度讲,相当于设置参数'-a X'来修剪adapter). 低质量值修剪会正常执行。默认:NOT SELECTED (运行自动检测接头程序)。                     
--max_length <INT>      高于此<INT> bp的read会被裁剪,仅建对smallRNA测序去除non-small RNA序列设置。
-s/--stringency <INT>      设定可以忍受的前后adapter重叠的碱基数,默认为1(非常苛刻)。即使一个碱基重叠也将修剪掉read的3'端,可以适当放宽,因为后一个adapter几乎不可能被测序仪读到。read.
-e <ERROR RATE>         运行的最大错误率(错配数除以匹配区域的长度),默认为0.1
--gzip                  将输出结果压缩为GZIP格式。
--dont_gzip             不压缩文件,会覆盖--gzip选项。

--length <INT> 			保留read的最小长度,由于质量或接头修剪后序列长度小于设定值会被修剪。设置成'0'表示关闭该参数。默认:20bp。对于双末端read,两条序列必须都大于 <INT> 才符合要求(详见--paired)。如果只有双末端序列中只有一条长度太短,可以保存成单末端的read(详见 --retain_unpaired)。

--max_n COUNT           read包含N的最大数,在双末端中,read中N数目超过此限制的一对read都会被去除。

--trim-n                去除read中的N,不用于RRBS模式。

-o/--output_dir <DIR>   指定输出目录。需要提前建立目录,否则运行会报错。

--no_report_file        不输出报告。
--suppress_warn         不输出 STDOUT  STDERR 信息。
--clip_R1 <int>         从read(双末端的read 1或单末端read) 5'端切除<int> bp碱基,对于5'端质量值差的read,或需要切除5'端的碱基的情况设置此参数。默认: OFF

--clip_R2 <int>         从双末端的read2 5'端切除<int> bp碱基,对于双末端read 2 5'端质量值差的read,或需要切除5'端的碱基的情况设置此参数,只对双末端read生效。。对于双末端BS-Seq, 建议删除5'端前几个bp,因为末端修复反应可能产生低甲基化。请参考Bismark用户指南中的M-bias图部分。默认: OFF

--three_prime_clip_R1 <int>     在去除Adaptor序列和低质量序列后,从read(双末端的read 1或单末端read) 3'端切除<int> bp碱基,对于3'端质量值差的read,或需要切除3'端的碱基的情况(这些序列与接头序列、质控无关)设置此参数。默认: OFF

--three_prime_clip_R2 <int>     在去除Adaptor序列和低质量序列后,从双末端read 2 3'端切除<int> bp碱基,对于3'端质量值差的read,或需要切除3'端的碱基的情况(这些序列与接头序列、质控无关)设置此参数。默认: OFF

--2colour/--nextseq INT 		该选项激活Cutadapt的 '--nextseq-trim=3'CUTOFF',设置一个质量控制的阈值(一般使用 -q 设置), 但忽略G碱基的质量值。该参数适用于NextSeq和NovaSeq平台, 因为在这些平台中,没有任何信号值的见解被称为高质量G碱基。该参数与 '-q INT'的功能不同,注意区分。

--path_to_cutadapt </path/to/cutadapt>     指定Cutadapt软件的安装路径,如/my/home/cutadapt-1.7.1/bin/cutadapt,如果不指定,默认Cutadapt存在与系统环境变量PATH中。

--basename <PREFERRED_NAME>		输出文件的前缀名PREFERRED_NAME,单末端数据的输出名类似PREFERRED_NAME_trimmed.fq(.gz), 双末端数据的输出名类似PREFERRED_NAME_val_1.fq(.gz)  PREFERRED_NAME_val_2.fq(.gz) for paired-end data. --basename只针对1个文件(单末端)或2个文件(双末端)有效,不能用于过长的列表。

-j/--cores INT          用于质控的线程数 [默认: 1] 需要安装Python 3和多线程压缩命令pigz(parallel gzip)。Trim Galore可以从Cutadapt命令的配置信息中获得这些信息(Cutadapt命令所在路径可从环境变量或--path_to_cutadapt参数中获得)。如果检测到Python 2,--cores参数设置为1。如果没有找到pigz命令,Trim Galore使用gzip命令压缩,gzip压缩会减慢多线程的速度(https://github.com/FelixKrueger/TrimGalore/issues/16#issuecomment-458557103 for more info)。
实际使用的线程数:如果安装了Python 3和pigz, --cores 2表示使用2个线程读取输入文件,2个线程输出结果,Cutadapter使用2个线程并占用额外的两个线程,最后Trim Galore使用1个线程,可以达到9个线程,大部分时间并不同时使用9个线程。双末端序列在输出时使用两倍的线程数,--cores 4表示:4 (read) + 4 (write) + 4 (Cutadapt) + 2 (extra Cutadapt) +	1 (Trim Galore) = 15

# 其他参数 - 不使用adapter/质量值进行质控

--hardtrim5 <int>      	切除5'端的指定的 <int> bp碱基序列,一旦完成,Trim Galore会退出。输出文件名类似 .<int>_5prime.fq(.gz)。使用示例:
                        before:         CCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAAT
                        --hardtrim5 20:CCTAAGGAAACAAGTACACT

--hardtrim3 <int>       切除3'端的指定的 <int> bp碱基序列,一旦完成,Trim Galore会退出。输出文件名类似 .<int>_3prime.fq(.gz)。使用示例:
												before:         CCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAAT
												--hardtrim3 20:TTTTTAAGAAAATGGAAAAT

--clock                 该模式下,会使用Mouse Epigenetic Clock (Multi-tissue DNA methylation age predictor in mouse, Stubbs et al.,Genome Biology, 2017 18:68 https://doi.org/10.1186/s13059-017-1203-5)方法来修剪reads。
dual-UMI RRBS reads的格式如下:
Read 1  5' UUUUUUUU CAGTA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF TACTG UUUUUUUU 3'
Read 2  3' UUUUUUUU GTCAT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF ATGAC UUUUUUUU 5'
UUUUUUUU是一个8-mer的UMI(unique molecular identifier),CAGTA是一个固定碱基序列,FFFFFFF...是RRBS片段测序得到的实际序列。Read 1 (R1)和Read 2 (R2)的UMI,以及fixed sequences (F1 or F2),写入read的ID,并从实际序列中去除。例如:
R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT ATCTAGTTCAGTACGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG
R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT    CAATTTTGCAGTACAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA

R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT:R1:ATCTAGTT:R2:CAATTTTG:F1:CAGT:F2:CAGT
CGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG
R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT:R1:ATCTAGTT:R2:CAATTTTG:F1:CAGT:F2:CAGT
CAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA
结果写入.clock_UMI.R1.fq(.gz)  .clock_UMI.R2.fq(.gz)。如果需要切除read 3'端潜在可能包含UMI和fixed sequence大约15bp的序列,命令:trim_galore --paired --three_prime_clip_R1 15 --three_prime_clip_R2 15 *.clock_UMI.R1.fq.gz *.clock_UMI.R2.fq.gz
使用Bismark进行比对,UmiBam的'--dual_index'来去除重复 (https://github.com/FelixKrueger/Umi-Grinder)。
 UmiBam UMI R1:(ATCTAGTT);UMI R2:(CAATTTTG)。

--polyA                 实验功能:识别并去除序列中的poly-A序列。设置该参数,软件自动识别包含'AAAAAAAAAA'或'TTTTTTTTTT'的read。对于双末端Read 1或单末端文件,去除PolyA或PolyT。对于双末端Read2,使用互补碱基进行修剪。默认使用A{20}对Read1进行3'端修剪,T{150}对Read2进行5'端进行修剪。这些参数可以通过-a和-a2修改。使用 _ 代替空格,修剪后的信息储存在read ID中,用于后续识别PolyA修剪过的序列。在read ID前后分别添加"32:A:"和"_PolyA:32"标签,例如:
@READ-ID:1:1102:22039:36996 1:N:0:CCTAATCC
GCCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAATAAAAACTTTATAAACACCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
@32:A:READ-ID:1:1102:22039:36996_1:N:0:CCTAATCC_PolyA:32
GCCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAATAAAAACTTTATAAACACC
注意:在进行poly-A修剪之前,首先要进行adapter和质量值控制,建议的分析流程如下:
1) trim_galore file.fastq.gz
2) trim_galore --polyA file_trimmed.fq.gz
3) zcat file_trimmed_trimmed.fq.gz | grep -A 3 PolyA | grep -v ^-- > PolyA_trimmed.fastq
1) 第一步去除质量值和Illumina adapter的污染,2) 查找并去除PolyA污染,3) 如果需要的话,可以输出PolyA修剪的序列到特定的FastQ文件。


# RRBS选项 (针对MspI处理的样本):

--rrbs                  指定输入文件类型为MspI降解的RRBS样本(识别位点: CCGG)。 裁剪单末端或双末端Read 1序列3'端的2bp序列。经过低质量裁剪的read不再处理。双末端Read 2文库会额外从5'端的前2bp的碱基(通过 '--clip_r2 2'选项设置)。避免测序片段中靠近3' MspI位点填充的胞嘧啶位点对后续甲基化分析的影响。不使用NuGEN ovation RRBS System 1-16 kit处理的样本(说明见下文)。

--non_directional       针对non-directional RRBS文库设置,修剪以 'CAA'  'CGA' 开头的read,去除这类read的前两个碱基。消除末端修复反应中使用胞嘧啶位点对后续甲基化分析的影响(参考'--rrbs' 选项)。'--non_directional' 需要和 '--rrbs'一起使用,注意在双末端模式时,不要设置 '--clip_r2 2'
--keep                  保留修剪过程中的中间文件。默认: off。不设置此参数,这些中间文件会在软件运行完成后删除。只对RRBS样本有效。
对于使用NuGEN Ovation RRBS System 1-16 kit的RRBS:由于NuGEN Ovation kit在每个MspI位点后添加了可变长度(0-3)的核苷酸序列,该模式不需要设置 --rrbs。这一步的去除在后续步骤中进行(查看NuGEN Ovation kit的参考手册)MseI RRBS:使用MseI(识别位点:TTAA)代替MspI对DNA进行降解的数据,不需要设置 --rrbs  --non_directional,基本所有的序列都使用'TAA'开头, 'TAA' 限制性位点的末端修复反应不涉及任何胞嘧啶cytosines,因此不需要特殊处理,以标准模式运行即可。


双末端参数:						
--paired                指定输入文件是双末端测序序列(质量值/adapter/RRBS)。 双末端序列需要一个的最小长度作为阈值,通过 --length 选项指定。如果只有一个read长度达到阈值,参考 --retain_unpaired 选项,使用此参数去除较短的read对,不影响FastQ文件中的序列排序,对于双端测序结果,一对reads中,如果有一个被剔除,那么另一个会被同样修剪,而不管是否达到标准。许多比对软件对FastQ的顺序有要求。Trim Galore输出成对的文件名,如 file1_1.fq file1_2.fq SRR2_1.fq.gz SRR2_2.fq.gz ... .

-t/--trim1              从每个read的3'端切除1bp。如果计划使用Bowtie1将FastQ文件作为双末端数据比对。需要设置此参数。
Bowtie (1) 比对如下图:
R1 --------------------------->
R2 <---------------------------
# 或者:
R1 ----------------------->
R2       <-----------------
起始/终止位点包含在其他read中。
注意: 如果你计划使用Bowtie2, BWA 等,不需要设置此参数

--retain_unpaired       对于双端测序结果,一对reads中,如果一个read变得很短,较长的read会被单独保存以'.unpaired_1.fq'或'.unpaired_2.fq'为后缀名的文件内。长度的阈值通过 -r1/--length_1 和 -r2/--length_2。 默认:OFF

-r1/--length_1 <INT>    保留Unpaired single-end read 1长度的阈值写入以'.unpaired_1.fq' 为后缀名的输出文件。这些序列可以以单末端形式比对到参考基因组上。默认: 35 bp.
-r2/--length_2 <INT>    保留Unpaired single-end read 2长度的阈值写入以 '.unpaired_2.fq' 为后缀名的输出文件。这些序列可以以单末端形式比对到参考基因组上。默认: 35 bp.

trim_galore的使用示例 编辑

# 创建4.trim_galore的文件夹来存储修剪过的
mkdir ~/trim_galore/
# 切换到4.trim_galore的文件夹
cd ~/trim_galore/
# 将2.raw_fq目录下的fastq文件链接到4.trim_galore目录下
ln -s ~/raw_fq/*.fastq.gz ~/trim_galore/

# 使用for循环批量生成命令列表,并写入到trim_galore.command文件
for i in `ls *_1.fastq.gz`
do
i=${i/_1.fastq.gz/}
echo "trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired -o ~/4.trim_galore ${i}_1.fastq.gz ${i}_2.fastq.gz"
done > trim_galore.command
# 完成fastq文件质量控制后,整合修剪后的质量报告
echo "multiqc ~/trim_galore/*zip -o ~/trim_galore/" >> trim_galore.command
# 检查生成的命令列表是否正确
cat trim_galore.command
# 批量运行质量控制
sh trim_galore.command
# 查看生成的结果
ls -l ~/4.trim_galore/