1 为什么需要进行QC ?

测序质量问题包括[1]：

• low-confidence bases

• sequence-specific bias

• 3′/5′ positional bias

• polymerase chain reaction (PCR) artifacts

• untrimmed adapters

• sequence contamination

这些问题的存在会影响reads mapping到参考基因组、组装、表达量估计等。有的问题可以通过过滤、剪切、错误或偏差纠正来解决。而有的问题我们是没办法解决的。但至少，我们在解释分析结果的时候，能意识到这些问题的存在[1]。

导致质量问题的原因：

• 在文库构建的时候可能引入污染，比如取水培拟南芥根组织的时候，可能会把无法清除的部分绿藻混入样品中。

• 部分文库在构建时，会使用PCR进行扩增，PCR过度扩增会导致产生过多的一致序列（duplicates），致使后续数据处理（比如 RNAseq表达量计算、BSA计算snp频率等）时，出现误差。

• 在RNAseq文库构建时，一般会加入随机六聚体引物（random hexamers）。以华大的建库为例，华大在进行转录组测序建库的时候，默认使用的是随机引物进行反转录（图一）。随机六聚体引物的使用会引入碱基组成偏离。

• 测序本身有一定的错误率。

• 其他

网页链接：BGI常见问题解答：http://www.bgitechsolutions.com/science/cat-91-95-409.html

质控并没有一个确定标准。在质控过程中，会去除一些低质量reads、切除质量差的碱基、在切除质量差的碱基后还可能会去除一些长度过低的reads……而这些处理，其实相当于我们损失了一些信息。而由于信息的损失，也可能会造成最后结果有偏差。

质控是信息量和信息准确度之间的一个平衡。所以，质控并见得是越严格越好。应该根据自己实际情况的需要来选择质控的严格程度。一般来说，对基因组知道的越少，质控应该越严格[2]。

质控分两个阶段[2]：

• 比对前（Pre-alignment）: 处理“原始数据（raw data）”

• 比对后（Post-alignment）: “数据过滤（data filtering）”

2. 序列质量检测及可视化软件

• FastQC: 只能用于检查测序质量，不能用于质控处理。

• PRINSEQ: 除了可以检测测序质量，还能用于质控处理。

2.1 FastQC的安装和使用

FastQC简介[1]

FastQC是一个Java程序。既有图形操作界面（graphical user interface，GUI），也可以使用命令行操作。图形界面的高通量数据分析软件Galaxy[3]和Chipster[4]也将其整合其中。

FastQC运行速度比较快。支持的输入文件格式有[5]：

• FastQ (所有质量编码系统的fasqc)

• Casava FastQ文件

• Colorspace FastQ（SOLiD平台的测序数据）[6]

• GZip压缩的FastQ

• SAM

• BAM

• SAM/BAM Mapped only (normally used for colorspace data)

FastQC安装

以下安装流程仅供参考，读者可根据自己情况选择别的安装方式。

# 进入你软件源码存放的文件夹。
cd ~/software
# 下载软件
curl -O http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
# 或使用wget：wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
# 用Linux的: 用sudo apt-get 安装和解压。
unzip fastqc_v0.11.5.zip
cd FastQC
ls -l
# 使fastqc执行。
chmod +x fastqc
# 把可执行的fastqc文件链接到入~/bin目录下
ln -s ~/software/FastQC/fastqc ~/bin/fastqc
# 或者把可执行的fastqc文件直接复制到到入~/bin目录下
# cp ~/software/FastQC/fastqc ~/bin/fastqc
# 测试软件是否安装成功
fastqc -h
# 出现帮助文档则为安装成功。

在下边这一步中，需要注意的是，你的~/bin目录的路径被添加到你的~/.profile (MAC)或者~/.bashrc (Linux)文件中。

# 把可执行的fastqc文件放入~/bin目录下
ln -s ~/software/FastQC/fastqc ~/bin/fastqc

如果没有添加，可以按下边命令进行添加

# 将~/bin 文件夹加到PATH
# Mac下:
echo 'export PATH=~/bin:$PATH' >> ~/.profile
# 使设置生效
source ~/.profile

# 在Linux 下这么做:
echo 'export PATH=~/bin:$PATH' >> ~/.bashrc
source ~/.bashrc

添加完了之后再按照上边的命令进行软件安装。

此外，也可以直接含有可执行命令的目录~/software/FastQC/添加到环境变量PATH中

当然，关于软件安装还有许多其他的方法，比如使用conda等工具，这里不赘述。

关于如何使用conda，请参考洲更同学的

• 生信软件安装神器

• 生信软件的好帮手-bioconda

• 什么，你嫌bioconda下载速度太慢？

FastQC命令使用

使用下边的命令对reads.fastq.gz这个压缩格式的fastq文件进行质量检测，并且对每个测序文件生成两个结果文件: html格式的报告文件以及一个zip压缩文件。

