fastqc和Trimmomatic的使用

news/发布时间2024/8/25 17:05:16

1.FastQC分析检测报告

在先前的记录中，我们已经得到了我们的QC报告，现在要针对我们的报告对原始数据进行过滤
其中
和
都表明该数据需要去接头，并对序列进行处理

2.Trimmomatic的下载

首先，使用conda安装Trimmomatic
conda install Trimmomatic
就可以安装完毕了，安装完使用trimmomatic -h检测是否安装成功

成功了是这样的。
注意：安装前记得先安Java
https://www.jianshu.com/p/43b564783e32

3.开始过滤

可以参考
https://blog.csdn.net/I_LiYY/article/details/105533946
在开始过滤之前，先准备好几个东西：
①.确认好数据是phred33还是phred64，具体的讲解内容在这个帖子里：
https://www.jianshu.com/p/248308513e2e
但是这里只要区分清楚这两种是有区别的就可以，不要搞错了，可以使用脚本来判断，原理是解压1000条出来看看是哪种编码

使用的脚本在这里，非常方便：https://www.jianshu.com/p/9ceabb21be12
②.提前写好的命令：按照自己的需求提前写好，现场一个一个输很麻烦的，而且写好以后可以批量处理
使用命令
trimmomatic PE -threads 8 -phred33 SRR13810477_1.fastq.gz SRR13810477_2.fastq.gz paired_1_R1_paired.fq.gz unpaired_1_R1_unpaired.fq.gz paired_1_R2_paired.fq.gz unpaired_1_R2_unpaired.fq.gz ILLUMINACLIP:/root/anaconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:30 LEADING:3 TRAILING:3 MINLEN:30
主要分为三部分，这里解释一下其中各参数都代表什么：
第一部分：
PE：是双端模式，给双端测序数据用的，如果是单端测序，用SE
-threads：这是线程数，你可以选择合适的线程加快进度，不然默认单核，超慢
-phred33：这就是之前提到的，不多赘述
SRR13810477_1.fastq.gz SRR13810477_2.fastq.gz：这是你要处理的数据，因为是双端测序所以有两个
paired_1_R1_paired.fq.gz unpaired_1_R1_unpaired.fq.gz paired_1_R2_paired.fq.gz unpaired_1_R2_unpaired.fq.gz：这是要输出的四个文件，输出文件有四个，使用 -baseout 参数指定输出文件的 basename，软件会自动为四个输出文件命名，过滤之后双端序列都保留的就是 paired，反之如果其中一端序列过滤之后被丢弃了另一端序列保留下来了就是 unpaired(即成对的clean reads, 未成对的正向序列以及未成对的反向序列)，如图：
接下来的参数主要参与去接头的第二部分：
ILLUMINACLIP:/root/anaconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10
/root/anaconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa：参数后面分别接adapter序列的fasta文件：第一步 seed 搜索时允许的最大错配碱基个数2：palindrome模式下匹配碱基数阈值30：simple模式下的匹配碱基数阈值10(7-15之间）：palindrome 模式允许切除的最短接头序列为 8bp（默认值）：palindrome 模式去除与 R1 完全反向互补的 R2（默认去除false），但在有些情况下，例如需要用到 paired reads 的 bowtie2 流程，就要将这个参数改为 true，否则会损失一部分 paired reads。

按照规定顺序，ILLUMINACLIP 各个参数之间以冒号分开，PE测序需要注意最后一个参数。对于SE测序最后两个参数可以不设置

其中，导入adapter序列的fasta文件时，记得使用绝对路径，不明白的朋友可以直接find来找，如find / -name TruSeq3-PE-2.fa，找到路径之后可以去看一下，是这个样子的：
一般测序公司在给你数据的时候这个也会给你，不过Trimmomatic自带这些Illumina 平台的接头，所以直接用也可以
第三部分，就是关于过滤剪切的参数了，强烈建议一步一步来，先去接头，再过滤剪切。看自己的需要调整选择：
LEADING:3 切除首端碱基质量小于3的碱基
#Illumina平台有些低质量的碱基质量值被标记为 2 ，所以设置为3可以过滤掉这部分低质量碱基。
TRAILING:3 切除尾端碱基质量小于3的碱基
SLIDINGWINDOW:15:20
滑窗质量过滤，一般一个read的低质量序列都是集中在末端，也有很少部分在开头。从5'端开始进行滑动，当滑动位点周围一段序列(window)的平均碱基低于阈值，则从该处进行切除。Windows的size是15个碱基（一般设置在10-30之间），其平均碱基质量小于20，则切除
MINLEN:50 可被保留的最短reads长度，应根据原始序列的长度而定
HEADCROP: 在reads的首端切除指定的长度
CROP: 保留reads到指定的长度
TOPHRED33 将碱基质量转换为pred33格式
TOPHRED64 将碱基质量转换为pred64格式