はじめに

FASTX toolkitは、ショートリードのfastqファイルの前処理に使用されるコマンドラインツールの集合です。
低クオリティーのリードを除去したい場合や、クオリティーを基準に塩基をトリミングしたい際等に使用されます。

似たようなツールとしては、trimmomaticや、seqtk、cutadaptなんかが挙げられます。

以下では最良のトリミングツールについて議論が繰り広げられていました。 https://www.biostars.org/p/178647/

Trimmomaticはシンプルで速い
seqtkは並列処理が可能だがアダプター除去ができない。
最近ではシーケンサーの精度が上がっていてそもそもトリミングする必要がない。
アダプターのトリミングはしたほうが良いよ。
trim_galore はcutadaptのラッパー
cutadaptは速い使いやすいドキュメントわかりやすい。
BBDuk とseqtkはアルゴリズム同じ。

なんか俺のツールめっちゃいいよ！って人と根拠もなくこれ使ってます。イイヨ。って人がいますね。

結論としては「アダプター除去を正確に行えるツールでアダプター除去だけしましょう」で良いのかな。

そんなわけでfastx-toolkitで

クオリティーチェック
リードのフィルタリング
塩基のトリミング
アダプター除去

するやり方を見ていきましょう（強引）。

インストール

コンパイル済みのバイナリが配布されてます。以下の通りに行いましょう。

