转录组入门(3):了解fastq测序数据

作业要求

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


实验过程

1.fastq-dump将sra数据转换成fastq格式

# 需要将作业2(http://www.jianshu.com/p/da377252ee96)中下载的测序数据用工具sratoolkit转换成fastq的格式。
$ for ((i=56;i<=62;i++));do fastq-dump --gzip --split-3 -A ~/disk2/sra/SRR35899$i.sra -O ~/disk2/data/rna-seq;done

fastq-dump 用法:

--gzip 使得输出的结果是.gz 的格式
--split-3 对于PE测序,输出的结果是两个
_1.fastq.gz
-A| --accession 输入你的sra文件可以是绝对路径,我的数据来源是~disk2/sra/SRR35899$i.sra ( 如果你直接写accession,那么fastq-dump 会默认重新下载数据,并且会放在~/ncbi/public/sra目录下
-O 是设置输出的目录

fastq格式:

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

fastq数据格式

每条read由4行字符构成:
第一行:必须以@开头,后面跟着序列的唯一ID以及相关说明内容。
第二行:核酸序列,是有ATCGN字符组成。
第三行:“+”开头,内容和第一行@后面的一样。
第四行:每个测序碱基质量,是用ASCII码来表示的,与第二行的字符数一致。
碱基质量得分与错误率的换算关系:
Q = -10log10p(p表示测序的错误率,Q表示碱基质量分数)
ASCII值与碱基质量得分之间的关系:
Phred64 Q=ASCII转换后的数值-64
Phred33 Q=ASCII转换后的数值-33

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

2.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
质控结果文件

3.质控结果查看

质控结果有14个html文件,你可以选择用浏览器打开查看最终的QC reports。

  • 首先来大概看一下QC结果报告。


    QC可视化结果——双击html文件,在浏览器中直接打开
  • 左边是目录概要,可以点击想要看的结果,右边会跳转到特定详细的可视化结果。绿色代表“通过”,黄色代表“警告”,红色代表“不通过,失败”。


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


    基本数据统计
  • Per base sequence quality,每个read各位置碱基的测序质量。横轴碱基的位置,纵轴是质量分数,Quality score=-10log10p(p代表错误率),所以当质量分数为40的时候,p就是0.0001,质量算高了。红色线代表中位数,蓝色代表平均数,黄色是25%-75%区间,触须是10%-90%区间(黄色和触须我不是特别明白)。若任一位置的下四分位数低于10或者中位数低于25,出现“警告”;若任一位置的下四分位数低于5或者中位数低于20,出现“失败,Fail”
    各位置碱基质量
  • Per tile sequence quality,检查reads中每一个碱基位置在不同的测序小孔之间的偏离度,蓝色代表偏离度小,质量好,越红代表偏离度越大,质量越差。


    偏离度
  • Per sequence quality scores,reads质量的分布,当峰值小于27时,警告;当峰值小于20时,fail。我的报告峰值在38。
    reads质量分布
  • Per base sequence content,对所有reads的每一个位置,统计ATCG四种碱基的分布,横轴为位置,纵轴为碱基含量,正常情况下每个位置每种碱基出现的概率是相近的,四条线应该平行且相近。当部分位置碱基的比例出现bias时,即四条线在某些位置纷乱交织,往往提示我们有overrepresented sequence的污染。本结果前10个位置,每种碱基频率有明显的差别,说明有污染。当任一位置的A/T比例与G/C比例相差超过10%,报"WARN";当任一位置的A/T比例与G/C比例相差超过20%,报"FAIL"。
    碱基分布
  • Per Sequence GC Content,统计reads的平均GC含量的分布。红线是实际情况,蓝线是理论分布(正态分布,均值不一定在50%,而是由平均GC含量推断的)。 曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差(overrepresented reads)。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时,报"WARN";偏离理论分布的reads超过30%时,报"FAIL"。
    reads 平均GC含量分布
  • Per base N content,当测序仪器不能辨别某条reads的某个位置到底是什么碱基时,就会产生“N”,统计N的比率。正常情况下,N值非常小。当任意位置的N的比例超过5%,报"WARN";当任意位置的N的比例超过20%,报"FAIL"。
    各位置N的reads比率
  • Sequence Length Distribution,reads长度分布,当reads长度不一致时报"WARN";当有长度为0的read时报“FAIL”。
    reads 长度分布
  • Sequence Duplication Levels,统计不同拷贝数的reads的频率。测序深度越高,越容易产生一定程度的duplication,这是正常的现象,但如果duplication的程度很高,就提示我们可能有bias的存在。横坐标是duplication的次数,纵坐标是duplicated reads的数目,以unique reads的总数作为100%。下图中,大于10个重复的reads占总序列的20%以上,其他依次类推。当非unique的reads占总数的比例大于20%时,报"WARN";当非unique的reads占总数的比例大于50%时,报"FAIL“。
    统计不同拷贝数的reads的频率
  • Overrepresented sequences,一条序列的重复数,因为一个转录组中有非常多的转录本,一条序列再怎么多也不太会占整个转录组的一小部分(比如1%),如果出现这种情况,不是这种转录本巨量表达,就是样品被污染。这个模块列出来大于全部转录组1%的reads序列,但是因为用的是前200,000条,所以其实参考意义不大,完全可以忽略。


    一条序列的重复数
  • Adapter content,接头含量


    接头含量
  • Kmer content


    Kmer含量

参考资料:https://www.plob.org/article/5987.htmlhttp://blog.sina.com.cn/s/blog_1319a10ee0102vfbx.html

4.质控结果批量查看神器——MultiQC

知乎青山屋主写的为知笔记介绍了multiQC软件--批量显示QC结果。

# 利用Anaconda来安装MultiQC非常方便,首先得安装Anaconda,用清华源下载,特别快,而官网实则难以接受。
# 清华源地址:https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
# 官网:https://www.continuum.io/downloads/
$ cd ~/src
$ wget  https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive/
# 下载到的是shell脚本文件,直接运行,安装完成
$ bash Anaconda2-4.4.0-Linux-x86_64.sh
# 添加 Anaconda Python 免费仓库
$ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
$ conda config --set show_channel_urls yes
# 然后直接安装MultiQC
$ conda install -c bioconda multiqc
# 测试
$ multiqc --help
# 进入存放QC结果的文件夹,并执行multiqc
$ cd ~/disk2/data/QC
# 扫描结果文件,忽略html文件
$ multiqc /data/*fastqc.zip --ignore *.html
# 最后会默认生成一个名为multiqc_report.html文件,用浏览器查看,具体看青山屋主的介绍。

参考资料:https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
https://www.continuum.io/downloads/
http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ1iRTvV2GwkwL2AaxYi2fXHP7

推荐阅读更多精彩内容