搜索
查看: 5195|回复: 1

一行脚本判断你的fastq测序数据的质量编码方式

[复制链接]

365

主题

512

帖子

1713

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1713
发表于 2017-8-26 21:10:46 | 显示全部楼层 |阅读模式
这个有时候很重要,尤其是处理一些古代的数据。

代码如下:
[Shell] 纯文本查看 复制代码
head input.fastq | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | 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<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}'


当然,现代人用的是数据都是phred33啦,质量值都是大写字母哦。古代人才用phred64,一堆小写字母和字符



上一篇:bam文件里面的测序质量值问题 phred64 vs phred33
下一篇:转录组入门分析5_序列比对
回复

使用道具 举报

365

主题

512

帖子

1713

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1713
 楼主| 发表于 2017-8-26 21:12:24 | 显示全部楼层
本帖最后由 ydchen 于 2017-8-26 21:14 编辑

最近在学习质控知识时, 对于质量值体系及转换产生了一些疑问, 作了一些尝试, 趁集群故障, 在此总结一下
质量值体系
相比之前培训时所学的质控内容, (我拿到的) 流程中还多了一步 phred33to64, 也就是把 .fastq 格式的数据从 Phred33 质量值体系转换为 Phred64 质量体系, 于是先补充学习了下质量值体系:
首先要从质量值说起, 测序仪器下机数据的 fastq 文件中, 每条序列都对应了相同长度的质量值, 反映出每个碱基的准确性和可靠性, (现在主流用的) 计算公式为:
Q = -10log10p
而这个 p 值就是 Phred 计算出来的, 表示一个碱基被识别错误的可能性, Phred 一开始是一个软件 (或者说计算方法), 对测序仪器识别到的荧光强度 (三代的不了解) 进行评估, 针对不同仪器有不同的标准表, 然后根据表中荧光强度的范围和分辨率分析得出碱基的 p 值, 由于比较可靠, 逐渐被各大公司采纳 (以上脑补翻译自 Phred quality score 和 Phred base calling)
反正, Q 值为 10 就表示这个碱基有 90% 的概率是正确的, 20 就是 99%, 40 就是 4 个 9, 很好记, 相信大家也都很熟悉
但是这时候问题就来了, 因为一个碱基对应一个质量值, 可 Q 值可以是一位数也可以是两位数, 连在一起的话就分不出哪个对应哪个碱基 (比如某两个碱基的质量值序列为 123 , 则可能为 1|23 或 12 | 3, 当然这是个极端的例子), 此外也浪费存储空间, 因为一个数字只有 10 种可能性, 却要占去一个字符位, 而一个字节有 8 位, 理论上已经可以代表 256 种状态, 这还没换算一个字符要占多少字节, 因此就会把碱基质量值转换为相应的 ASCII 码, 这是计算机中最基本的一套字符体系, 用来把常用符号 (共 127 个) 转为二进制以便于机器使用, ASCII 码表见 ASCII code chart, 这样一个字符就可以搞定了, 很方便很省事
但是这时候问题又来了, 最理想的情况当然是直接把质量值作为序号找出对应的那个 ASCII 码, 比如质量值为 40 就换成十进制 40 对应的 ASCII 码,可惜质量值根据测序仪公司的标准不同, 范围也各不相同, 基本都包括了 0 至 40 的区间, 甚至还可能是负值, 这就没法愉快地玩耍了, 而且人家 ASCII 码表也不配合, 0 到 31 对应的都是些控制字符 (比如回车, 退格), 根本不适合打印和保存, 可打印的都得从 32 号排起 (参见 ASCII printable code chart ), 所以各家测序仪器公司就把质量值再加上某个固定值作为 ASCII 码转换成了可打印字符从而保存在 FASTQ 文件中
可是问题还没有完, 这个固定值是多少好呢, 各家公司是竞争对手, 怎么可能你用什么我也用什么, 所以 Sanger 公司加了 33, 也就是质量值为 0 就转换成 ASCII 码 33, 查表可知为 !, 也即从可打印的字符开始 (排除了空格), 这就是现在所谓的 Phred33 体系, 当时的 Solexa (后来被 Illumina 收购) 公司就偏不用 33 (此处为个人脑补), 偏要加个 64, 这样质量值为 0 就用 @ 表示, 后面从 1 开始的就依次对应了 ABCD, 于是就成了 Phred64 体系, 至于当时三巨头的另一家测序仪公司 454 Life Sciences (后被 Roche 收购) 就更绝, 人家从碱基开始就不用 ACTG 表示, 直接整了个 ColorSpace 体系出来, 根本不和你们玩, (话说 Color Space 这玩意儿曾经把我狠狠地坑了好久), 当然后来大家也不跟 454 玩了, 最后他也就没得玩了
回到质量值体系, 这样就由 Sanger 公司和 Illumina 公司产生了 Phred33 体系和 Phred64 体系, 两家互相拗着, 这就辛苦了写生物信息分析软件的人, 两种质量体系都要考虑, 当然好一点的软件都是有参数接受体系类型的, 更好一点的软件就会自动判断体系类型进行对应转换
本来就这样结束的话也算是个圆满的故事了, 可是墨菲他老人家不高兴了 (Murphy's law), 所以在 2011 年, Illumina 公司表示他们又要改成 Phred33 体系了 (Upcoming changes in CASAVA), 真是大(wo)快(le)人(ge)心(qu)啊! 这么来回一折腾, 结果还是回到了最初的起点, 与老基友相拥而泣了, Phred33 一统江湖! (当然实际上现在 Ilumina 的 Phred33 和最初 Sanger 的 Phred33 还是有点区别的, 详见后文)
[AppleScript] 纯文本查看 复制代码
fqtype () {
        less $1 | head -n 999 | 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<=74 && min<59) print "Phred+33"; \
        else if(max>73 && min>=64) print "Phred+64"; \
        else if(min>=59 && min<64 && max>73) print "Solexa+64"; \
        else print "Unknown score encoding"; \
        print "( " min ", " max, ")";}'
}