$ wget http://hannonlab.cshl.edu/fastx_toolkit/fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2
$ tar -xjf fastx_toolkit_0.0.13_binaries_Linux_2.6_amd64.tar.bz2
$ sudo cp ./bin/* /usr/local/bin

使用法

それでは実際の使用法を見ていきましょう。
サンプルは適当に探した以下のfastqを使用することにします。
http://ddbj.nig.ac.jp/DRASearch/study?acc=DRP000011

クオリティーチェック

fastx_quality_statsはリードごとのクオリティーを一覧表示します。

fastx_quality_stats -i DRR000031.fastq -o DRR000031_stats.txt

出力は以下

column  count   min     max     sum     mean    Q1      med     Q3      IQR     lW      rW      A_Count C_Count G_Count T_Count N_Count Max_count
1       4589774 0       40      178945470       38.99   40      40      40      0       40      40      1281908 1041157 1168850 1096774 1085    4589774
2       4589774 0       40      177402660       38.65   40      40      40      0       40      40      1317538 973213  1159883 1112625 26515   4589774
3       4589774 0       40      176764302       38.51   40      40      40      0       40      40      1320776 1020643 1118013 1101889 28453   4589774
4       4589774 0       40      176334634       38.42   40      40      40      0       40      40      1312689 1011467 1160762 1080964 23892   4589774
・
・
・

このファイルを使って、クオリティーの可視化を行うこともできます。

↓クオリティーのボックスプロット

$ fastq_quality_boxplot_graph.sh -i DRR000031_stats.txt -o DRR000031_quality.png -t "My Library"

f:id:kimoppy126:20181022232350p:plain

↓塩基数の分布

$ fastx_nucleotide_distribution_graph.sh -i DRR000031_stats.txt -o DRR000031_nuc.png -t "My Library"

f:id:kimoppy126:20181022232430p:plain

リードのフィルタリング

リードのフィルタリングを行う際には、fastq_quality_filterを使用します。
基本的にヘルプメッセージを見れば使い方がなんとなくわかるようになっています。

$ fastq_quality_filter -h
usage: fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

   [-h]         = This helpful help screen.
   [-q N]       = Minimum quality score to keep.
   [-p N]       = Minimum percent of bases that must have [-q] quality.
   [-z]         = Compress output with GZIP.
   [-i INFILE]  = FASTA/Q input file. default is STDIN.
   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.
   [-v]         = Verbose - report number of sequences.
                  If [-o] is specified,  report will be printed to STDOUT.
                  If [-o] is not specified (and output goes to STDOUT),
                  report will be printed to STDERR.

入力ファイルを-i 、出力ファイルを-oで与えます。
-qで残すクオリティーの下限、-pでクオリティーを持つ塩基数の下限を指定します。

例えば、クオリティーを持った値が90%以上、クオリティースコアが20以上のものだけを出力したい場合、

fastq_quality_filter -i DRR000031.fastq -o DRR000031.fastq.filtered -p 90 -q 20

fastqファイルのクオリティ

fastqファイルには、Sanger形式と、llumina 形式の2種類が存在します。

Sanger形式ではPhredクオリティスコアの0から93の値は、ASCIIコードでは33から126の間の文字としてエンコードされる。 Fastq - Wikipedia

Sanger形式のfastqファイルでは、offset値として33が使用されています。

llumina 1.3+ 形式ではPhredクオリティスコアの 0 から 62 を ASCII の64 から 126でエンコードする。 Fastq - Wikipedia

一方、Illumina形式のfastqファイルでは、offset値として64が使用されています。

これらの違いがあるということは、クオリティーのthreshold値をfastx-toolkitで指定する際に、どちらのエンコーディングを使用するか、意識する必要があります。

Sanger形式が一般に用いられているが、fastx-toolkitのデフォルト値はoffset値64です。fastx-toolkitは2009年で開発が止まっているから仕方のないことなのです。

offset値を現在出回っているSanger形式のエンコーディングに対応させるには、 -Q33をオプションで与えてあげる必要があります。

fastq_quality_filter -i DRR000031.fastq -o DRR000031.fastq.filtered -p 90 -q 20 -Q33

フィルタリング前後のリード数を比較すると、確かに減っていることがわかります。

$ wc -l DRR000031.fastq
18359096 DRR000031.fastq
$ wc -l DRR000031.fastq.filtered
13571236 DRR000031.fastq.filtered

クオリティーを指定した塩基のトリミング

クオリティーに応じた塩基のトリミングを行うには、fastq_quality_trimmerを使用します。

$ fastq_quality_trimmer -h
usage: fastq_quality_trimmer [-h] [-v] [-t N] [-l N] [-z] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

   [-h]         = This helpful help screen.
   [-t N]       = Quality threshold - nucleotides with lower
                  quality will be trimmed (from the end of the sequence).
   [-l N]       = Minimum length - sequences shorter than this (after trimming)
                  will be discarded. Default = 0 = no minimum length.
   [-z]         = Compress output with GZIP.
   [-i INFILE]  = FASTQ input file. default is STDIN.
   [-o OUTFILE] = FASTQ output file. default is STDOUT.
   [-v]         = Verbose - report number of sequences.
                  If [-o] is specified,  report will be printed to STDOUT.
                  If [-o] is not specified (and output goes to STDOUT),
                  report will be printed to STDERR.

入力ファイルを-i 、出力ファイルを-oで与えます。
-tでクオリティーの閾値を指定します。

例えば、クオリティーが20 以上の塩基配列のみ残したい場合、以下のように指定します。

fastx_trimmer -t 20 -i DRR000031.fastq-o DRR000031.fastq.q20

場所を指定した塩基のトリミング

クオリティーに関わらず、リードの決まった箇所を取り除くには、fastx_trimmerを使用します。

$ fastx_trimmer -h
usage: fastx_trimmer [-h] [-f N] [-l N] [-t N] [-m MINLEN] [-z] [-v] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

   [-h]         = This helpful help screen.
   [-f N]       = First base to keep. Default is 1 (=first base).
   [-l N]       = Last base to keep. Default is entire read.
   [-t N]       = Trim N nucleotides from the end of the read.
                  '-t'  can not be used with '-l' and '-f'.
   [-m MINLEN]  = With [-t], discard reads shorter than MINLEN.
   [-z]         = Compress output with GZIP.
   [-i INFILE]  = FASTA/Q input file. default is STDIN.
   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.

入力ファイルを-i 、出力ファイルを-oで与えます。 -fで何塩基から残すかを、-lで何塩基まで残すかを指定します。

例えば1塩基目から32塩基までを残したい場合、以下のように指定する。

fastx_trimmer -f 1 -l 32 -i DRR000031.fastq.filtered -o DRR000031.fastq.out32

アダプターのトリミング

アダプターのトリミングを行うには、fastx_clipperを使用します。

$ fastx_clipper -h
usage: fastx_clipper [-h] [-a ADAPTER] [-D] [-l N] [-n] [-d N] [-c] [-C] [-o] [-v] [-z] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.13 by A. Gordon (gordon@cshl.edu)

   [-h]         = This helpful help screen.
   [-a ADAPTER] = ADAPTER string. default is CCTTAAGG (dummy adapter).
   [-l N]       = discard sequences shorter than N nucleotides. default is 5.
   [-d N]       = Keep the adapter and N bases after it.
                  (using '-d 0' is the same as not using '-d' at all. which is the default).
   [-c]         = Discard non-clipped sequences (i.e. - keep only sequences which contained the adapter).
   [-C]         = Discard clipped sequences (i.e. - keep only sequences which did not contained the adapter).
   [-k]         = Report Adapter-Only sequences.
   [-n]         = keep sequences with unknown (N) nucleotides. default is to discard such sequences.
   [-v]         = Verbose - report number of sequences.
                  If [-o] is specified,  report will be printed to STDOUT.
                  If [-o] is not specified (and output goes to STDOUT),
                  report will be printed to STDERR.
   [-z]         = Compress output with GZIP.
   [-D]  = DEBUG output.
   [-M N]       = require minimum adapter alignment length of N.
                  If less than N nucleotides aligned with the adapter - don't clip it.   [-i INFILE]  = FASTA/Q input file. default is STDIN.
   [-o OUTFILE] = FASTA/Q output file. default is STDOUT.

色んなオプションがありますが、実際使うのは-aくらい。

例えば、CTAATCGAというアダプター配列を除去したい場合、以下のように指定します。

fastx_clipper　-a CTAATCGA　-i DRR000031.fastq -o DRR000031.noadapter.fastq

ばいばいバイオ

配列のクオリティーコントロール - FASTX-toolkit