抱歉拖了这么久，主要是最后分析出来的数据跟某平台对比后发现差了很多。无法解释这个结果。用了几个软件测试，最后的结果都跟之前平台的Exp结果对不上。最后猜测可能还是这套数据有问题……之前查看数据库的时候也这套数据明明是单端结果但是NCBI和ENA都记录的是双端测序。
原因找到了，数据可能是把双端测序的结果保存在同一个fastq文件里了（感兴趣的可以搜索interleaved fastq file）。做seq之前需要把数据分开，具体的指令为：
cat input.fastq | paste - - - - - - - - | tee | cut -f 1-4 | tr "\t" "\n" | egrep -v '^$' > R1.fastq
cat input.fastq | paste - - - - - - - - | tee | cut -f 5-8 | tr "\t" "\n" | egrep -v '^$' > R2.fastq
懒得弄了，这个数据也没查到对应数据的文章，不折腾了……
接上一节，我们发现了下载的数据中有一些短重复片段，保险起见我们还是需要去除一下接头和一些低质量的读段。在这里我主要介绍两个软件：trimmomatic和fastp。同样，我也会介绍TBtools的trimomatic wrap插件的使用方法。
1. Trimmomatic
老板说图多才有人看
Trimmomatic基本上是illumina专用的去接头软件了。由于Trimmomatic也是基于java的程序，所以TBtools兼容这块只需要做个接口，然后留几个重要参数就好了。
安装Trimmomatic的方法很简单，在随便一个位置打开terminal，然后输入：
conda install -c bioconda trimmomatic
相信大家已经发现了，目前用到的生物信息软件基本都是通过conda安装的。总之，遇事不决就conda，conda不行就再加个上bioconda。
首先Trimmomatic会修剪掉下列四种情况的接头序列或者是人为设定的technical sequence。
A：完整匹配的technical sequence
B：部分完整匹配technical sequence的3'末端
C：长读段中的接头序列
D：部分匹配的3'末端短接头
接下来是一些关键参数：
-
ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the read.
按照你的数据选择接头文件列表TruSeq3对应HiSeq和MiSeq
TruSeq2 (as used in GAII machines)
TruSeq3 (as used by HiSeq and MiSeq machines), -
SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average quality within the window falls below a threshold.
SLIDINGWINDOW:<windowSize>:<requiredQuality>
对应两个参数窗口大小（碱基数）和对应碱基序列的质量。一般就是4和15，没必要乱改。除非数据质量实在是很差。 -
LEADING: Cut bases off the start of a read, if below a threshold quality
因为机器对初始几个序列检测不太准，一般默认依次把质量低于3的碱基切掉 -
TRAILING: Cut bases off the end of a read, if below a threshold quality
同理，尾部也能切掉，不过没必要。尤其是当你数据是双端测序结果的时候 -
CROP: Cut the read to a specified length
直接从中间切断丢弃尾部序列，慎用 -
HEADCROP: Cut the specified number of bases from the start of the read
切掉头部对应碱基数并丢弃，同样，慎用 -
MINLEN: Drop the read if it is below a specified length
这个参数重要也不重要，你需要看一眼你的FastQC结果，一般读段都在100 bp左右，这个时候默认36就好。如果你的读段是50 bp甚至更短，你就需要修改这个参数。改的越低，结果里就有越多的错误读段。
综上，trimmomatic的代码是：
mkdir trimmomatic;
for i in `seq 56 79`;
do
trimmomatic SE -phred33 \
SRR51316${i}.fastq \
trimmomatic/SRR51316${i}_clean.fastq.gz \
ILLUMINACLIP:/home/barnett/miniconda3/share/trimmomatic-0.39-1/adapters/TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:0 SLIDINGWINDOW:4:15 MINLEN:36;
done
这里需要注意一下ILLUMINACLIP的位置，由于版本，平台等问题，接头文件的位置不一定相同，最好用Everything这个软件找一下，填上正确的文件路径。
另外这里介绍的是SE的代码，如果是双端结果那么代码就是：
mkdir trimmomatic;
for i in `seq 56 79`;
do
trimmomatic PE -threads 8 -phred33 \
#需要剪切的双端fastq文件名
SRR51316${i}_1.fastq SRR51316${i}_2.fastq \
#双端保存的位置
trimmomatic/SRR51316${i}_paired_clean_R1.fastq.gz \
trimmomatic/SRR51316${i}_unpaired_clean_R1.fastq.gz \
trimmomatic/SRR51316${i}_paired_clean_R2.fastq.gz \
trimmomatic/SRR51316${i}_unpaired_clean_R2.fastq.gz \
ILLUMINACLIP:/home/barnett/miniconda3/share/trimmomatic-0.39-1/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:0 SLIDINGWINDOW:4:15 MINLEN:36;
done
这里注意一下你的双端fastq文件名需要修改成自己文件的名字。不是很建议把raw data和clean data放在同一位置。个人推荐重新建立一个02_clean_data保存剪切后的结果。
../02_clean_data/filename_paired(or unpaired)_clean_R1.fastq.gz）
2. Fastp
接下来介绍的是另一款质控、剪切一条龙的软件fastp。
Fastp的优点是集成度高，只要一行代码就能完成Fastqc→trimmomatic→Fastqc的流程。
当然作者还表示这软件速度比trimmomatic+fastqc快三倍，但实际上因为本人硬盘的I/O不是很高，实际体验下来区别没那么大。具体代码跟trimmomatic类似：
#fastp做质控和接头剪切
for i in `seq 56 79`;
do
mkdir SRR51316${i};
cd SRR51316${i};
fastp -i ../SRR51316${i}.fastq -o SRR51316${i}.fastq.gz -Q --length_required=25 --cut_front --cut_window_size 4 --cut_mean_quality 15 --compression=6;
cd ..;
done
注意：这里还是只介绍单端的用法。如果你的数据是双端，请参考：
fastp
-i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq
Fastp每次操作会生成一个网页报告，为了让每个网页都能保存下这里选择将每一组数据保存单独的文件夹里。
最后如果你对其它参数和功能感兴趣的可以看看这里：
https://www.jianshu.com/p/6f492058da5b
3. TBtools的Trimmomatic插件
没啥好讲的……枯燥。
点这里查看接头文件
不确定哪种数据就用Merged.Adapter.fa
单端数据是这样
单端格式
双端数据格式就是这样
双端格式
想改参数看这里
操作太简单了不知道该讲啥，用就完事了。
网友评论