Fastq 格式说明 & (Phred33 or Phred64)

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

image

每条序列由4行字符表示,上述样例显示有两条序列:

第一行:必须以“@”开头,后面跟着唯一的序列ID标识符,然后跟着可选的序列描述内容,标识符与描述内容用空格分开。

第二行:序列字符(核酸为[AGCTN]+,蛋白为氨基酸字符)。

第三行:必须以“+”开头,后面跟着可选的ID标识符和可选的描述内容,如果“+”后面有内容,该内容必须与第一行“@”后的内容相同。

第四行:碱基质量字符,每个字符对应第二行相应位置碱基或氨基酸的质量,该字符可以按一定规则转换为碱基质量得分,碱基质量得分可以反映该碱基的错误率。这行的字符数与第二行中的字符数必须相同。字符与错误率的具体关系见下文介绍。

在满足上述要求的前提下,不同的测序仪厂商或数据存储商对第一行和第四行的定义有些差别。

第一行,即标识行在Illumina和NCBI SRA中的样式如下:

Illumina casava 1.8+(详细的解释可参考wiki):

@HWI-ST1276:97:D1DCYACXX:7:1101:1406:2170 1:N:0:CGACGT

NCBI SRA:

@SRR387514.1 ILLUMINA-C4D679_0049_FC:1:12:3317:1141 length=40

对于第四行的编码,最初由Phred程序的开发者定义,一般称为Phred qualitiy. 在Illumina早起版本(v1.3,v1.4)中,因为对quality的定义与Phred的不同,这行应该称为 Solexa quality。但从Illumina v1.5以后,也开始采用Phred的定义。

碱基质量得分是怎么来的?

Phred最初是一个从测序仪中产生的荧光记录数据中识别碱基的程序。在早起的荧光染料测序中,每次发生碱基合成时会释放出荧光信号,该信号被CCD图像传感器捕获。记录下荧光信号的峰值,生成一个实时的轨迹数据(chromatogram)。因为不同的碱基用不用的颜色标记,检测这些峰值即可判断出对应的碱基。但由于这些信号的波峰、密度、形状和位置等是不连续或模糊的,有时很难根据波峰判断出正确的碱基。

image

图1 chromatogram样图

Phred计算许多与波峰大小和分辨率相关的参数,根据这些参数,从一个巨大的查询表中找出碱基质量得分。这个查询表是根据对已知序列的测序数据分析得到的(应该是分析得到波峰参数与碱基错误率的关系,再通过公式把错误率转换成质量得分,得到波峰参数与质量得分的直接对应表)。不同的测序试剂和机器用不同的查询表。为了节约磁盘空间,质量得分(可能占用两个字符)按一定规则(Phred+33或Phred+64)被转换为单个字符表示。

碱基错误率与质量得分的关系有如下两种:

Qphred = -10log10p

Qillumina-prior to v.1.4 = -10log10(p/(1-p))

image

图 2 质量得分Q和错误率p的关系,红色的为phred,黑色的为Illumina早期版本,虚线表明p=0.05,对应的质量得分为Q≈13

在不同版本的编码中,除了质量得分与错误率有所差别外,在字符与得分的转换上也有差别。

image

图3 不同版本质量得分与质量字符ASCII值的关系

质量字符的ASCII值和质量得分的关系有如下两种:

Phred+64 质量字符的ASCII值 - 64

Phred+33: 质量字符的ASCII值 - 33

可以粗略分为 Phred+33和Phred+64,这里的33和64就是指ASCII值转换为得分该减去的数值。

在处理测序数据时,因为一些软件会根据碱基质量得分的不同做不同的处理,常要指定正确的编码方式,有必要对质量字符与质量得分的关系(Phred+33或Phred+64)作出正确的判断。当然,如果处理的是最近两年产生的测序数据,基本上都是Phred+33的,但从NCBI SRA数据库下载的旧数据就不一定了。

根据图3中Phred+33与Phred+64所使用的质量字符范围的不同,可以对fastq文件中质量得分的编码方式做出判断。图3中显示,ASCII值小于等于58(相应的质量得分小于等于25)对应的字符只有在Phred+33的编码中被使用,所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情况下,ASCII值大于等于74的字符只出现在Phred+64中。利用这些信息即可在程序中进行判断。

文章末尾是一个对Phred+33或Phred+64做区分的perl脚本。

该脚本的判断思想如下:

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

1. 如果有2个以上的质量字符ASCII值小于等于58(即有两个碱基的得分小于等于25),同时没有任何质量字符的ASCII值大于等于75,即判断是Phred+33。