# 使用prefetch从NCBI下载测试数据集。下载的数据会放在目录~/ncbi/public/sra/
prefetch SRR1553607
# 创建文件夹并把sra文件拆分为两个测序文件（双端测序）
mkdir QC
cd QC
fastq-dump --split-files SRR1553607
# 获得2个数据文件 SRR1553607_1.fastq 和 SRR1553607_2.fastq
ls # 查看文件
# 使用fastqc进行质量检查
fastqc *fastq
# 此时会在当前文件夹生成结果文件。

# 假如你想把所有的结果文件统一放到一个文件夹：
mkdir results
fastqc -o results *fastq

可以从百度云下载该测试数据：链接: http://pan.baidu.com/s/1o8r8Lz8  密码:xm9v

关于从NCBI上下载数据，更多信息请参考我们公众号的文章《从NCBI下载测序数据 | 也许是目前最详细的版本》 。

FastQC给出的报告结果包括以下几个内容：

• 基本的统计信息（Basic Statistics）

• 每个位置的碱基的质量情况（Per base sequence quality）

• reads每个碱基的平均质量的频率分布图（Per sequence quality scores）

• 每个位置上的碱基的比例分布（Per base sequence content）

• read的GC含量的频率分布图（Per sequence GC content）

• 每个位置的碱基的N（测序时，无法判断是哪个碱基，则记为N）含量（Per base N content）

• 序列长度频率分布（Sequence Length Distribution）

• read重复的频率分布（Sequence Duplication Levels）

• Overrepresented sequences

• 接头含量（Adapter Content)

• Kmer含量（Kmer Content）

并且给出一个判断：通过（pass）、警告（warn）和不通过（fail）。不过，软件设置的判断阈值未必适合你的数据。

所以即便有一些项目出现警告或者不通过也不用过于担心。

比如，RNAseq的数据在序列重复水平（Sequence Duplication Levels）这一项不通过就很正常。

而RNAseq数据由于使用了随机六聚体引物，在reads开头会出现碱基成分偏离，所以Per base sequence content一项也会出现警告或者不通过的信号，而且在5’端会产生比较多的K-mer。这是一个真正的技术偏差,也无法通过修剪reads来纠正,不过在大多数情况下似乎不会对下游分析产生负面影响[15]。

更多fastqc结果文件的解读可以参见微信公众号“生信菜鸟团”的文章《测序数据质量控制之FastQC》或者FastQC官方帮助文档[14]。这里不再赘述。

2.2 PRINSEQ的安装和使用

PRINSEQ有网页版本, 也有命令行版本。Chipster GUI中也有PRINSEQ。

PRINSEQ输出的报告内容包括[1]：

• reads数目以及长度频率分布

• 碱基质量频率分布

• 序列复杂度

• GC 含量

• Ns的含量

• polyA/T 尾巴

• duplicates

• 接头

PRINSEQ可以处理的文件类型为非压缩的FASTQ、FASTA以及QUAL[1]。不支持压缩格式的fastq。

安装

# 进入你软件源码存放的文件夹。
cd ~/software
# 下载软件
curl -OL http://downloads.sourceforge.net/project/prinseq/standalone/prinseq-lite-0.20.4.tar.gz
# 解压
tar zxvf prinseq-lite-0.20.4.tar.gz
# 进入文件夹
cd prinseq-lite-0.20.4
ls # 查看目录里的内容
# 把该目录下的三个perl脚本(prinseq-lite.pl、prinseq-graphs.pl和prinseq-graphs-noPCA.pl )变得可执行。
chmod +x *.pl

# 把目录 ~/software/prinseq-lite-0.20.4/ 添加到环境变量

# Mac下:
echo 'export PATH=$PATH:~/software/prinseq-lite-0.20.4/' >> ~/.profile
# 使设置生效
source ~/.profile

# 在Linux 下这么做:
echo 'export PATH=$PATH:~/software/prinseq-lite-0.20.4/' >> ~/.bashrc
source ~/.bashrc

# 测试软件是否安装成功
prinseq-lite.pl -h
# 出现帮助文档则为安装成功。

通过上面的安装，prinseq-lite.pl是可以正常使用的。但是由于画图还需要别的perl模块，所以需要进行模块安装才能正常使用。

非root权限下的perl模块安装：

查看软件目录下的README文档可以知道：

• prinseq-lite.pl会用到的模块包括：

• Getopt::Long

• Pod::Usage

• File::Temp

• Fcntl

• Digest::MD5

• Cwd

• List::Util

prinseq-lite.pl无需安装perl模块。

• prinseq-graphs.pl会用到的模块包括：

• Getopt::Long

• Pod::Usage

• File::Temp

• Fcntl

• Cwd

• JSON

• Cairo

• Statistics::PCA

• MIME::Base64
  由于有的平台上Statistics::PCA模块会有一些问题，所以作者写了prinseq-graphs-noPCA.pl作为替代。

• 假如报告需要为网页版，需要用到以下模块:

• CGI

• File::Path

• IO::Uncompress::AnyUncompress

