【RSS】[3] Filtering FASTQ files & cutadapt用法

在分析了QA报告之后，人们可能想放弃some of the reads based on several criteria;
我们将使用两个不同的工具来执行这些过滤步骤：PRINSEQ 和 fastq mcf;

PRINSEQ提供了广泛的过滤选项，我们可以在manual中了解它们：

PRINSEQ lite-h

根据我们从质量保证报告中了解到的情况，我们可以决定应用following filters：

zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt1 \
    -out_bad null \
    -trim_qual_right 30 -min_len 32 \
    -ns_max_p 5 \
    -lc_method dust -lc_threshold 10
    
zcat SRR031714_2.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_2_filt1 \
    -out_bad null \
    -trim_qual_right 30 -min_len 32 \
    -ns_max_p 5 \
    -lc_method dust -lc_threshold 10

练习：我们使用哪些标准来丢弃reads？
解决方案：用prinseq-lite -h查看参数，指定标准，去除不符合标准参数的fastq；

练习：如果我们的文件是phred 64格式的，我们应该使用哪个option？
答：我们应该使用-phred64选项。正确选择此选项非常重要，否则软件将无法正确解释质量字符串，这将影响我们的输出。

练习：过滤完fastq文件后，再次fastQC QA报告来决定我们是否对结果满意，你认为我们能解决最初发现的主要问题吗？
答案：我们可以使用以下命令在筛选的fastq文件上运行FastQC：

fastqc SRR031714_1_filt1.fastq SRR031714_2_filt1.fastq

然后我们可以打开并检查报告。

练习：通常在过滤步骤之前和之后，跟踪fastq文件中可用的读取次数也是有用的。你能从FastQC报告中收集这些信息吗？现在想象一下我们正在处理更多的文件；你能想出一个更自动化的方法来检查吗？
答案：我们可以从FastQC报告的第一个表中获取每个fastq文件的读取次数，如果我们希望自动化流程，还可以执行以下命令：

grep '^@SRR' SRR031714_1_filt1.fastq | wc -l
# find the lines starting with "@SRR" and count how many there are
# 5,150,415
# 利用wc指令我们可以计算文件的Byte数、字数、或是列数
# -l或--lines 只显示行数

请注意，在某些筛选步骤中过于严格可能会导致信息的严重丢失。从一开始就丢弃太多的读取通常不是一个好的选择，而且还可以依赖mapping step来丢弃低质量的读取。

如果使用paired-end data，筛选步骤可能会变得复杂，因为必须确保两个fastq文件（每个读取对一个）以相同的顺序包含相同数量的读取。有一些工具可以简化此任务，例如fastq mcf。我们只需输入工具名称即可打印可用选项列表：

fastq-mcf

我们注意到fastq mcf的主要功能是从fastq文件中删除adapter。因此，我们需要提供一个带有adapter序列的fasta文件。在我们的案例中，我们决定只检查标准Illumina paired-end adapter：

cat adapters.fa

现在我们已经很好地理解了所需的输入，我们可以继续执行该工具，尝试匹配我们之前使用的PRINSEQ：

fastq-mcf adapters.fa SRR031714_1.fastq.gz SRR031714_2.fastq.gz \
    -o SRR031714_1_filt2.fastq -o SRR031714_2_filt2.fastq \
    -q 30 -P 33 -l 32 --max-ns 1

根据所使用的filtering options，我们可能最终得到一组具有不同长度的reads，导致我们可能会在下游分析中遇到一些困难；

因此，为了方便和safe side，我们可以使用针对所有reads的filtering options；

zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt3 \
    -out_bad null \
    -trim_right 5
    
zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt3 \
    -out_bad null \
    -trim_right 5

# zcat(选项)(参数)
# zcat命令 用于不真正解压缩文件，就能显示压缩包中文件的内容的场合。

总之，您可以应用许多filtering combinations，具体选项在很大程度上取决于type of data和后续分析。我们建议查看PRINSEQ manual以获得关于该主题的详细概述。

使用cutadapt去除adapter

在生物信息学100个基础问题 —— 第10题读懂FastQC报告之adapter与kmer中，我们知道，测序结果中可能会有若干条序列存在adapter的信息，而adapter的信息一般是不在基因组上存在的。所以，在比对之前如果不把adapter去干净，我相信你会得到1个非常非常低的mapping rate。

通常情况下，我们都是使用cutadapt这个软件进行adapter（接头）序列的去除。cutadapt这个软件不但支持单端序列，还支持双端序列的切除，同时还支持gz格式的自动压缩与解压缩。1个常用的切除命令类似：

# 在linux 命令行模式下
cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq

# -a是第1个文件的adapter序列
# -A是第2个文件的adapter序列
# -o是第1个输出文件
# -p是第2个输出文件
# reads.1.fastq 是第1个输入文件，也就是双端测序中的read-1
# reads.2.fastq 是第2个输入文件，也就是双端测序中的read-2

cutadapt中-a/-A 参数与-g/-G参数分别代表什么意思？Illumina测序过程中，一般不会用到哪个参数？
cutadapt可以过滤一些非常短的reads，请解释其中-m 参数是什么意思？为什么要过滤一些非常短的reads？
在测序的过程中，我们经常发现一些序列的3'端的测序质量不太好（如图2所示），即使去掉adapter以后还是需要把低质量的序列再去除1次，从而保证后续的mapping质量。cutadapt可以使用一些办法来去除3'端质量不太好的序列。请说明用哪个参数来设置相关的cutoff，并简要说明cutadapt对read质量判断的策略与方法。

cutadapt软件的介绍

这个cutadapt软件是最常用的去adapter的工具。它是基于Python编写的一个Python包，还是用bioconda安装：

conda install -c bioconda cutadapt

关于adapter的使用汇总

image.png

其实主体还是利用fastqc进行质量评估，然后根据质量评估的信息。进行相应的质控。利用到的软件是cutadapter 和 fastx-toolkit，其实还有一些别的软件可以用，比如说trimmomatic或者是Trim Galore；cutadapter序列过滤参考

illumina的通用adapter序列是什么？

Illumina Paired End Adapters(cannot be used for multiplexing)
Top adapter
    5' ACACTCTTTCCCTACACGACGCTCTTCCGATC*T 3'
Bottom adapter
    5' P-GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG 3'
    
**来源**
http://bioinformatics.cvr.ac.uk/blog/illumina-adapter-and-primer-sequences/
不同颜色的图例代表不同的测序通用的adapter；
如果在当时fastqc分析的时候-a选项没有内容，则默认使用图例中的通用adapter序列进行统计。