使用trim_galore对数据进行质量控制-过滤

写在前面
参考1
参考2

cutadapt 软件可以对NGS数据进行质量过滤
FastQC 软件可以查看NGS数据的质量分布
trim_galore 将这两个软件封装到一起，使用起来更加方便

• Trim galore，是可以自动检测adapter
• 去除reads 3’端的低质量碱基 # 自动调用cutadapt
• 去除adapter序列 # 自动去除3’端的adapter,可以通过设定--Illumina，--small_rna，--nextera参数来指定对应的adapter类型
• 去除长度太短的序列 #通过设定--length参数，小于设定值被去除
• 其它过滤

说明

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

• Illumina: AGATCGGAAGAGC # 如果输入参数 --Illumina,就会默认trimmed前13bp的adapter
• Small RNA: TGGAATTCTCGG # 同上 如果输入参数 --small_rna,就会默认trimmed前12bp的adapter
• Nextera: CTGTCTCTTATA # 同上 如果输入参数 --nextera,就会默认trimmed前12bp的adapter

对于单端测序数据，基本用法如下
trim_galore --quality 20 -a AGATCGGAAGAGC --length 20 -o out_dir input.fq
其中-a后面可以跟着序列（-a AGATCGGAAGAGC）
对于双端测序数据，基本用法如下
trim_galore --paired --quality 20 -a AGATCGGAAGAGC -a2 AGATCGGAAGAGC --length 20 -o out_dir R1.fq.gz R2.fq.gz

参数说明

--quality/-q<INT>：设定Phred quality score阈值，默认为Phred 20 切除质量得分低于设定值的序列
--phred33/64：使用ASCII+33/64质量得分作为Phred得分选择-phred33或者-phred64，表示-测序平台使用的Phred quality score。
（需要确认:anger/Illumina 1.9+ encoding为33; Illumina 1.5 encoding为64）
--adapter/-a ：输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
-s/--stringency<INT>：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length<INT>：设定输出reads长度阈值小于设定值会被抛弃，默认值20bp；即小于20bp的被去除。注意，在pe150下，不要设置太高，可以50或36（默认20）
--max_length : 设置长度大于此值被丢弃
-e <ERROR RATE> ：默认0.1
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
-o/--output_dir：输入目录 [需要提前建立目录，否则运行会报错]。
-- trim-n : 移除read一端的reads
-j :使用线程数, 注意假设已使用Python 3并安装了Pigz，那么内核设置为4，实际使用内核是15，因此，最高设为4.
--fastqc #当分析结束后，使用默认选项对结果文件进行fastqc分析

ls *.fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done
## Sequence length 51
## %GC 39
## Adapter Content passed
mkdir QC_results
mv *zip *html QC_results/

https://www.cnblogs.com/weipeng-loaded/p/10601285.html

FastQC（测序质量分析）

$: fastqc -q -t 4 -o ./fastqc_result/ *.fastq.gz &
> -t：调用核心数目
> -q：安静运行，运行过程中不会生成报告，在结束时将报告生成一个文件
> -o ../path/to/file ：文件输出位置
> *. fq.gz：输入文件，**当前目录下**所有名字中有“ .fq.gz ”的文件

查看QC结果
单个查看：鼠标双击打开html文件查看
批量查看：使用MultiQC软件: multiqc *fastqc.zip

Fastqc结果报告关注重点：

1）.basic statistics
2）.per base sequence quality
3）.per base sequcence content
4）.adaptor content
5）.sequence duplication levels
最主要的几个指标是

• GC含量
• Q20和Q30的比例以及是否存在接头（adaptor）
• index以及其他物种序列的污染等。

先看利用fastqc检测原始序列的质量：

在*.gz所在目录：
ls *.gz | while read id ; do /home/yangjy/miniconda3/envs/chipseq/bin/fastqc $id;done

从页面左侧的的【summary】中可以看出有哪些选项没有通过，可以看出此数据的测序质量。从【Basic Statics】看出数据的序列数量，测序平台以及GC含量等相关信息

SRR391033_fastqc.html

Per base sequence quality : 每个read各位置碱基的测序质量。
横轴碱基的位置，纵轴是质量分数，
Quality score= -10log10p
（p代表错误率），所以当质量分数为40的时候，p就是0.0001，质量算高了。
[红色线]代表中位数，
[蓝色线]代表平均数
[黄色]是25%-75%区间，
[权]是10%-90%区间。
若任一位置的下 [四分位数] 低于10或者 [中位数] 低于25---->出现 “警告”
若任一位置的下 [四分位数] 低于5或者 [中位数] 低于20，出现“失败，Fail”

Per base sequence quality

 此图中
- 横轴是测序序列第1个碱基到第51个碱基
- 纵轴是质量得分，Q = -10*log10（error P）即20表示1%的错误率，30表示0.1%
- 每1个boxplot，都是该位置的所有序列的测序质量的一个统计:
*上面的bar是90%分位数
*下面的bar是10%分位数，
*箱子的中间的横线是50%分位数(上边是75%分位数，下边是25%分位数)
- 图中蓝色的细线是各个位置的平均值的连线(一般要求此图中，所有位置的10%分位数大于20,即Q20过滤)
。。上面的这个测序结果，需要把后面的27bp以后的序列切除，从而保证后续分析的正确性
【Failure 报错 如果任何碱基质量低于5,或者是任何中位数低于20】

较好

Per tile sequence quality 每tile的碱基质量
【蓝色】代表测序质量很高，【暖色】代表测序质量不高

Per tile sequence quality

此图中
- 横轴代表51个碱基的每个不同位置
- 纵轴是tail的Index编号
- 这个图主要是为了防止，在测序过程中，某些tile受到不可控因素的影响而出现测序质量偏低
。。 某些tail出现暖色可以在后续分析中把该tail测序的结果全部都去除