• LWP::Simple

• File::Copy

• File::Basename

非root用户运行下面的代码获取自己的私人cpan下载器cpanm。

wget -O- http://cpanmin.us | perl - -l ~/perl5 App::cpanminus local::lib
# 不知为何，在这一步我需要运行两次才能安装好。
eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`
echo 'eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`' >> ~/.bashrc
echo 'export MANPATH=$HOME/perl5/man:$MANPATH' >>  ~/.bashrc
echo 'export PATH=$PATH:~/perl5/bin' >>  ~/.bashrc
source ~/.bashrc

# ~/.bashrc可能需要更改为.bash_profile, .profile, etc等等，取决于你的linux系统！ 
# 然后你直接运行cpanm Module::Name，就跟root用户一样的可以下载模块啦！

注：安装cpanm的代码根据perl模块安装大全略有修改。

用cpanm基本能安装上大部分的模块。但是prinseq需要的模块，有的安装起来真心心累。比如在安装Cairo的时候，需要系统安装了cairo这个画图的工具库。而cairo画图工具库的安装需要依赖另外两个库：libpng和Pixman。没有root权限，安装很难。而且看质量情况一般用FastQC也比较足够，所以果断弃坑。

命令使用

PRINSEQ的质控报告可以是文本格式或者html格式，如果想输出html格式，需要把前面提到的模块安装好。

输入html格式需要分为两步[1]：

• 第一步：生成临时的图文件

  prinseq-lite.pl -fastq SRR1553607_1.fastq -out_good null -out_bad null -graph_data graph

  如果是质量编码体系是Illumina 1.8之前的版本，需要多加一个参数-phred64：

  prinseq-lite.pl -fastq SRR1553607_1.fastq -phred64 -out_good null -out_bad null -graph_data graph

  由于没有处理我们的数据，所以没有需要被接受或者被扔掉的reads，我们把这些输出文件(-out_ good and -out_bad) 设为null。

• 第二步：生成html格式或是一个图一张的png。

# 一张张的png图格式
prinseq-graphs.pl -i graph -png_all -o QCreport
# html格式
prinseq-graphs.pl -i graph -html_all -o QCreport

3 reads数据的质控

这部分内容主要来自：RNAseq Data Analysis: A Pratical Approach

3.1 碱基质量

碱基质量一般用Phred值表示。碱基错误的可能性，用log~10~进行变换，再乘以-10。举例来说，某个碱基有1%的概率是错误的，那么它的质量值为q = −10*log~10~(0.01) = 20。质量值的范围一般是0到40。在fastq文件里，为了节省存储，质量值是用ASCII字符表示的。目前的fastq文件都用的Sanger编码方式：第33个ASCCII字符表示质量为0。但是在1.8版本之前的fastq文件是以第64个ASCCII字符表示质量为0。如果你不知道你的数据是哪种编码系统，用FastQC可以检测出来，如图。而Trimmonatic软件可以帮你把fastq的质量编码方式进行转换[1]。

关于测序质量编码方式的讲解，可以查看《小白生信学习记3》。

• 通常一条read，后面的测序质量会变差，如下图fastqc的检测结果。RRINSEQ也有该项检测功能[1]。

  
  

• 此外，从上图也可以看出来，双端测序时，左端的测序质量会比右端的测序质量更好（原因我没有查到，有知道的朋友望能不吝赐教，感谢！）。

• 查看reads的平均质量的频率分布图可以更好体现总体的测序质量。理想的情况是，大部分reads的平均碱基质量高于25[1]。

  
  

3.2 低质量碱基的reads的处理方式

含有低质量碱基的reads的处理方式包括过滤和切除。

如果你决定要过滤或者切除双端测序的reads，则需要选择一个能保持双端序列顺序对应的工具。因为当把reads mapping到参考序列时，比对软件会以相同的顺序在两个测序文件里寻找配对的reads。注意，不单是过滤，切除也可能会对顺序造成影响：当一些reads被完全切除或者切除得太短而被过滤掉，从而导致顺序的变化[1]。

3.2.1. 过滤：去除整条reads

Trimmomatic, FastX,以及PRINSEQ都可以通过质量过滤reads。

• FastX: 允许用户设定一个最低质量值以及有多少比例的碱基需要高于这个最低值。

• PRINSEQ and Trimmomatic :

• 通过read的碱基平均质量来过滤

• 支持双端文件处理

Trimmomatic的安装

安装方法参考了《Biostar课程7、8-二代测序质控：fastqc、prinseq、trimmomatic、seqtk》

# 安装trimmomatic
cd ~/software
curl -O http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-Src-0.36.zip
unzip Trimmomatic-Src-0.36.zip
cd Trimmomatic-0.36

# 可以到这查看一下手册：http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf

# 运行的命令比较长
java -jar trimmomatic-0.36.jar

