之前分析过的测序数据,数据质量都很好,给了我一个错觉质控的前后差别不是很大,内心里对质控这一步也就不是很重视,跑完质控有时也懒得看结果,拿到一批测序数据后,也总是忽略了去看测序策略是什么,按照固定的流程用fastqc做质控和trimmomatic去除adpater和低质量的碱基,结果今天翻车了!!
批量跑完了20个数据的ATAC-seq流程,结果发现样本的比对率参差不齐,有的能达到80-90%,有的却低于50%,回过头来找原因发现cleandata中的adpter并没有去除,也就是说我用的adapter并不是建库所用的,那么怎么知道这批数据建库时用的什么adapter呢?
求助了健明师兄,他推荐我使用Trim galore做质控,并且一眼看出我这个测序策略是nextseq(不知道他怎么看出来的)。
查了一下Trim galore,是可以自动检测adapter,也发现了自己的错误,trimmomatic只是针对Illumina高通量测序平台设计的接头去除和低质量reads清洗软件,Nextera的接头和它是不一样的(基础知识很重要!!!)。
下面就对Trim galore的下载,安装和使用做一个简要介绍,并总结了另外两个质控软件Trimmomatic 和cutadapter 的使用。
1. Trim galore简介
Trim Galore是对FastQC和Cutadapt的包装。适用于所有高通量测序,包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera 和smallRNA测序平台的双端和单端数据。主要功能包括两步:
第一步首先去除低质量碱基,然后去除3' 末端的adapter, 如果没有指定具体的adapter,程序会自动检测前1million的序列,然后对比前12-13bp的序列是否符合以下类型的adapter:
- Illumina: AGATCGGAAGAGC
- Small RNA: TGGAATTCTCGG
- Nextera: CTGTCTCTTATA
2. 下载安装软件 trim galore
conda安装:
conda install trim-galore
conda安装时可以看出依赖的环境很多,我们的大机环境很复杂,我并没有安装成功。
下载安装包安装:
下载安装包安装很简单,下载后解压,配置下环境变量就可以使用。
## 需先安装fastqc和cutadapt
wget https://github.com/FelixKrueger/TrimGalore/archive/0.4.5.tar.gz
tar zxvf 0.4.5.tar.gz
3.使用
# 处理双端测序结果
echo " trim_galore cut adapters started at $(date)"
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 \
--paired $dir/cmp/01raw_data/$fq1 $dir/cmp/01raw_data/$fq2 \
--gzip -o $input_data
echo "trim_galore cut adapters finished at $(date)"
参数说明:
--quality:设定Phred quality score阈值,默认为20。
--phred33::选择-phred33或者-phred64,表示测序平台使用的Phred quality score。
--adapter:输入adapter序列。也可以不输入,Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个,也直接显式输入这三种平台,即--illumina、--nextera和--small_rna。
--stringency:设定可以忍受的前后adapter重叠的碱基数,默认为1(非常苛刻)。可以适度放宽,因为后一个adapter几乎不可能被测序仪读到。
--length:设定输出reads长度阈值,小于设定值会被抛弃。
--paired:对于双端测序结果,一对reads中,如果有一个被剔除,那么另一个会被同样抛弃,而不管是否达到标准。
--retain_unpaired:对于双端测序结果,一对reads中,如果一个read达到标准,但是对应的另一个要被抛弃,达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip:清洗后的数据zip打包或者不打包。
--output_dir:输入目录。需要提前建立目录,否则运行会报错。
-- trim-n : 移除read一端的reads
其他质控方法:
Trimmomatic
Trimmomatic是针对Illumina高通量测序平台设计的接头去除和低质量reads清洗软件。软件中包括有Illumina平台常见接头序列,可以很方便处理单端和双端RNA-seq数据。Trimmomatic也支持自己设计要去除的接头序列文件。
运行代码:
## 双端测序
echo "trimmomatic cut adapters started at $(date)"
java -jar /software/biosoft/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 8 $dir/cmp/01raw_data/$fq1 $dir/cmp/01raw_data/$fq2 \
$input_data/$sample\_paired_clean_R1.fastq.gz \
$input_data/$sample\_unpair_clean_R1.fastq.gz \
$input_data/$sample\_paired_clean_R2.fastq.gz \
$input_data/$sample\_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/software/biosoft/software/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
echo "trimmomatic cut adapters finished at $(date)"
重要参数解释:
-threads
:设置线程数目。-phred33
:选择-phred33
或者-phred64
,表示测序平台使用的Phred quality score。查询方法:首先,运行FastQC,在结果报告第一项会“猜出”测序平台。之后,查询平台对应Phred列表。-trimlog
:输出运行日志。日志中包括对每一个read具体选择数据,所以文件会比较大。ILLUMINACLIP
:跟随四个参数,分别是:<fastaWithAdaptersEtc>
为adaptesr文件完整路径(在Trimmomatic的默认安装目录下的adapter
,有整理好的);<fastaWithAdaptersEtc>
为seed matches(16bases)在匹配时的最大错配数目;<palindrome clip threshold>
对于一对reads当得分超过30(约50 bases),seeds会被延伸和固定;<simple clip threshold>
,对于单端reads当得分超过10(约17 bases),seeds会被延伸和固定。LEADING
和TRAILING
:分别为去除read头部和尾部的低质量(低于quality3)碱基数目。SLIDINGWINDOW
:跟随两个参数,分别是<windowSize>
为扫描“窗口”长度;<requiredQuality>
为窗口碱基质量的平均阈值,低于此会被删除。MINLEN
:设置最短reads数目。需要根据下游alignment软件设定,比如Bowtie适用于短序列,比如50bp以下;而Bowtie2适用于50bp以上。TopHat 则根据实际使用Bowtie或者Bowtie2选择。
FastQC
FastQC是用于对二代测序数据质量快速检验的工具,可以输入fastq(fastq.gz)、sam或者bam文件。查看输出结果解释。可以联合multiqc使用,查看多个qc的报告。
代码示例:
echo "fastqc started at $(date)"
###method (3)
fastqc -o $dir/cmp/qc $dir/cmp/01raw_data/*gz
#multiqc *fastqc.zip --ignore *.html
echo "fastqc finished at $(date)"
cutadapter
cutadapt -q 30 -b CTGTCTCTTATACACATCTGACGCTGCCGACGA --minimum-length 20 --overlap=5 -o tmpl1.1.fastq --paired-output tmpl1.2.fastq SRR2920469_1.fastq.gz SRR2920469_2.fastq.gz