FASTQ形式の解釈と操作（1）

FASTQとは

FASTQ形式はテキストベースの形式で、DNAなどの塩基配列とそのクオリティスコアを1つのファイルに一緒に保存する際に用いられる。 塩基配列とクオリティスコアは各1文字のASCII文字で表され、これにより塩基とクオリティの対応関係が分かりやすくなっている。この形式はWellcome Trust Sanger Instituteで開発されたもので、FASTA形式のデータにクオリティ情報を追加するためのものであった。しかし今日ではIllumina Genome Analyzerのような次世代シーケンサー等から出力された塩基配列のデータを保存する際のフォーマットとしてデファクトスタンダードとなっている。
（Wikipedia の Fastq より）

ちなみに、

• 英語版Wikipediaの方が詳しい：FASTQ format
• 塩基配列だけのものをFASTA形式と呼ぶ：FASTA format

こんな感じ

ニホンザルの腸内細菌叢の研究から16S rRNA配列データを拝借（DDBJアクセッション番号：DRA006677）。

@HWI-M01574:19:000000000-A9784:1:1101:13519:1978
GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGGGTAGAAATACCCTTAGTGGCGGACGGGTGAGTAACGCGTGAGCAATCTGCCCTTGTCAGGGGAATAACACATAGAAATGTGTGCTAATACCGCATGACGTTACGAAGTGGCATCACTTTGTAACCAAAGAATTATCGGACAAGGATGAGCTCGCGTCCCATTAGGTAGTTGGTGAGTTAACGGCCCACAAAGCCGACGATGGGTAGCCGACCTGAGAGGGTGATCGGCGACATTG