# 我们写一个脚本来替代我们运行。
# 你也可以通过命令行来达成
echo '#!/bin/bash' > ~/bin/trimmomatic
echo 'java -jar ~/software/Trimmomatic-0.36/trimmomatic-0.36.jar $@' >> ~/bin/trimmomatic
chmod +x ~/bin/trimmomatic

# 安装成功。
trimmomatic

# 也可以使用conda进行安装。
conda install trimmomatic

Trimmomatic[1]

Trimmomatic是一个多功能的基于JAVA的reads处理工具。Trimmomatic既有命令版本，也被Galaxy 和 Chipster图形界面软件整合了。Trimmomatic的输入和输出文件都是FASTQ（压缩或者非压缩的都可以）Trimmomatic是多线程的，所以运行起来很快。

Trimmomatic的功能：

• 可以以不同方式，基于质量去除接头以及剪切reads

• 可以基于reads质量和长度过滤reads

• 把碱基质量值从一种编码体系转换成另一种

Trimmomatic可以一条命令，按着你给的顺序把多个步骤完成。

下边这条命令是用Trimmomatic过滤双端reads，设定read的碱基平均质量阈值为20 (AVGQUAL:20)，输入和输出文件都为压缩的[1]。

trimmomatic PE -phred33 SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20
# 如果没有按照上边那样写一个脚本进行命令名称替代，则按下边写:
java -jar trimmomatic-0.36.jar -phred33 SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20

-phred33 表质量编码体系是示Illumina 1.8之后的版本。质量编码体系是Illumina 1.8之前的版本，需要改为-phred64。不过，从0.32版本，软件会自动识别是33还是64。所以 -phred参数可以不加：

trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20

用PRINSEQ进行同样的质量过滤[1] ：

prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -min_qual_mean 20 -out_good qual_filtered -out_bad null –no_qual_header –log -verbose

# 同样，如果质量编码体系是Illumina 1.8之前的版本需要加上-phred64：
prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -phred64 -min_qual_mean 20 -out_good qual_filtered -out_bad null –no_qual_header –log -verbose

• 过滤剩余的reads在文件 qual_filtered_1.fastq和 qual_filtered_2.fastq里

• 过滤掉的则在 qual_filtered_1_singletons.fastq 和 qual_filtered_2_singletons.fastq 里

• –verbose 使得可以再运行完之后查看数据统计。

• –no_qual_header 告诉PRINSEQ只写“ + ”，而不是写“+ read 名称”作为每个read的质量行的header，从而减少fastq文件的大小。

3.2.2. 切除：从read一端切除低质量碱基

切除方案1：设定切除长度或者是切剩余的长度。

假如在reads末端发现低质量碱基，最简单的方法就是从某端开始切除，并设定切除reads的剩余长度或者给一个需要切除的长度。

但是这个方法会把质量好的那部分序列也去除了。如何减少高质量序列的损失？

切除方案2：切除低于质量阈值的碱基（从3’或者5’端开始，逐个切除低于你设定质量阈值的碱基。）

• FastX 只能从3’端开始切除

• PRINSEQ 和 Trimmomatic 可以从任意一端开始切除。

• PRINSEQ, Trimmomatic和Cutadapt支持双端文件，当reads被切到低于用户定义的最小长度而被过滤时，可以保证reads配对的一致性。

用Trimmomatic从3’端开始切除低质量碱基（TRAILING表示3’端开始，如果想从5’开始，则改为LEADING），碱基质量阈值为20 （TRAILING:20），过滤切除后长度小于50个碱基的reads（MINLEN:50）。切除和过滤参数需要按正确顺序排。

trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz trim_paired1.fq.gz trim_unpaired1.fq.gz trim_paired2.fq.gz trim_unpaired2.fq.gz TRAILING:20 MINLEN:50

切除方案3：使用滑窗的方法（设定一个滑窗，将窗口的碱基质量与设定的阈值进行比较。）

• Trimmomatic从reads的5’端开始滑窗

• PRINSEQ允许用户选择从哪个端开始

注意：

• 从5’端开始滑窗：是从5’端开始保留高质量的序列，直到到达质量值低于阈值的地方停止。

• 从3’端开始滑窗：是从3’端开始切除低质量的序列，直到到达质量值高于阈值的地方停止。

由于有的reads有可能会出现中间质量很低，而两端比较高的情况。而测序质量一般都是从5’高而3’低，所以使用5’端开始滑窗的方法通常会导致质控后产生的reads更短。导致低于read长度阈值的序列更多，使得被过滤掉reads也随之更多。而这个问题对于双端测序影响会更大，因为其中一端测序的read被过滤，在另一端测序文件里对应配对的序列也会被去除。

除了可以选择reads扫描方向，PRINSEQ还有其他灵活的设置：Trimmomatic允许用户设置窗口大小，但只能使用窗口内碱基质量的平均值作为窗口的质量值，而PRINSEQ还允许设定滑窗的步长、允许选择使用平均值或最小值来跟阈值作比较。