Per sequence quality scores 各个序列质量的分布情况

Per sequence quality scores

假如1条序列长度为51bp,则这51个位置的 每个位置Q之的平均值就是这条reads的质量值
- 该图横轴是0-33，表示Q值
- 纵轴 对应的reads数目
数据中，测序结果主要集中在高分中，证明测序质量良好！

Per base sequence content，碱基分布 ; 对所有reads的每一个位置，统计ATCG四种碱基的分布，理论上来说，A和T应该相等，G和C应该相等，但是一般测序的时候，刚开始测序仪状态不稳定，很可能出现严重分离的情况。像这种情况，即使测序的得分很高，也需要cut开始部分的序列信息，一般像这种情况，会cut前面5-10bp。

当任一位置的A/T比例与G/C比例相差超过10%，报'WARN'；当任一位置的A/T比例与G/C比例相差超过20%，报'FAIL'

Per base sequence content

- 横轴为位置，
- 纵轴为碱基含量，
* 正常情况下每个位置每种碱基出现的概率是相近的，四条线应该平行且相近。
* ！当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。

Per sequence GC content 序列平均GC含量分布图

Per sequence GC content

- 横轴是0 - 100%； 
- 纵轴是每条序列GC含量对应的数量
- 蓝色的线是程序根据经验分布给出的理论值
- 红色是真实值，两个应该比较接近才比较好
- 当红色的线出现双峰，基本是混入了其他物种的DNA序列
- 这张图中的信息良好

Per base N content 序列中各个位点的N含量，越小越好

Per base N content

sequence_length_distribution 序列测序长度统计，从图中可以看出序列的平均长度为51。
每次测序仪测出来的长度在理论上应该是完全相等的

sequence_length_distribution

当测序的长度不同时，如果很严重，则表明测序仪在此次测序过程中产生的数据不可信

Sequence Duplication Levels 指在测序前建库PCR过程中 由一些序列扩增次数过多导致的。若重复较高则需要进行处理这些dup。

Sequence Duplication Levels
较好

Overrepresented sequences
有某个序列大量出现被称为(over-represented）；fastqc的标准是占全部reads的0.1%以上。
*为了计算方便，只取了fq数据的前200,000条reads进行统计，所以有可能over-represented reads不在里面。而且大于75bp的reads也是只取50bp。
当发现超过总reads数0.1%的reads时报"黄色!"
当发现超过总reads数1%的reads时报"红色×"
正常显示No overrepresented sequences

Overrepresented sequences

Adapter Content
衡量 序列中两端adapter的情况
*如果在 fastqc分析时候 -a选项没有内容，则默认使用图例中的四种通用adapter序列进行统计

Adapter Content

有adapter序列没有去除干净的情况
在后续分析的时候需要先使用cutadapt软件进行去接头也可以用 trimmomatic来去除接头

使用trim_galore对数据进行质量控制-过滤

trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq --gzip -o ../clean

重新用fastqc检测进行过滤后的reads质量

fastqc -o out_dir *fq.gz
multiqc *fastqc.zip --ignore *.html

附注，处理信息

Multicore support not enabled. Proceeding with single-core trimming.  [未启用多核支持]
Path to Cutadapt set as: 'cutadapt' (default)
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt version: 1.18                             Cutadapt版本： 1.18 
single-core operation.
Output will be written into the directory: /home/yangjy/chipseq_test/clean/  

AUTO-DETECTING ADAPTER TYPE
===========================
Attempting to auto-detect adapter type from the first 1 million sequences of the first file (>> SRR391032.fastq.gz <<)
Found perfect matches for the following adapter sequences:
Adapter type Count Sequence Sequences analysed Percentage
Illumina 1926 AGATCGGAAGAGC 1000000 0.19
smallRNA 1 TGGAATTCTCGG 1000000 0.00
Nextera 0 CTGTCTCTTATA 1000000 0.00
Using Illumina adapter for trimming (count: 1926). Second best hit was smallRNA (count: 1)
Writing report to '/home/yangjy/chipseq_test/clean/SRR391032.fastq.gz_trimming_report.txt'
SUMMARISING RUN PARAMETERS
==========================
Input filename: SRR391032.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.6
Cutadapt version: 1.18
Number of cores used for trimming: 1
Quality Phred score cutoff: 20
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 20 bp
Output file(s) will be GZIP compressed
Cutadapt seems to be reasonably up-to-date. Setting -j 1
Writing final adapter and quality trimmed output to SRR391032_trimmed.fq.gz

输入文件名：SRR391032.fastq.gz
整理模式：配对端
Trim Galore版本：0.6.6
Cutadapt版本：1.18
用于微调的内核数：1
Quality Phred分数截止：20
已选择质量编码类型：ASCII + 33 适配器序列：“ AGATCGGAAGAGC”（Illumina TruSeq，Sanger iPCR；自动检测）
最大修整错误率：0.1（默认值）
适配器所需的最小重叠量（严格度）：3 bp
最小删除序列对之前，两次读取都需要序列长度：20 bp
输出文件将经过GZIP压缩
设置-j 1
将最终适配器和质量调整后的输出写入SRR391032_trimmed.fq.gz

再回顾一下我们的命令：
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired fq --gzip -o ../clean

##细节：

出现了一个报错

first_line = file.readline()
AttributeError: 'str' object has no attribute 'readline'

与python有关；解决方案
https://www.cnblogs.com/lovebing/p/13048948.html
https://blog.csdn.net/qq_41185868/article/details/82079079
########暂时没有解决!

其他软件：
数据过滤之 Trimmomatic 详细说明