2. 如果有2个以上的质量字符ASCII值大于等于75(即有两个碱基的得分大于等于10),同时没有任何质量字符的ASCII值小于等于58,即判断是Phred+64。

3. 如果所有质量字符的ASCII值介于59到74之间,即判断可能是Phred+33,但建议使用更多的序列做进一步测试(出现这种结果可能有两种情况:1, Phred+33编码,所有碱基质量得分介于26到42之间;2,Phred+64编码,所有碱基质量得分介于-5到10;是前者的可能性大)。

4. 如果出现上述3种以外的情况,建议打印出质量字符的ASCII值人工判断。

理解错误的地方欢迎指正。

附录检测格式的perl:


#!/usr/bin/perl -w
use strict;
use Getopt::Long;
#
# fastq_phred.pl - Script for judge the fastq's encoding, whether it is phred33 or phred64.
#
# Version: 0.3 ( May 19, 2014)
# Author: Wencai Jie (jiewencai<@>qq.com), NJAU, China.
#
# Permission is granted to anyone to use this software for any purpose, without
# any express or implied warranty. In no event will the authors be held liable 
# for any damages arising from the use of this software.
#

#Get options.
my ($help, $print_score, $detail, $print_ascii, $reads_num, $reads_start_arg, $reads_end_arg);
my $reads_end_turn;
GetOptions(
    'help|h!' => \$help,
    'score|s!' => \$print_score,
    'detail|d!' => \$detail,
    'ascii|a!' => \$print_ascii,
    'reads_num|n=i' => \$reads_num,
    'reads_start|b=i' => \$reads_start_arg,
    'reads_end|e=i' => \$reads_end_arg,
);

my $usage = "
fastq_phred.pl:
This program can print fastq file's reads quality scores, ASCII value, and help to judge it's 
encoding by the ASCII value range, whether it is phred33 or phred64.

Usage:
    perl fastq_phred.pl [options] <file1.fq [file2.fq ...]>
Options:
    -h|--help         print this help message.
    -s|--score        print scores.                                 [default: Do not print scores] 
    -d|--detail       print detail scores or ASCII value when       [default: Do not print scores in detail] 
                          --score or --ascii set.
    -a|--ascii        print quality character's ASCII value. if     [default: Do not print ASCII vaule] 
                          this option set, the --score will disabled.
    -n|--reads_num    reads number used to test phred encoding      [default: 1000]
                          and print scores. It's advised to use more
                          than 100 reads to do the test.               
    -b|--reads_start  reads start position used to test phred       [default: 1]
                          encoding and print scores.
    -e|--reads_end    reads end position used to test phred         [default: the length of the read]
                          encoding and print scores.

";
    
if ($#ARGV < 0 or $help){ 
    print "$usage";
    exit;
}

#Check parameters.
unless ($reads_num){
    $reads_num = 1000;
}
if ($reads_start_arg && $reads_start_arg  < 0){
    print STDERR "ERROR:The reads start position should great than 0.\n\n";
    exit;
}
if ($reads_end_arg && $reads_end_arg  < 0){
    print STDERR "ERROR:The reads end position should great than 0.\n\n";
    exit;
}
if ($reads_start_arg && $reads_end_arg && $reads_end_arg < $reads_start_arg){
    print STDERR "ERROR:The reads start position should great than end position.\n\n";
    exit;
}

&main;