下边的命令是使用Trimmomatic软件，设置窗口大小为3碱基，从5’开始切除reads，read的平均质量低于20（SLIDINGWINDOW:3:20），开始切除后边的碱基序列。过滤长度小于50碱基的reads（MINLEN:50）。使用更大的窗口可以减轻切除的力度。

trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz win_paired1.fq.gz win_unpaired1.fq.gz win_paired2.fq.gz win_unpaired2.fq.gz SLIDINGWINDOW:3:20 MINLEN:50

下边的命令是使用PRINSEQ软件，设置窗口大小为3碱基，从3’开始切除平均质量低于20的序列（高于20时停止）。同样过滤长度小于50碱基的reads。
保留的reads放在window_1.fastq和window_2.fastq，因配对read被去除而同时被去除的reads放在window_1_singletons.fastq和window_2_singletons.fastq。

prinseq-lite.pl -trim_qual_window 3 -trim_qual_type mean -trim_qual_right 20 -trim_qual_rule lt -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -out_good window -out_bad null -verbose -min_len 50 -no_qual_header

切除方案4：另一种滑窗的方法是“running sum”，最初是在BWA比对软件中采用，所以通常也被称为“BWA quality trimming” 。

从右端（3’端）开始扫描reads，将每个碱基的质量值与给定阈值进行比较，并计算累计差异的加和。在累计“badness”最高的read位置切除。这个方法被Cutadapt软件采用。

最后，Trimmomatic提供了一种适应性的切除方法，叫MAXINFO。其目标是在尽可能在保留reads长度和去除错误碱基序列之间取得一个平衡。

它包括两个参数：

• 目标长度（targetLength）

• 严格度（strictness）

可以通过严格度参数来控制这个平衡，严格度是从0到1的值，数值越高越严格。参数值低 (<0.2)有利于保留更长的reads，而高(>0.8) 则保证reads的准确性。

MAXINFO过滤的原理参见：

• Trimmomatic的使用手册《Trimmomatic Manual: V0.32》：http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf

• 这里有篇中文版的：《NGS 数据过滤之 Trimmomatic 详细说明》：http://www.jianshu.com/p/a8935adebaae

下边的MAXINFO切除命令是设置目标长度为50，严格度为0.7[1]。

trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz MI_paired1.fq.gz MI_unpaired1.fq.gz MI_paired2.fq.gz MI_unpaired2.fq.gz MAXINFO:50:0.7 MINLEN:50

3.2.3 不确定碱基（Ambiguous Bases）

在测序时，如果一个碱基无法判定是哪种，则记为N。不同的组装软件和比对软件对于这种不确定碱基的处理方式不一样：有的把N用4种碱基进行随机替代，有的则是用固定某个碱基替代N。由于Ns会导致错误的组装和mapping，含N多的reads需要被去掉。

PRINSEQ的质控报告会给出一个Ns含量的图，可以看reads含N比例的一个频率分布情况。

使用PRINSEQ过滤含Ns过多的reads：

prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -ns_max_n 2 -out_good nfiltered -out_bad null -no_qual_header -log -verbose

-ns_max_n 2表示过滤含超过2个N的reads。可以用-ns_max_p来过滤含超过设定比例Ns的reads。

保留的reads放在nfiltered_1.fastq 和 nfiltered_2.fastq中，质控中丢失配对read的reads放在nfiltered_1_singletons.fastq和nfiltered_2 _singletons.fastq

3.2.4 接头（Adapters）

测序时候，一般会使用测序接头序列，而接头序列的存在可能会导致基因组组装和转录本组装出现问题，所以需要在分析数据之前去除掉。此外，其他的tags，比如multiplexing identifiers以及引物也需要被去除[1,2]。

去接头是很有挑战性的[1]：

• 首先，这些tags会有测序错误，所以切除软件需要能应对错配、indels以及不确定碱基（N）

• 其次，在对小RNA进行测序时，reads可能会延伸到3’端的接头序列。而这种“读穿”的情况会导致read中含有的3’端的接头序列只是一部分序列而无法被识别出来。

• 最后，假如数据是从公共数据库获取的，可能根本无法知道接头序列信息。

fastqc支持你把接头序列加到adapter文件中，进行接头分析。

查看接头序列。

cat ~/software/FastQC/Configuration/adapter_list.txt

假如接头序列未知，可以用TagCleaner软件预测。而FastQC的k-mer overrepresentation图和PRINSEQ的Tag Sequence Check图也可以检测出是否有接头。

Trimmomatic, FastX, TagCleaner以及Cutadapt都可以去除接头序列。

其中Trimmomatic, TagCleaner和 Cutadapt：

• 可以应对错配

• 两端都可以切除接头

• 允许用户指定read和tag序列的最小重叠

• TagCleaner还可以应对indels以及不确定碱基

• Trimmomatic运行很快，因为它先用短的种子序列（seed sequences）扫描reads并只对那些跟种子序列能很好匹配的reads进行完全比对。

