生物信息學/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/