GGGGGGGGGGCGGGGGGGGGGGGGGGGGGGGGGGGGGCFFGGGGCGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGCFEGFGGGGGGGGGGGFGAFGFGGGGGEGFFGGGGCFFDGAAFGGGGGCGGGGGG54>EEGGGGGGF4?DFGFFFGGFGFA=FFFFFFA6CFGFFFF>)1:BFGFG4;@3817(47;596ADF<>4FAFB>5223:?.6=?>(1096(2(3:8,-1291:BF:2?)410(397?().2::366<:04(310(66 @HWI-M01574:19:000000000-A9784:1:1101:12858:1993 GATGAACGCTAGCTACAGGCTTAACACATGCAAGTCGAGGGGAAACGACATTGGAAGCTTGCTTCCGATGGGCGTCGACCGGCGCACGGGTGAGTAACGCGTATCCAACCTGCCCACCACTTGGGGATAACCTTGCGAAAGTAAGACTAATACCCAATGAGATTCAAAGCGGACATCCAAGTTGAATTAAAGATTTATCGGTGATGGATGGGGATGCGTCTGATTAGCTTGTTGGCGGGGTAACGGCCCCCCAAGGCGACGATCAGTAGGGGTTCTGAGA + GGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGFGGGGGGGGGGGEFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGGGDGGGGGFGGGGGGGGGGGGGGGGCCGGFGFFFGGFGFF6BFFGGGEEGFGFEFA=FCFFEFFCGFEFFFF3CD>5(3AA<<<0>:?F?6:B:06B>F99?>?69>1-(-(48:?239>6(3:?B<9)8:041324)).

1本の配列は4行で記述される。1行目は@ではじまって配列名（ヘッダー）を表す。2行目は塩基配列。3行目は+で区切る。4行目は各塩基のクオリティ値を表す。原理上、2行目と4行目は同じ文字数になる。2行目と塩基配列や4行目のクオリティ配列は、原則、途中で折り返さない（※FASTA形式は折り返すことが多い）。クオリティ値が、1行目の文頭の@や、3行目の+と紛らわしくなるため。

上の例では2本の配列（8行）が並んでいる。@と+を目印にすることで、複数配列の記述が可能（Multi-FASTQ形式）。

クオリティ値の解釈

• 塩基配列は（縮重表現を使わなければ）A/T/G/C/Nの5種類の文字だけ使う。NはA/T/G/Cのどれかわからない。
• クオリティ値はASCII文字で表現する。イルミナのシークエンサーではPhredクオリティスコアが使われている。

出力された塩基（ベースコール）が**正しくない確率（probability of incorrect base call）**を p としたときに、クオリティ値 Q は、

Q = −10 * log_10(p)　⇔　p = 10^(−Q / 10)**

で計算される。例えば、クオリティ値 Q が 10 のときは、0.1 の確率でそのベースコールは正しくない。つまり90%確からしい（base call accuracy）。まとめてみると、

• Q = 0 ならば、確からしさは100% − 100% = 0%である。
• Q = 10 ならば、確からしさは100% − 10% = 90%である。
• Q = 20 ならば、確からしさは100% − 1% = 99%である。
• Q = 30 ならば、確からしさは100% − 0.1% = 99.9%である。
• Q = 40 ならば、確からしさは100% − 0.01% = 99.99%である。

90%の確からしさではさすがに信用できないので、イルミナのシークエンサーの結果をカットオフするときは、ベースコールが Q ≥ 20 とか Q ≥ 30 くらいの閾値を設定するのが普通である。

Phredクオリティスコア

• 現在のイルミナのシークエンサーは基本的にASCIIコードの10進数表現から、33を引いた値を Q とする（Phred+33表現）。
• 古いバージョン（Illumina 1.3+など）はPhred+64を使ってるので注意。

例えば!のASCIIコードは33なので、Q = 0。確からしさは0%。まとめてみると、

• クオリティ値が! (ASCII = 33)ならば、Q = 0 で、確からしさは0%。
• クオリティ値が+ (ASCII = 43)ならば、Q = 10 で、確からしさは90%。
• クオリティ値が5 (ASCII = 53)ならば、Q = 20 で、確からしさは99%。
• クオリティ値が? (ASCII = 63)ならば、Q = 30 で、確からしさは99.9%。
• クオリティ値がI (ASCII = 73)ならば、Q = 40 で、確からしさは99.99%。

FASTQデータを解読してみよう

前述のニホンザルの腸内細菌叢の研究から。

@HWI-M01574:19:000000000-A9784:1:1101:13519:1978
GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGAAGGGTAGAAATACCCTTAGTGGCGGACGGGTGAGTAACGCGTGAGCAATCTGCCCTTGTCAGGGGAATAACACATAGAAATGTGTGCTAATACCGCATGACGTTACGAAGTGGCATCACTTTGTAACCAAAGAATTATCGGACAAGGATGAGCTCGCGTCCCATTAGGTAGTTGGTGAGTTAACGGCCCACAAAGCCGACGATGGGTAGCCGACCTGAGAGGGTGATCGGCGACATTG

GGGGGGGGGGCGGGGGGGGGGGGGGGGGGGGGGGGGGCFFGGGGCGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGCFEGFGGGGGGGGGGGFGAFGFGGGGGEGFFGGGGCFFDGAAFGGGGGCGGGGGG54>EEGGGGGGF4?DFGFFFGGFGFA=FFFFFFA6CFGFFFF>)1:BFGFG4;@3817(47;596ADF<>4FAFB>5223:?.6=?>(1096<?(2(3:8,-1291:BF:2?)410(397?().2::366<:04(310(66

• 配列名はHWI-M01574:19:000000000-A9784:1:1101:13519:1978。
• A/T/G/Cが280個あるので、この配列の塩基長は280 bp。
• 塩基の1文字目のクオリティ値はG。GのASCIIコードは71なので、Q = 71 − 33 = 38。
• したがって、塩基の1文字目が正しくない確率は P = 10^(−3.8) = 0.00016。つまり、99.984%確からしい。
• クオリティ値Gがしばらく続くので、この配列の前半部は非常に確からしい。
• 140文字目に5 (ASCII = 53)が登場。これは Q = 20 で、確からしさは99%なので、140文字目以降からクオリティが下がっていくのがわかります。

クオリティ値の早見表

• クオリティ値がアルファベットのとき

  • クオリティ値が[A,B,...,I,J]なら Q 値は[32,33,...,40,41]。確からしさは99.9%から99.99%のオーダーなので非常に信頼できる。

• クオリティ値が? @のとき

  • それぞれ Q 値は30 31。確からしさは99.9%以下。

• クオリティ値が: ; < = >のとき

  • それぞれ Q 値は25 26 27 28 29。確からしさは99.7%～99.9%。

• クオリティ値が数字のとき

  • クオリティ値が[0,1,...,7,8,9]なら Q 値は[15,16,...,22,23,24]。確からしさは96.8%～99.6%。信頼しすぎない方がいい。

• クオリティ値がその他の記号! " # $ % & ' ( ) * + , - . /のとき

  • 確からしさは96.0%以下。適切な閾値を決めてカットするのが無難。

終わりに

• 今日はここまで。（2019年10月25日、了。）