• Trimmomatic支持双端测序文件，并且可以利用双端信息来检测接头。这种所谓的回文方法是基于“读穿”会在两端的测序中都发生，而插入片段会被完全测通。所以reads才能被对上，允许只测到一部分的接头序列被检测到，即便少到只有1个碱基。

注意，如果你要把general、quality以及adapter trimming 同时放到Trimmomatic的一个命令中，你要把切除接头放到第一步，因为识别完整的接头序列比识别部分接头序列要容易得多。

我们也可以自己创建接头序列文件[2]，或使用Trimmomatic自带的：

创建接头序列文件:

echo ">illumina" > adapter.fa
echo "AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC" >> adapter.fa

使用trimmomatic去除另一个sra文件SRR519926的接头:

# 单端去除接头
trimmomatic SE SRR519926_1.fastq output.fq ILLUMINACLIP:adapter.fa:2:30:5
# 查看质量情况
fastqc SRR519926_1.fastq
fastqc output.fq

去除接头后的reads存放在output.fq里。

左边的图示去除接头之前，右边的图是去除接头之后。

使用bbduk进行切除接头:

bbduk.sh in=SRR519926_1.fastq  ref=adapter.fa out=bbduk.fq

大部分工具都允许多个步骤放到同一条命令运行。

# 双端模式
trimmomatic PE SRR519926_1.fastq SRR519926_2.fastq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:adapter.fa:2:30:5

先后顺序一般来说很重要，因为后面的处理依赖于前面的结果。

先以质量值切除，再切adapter，还是反过来更好？由于adapters有重叠区，如果我们先以质量切除，那么有可能因为切除了部分碱基而导致无法识别出adapter，另一方面，如果碱基质量很差，也可能无法比对上adapter而导致adapter无法被检测出来。通常，这个顺序不会导致特别大的差异。不过，对于有问题的数据，最好还是先试验一下[2]。

3.2.5 随机六聚体引物导致的sequence-Specific Bias以及错配

在文库构建时，RNA是片段化的，随机引物六聚体会被作为引物来反转录产生cDNA，然后再从末端进行测序。所以，理论上说，我们期待reads从转录本的随机位置开始，因此reads中不应该会有碱基组成偏离的情况。但是，实际情况表明，随机六聚体会在RNAseq的reads开头引入碱基组成偏离[7]。这种序列特异的偏差会影响基因以及转录本表达量的估算，因为测序reads对转录本的覆盖是不均一的。在不同碱基位置的严重偏差也可能是未切除接头或者是文库测序过度的信号[1]。

FastQC会给一个reads每个位置的碱基组成图（如下图），理论上应该每种碱基都是一条直线，该值应近似等于物种中该碱基的含量。在任一碱基位置上，如果A和T的差异或者C和G的差异大于10%，则FastQC会显示警告[1]。

序列特异的偏离无法通过过滤或者切除来去除，不过 Cuffinks 和 eXpress包提供了一个校正方法，这个方法学习从数据中挑选出来的序列并把这些信息包含到丰度估计中[1]。

如果只是前面几个碱基出现碱基成分偏差，也可以直接切除。

注意，除了造成序列特异性偏差，随机六聚体引物还会导致Illumina RNAseq reads开头的错配[8]。错配率最高的是第一个核苷酸，虽然有高达7个核苷酸会被影响。错配的碱基通常有很好的质量值，因此不同通过质量来检测它。为了避免这种情况，第一个碱基应当用Trimmomatic 或 PRINSEQ给切除。

3.2.6 GC 含量

reads的GC含量通常服从正态分布，中间的GC值应该近似为物种的GC含量。异常的分布性状或跟源基因组的GC含量有很大偏差，可能意味着文库污染。

FastQC和PRINSEQ都会将reads的平均GC含量的频率分布给出来。

注意，标准的文库构建流程会使用PCR扩增，而PCR对于GC含量过高或过低的区域扩增比较难。这些GC含量是样本特异的，从而导致差异表达分析的复杂度增加。GC偏差无法用reads预处理的方式消除。不过，在后续分析阶段，各种各样的标准化和校正方法会考虑这一点[9,10]。

3.2.7 重复序列（Duplicates）

duplicates有两种可能：

• 在样品中真实存在的，比如RNAseq数据

• 是制备文库等过程中人为引入的，比如PCR过度扩增

duplicates检测方法：

• 序列一致性：忽略掉了由于测序错误而序列不一致的那部分duplicates

• 比对一致性：只有在有参考序列的前提下才可行。

• k-mer法

是否该去除duplicates？

• 在有的高通量测序中，一致的序列的水平太高意味着存在PCR的过度扩增（比如基因组重测序）

• 对于RNAseq，duplicates常常是高表达转录本的一种自然结果。对于差异表达分析，不建议去除duplicates。

• 如果一个稀疏覆盖的转录本在一个位置上有堆积了特别多的reads，那么这很可能就是PCR导致的[1]。

