利用Trimmomatic进行转录组原始数据过滤：小果吐血整理！
生信人R语言学习必备
立刻拥有一个Rstudio账号
开启升级模式吧
（56线程，256G内存，个人存储1T）
随着高通量测序技术的发展，转录组测序已经成为了研究基因表达的重要手段。但是，由于实验条件、测序平台等原因，转录组测序数据中常常存在一些噪声，这些噪声会影响我们对基因表达情况的判断。
因此，在进行进一步的差异分析、功能注释和生信分析之前，我们需要对转录组原始数据进行过滤，去除掉低质量的序列和噪声，提高数据的可靠性和准确性。
一、什么是转录组
转录组是指在一个时刻内细胞内转录生成的RNA的总和。与基因组相比，转录组更反映了细胞内基因的表达情况，能够为我们揭示细胞内部的基因调控网络以及基因功能的变化。因此，转录组测序已经成为研究基因功能和基因调控机制的重要手段。
二、为什么要过滤转录组原始数据？
转录组测序产生的原始数据往往具有以下问题：
低质量的序列：由于测序仪器或反应体系等因素的影响，会产生一些低质量的测序序列，这些序列不仅会占用服务器空间和计算资源，也会影响后续的数据分析结果。
读长偏差：由于测序仪器等因素的限制，测序数据中可能存在长度偏差，如头尾序列较短或较长等，这些偏差会对基因表达水平的计算造成影响。
噪声：由于实验条件等原因，测序数据中常常会出现一些异质性，比如SNP（单核苷酸多态性）、INDEL（插入/缺失）等，这些异质性会影响后续的差异分析和功能注释。
因此，我们需要对原始数据进行过滤，去除掉低质量的序列和噪声，提高数据的可靠性和准确性。
三、如何对转录组原始数据进行过滤？
目前，有很多工具可以对转录组原始数据进行过滤。其中比较常用的工具包括Trimmomatic、Fastp和Cutadapt等。下面小果以Trimmomatic为例，介绍如何对转录组原始数据进行过滤。
1.准备工作：仍然是用miniconda安装Trimmomatic
conda create –name trimmomaticconda activate trimmomaticconda install -c bioconda trimmomatictrimmomatic -version
可以看到安装成功，版本为0.39
2. 运行Trimmomatic
对于双末端测序来说，可以运行以下命令：
java -jar trimmomatic-0.39.jar PE input_forward.fq.gz input_reverse.fq.gzoutput_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:True LEADING:3 TRAILING:3 MINLEN:36
为了方便理解小果将解释放在下面：
input_forward.fq.gz 和 input_reverse.fq.gz是指输入的双端测序数据，包括 forward 和reverse reads。output_forward_paired.fq.gz 和 output_reverse_paired.fq.gz是指输出的过滤后的双端测序数据，包括 forward 和 reverse reads，且两端 reads 均通过质量控制和过滤。output_forward_unpaired.fq.gz 和 output_reverse_unpaired.fq.gz指输出的未匹配的reads，其中output_forward_unpaired.fq.gz 包括未匹配的 forward reads， output_reverse_unpaired.fq.gz 包括未匹配的 reverse reads。ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:True: 使用 TruSeq3-PE.fa 作为 adapter 序列文件，设定的参数依次为：seed mismatches (2), palindrome clip threshold (30), simple clip threshold (10), keep both reads (True)。LEADING:3: 将 reads 开头的质量值低于 3 的部分切除。TRAILING:3: 将 reads 末尾的质量值低于 3 的部分切除。MINLEN:36: 丢弃长度小于 36 bp 的 reads。
!!!注意：小果高亮的两个文件，在执行命令时要输入正确的文件路径和文件名称。另外因为在实际工作中我们往往需要处理很多数据，比如这次小果准备过滤9个文件，每个.sra文件又转化为了read1，read2两个文件，所以小果为了偷懒写了个脚本，一次性过滤所有数据，脚本支持双端测序，过滤的数据放在同一个文件夹下，其中的路径根据自己的修改即可。
脚本如下：
#!/bin/bash# 定义输入和输出目录input_dir=~/my_folder/trimmoutput_dir=~/my_folder/trimm/ninanjie_output# 遍历输入目录中以 _1.fastq 结尾的所有文件.gzfor file1 in ${input_dir}/*_1.fastq.gzdo# 获取基本文件名（删除 _1.fastq.gz 扩展名）base=$(basename ${file1} _1.fastq.gz)# 定义相应文件2的路径file2=${input_dir}/${base}_2.fastq.gz# 定义输出文件名out_forward_paired=${output_dir}/${base}_forward_paired.fq.gzout_forward_unpaired=${output_dir}/${base}_forward_unpaired.fq.gzout_reverse_paired=${output_dir}/${base}_reverse_paired.fq.gzout_reverse_unpaired=${output_dir}/${base}_reverse_unpaired.fq.gz# 使用nohup在后台运行Trimmomaticnohup java -jar /media/desk16/iyun003/miniconda3/envs/trimmomatic/share/trimmomatic-0.39-2/trimmomatic.jarPE${file1}${file2}${out_forward_paired}${out_forward_unpaired}${out_reverse_paired}${out_reverse_unpaired}ILLUMINACLIP:/media/desk16/iyun003/miniconda3/envs/trimmomatic/share/trimmomatic-0.39-2/adapters/TruSeq3-PE.fa:2:30:10:2:TrueLEADING:3TRAILING:3MINLEN:36 > ${output_dir}/${base}_trimmomatic.log 2>&1 &Done
查看其中一个文件的日志，最后一行显示成功！小果泪目！
部分输出结果
单末端测序运行以下命令：
java -jar trimmomatic-0.35.jar SE -phred33 input.fq.gz output.fq.gzILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
还是一样注意修改文件路径和文件名称！
好啦，不多说了，今天的内容就到这里了，你学会了吗！
欢迎使用：云生信 – 学生物信息学 (biocloudservice.com)
如果想用服务器可以联系微信：18502195490（快来联系我们使用吧！）
扫码加小果
领取生信大礼包
点击“阅读原文”立刻拥有
↓↓↓