再来说说这两个质量体系在实际使用中的影响, 其实按道理自 Illumina 改回 Phred33 后就省事多了, 各种软件都默认支持 33 体系, 再保留个参数接受 64 体系的就行了, 事实上, 现在好多软件也就是这么做的, 不过仔细看前文的各体系示例图, 还是可以发现 Illumina 现在的 Phred33 已经不再是以前的那个 Phred33 了, 质量值的范围是 0 至 41, 而且还特地标了个 2 出来, 这又是 Illumina 自己搞出来的新花样, 官方解释是: 如果一条 read 结尾处有大量低质量碱基, 那么这整条 read 都会被标为 2, 所以其实就是提前给过滤了下
总结一下, 想要快速判断一个 fastq 文件是什么质量体系, 最简单的方法就是直接用肉眼看文件里的质量值, 如果大多数都是小写字母, 那就是 Phred64, 如果大多数是大写字母, 数字和字符, 那就是 33 体系了, 此外, 还要注意的是, 大多数情境下或者软件参数中, Illumina 就是指 Phred64 体系, Sanger 指 Phred33 体系, 同时, 这两种体系一般也会写作 Phred+64 和 Phred+33
实际使用中的作用
现在来说说在流程中, 从原始数据到各自分析模块所需要的数据之间发生了什么, 主要是过滤和比对这两个环节, 而质量控制发生在这两步前后, 也就是过滤前, 过滤后比对前, 比对后, 共三次质控
单讲质量分析, 首先, 我拿到的流程里有一步 phred33to64 , 用 C 脚本来把 phred33 体系的 fastq 文件转换成 phred64 体系, 我看了下转换前后的文件, 好像就真是单纯的转换而已 (其实 read 名称也稍微动了下, 后面会讲到), 然后会用 check_gz_fastq 和 distribute_fqcheck.pl 这两个脚本分别进行数据质量统计和作图, 从而完成了一次简单的质量报告, 根据几个关键数据, 碱基含量分布图和碱基质量分布图来判断数据是否符合要求, 能否继续下一步分析, 如果不符合要求, 问题出在哪里, 要如何调整等等
现有流程中各脚本作用分析
根据我的理解, 之所以要把好好的 33 体系转成 64 体系, 应该只是因为 check_gz_fastq 这个脚本没有办法接收 33 体系的数据 (至少我不知道), 而 distribute_fqcheck.pl 这个画图脚本要用到 check_gz_fastq 这个脚本输出的信息, 仔细分析下: check_gz_fastq 获得的 fqcheck 统计文件中, 有用信息为总 reads 数, 序列长度总长, error rate, %GC, Q20, Q40, %G-%C, %A-%T) , distribute_fqcheck.pl 作出随读长的 GC/AT 含量分布图和碱基质量值分布图
如果原始数据质量没有问题, 就用 SOAPnuke 软件进行过滤, 而 SOAPnuke 命令行参数中有 -Q 和 -G 来指定接收和 (或) 输出不同的数据质量值编码类型, 过滤完了再作一次质控, 如果也能用, 就用 BWA 软件进行比对 (无参的不了解), 而 BWA 软件文档中未明确提及接收的是什么质量值体系, 不过根据搜索作者的回复可知, 其默认是接受 Phred33 体系的, 同时如果在比对时不考虑质量值的影响, 则输出结果中的质量值应该是与输入相同的
在下游的分析软件或脚本中, 则需要根据文档说明和具体分析内容, 来决定是否能够修改 (或是否需要修改) 参数设定和脚本内涉及质量值计算的部分
现存问题和解决
从上面的分析可以看到, 我们现有流程中的 phred33to64 主要是为了符合质量报告环节中现有工具的输入需要, 也即 check_gz_fastq 和 distribute_fqcheck.pl 这两个脚本, 那么这一步额外的转换是否有弊端, 又能否改进, 改进的话会有什么影响呢, 下面就来分析一下
首先, 多一步转换自然需要多一些时间和运算资源, 并且我观察到 check_gz_fastq 输出的 fqcheck 文件中质量值的分布是个稀疏矩阵, 也即大多数质量值都不存在, 似乎与常识有异, 同时 phred33to64这个 C++ 脚本并没有人在维护, 也没有明确的文档, 如果出了问题, 也很难定位和修正
然后, 既然我们只是为了一份质量报告, 那么现在同样有很多 fastq 文件的分析软件, 是不是可以考虑用其它工具来实现呢, 我测试了用 java 写的 FastQC 软件和 R 包 Rqc, 简要介绍和比较如下
FastQC 软件支持 fastq, sam, bam 等多种输入格式, 支持多进程, 文档详细: Analysis Modules, 输出 html 格式的报告, 同时附有数据的 summary 文件, 在报告中会对各项指标作一个简单的评估, 方便快速定位问题, 并且能在对应文档说明中找到可能造成这些问题的原因
此外, 如果输入文件是未经转换的原始 Phred33 数据, FastQC 还可以展示测序过程中, 测序仪样品池中各部分的质量值分布 (如果是 phred33to64 得到的数据就没这项分析内容), 从而帮助分析究竟是测序过程中存在问题, 还是建库 PCR 中存在问题
Rqc 相对要简单些, 只接受 fastq 文件, 只输出 html 文件, 但是可以多个文件一起分析, 在报告中就把同一指标下多个样本的结果直接画在同一张图上, 另外还有个好处就是因为是开源的, 内部函数都知道, 可以根据自己的需求调整改进, 或者输出中间结果, 另外 Rqc 画出来的图也要好看一些
当然, 这种替代方案也有其不足之处: 比如 FastQC 出于性能考虑, 对部分指标计算时只取输入文件中的前 100,000 条 reads (不过已经有足够的代表性), 另外更关键的一点就是这两种方法都无法得到 %Q20 和 %Q40 的值, 而这两个指标似乎是国内比较看重的, 至于这个指标是否不可替代, 此处不作更多讨论
初步结论
由于目前我们得到的数据基本都是 Illumina 1.5+ 的格式, 也就是 Phred+33 质量值体系的, 而主流软件也都对 Phred33 和 Phred64 有很好的支持, 因此对于我们的现有流程而言, 从原始下机数据开始, 质量控制可以用 FastQCRqc 实现, 过滤 (SOAPnuke), 比对 (BWA) 及后续分析中各软件也都是可以调节的, 那么, 出于减少复杂度, 和更好地适应以后的数据格式 (基本都是 Phred33 了, 很可能以后的软件会逐渐放弃支持 Phred64 体系) 的目的, 应该尽量转移到直接使用 Phred33 体系数据这种方案上来, 而唯一的不足恐怕就是上文提到的, 质量报告中没有 %Q 这一指标
最后, 如果一定要转换一次质量值, 也推荐用 Li Heng 写的 seqtk 软件来实现


回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-24 09:16 , Processed in 0.036817 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.