• 在SNP calling和基因组变异检测时，应当去除duplicates。因为变异的检测是通过reads在变异发生的位置的覆盖度来推算的。所以，可能会由于duplicates的缘故而错误认为某个SNP或者InDel是可靠的[2]。

duplicates的检测和去除

FastQC 和 PRINSEQ 都能分析duplicates。FastQC只检测精确的duplicates，而PRINSEQ还会检测5’和3’端的duplicates。

• PRINSEQ的过滤功能允许用户设定一条reads允许有多少的duplicates。

• FastX Collapser会把一致的序列合为一条，并对有多少条一致序列进行计数。

这些工具用于原始reads（raw reads），基于序列识别duplicates。注意，实际中，PCR同样片段产生的duplicates的测序reads未必会完全一致，因为有测序错误。所以，基于序列的方法会低估了duplicates的数目。

当reads被比对到参考基因组，duplicates可以通过mapping位置的一致性（而非序列内容）来识别。

使用PRINSEQ去除原始reads中的duplicates。下边的命令是去除出现超过100次(-derep _min 101)

prinseq-lite.pl -fastq reads1.fastq -fastq2 reads2.fastq -derep 1 -derep_min 101 -log -verbose -out_good dupfiltered -out_bad null -no_qual_header

FastQC的duplicate图怎么看？[2,11]

看懂这个图其实有点难。首先，红色圈里的值很重要，她表示有多少比例的reads是不一样的。

有两种duplicate的计算方式：

• 对所有序列（all sequences）

• 对不一样的序列（distinct sequences）

大神Istvan Albert在Tutorial: Revisiting the FastQC read duplication report（https://www.biostars.org/p/107402/）对这作了更详细的讲解：

首先，定义什么是distinct sequences：

distinct sequences = singletons（重复值为1的序列）的数目 + doubles（重复值为2的序列）的数目 + triplets（重复值为3的序列）的数目+ …

这里，每种序列不论重复了几次，都只算一种distinct sequences。

图标题上的值 = distinct/总序列数 * 100%

这就是蓝线和红线想要展示的：

• 蓝线表示在给定的重复值（X轴），所有发生duplicate的序列的比例（Y轴）。这个比例是对于总reads数来计算的。

• 红线表示在给定的重复值，发生了duplicate的distinct sequences的比例。这个比例是对于总distinct sequences数来计算的。

举例来说，有两种情况，都含有20条reads。
Case 1: 10 条unique reads +  5 种重复2次的reads

A B C D E F G H I J K K L L M M N N O O

Case 2: 10 unique reads + 1种重复了10次的read

A B C D E F G H I J K K K K K K K K K K

Case 1 （在上面的图）展示了15种distinct reads，即比例为15/20=75%

对于all sequences的计算方法（蓝线）：

总数为20的总reads数中：

• singletons数目为1x10 =10，比例为10/20 = 50%

• doubles数目为5x2 =10，比例为10/20 = 50%

所以蓝线是一个平台，因为二者比例都为50%。

对于distinct sequences的计算方法（红线）：

对于总共为15种distinct sequences有：

• 10个singletons，比例为 10/15=66%

• 5个重复两次的重复序列（duplicates），比例为 5/15=33%

所以红线是一条斜线。

Case 2 （在下面的图）则有11种distinct reads，所以比例为11/20=55%。

对于all sequences的计算方法（蓝线）：

总数为20的总reads数中：

• singletons数目为1x10 =10，比例为10/20 = 50%

• 重复值为10的read类型的read数目为1x10 =10，比例为10/20 = 9%

所以，蓝线在重复值为1和10的地方，比例都为50%。

对于distinct sequences的计算方法（红线）：

对于总共为11种distinct sequences有：

• singletons类型为10种，比例为10/11 = 91%

• 重复值为10的read类型的read类型为1种，比例为1/11 = 9%

所以红线在重复值为1时，比例为91%，而在重复值为10时，比例为9%。

3.2.8 序列污染（Sequence Contamination）

有的时候reads可能会被别的生物体或者载体的序列污染。这会在GC含量中体现。PRINSEQ的dinucleotide frequency 图也可以给出一些样品污染的线索。但也许最直接的方法就是把reads比对到可能的污染源序列上。FastQ Screen工具可以把reads mapping到用户怀疑的污染源序列上，用的Bowtie比对工具，并且会把结果总结成文本和图形的形式[12]。或者，你可以随机抽取一部分reads，使用BLAST对常规核苷酸数据库进行比对[1]。

3.2.9 低复杂度序列和polyA尾巴（Low-Complexity Sequences and PolyA Tails）

低复杂度的序列

低复杂度的序列由于缺乏信息量，所以很难可靠地map到参考序列上。比如，单一（homopolymer）重复序列 比如：AAAAAAAAAA，二核苷酸（dinucleotide）重复序列如： CACACACACA，三核苷酸（ trinucleotide）重复序列如：CATCATCATCAT。

PRINSEQ通过两种方法计算reads的序列复杂度：DUST 和 Entropy。