sub main{
    my $filename = '';
    while ($filename = shift @ARGV){
    my @FQ = ();
    my @all_ascii = ();
    my ($file_end, $phred_result) = ('','');
    my ($Q, $count, $lt_58, $gt_75) = (0, 0, 0, 0);
    open FQ,"<$filename" or die "Can not open $filename:$!\n";
    #Read sequences.
    while($count < $reads_num){
        $count++;
        @FQ=(); 
        #read four lines from fastq file.
        for(my $i=0; $i<=3; $i++){
            if (eof(FQ)){
                $file_end = 'yes';
                last;
            }
            $FQ[$i]=<FQ>;
            if ($FQ[0] !~ m/^@/){
                my $line = $count*4-3;
                print STDERR "ERROR:\n$filename: It's not a correct fastq format.\nline '$line': $FQ[0]\n";
                exit;
            }
        }
        if ( $file_end eq 'yes'){
            next;
        }
        my @ascii_ref = &cal_ascii($FQ[3], $reads_start_arg, $reads_end_arg);
        push @all_ascii, [@ascii_ref];
    }

    #print ASCII.
    if ($print_ascii){
        print "\n","."x50," ASCII Value: $filename ","."x50,"\n";
        &print_array_of_array(\@all_ascii, 0, $detail);
        next;
    }

    #Stastic ASCII value range.
    foreach my $ascii_ref (@all_ascii){
        $lt_58 += (grep { $_ <= 58} @{$ascii_ref});
        $gt_75 += (grep { $_ >= 75} @{$ascii_ref});
    }

    #Guess the Phred with ASCII value range. 
    if ($lt_58 > 1 && $gt_75 == 0 ){
        $Q = 33;
        $phred_result = "$filename: The encoding should be Phred33.\nThe quality score character number that ASCII value less than 58 : $lt_58\nThe quality score character number that ASCII value great than 75: $gt_75";
    }elsif($lt_58 == 0 && $gt_75 > 1){
        $Q = 64;
        $phred_result = "$filename: The encoding should be Phred64.\nThe quality score character number that ASCII value less than 58 : $lt_58\nThe quality score character number that ASCII value great than 75: $gt_75";
    }elsif($lt_58 == 0 && $gt_75 == 0){
        print STDERR "$filename: The encoding should be Phred33 that all of the nucleotide quality score great than 25 and less than 41, but it's advised to send more reads to be tested with '-n <int>' options.\n";
        exit;
    }else{
        print STDERR "$filename\nWarning: Abnormal endoding, Please test again with more reads or make a judgement by yourself with ASCII value by '-ascii' options.\n"; 
        exit;
    }

    #print score.
    if ($print_score){
        print "\n","."x50," Quality Score: $filename ","."x50,"\n";
        &print_array_of_array(\@all_ascii, $Q, $detail);
    }

    #Print the phred encoding result.
    print STDERR "$phred_result\n\n";
    }
}

#Print Score or ASCII value.
sub print_array_of_array{
    my ($array_of_array_ref, $Q, $detail) = @_;
    my ($average_value, $total_value, $value_num) = (0, 0, 0);
    my @array_of_array = @{$array_of_array_ref};
    my %value_h;
    foreach my $array_ref (@array_of_array){
        for (my $i=0;$i<=$#{$array_ref};$i++){
            my $out_value  = ${$array_ref}[$i] - $Q;
            $value_h{$out_value}++;
            $total_value += $out_value;
            $value_num ++;
            print "$out_value " if ($detail);
        }
        print "\n" if ($detail);
    }
    unless ($detail){
        foreach my $out_value (sort {$a <=> $b} keys %value_h){
            print "$out_value\t$value_h{$out_value}\n";
        }
    }
    $average_value = (int ($total_value/$value_num)*100) /100;
    print "Average: $average_value\n";
}

#Calculate phred score.
sub cal_ascii{
    my ($read,$reads_start, $reads_end) = @_;
    my @all_ascii = ();
    my $ascii = 0;
    #The $read string's end is a "\n";
    my $reads_len = length($read) - 1;
    #The $reads_end should be less than the read length.
    if( $reads_end_arg && $reads_end_arg <= $reads_len){
        $reads_end = $reads_end_arg;
    }else{
        $reads_end = $reads_len;
    }
    #If the the reads start position set, the $reads_start equal to it,
    #else the reads start position set to 1.
    if ( $reads_start_arg && $reads_start_arg <= $reads_end){
        $reads_start = $reads_start_arg;
    }else{
        $reads_start = 1;
    }
    #Convert 1 base coordinate system to 0 base coordinate system.
    for(my $j=$reads_start-1; $j<=$reads_end-1; $j++){
        $ascii = ord(substr($read,$j,1));
        push @all_ascii, $ascii;
    }
    return @all_ascii;
}



```




参考资料:

1. [https://en.wikipedia.org/wiki/FASTQ_format](https://en.wikipedia.org/wiki/FASTQ_format)

2. [https://en.wikipedia.org/wiki/Phred_quality_score](https://en.wikipedia.org/wiki/Phred_quality_score)

3. [https://en.wikipedia.org/wiki/Phred_base_calling](https://en.wikipedia.org/wiki/Phred_base_calling)

4. [http://maq.sourceforge.net/fastq.shtml](http://maq.sourceforge.net/fastq.shtml)

5. [http://maq.sourceforge.net/qual.shtml](http://maq.sourceforge.net/qual.shtml)

6. [http://supportres.illumina.com/documents/myillumina/a557afc4-bf0e-4dad-9e59-9c740dd1e751/casava_userguide_15011196d.pdf](http://supportres.illumina.com/documents/myillumina/a557afc4-bf0e-4dad-9e59-9c740dd1e751/casava_userguide_15011196d.pdf)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容