转录组（三）：了解fastq测序数据

时间：2024-03-10 12:58:33

通过fastq-dump将sra文件转换为fastq格式

#将之前的SRR3589956~SRR3589962转换为fastq格式
#批量生成脚本，并让每条命令后台执行，速度更快
for i in `seq 56 62`;do echo fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR35899$i -O /home/u3893/rna_seq_AKAP95/data/rna_seq \&;done >>raw_data.sh

fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589956 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589957 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589958 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589959 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589960 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589961 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &
fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR3589962 -O /home/u3893/rna_seq_AKAP95/data/rna_seq &

#也可以循环执行命令，但时间较长
for i in `seq 56 62`;do fastq-dump --gzip --split-3 /home/u3893/rna_seq_AKAP95/sra/SRR35899$i -O /home/u3893/rna_seq_AKAP95/data/rna_seq;done

fastq-dump用法

默认情况下fastq-dump不对reads进行拆分, 对于很早之前的单端测序没有出现问题.但是对于双端测序而言,就会把原本的两条reads合并成一个,后续分析必然会出错.

常见的参数有三类:

1.--split-spot: 将双端测序分为两份,但是都放在同一个文件中
2.--split-files: 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads直接丢弃
3.--split-3 : 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads会单独放在一个文件夹里

关于遇到的Rejected 35403447 READS because of filtering out non-biological READS就是因为原来是SE数据,但是用--split-3当作PE数据处理,出现的问题. 看起来好像有问题,但是对后续结果分析没有太多影响.
因此,对于一个你不知道到底是单端还是双端的SRA文件,一律用--split-3.

--gzip 输出的压缩格式
-O|--outdir : 输出到指定文件夹
-A参数，转换文件名称

fastq格式

fastq文件
第一行以@开头，之后为序列的标识符以及描述信息（与FASTA格式的描述行类似）
第二行为序列信息
第三行以+开头，之后可以再次加上序列的标识及描述信息（可选）
第四行为质量得分信息，与第二行的序列相对应，长度必须与第二行相同
碱基质量得分与错误率的换算关系：
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

Fastqc进行质控

#先生成批量脚本
#do后面接两条命令也同样用；连接，和多条命令写在一行无异
for i in `seq 56 62`;do echo fastqc -o ~/rna_seq_AKAP95/data/QC/ ~/rna_seq_AKAP95/data/rna_seq/SRR35899$i\_1.fastq.gz \&;echo fastqc -o ~/rna_seq_AKAP95/data/QC/ ~/rna_seq_AKAP95/data/rna_seq/SRR35899$i\_2.fastq.gz \&;done >>raw_data.sh

fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589956_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589956_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589957_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589957_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589958_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589958_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589959_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589959_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589960_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589960_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589961_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589961_2.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589962_1.fastq.gz &
fastqc -o /home/u3893/rna_seq_AKAP95/data/QC/ /home/u3893/rna_seq_AKAP95/data/rna_seq/SRR3589962_2.fastq.gz &

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 序列。

质控结果

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

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

    image

  • 左边是目录概要，可以点击想要看的结果，右边会跳转到特定详细的可视化结果。绿色代表“通过”，黄色代表“警告”，红色代表“不通过，失败”。

    image

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

    image

  • Per base sequence quality，每个read各位置碱基的测序质量。横轴碱基的位置，纵轴是质量分数，Quality score=-10log10p（p代表错误率），所以当质量分数为40的时候，p就是0.0001，质量算高了。红色线代表中位数，蓝色代表平均数，黄色是25%-75%区间，触须是10%-90%区间（黄色和触须我不是特别明白）。若任一位置的下四分位数低于10或者中位数低于25，出现“警告”；若任一位置的下四分位数低于5或者中位数低于20，出现“失败，Fail”。

    image

  • Per tile sequence quality，检查reads中每一个碱基位置在不同的测序小孔之间的偏离度，蓝色代表偏离度小，质量好，越红代表偏离度越大，质量越差。

    image

  • Per sequence quality scores，reads质量的分布，当峰值小于27时，警告；当峰值小于20时，fail。我的报告峰值在38。

    image

  • Per base sequence content，对所有reads的每一个位置，统计ATCG四种碱基的分布，横轴为位置，纵轴为碱基含量，正常情况下每个位置每种碱基出现的概率是相近的，四条线应该平行且相近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。本结果前10个位置，每种碱基频率有明显的差别，说明有污染。当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"。

    image

  • Per Sequence GC Content，统计reads的平均GC含量的分布。红线是实际情况，蓝线是理论分布（正态分布，均值不一定在50%，而是由平均GC含量推断的）。 曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"。

    image

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

    image

  • Sequence Length Distribution，reads长度分布，当reads长度不一致时报"WARN"；当有长度为0的read时报“FAIL”。

    image

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

    image

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

    image

  • Adapter content，接头含量

    image

  • Kmer content

    image

参考资料：https://www.plob.org/article/5987.html；http://blog.sina.com.cn[图片上传中...(image-f24716-1580318982095-14)] /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文件，用浏览器查看，具体看青山屋主的介绍。

参考
http://www.biotrainee.com/thread-1744-1-1.html

https://www.jianshu.com/p/14fd4de54402

相关文章