• DUST分值从0到100，单一重复序列拥有最高的分值。一条read，假如DUST分值高于7，则被认为是低复杂度。

• Entropy打分则相反，单一重复序列分值最低，为0，而分值高于50则认为复杂度低。

你可以通过选项–lc-threshold来过滤低质量的reads，用–lc-method (dust/entropy)来选择使用哪种方法。

PolyA/T尾巴

PolyA/T尾巴是reads末端的A或T重复序列。

PRINSEQ会报告有多少reads含有大于等于5个的碱基长度的polyA/T尾巴，以及尾巴长度的频率分布。你可以通过设置最小尾巴长度选项–trim_tail_right (或者 –trim_tail_left) 来去除尾巴。

4 Reads过滤和切除软件推荐

Reads过滤和切除软件那么多，应该如何选择？每个人都有自己的偏好。这里笔者选了几个好友以及书籍《Biostar Handbook》和《RNAseq Data Analysis: A Practical Approach》的推荐罗列如下供大家参考：

4.1 健明兄

• trim_galore：该软件需要调用FastQC和cutadapt ，需要提前装好。

4.2 天地本无心

三种软件组合使用：

• cutadapt

• lengthsort (SolexaQA)

• DynamicTrim (SolexaQA)

4.3 Viswanathan Satheesh

• BBmap里的BBDuk

4.4 Biostar Handbook (January 2017)[2]

在Biostar Handbook里罗列了许多QC软件，如下：

• BBDuk： BBMap里的一个包

• BioPieces：序列处理的一个套件（含多个软件）

• CutAdapt：2011年发表在Embnet Journal的application note部分。

• fastq-mcf： 2013年发表在The Open Bioinformatics Journal上

• Fastx Toolkit: 短reads FASTA/FASTQ文件处理的命令行软件集合

• FlexBar：灵活的barcode和接头去除软件。2012年发表在Biology上

• NGS Toolkit：2012年发表在Plos One上

• PrinSeq： 2011年发表在Bioinformatics的application note部分。

• Scythe：一个贝叶斯接头切除软件

• SeqPrep - 一个可以去除接头以及/或者把有重叠区域的配对reads合并为一条read

• Skewer: 一个能快速且准确地处理NGS双端reads的接头去除软件

• TagCleaner：2010年发表在BMC Bioinformatics上

• TagDust： 2009年发表在Bioinformatics上

• Trim Galore - 基于Cutadapt和FastQC的一个工具。能对FastQ文件进行质控和接头去除。同时，还能用于MspI-酶切 RRBS-型的（简化的重亚硫酸盐测序（Bisufite-Seq）） 文库

• Trimmomatic：2012年发表在Nucleic Acid Research的application note上

用于QC的R（Bioconductor）包: PIQA、ShortRead

其中，作者推荐使用以下几种：

• Trimmomatic

• BBDuk

• flexbar

• cutadapt

4.5 RNAseq Data Analysis: A Practical Approach (June 2014)[1]

在这本书里，作者主要介绍了如下几个软件：

• PRINSEQ

• Trimmomatic

• Cutadapt

• FASTX-toolkit

最后

质控软件太多，这里不可能一一予以介绍。以后如果有时间和经历可能会把软件安装和使用继续写写。

注：本文大部分内容参考了《RNAseq Data Analysis: A Pratical Approach》一书。

参考资料

1. RNAseq Data Analysis: A Pratical Approach

2. István Albert. 2017. Biostar Handbook

3. Goecks, J., Nekrutenko, A., and Taylor, J. Galaxy: A comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol, 11(8):R86, 2010.

4. Kallio, M.A., Tuimala, J.T., Hupponen, T. et al. Chipster: User-friendly analy- sis so ware for microarray and other high-throughput data. BMC Genomics, 12:507, 2011.

5. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/2%20Basic%20Operations/2.1%20Opening%20a%20sequence%20file.html

6. https://en.wikipedia.org/wiki/FASTQ_format

7. Hansen, K.D., Brenner, S.E., and Dudoit, S. Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res, 38(12):e131, 2010.

8. van Gurp, T.P., McIntyre, L.M., and Verhoeven, K.J. Consistent errors in  rst strand cDNA due to random hexamer mispriming. PLoS ONE, 8(12):e85583, 2013.

9. Roberts, A., Trapnell, C., Donaghey, J., Rinn, J.L., and Pachter, L. Improving RNA-seq expression estimates by correcting for fragment bias. Genome Biol, 12(3):R22, 2011.

10. Roberts, A., Trapnell, C., Donaghey, J., Rinn, J.L., and Pachter, L. Improving RNA-seq expression estimates by correcting for fragment bias. Genome Biol, 12(3):R22, 2011.

11. Revisiting the FastQC read duplication report

12. FastQ Screen. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/.

13. MacManes, M.D. On the optimal trimming of high-throughput mRNA sequence data. Front Genet, 5:13, 2014.

14. FastQC官方帮助文档(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)