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

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

ちなみに、

- 英語版Wikipediaの方が詳しい：[FASTQ format](https://en.wikipedia.org/wiki/FASTQ_format "FASTQ format")
- 塩基配列だけのものをFASTA形式と呼ぶ：[FASTA format](https://en.wikipedia.org/wiki/FASTA_format "FASTA format")

# こんな感じ
[ニホンザルの腸内細菌叢の研究](https://rdcu.be/2KC6)から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<FFBFFFF?(9?BF>>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形式）。

# クオリティ値の解釈

- 塩基配列は（[縮重表現](https://en.wikipedia.org/wiki/FASTA_format#Sequence_representation)を使わなければ）`A/T/G/C/N`の5種類の文字だけ使う。`N`は`A/T/G/C`のどれかわからない。
- クオリティ値はASCII文字で表現する。イルミナのシークエンサーでは[Phredクオリティスコア](https://en.wikipedia.org/wiki/Phred_quality_score)が使われている。
 
出力された塩基（ベースコール）が**正しくない確率（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コード](https://en.wikipedia.org/wiki/ASCII#Printable_characters)の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データを解読してみよう
前述の[ニホンザルの腸内細菌叢の研究](https://rdcu.be/2KC6)から。

>@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日、了。）