专栏首页Y大宽fastq格式文件及phred33的判断

fastq格式文件及phred33的判断

Fastq格式文件储存了生物序列的信息及其质量信息。以电脑中的一个文件的为例

$ zcat SRR8518347.fq.gz |head
@SRR8518347.1 1 length=150
GTATTTTCTAAGACGTAATGCTACTTTACACTTAACAGACTACAGTATAGTGTAAACGTAACGGTTATGTCCACTAGGAAACCAAAATATTCGTGTTACTTGCTTTACTGTACTTGTTTTATTTCAGTAGTCGGGAACTAAACTCACAGT
+SRR8518347.1 1 length=150
AAA,,F<F,FKKAFFF<FFFKF<FA,KKFFKKKKKKKKKF,AKFK,A7<FKK,AFFKKFKKKA,AFF,FKAKKKA<FFKAFKKKKKA,F<KKKKFFKKKKK,FKKAKAKKKKF7FFFAKFAFFKK7,,FF<F,,<,AA<KAFKKKKKAF,
@SRR8518347.2 2 length=150
CCAGCTAATTTTTGTATTACTAGTAGAGACGGGGTTTCACCAGGTTGGTCAGGATGGTCCTGAACTCCTGAGCTCAGGCGATCCACCCACCTCAGCCTCCCAAAGTGGTAGGATTACAGGCGTGAGATCGGAAGAGCACACGTCTGAACT
+SRR8518347.2 2 length=150
AAAFAAFFFKKFFKKKKKKKKKKKKFKKA7FKKK(FFFKKKFFFKK,A,7F,FFKK,FK7AKKFKKK,FKA,AFF7<AA<AAFFFKFK<KKK,FKFFKKKKKKKKA,,7FKKKFKKK7FAFKAA7,A7FF,FFAK7<FAAFKKKK7FFF7

1 格式说明

第一行:必须@开头,紧跟唯一的序列的ID标识符,后面可跟其他描述性内容,但序列ID与描述部分空格分开,这点很重要,是很多分析的基础 第二行:序列的顺序,基本不可能出现重复。这里是核酸。 第三行:必须+开头,后面可以跟描述信息,此处是跟了第一行的信息。但为了节约空间,很多时候这一行只有一个+,而不跟任何内容。 第四行:碱基或氨基酸(此处碱基)的质量字符,对应着第二行的碱基,反应的是该碱基的错误率,所以这一行的字符数和第二行要一一对应,否则就乱了。这里就引入了ASCII code。

fastq format 1

fastq format 2

所以引入下面,碱基质量值是什么,如何获得?怎么表示?如何转换?

2 碱基质量值Q值和ASCII码之间的关系

因为第四行的编码,开始由Phred程序开发者定义,所以一般称为Phred quality。那碱基质量得分怎么来的?

Chromas软件展示的一个DNA 序列质量结果

每合成一个碱基,即可发出一个荧光信号,该信号可以被捕捉到,并生成是是轨迹数据。不同的碱基用不同颜色标记,检测相应峰值即可判断碱基。 而Phred通过计算相应波峰参数,去查询通过已知序列测序分析得到的一个表,即可把错误率转换为质量得分。也就是把波峰参数和质量得分对应起来。 碱基错误率与质量得分的关系如下

Phred quality score

也就是说,质量值Q是测序错误率的对数*-10。假如错误率是0.01,则Q值为20。可见,错误率越低,其Q值越高。即Q值越高越可靠。 Q10——0.1 Q20——0.01 Q30——0.001 Q40——0.0001 那就直接以上面这种Q值表示碱基质量不就行了吗?为什么要用ASCII码表示。 如果用数字表示,数字和数字之间需要有间隔符号以区分,再者会浪费存储空间,所以可以把质量值转变为相应的ASCII码,这样就完成了把质量数向ASCII码的转换,那现在看下ASCII码

ASCII码

如果直接把Q值直接对应ASCII码,应该挺方便的,但是,Q值有时会有负值,再者,看ASCII码的0-31位都是控制字符,没法打印和保存,能打印的从要从32位的Space开始,所以就可以给实际的Q值加上一个固定值,这样就可以打印出来而保存在fastq文件中了。 所以下面就是固定值加多少的问题。phred33,就是加了33,也就是10变成43,查表是+。例子中很多F,Q值是70-33=37。而当时的solexa加的是64,这就是Phred64.他们与ASCII的对应关系如下

Q score vs Phred

数据处理时,有些软件会根据碱基质量得分的不同做不同处理,所以需要指定正确的编码方式,也就是需要对Phred33还是64正确判断。不过现在基本都33了,但如果下载以前的数据不一定。 下面是不同版本质量得分和质量字符ASCII的关系

不同测序标记中的Phred的使用

从上面可以看出,Phred33的字符使用33-73,而+64使用包括59(包括)-104之间的ASCII码。所以,只要ASCII小于59的仅仅在Phred+33中使用,而+64的都大于59,而从表可以看到大于等于74的旨在Phred+74中使用,所以,如果软件没有自动判断,根据以上特点,就可以自行判断是Phred33还是Phred64。

3 如何判断是Phred33还是Phred64

默认读取1000条序列,在这1000条序列中:

    1. 如果有2个以上的质量字符ASCII值小于等于58(即有两个碱基的得分小于等于25),同时没有任何质量字符的ASCII值大于等于75,即判断是Phred+33。
    1. 如果有2个以上的质量字符ASCII值大于等于75(即有两个碱基的得分大于等于10),同时没有任何质量字符的ASCII值小于等于58,即判断是Phred+64。
    1. 如果所有质量字符的ASCII值介于59到74之间,即判断可能是Phred+33,但建议使用更多的序列做进一步测试(出现这种结果可能有两种情况:1, Phred+33编码,所有碱基质量得分介于26到42之间;2,Phred+64编码,所有碱基质量得分介于-5到10;是前者的可能性大)。
    1. 如果出现上述3种以外的情况,建议打印出质量字符的ASCII值人工判断。 脚本如下
$ cat fq_qual_type.sh 
less $1 | head -n 1000 | awk '{if(NR%4==0) printf("%s",$0);}' \
| od -A n -t u1 -v \
| awk 'BEGIN{min=100;max=0;} \
{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END \
{if(max<=126 && min<59) print "Phred33"; \
else if(max>73 && min>=64) print "Phred64"; \
else if(min>=59 && min<64 && max>73) print "Solexa64"; \
else print "Unknown score encoding"; \

运行

sh fq_qual_type.sh example.fq
本文参与 腾讯云自媒体分享计划 ,欢迎热爱写作的你一起参与!
本文分享自作者个人站点/博客:https://www.jianshu.com/u/51a71446d509复制
如有侵权,请联系 cloudcommunity@tencent.com 删除。
登录 后参与评论
0 条评论

相关文章

  • lncRNA组装流程的软件介绍之trim-galore

    生信技能树
  • 了解fastq文件

    第一行:以‘@’开头,是这一条 read 的名字,这个字符串是根据测序时的状态信息转换过来的,中间不会有空格,它是每一条 read 的唯一标识符,同一...

    生信喵实验柴
  • Trimmomatic 数据过滤

    Trimmomatic 是一个很常用的 Illumina 平台数据过滤工具。支持 SE 和 PE 测序数据。主要用来去除 Illumina 平台的 fastq ...

    生信编程日常
  • 第2篇:原始数据的质控、比对和过滤

    首先对拿到的原始测序数据(fastq或fastq.gz格式)进行质控检测,直接用fastqc软件,再加上multiqc将多个检测结果一起展示。 如:

    生信技能树
  • 转录组分析 | 使用Trimmomatic过滤Fastq文件

    上一期,小编教大家使用FastQC评估了自己手中RNA-seq数据的质量,今天教大家使用Trimmomatic切除数据中的接头序列和低质量序列。

    生信小王子
  • 5 trim_galore去接头

    ls SRR85182.gz | while read id;do (nohup fastqc $id &);done find SRR851819.gz|x...

    Y大宽
  • 高通量测序数据质控神器Trimmomatic

    高通量测序下机的原始数据中存在一些低质量数据、接头以及barcode序列等,为消除其对后续分析准确性产生的影响,在数据下机以后对原始数据进行质控处理就成了至关重...

    kongxx
  • trim-galore并行处理时的几个问题

    config是需要进行处理的文件列表 trim_galore命令这里用的也比较简单,总结下处理时遇到的问题

    Y大宽
  • 批量与并行不一样

    入门级做好配置文件命令脚本文件提交至后台进阶级做好配置文件命令脚本文件提交至后台补充一个错误的例子

    生信技能树
  • RNA-seq入门实战(一):上游数据下载、格式转化和质控清洗

    连续两次求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,以及 生信技能树知识整理实习生招募,让我走大运结识了几位优秀小伙伴!大家开始根据我的ngs组学视频进行...

    生信技能树
  • ngs组学数据分析上下游分析都可以基于R语言吗?

    虽然有点难度,但其实确实是可以的,对生信工程师来说,就是整理流程(把Linux命令替换成为R语言代码)工作量比较大。如果大家感兴趣而且确实有需求,不妨看看这个文...

    生信技能树
  • 通过简单数据熟悉Linux下生物信息学各种操作

    关于trimmomatic相关内容请见https://zhuanlan.zhihu.com/p/28924793 https://www.jianshu.co...

    Y大宽
  • 生信软件 | Trimmomati (质量控制,修剪低质和接头序列)

    trimmomatic SE -phred33 input.fq.gz output.fq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10...

    白墨石
  • 宏基因组reads筛选:去除宿主序列

    基于环境的复杂性与研究对象的不同,宏基因组数据在组装之前常需要过滤掉一些序列以防干扰研究。例如要研究动植物组织或肠道的微生物组,往往需要去除宿主的DNA序列。假...

    SYSU星空
  • RNA-seq数据分析完全指北-05:去接头以及过滤

    一般来说,测序结果如下如所示,包括barcode和部分insert。而barcode部分在demultiplexing会被去除,剩下的就只有测到的一部分inse...

    生信菜鸟团
  • NGS基础 - FASTQ格式解释和质量评估

    FASTQ文件格式和命名 高通量测序之后用于下游分析的数据一般存储在FASTQ文件中。为了节省空间,又不影响下游使用,也一般用gzip压缩的格式。 单端测序每个...

    生信宝典
  • 3 wes测序质量的控制

    没有原视频中文件,我用了下载的三个文件做例子。乳腺癌的组织样本。所以原视频中命令我也用不上,但是还是列出来

    Y大宽
  • 全套MeRIP-seq文章图表复现代码

    这里面的MeDIP-seq指的是DNA,那么MeRIP-seq其实就是RNA水平的又叫做m6a测序,恰好看到了咱们的表观微信交流群我们的生信技能树优秀转录组讲师...

    生信技能树
  • 纳米孔数据处理

    背景 前面介绍了纳米孔测序的原理与碱基识别,本次带大家认识纳米孔测序数据的格式,以及怎么质控与处理。

    生信喵实验柴

扫码关注腾讯云开发者

领取腾讯云代金券