fastqファイルのクオリティチェック（fastq_illumina_filterを使ってみた）

illuminaのCASAVA-1.8 pipeline（Version 1.8）で生成されたfastqファイルは、クオリティの低いリードには「Y」、クオリティの良いリードには「N」が与えられています。（クオリティの良し悪しはCASAVA filteringにより落ちたリードかどうかで判断。フィルタリングの具体的な方法はよくわかっていない。）

fastqファイルは１つの配列の情報が４行で記述されており、CASAVA filterから落ちたリードかどうかはその１行目に書かれています。（fastqファイルのデータの読み方についてはwikiを参照。）

NCBIの「Gene Expression Omnibus（GEO）」か、DDBJの「Sequence Read Archive（SRA）」からダウンロードしたsraファイルを「SRA Toolskit」でfastqファイルに変換すると、illuminaのCASAVA pipelineから生成されたfastqファイルであっても、CASAVA filterの情報（NもしくはY）が抜け落ちています。（fastq-dumpコマンドを実行すると、自動的にクオリティの低いリードはフィルタリングされているよう。そのため、あまり気にする必要がないのかも。）

そのため、GEOやSRAのような「Public database」からダウンロードしてきたsraファイル（fastqファイル）については、以下の作業を行う必要はないと思います。

ただ、illuminaのCASAVA-1.8 pipelineを使って自前でfastqファイルを生成した場合では、以下のCASAVA filteringでクオリティが低いと判断されたリード（「Y」と記載）を排除したほうが良いと思います。

〈fastq_illumina_filter-0.1のインストール〉
（１）wgetコマンドで指定したサイトからファイルをダウンロードする。
$ wget http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/fastq_illumina_filter-0.1.tar.gz
⇒もしくは以下のサイトから、「Download」の「Source code」にあるファイルをダウンロード。
http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/

（２）ファイルを解凍する。
$ tar -xzf fastq_illumina_filter-0.1.tar.gz

（３）解凍したファイルのディレクトリに移動する。
$ cd fastq_illumina_filter-0.1

（４）コンパイルを実行する。
$ make
g++ -O2 -g -Wall -Wextra -Werror    fastq_illumina_filter.cpp   -o fastq_illumina_filter

（５）パスの通っているディレクトリへ、プログラム（fastq_illumina_filter）を移動する。
$ sudo cp fastq_illumina_filter /usr/local/bin
[sudo] password for imamachi:
⇒一時的に管理者権限を利用する。

（６）インストールされているか、実際にコマンドを実行してみる。
$ fastq_illumina_filter
fastq_illumina_filter (version 0.1) - Filters a FASTQ file generated by CASAVA 1.8
Copyright (C) 2011 - A. Gordon <gordon@cshl.edu> - Released under AGPLv3 or later.

Usage: fastq_illumina_filter  [--keep Y/N] [-NYhv] [-o OUTPUT] [INPUT]

   [INPUT]   = Input file. Reads from STDIN if no file is specified.
   [-o OUTPUT] = Output file. Default is STDOUT if no file is specified.
   --keep N  = Keep reads that were NOT filtered.
               (Reads that have 'N' in the read-ID line.)
   --keep Y  = Keep reads that were filtered.
               (Reads that have 'Y' in the read-ID line.)
   -N      = same as '--keep N'
   -Y      = save as '--keep Y'
   -v      = Report read counts to STDERR
             Use twice to show progress while procesing the file.
             When combined with '-o FILE', report goes to STDOUT.
   -h      = This helpful help screen

!! NOTE ABOUT Y/N !!:
  Reads that were filtered (have 'Y' in the read-ID)
  are the LOW QUALITY reads. You most likely DO NOT want them.
  (In previous CASAVA version, those are reads that have not passed filtering)

  Reads that were NOT filtered are the better-quality reads.
  So using '-N' or '--keep N' is probably the option you want to use.

FASTQ files from CASAVA-1.8 Should have the following READ-ID format:
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read>:<is filtered>:<control number>:<index sequence>

Example 1:
  $ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | head -n 4
  @EAS139:136:FC706VJ:2:5:1000:12850 1:N:18:ATCACG
  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
  +
  BBBBCCCC?<A?BC?7@@???????DBBA@@@@A@@

  $ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | fastq_illumina_filter -vvN | gzip > good_reads.fq.gz
  Processed 4,000,000 reads
  fastq_illumina_filter (--keep N) statistics:
  Input: 4,000,000 reads
  Output: 3,845,435 reads (96%)

Example 2 (The input and output files are not compressed):
  $ fastq_illumina_filter --keep N -v -v -o good_reads.fq NA10831_ATCACG_L002_R1_001.fastq
  Processed 4,000,000 reads
  fastq_illumina_filter (--keep N) statistics:
  Input: 4,000,000 reads
  Output: 3,845,435 reads (96%)

〈使用例〉
$ fastq_illumina_filter -N -v -o OUTPUT_file.fastq INPUT_file.fastq

・オプション
-N: CASAVA filterで落ちなかったリードをOUTPUT_file.fastqに出力。
-v: リード数の情報などを画面に出力。（STDERR）
など

〈参考文献〉
・Fastq - Wikipedia
http://ja.wikipedia.org/wiki/Fastq