转录组数据分析笔记（1）——如何用fastqc和trim-galore做测序数据质控

从文章末尾我们可以获得一些测序数据信息：

Data availability. ChIP-seq and RNA-seq data have been deposited with ArrayExpress database (www.ebi.ac.uk/arrayexpress), accession numbers E-MTAB-4680, E-MTAB-4684 and E-MTAB-5130.

可以看到作者将CHIP-seq和RNA-seq数据上传到ArrayExpress这个数据库中，这个数据库是欧洲生物信息研究所（European Bioinformatics Institute, EBI）旗下的公共数据库，主要用于存放芯片和高通量测序数据，我们可以直接从该数据库中下载我们需要的RNA-seq数据，自己动手分析。

EBI数据库可以直接下载fastq数据，不需要做SRA数据转换（NCBI数据库中下载sra数据则需要转换，需要用工具fastq-dump），这是EBI数据库下载高通量测序数据的优点，但是这个数据库经常网络连接不稳定，用aspera或者prefetch这种高速下载软件也不一定能稳定下载 最好的方法是科学上网。我们可以从ArrayExpress数据库中输入索引号E-MTAB-5130，直接获得样本信息和测序信息。

我们可以看到，下载的数据是双端测序产生的。我们不能直接用下载的raw data做后续分析，必须要进行质控查看测序质量如何。

2.1 使用fastqc对测序数据生成质控报告

下载好的fastq文件可以直接用fastqc工具做测序数据质控，输入以下命令一次生成所有qc报告：

$ fastqc *.fastq.gz -o ./ #在当前目录下对所有.fastq.gz文件生成qc报告，-o参数定义输出目录

运行结束后我们可以得到.html文件和.zip压缩包，这个就是质控报告。在虚拟机里，我们可以直接点开.html后缀的网页文件查看质控报告（和压缩包的内容是一致的）。

顺便介绍一个非常好用的工具multiqc，可以通过conda install直接安装，这个工具可以将批量生成的qc报告合并为一个，看起来更加直观。在生成qc报告的当前目录下，运行代码：

$ multiqc ./

2.2 质控报告解读

2.2.1 基本信息

绿色表示通过，黄色表示不太好，红色表示不通过。RNA-seq一般在Sequence Duplication Levels上结果会不好，一个基因可能会大量表达，测到好多遍。

2.2.2 核苷酸测序质量箱式图

这里测序质量（纵坐标）用Q值表示，p为出错率，Q值计算式为Q=-10*lg（p）。每一个核苷酸的测序质量可以从fastq文件第四行一一对应上，这里只是做了统计和可视化。我们可以看到每个位点的核苷酸测序质量Q值都在30以上，意味着每个位点的测序正确率都在99.9%以上，可以认为测序质量比较好。

箱式图解读：黄色箱子(25%和75%的分数线)，红色线(中位数)，蓝线是平均数，下面和上面的触须分别表示10%和90%的点。

2.2.3 测序泳道质量图

纵坐标为tile编号，这张图代表每次荧光扫描的质量。蓝色背景表明测序质量良好，白色和红色的背景表示测序过程中可能有小气泡或者测序泳道上有污染。直接的体现就是部分测序数据中出现连续的N，也就是不能读取，可能是任何一个核苷酸。

2.2.4 reads质量得分

可以看到平均质量在38，质量比较高。如果最高峰所对应的横坐标质量值小于27（错误率0.2%） 则会显示“警告”，如果最高峰的质量值小于20（错误率1%）则会显示“不合格”。

2.2.5 每条reads各个测序位点上各碱基出现概率

图上看得出比较稳定，测序刚开始的时候波动会大一点，这里的GC含量和AT含量不一致。如果任何一个位置上的A和T之间或者G和C之间的比例相差10%以上则报“警告”，任何一个位置上的A和T之间或者G和C之间的比例相差20%以上则报“不合格”。

2.2.6 GC含量和理论分布

可以看出GC含量在43%左右，与理论分布（也就是正态分布）比较吻合，中心峰值与所测转录组的GC含量一致。如果有不正常的尖峰，可能是测序文库有污染，接头的污染还会在过表达序列中体现。

2.2.7 每条reads的含N碱基数

不能识别的碱基会被读成N，这里没有N，测序质量非常好。横坐标表示reads的位置，纵坐标表示N的比例。如果任何一个位置N的比例大于5%则报“警告”，大于20%则报“失败”。

2.2.8 测序长度分布

这个测序仪一次测量长度是101bp。测序仪出来的原始reads通常是均一长度的，经过质控软件处理过的reads长度则不一样，这里说明测序结果较好。

2.2.9 重复序列水平

可以看到重复水平较低。图中横轴代表reads的重复次数，大于10次重复后则按不同的重复次数合并显示。纵坐标表示各重复次数下的reads数占总reads的百分比；蓝线展示所有reads的重复情况，红线表示在去掉重复以后，原重复水平下的reads占去重后reads总数的百分比；如果非unique的reads占总reads数的20%以上则报 ”警告“，占总read数的50%以上则报 ”不合格“。这项变黄是正常的。

2.2.10 过表达序列和接头序列

过表达的序列很可能是一些测序的接头序列，这里两种序列都看不到，说明质量良好。过表达序列是显示同一条reads出现次数超过总测序reads数的0.1%的统计情况，超过0.1%则报“警告”，超过1%则报“不合格”，会列出可能的接头序列。接头序列正常情况下含量接近于0。

2.3 trim-galore测序数据质控过滤

质控的目的使为了除去下机数据raw data中的接头序列和质量比较差的测序数据，Q<20，正确率小于99%，如果这样的核苷酸超过read长度的20%，则考虑将该read丢弃（只是建议，不是强制，根据需要可以自定义过滤条件）。

trim-galore也可以用conda install安装，非常方便，这是一个自动检测adaptor的软件，可以一个命令自动找出主流的测序接头并去除，还可以设置参数对测序数据质控。简单介绍一下trim-galore的一些参数：

因为是双端测序，16个样本每个都有_1和_2两个文件，可以写个脚本批量运行：

1
2
3
4
5
6
7
8

#!/bin/bash
for i in `seq 194 209`
    do 
    trim_galore -q 25 -phred33 -length 50 -stringency 3 -paired \
    -o /media/sf_/example/data/clean_data \
    /media/sf_/example/data/raw_data/ERR1698"$i"_1.fastq.gz \   #一端测序数据
    /media/sf_/example/data/raw_data/ERR1698"$i"_2.fastq.gz     #另一端测序数据
done

保存退出，运行，最后生成的_triming_report.txt文件就是生成的质控报告，_val_1.fq.gz就是过滤后瘦身的clean data，我们可以看到大小比原来小了10M左右，这个clean data才可以用于后续的分析流程

我截取了其中一个数据的质控结果，拉到最底下，可以看到两端测序数据中都有AGATCGGAAGAGC这个序列，在一个样本测序数据中出现240027次经过网上查找，AGATCGGAAGAGC这个序列确实是Illumina公司测序时的接头序列（点击这里查看），可以和上面fastqc质控报告中的测序平台Illumina相互验证。