转录组学习笔记03-了解fastq测序数据

任务

需要用安装好的sratoolkit把sra文件转换为fastq格式的测序文件，并且用fastqc软件测试测序文件的质量！作业，理解测序reads，GC含量，质量值，接头，index，fastqc的全部报告，搜索中文教程，并发在论坛上面。来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

将sra转换为fastq格式

sra格式

下载的原始数据是sra格式，SRA（Sequence ReadArchive）数据库是用于存储二代测序的原始数据，包括 454，Illumina，SOLiD，IonTorrent，Helicos 和 CompleteGenomics。除了原始序列数据外，SRA现在也存在raw reads在参考基因的比对信息。

fastq格式：

Fastq格式是一种基于文本的存储生物序列和对应碱基（或氨基酸）质量的文件格式。最初由桑格研究所（Wellcome Trust Sanger Institute）开发出来，现已成为存储高通量测序数据的事实标准。

fastq格式说明

FASTQ文件中每个序列通常有四行：

序列标识以及相关的描述信息，以‘@’开头；后面跟着序列的唯一ID以及相关说明内容。
第二行是核酸序列，是有ATCGN字符组成
第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加
第四行，每个测序碱基质量，是用ASCII码来表示的，与第二行的字符数一致。碱基质量得分与错误率的换算关系： Q = -10log10p（p表示测序的错误率，Q表示碱基质量分数） ASCII值与碱基质量得分之间的关系： Phred64 Q=ASCII转换后的数值-64 Phred33 Q=ASCII转换后的数值-33

一个典型的fastq文件如下所示：

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

如何判断是Phred64 还是 Phred33 ？ ASCII值小于等于58（相应的质量得分小于等于25）对应的字符只有在Phred+33的编码中被使用，所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情况下，ASCII值大于等于74的字符只出现在Phred+64中。如果是最近两年的测序数据，一般都是Phred33形式的。参考文章：http://blog.csdn.net/huyongfeijoe/article/details/51613827

转换工具：fastq-dump

将sra格式转换为fastq格式的工具是sratoolkit的fastq-dump，可以先用fastq-dump -h看一下帮助文件，分为如下几个部分：

输入： -A|—accession 序列号
处理中: Read Splitting, Full Spot Filters, Common Filters, Filters based on alignments, Filters for individual reads。基本都是些过滤参数。不太常用
输出： -O|—outdir 输出文件夹， -Z|—stdout 输出到标准输出， —gzip/—bzip2 输出为压缩格式
多文件选项：常用的就是—split-3，split-3表示如果是单端测序则一个sra文件出来一个fastq文件，如果是双末端，则一个sra问件对应两个fastq问件SRRXXXXXX_1.fastq,SRRXXXXXX_2.fastq
格式化：分为序列，质量等，不常用
—gzip 使得输出的结果是.gz 的格式

现在将disk2/sra下的sra文件转化为fastq格式，其中输入目录为disk2/sra，输出目录为/disk2/data/rna-seq，命令如下所示：

1
2
3

for ((i=56;i<=62;i++))
do fastq-dump --gzip --split-3 -A ~/disk2/sra/SRR35899$i.sra -O ~/disk2/data/rna-seq
done

转换后的结果如下所示：

对fastqc格式进行质控

工具：Fastqc

用法如下所示：

fastqc [-o output dir] [—(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

参数： -o 输出目录，需自己创建目录

—(no)extract 是否解压输出文件，默认是自动解压缩zip文件。加上—noextract不解压文件。

-f 指定输入文件的类型，支持fastq|bam|sam三种格式的文件，默认自动识别。

-t 同时处理的文件数目。

-c 是contaminant 文件，会从中搜索overpresent 序列。

代码如下所示：

mkdir -p ~/disk2/data/QC
cd ~/disk2/data/rna-seq
# 将所有的数据进行质控，得到zip的压缩文件和html文件
fastqc -o ~/disk2/data/QC *.fastq.gz

质控结果查看：

1. 质控后的文件如下所示：

质控结果有14个html文件，可以选择用浏览器打开查看最终的QC reports。随便打开一个质控结果文件，如下所示：

2. 质控报告的内容

左边是目录概要，一共12项内容，如下所示：

Basic Statistics
Per base sequence quality
Per tile sequence quality
Per sequence quality scores
Per base sequence content
Per sequence GC content
Per base N content
Sequence Length Distribution
Sequence Duplication Levels
Overrepresented sequences
Adapter Content
Kmer Content

可以点击想要看的结果，右边会跳转到特定详细的可视化结果。绿色代表“通过”，黄色代表“警告”，红色代表“不通过，失败”，如下所示：

2.1 Basic Statistics

基本的数据统计包括文件名，文件类型，编码形式，总的序列数，质量差的序列，序列平均长度，GC含量。

2.2 Per base sequence quality

图片太大，下部图片如下所示：

上图表示的是每个read各位置碱基的测序质量，具体的理解如下所示：

横轴表示碱基的位置，数字就是表示碱基，即1到51个碱基。
纵轴是质量分数，Quality score=-10log10p（p代表错误率），所以当质量分数为40的时候，p就是0.0001，质量算高了。
图中每1个boxplot，都是该位置的所有序列的测序质量的一个统计，上面的bar是90%分位数，下面的bar是10%分位数，箱子的中间的横线是50%分位数，箱子的上边是75%分位数，下边是25%分位数。
红色线代表中位数，蓝色代表平均数的连线，黄色是25%-75%区间，触须是10%-90%区间。若任一位置的下四分位数低于10或者中位数低于25，出现“警告”；若任一位置的下四分位数低于5或者中位数低于20，出现“失败，Fail”。

一般要求此图中，所有位置的10%分位数大于20,也就是我们常说的Q20过滤（即质量数低于20的reads清除掉）W。

2.3 Per tile sequence quality

检查reads中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色代表偏离度小，质量好，越红代表偏离度越大，质量越差。纵轴是tail的Index编号。
这个图主要是为了防止，在测序过程中，某些tail受到不可控因素的影响而出现测序质量偏低。蓝色代表测序质量很高，暖色代表测序质量不高，如果某些tail出现暖色，可以在后续分析中把该tail测序的结果全部都去除

2.4 Per sequence quality scores，reads

质量的分布，当峰值小于27时，警告；当峰值小于20时，fail。可以看出这个测序的报告峰值在38左右。
假如我测的1条序列长度为101bp，那么这101个位置每个位置Q之的平均值就是这条reads的质量值。该图横轴是0-40，表示Q值，纵轴是每个值对应的reads数目

2.5 Per base sequence content

对所有reads的每一个位置，统计ATCG四种碱基的分布，横轴为位置，这里是1-51bp，纵轴为碱基含量（百分数），正常情况下每个位置每种碱基出现的概率是相近的，四条线应该平行且相近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。本结果前10个位置，每种碱基频率有明显的差别，说明有污染。当任一位置的A/T比例与G/C比例相差超过10%，报”WARN”；当任一位置的A/T比例与G/C比例相差超过20%，报”FAIL”。
理论上来说，A和T应该相等，G和C应该相等，但是一般测序的时候，刚开始测序仪状态不稳定，很可能出现上图的情况。像这种情况，即使测序的得分很高，也需要cut开始部分的序列信息，一般遇到到这种情况，会cut前面5bp。

2.6 Per Sequence GC Content

统计reads的平均GC含量的分布。横轴是0 - 100%；纵轴是每条序列GC含量对应的数量。红线是实际情况，蓝线是理论分布（正态分布，均值不一定在50%，而是由平均GC含量推断的），两个应该比较接近才比较好。曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时，报”WARN”；偏离理论分布的reads超过30%时，报”FAIL”。当红色的线出现双峰，基本肯定是混入了其他物种的DNA序列。这张图中的信息良好。

2.7 Per base N content

当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”，统计N的比率。正常情况下，N值非常小。当任意位置的N的比例超过5%，报”WARN”；当任意位置的N的比例超过20%，报”FAIL”。

2.8 Sequence Length Distribution

reads长度分布，每次测序仪测出来的长度在理论上应该是完全相等的，但是总会有一些偏差。比如此图中，51bp是主要的，但是还是有少量的50和52的长度，不过数量比较少，不影响后续分析。当测序的长度不同时，如果很严重，则表明测序仪在此次测序过程中产生的数据不可信。当reads长度不一致时报”WARN”；当有长度为0的read时报“FAIL”。

2.9 Sequence Duplication Levels

统计不同拷贝数的reads的频率。测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在。横坐标是duplication的次数，纵坐标是duplicated reads的数目，以unique reads的总数作为100%。下图中，横坐标表示序列重复水平；纵坐标代表重复水平序列占所有序列的百分比。下图中大于10个重复的reads占总序列的20%以上，其他依次类推。当非unique的reads占总数的比例大于20%时，报”WARN”；当非unique的reads占总数的比例大于50%时，报”FAIL“。

2.10 Overrepresented sequences

一条序列的重复数，因为一个转录组中有非常多的转录本，一条序列再怎么多也不太会占整个转录组的一小部分（比如1%），如果出现这种情况，不是这种转录本巨量表达，就是样品被污染。这个模块列出来大于全部转录组1%的reads序列，但是因为用的是前200,000条，所以其实参考意义不大，完全可以忽略。

2.11 Adapter content

接头含量。①此图衡量的是序列中两端adapter的情况；②如果在当时fastqc分析的时候-a选项没有内容，则默认使用图例中的四种通用adapter序列进行统计；③本例中adapter都已经去除，如果有adapter序列没有去除干净的情况，在后续分析的时候需要先使用cutadapt软件进行去接头。

2.12 Kmer content

①这个图统计的是，在序列中某些特征的短序列重复出现的次数；② 我们可以看到1-8bp的时候图例中的几种短序列都出现了非常多的次数，一般来说，出现这种情况，要么是adapter没有去除干净，而又没有使用-a参数；要么就是序列本身可能重复度比较高，如建库PCR的时候出现了bias；④ 对于这种情况，我的办法是可以cut掉前面的一些长度，可以试着cut 5~8bp

批量质控结果查看MultiQC

# 安装multiqc
$ conda install -c bioconda multiqc
# 测试
$ multiqc --help
# 进入存放QC结果的文件夹，空白处右键“在终端打开”
# 扫描结果文件，忽略html文件
$ multiqc ./*fastqc.zip --ignore *.html
# 最后会默认生成一个名为multiqc_report.html文件，用浏览器查看

参考资料：

孟浩巍.20160410 测序分析——使用 FastQC 做质控