シーケンサーから得られた配列は何かしらの処理をすることが多い。
そこで結構便利かつ軽量なソフトであるseqkitを紹介する。
公式ページ（ドキュメント）はこちら
https://github.com/shenwei356/seqkit/releases
解凍してパスを通すだけ。
ちなみに、.tar.gzの解凍は tar -zxvf *.tar.gz でできます。
以下のような疑似FASTQを例とする。 以下をexample.fqとして保存し、以下を実行する。 シーケンスの長さが10〜14,クオリティは24〜28としている。
example.fqからもダウンロード可能
@seq0
ATGCATGCAT
+
9999999999
@seq1
ATGCATGCATG
+
:::::::::::
@seq2
ATGCATGCATGC
+
;;;;;;;;;;;;
@seq3
ATGCATGCATGCA
+
<<<<<<<<<<<<<
@seq4
ATGCATGCATGCAT
+
==============
@seq5
ATGCATGCATGCATG
+
>>>>>>>>>>>>>>>
基本的には seqkit ${コマンド} ${fastqファイル} で使用することができる。
デフォルトでは標準出力に出力されるため、ファイルに吐き出したい場合はリダイレクトする必要がある。
fastqの長さに関する簡単な統計値は seqkit stats *.fq で見ることができる。
タブ区切りで出力してくれる -T オプションが便利
$ seqkit stats -T example.fq
file format type num_seqs sum_len min_len avg_len max_len
example.fq FASTQ DNA 6 75 10 12.5 15
クオリティの統計値もみたい場合は、 -aオプションをつけることでクオリティの統計値なども見ることができる。
$ seqkit stats -a -T example.fq
file format type num_seqs sum_len min_len avg_len max_len Q1 Q2 Q3 sum_gap N50 Q20(%) Q30(%)
example.fq FASTQ DNA 6 75 10 12.5 15 11 12 14 0 13 100.00 0.00
seqkit fq2fa *.fq でfastqからfastaへ変換が可能。
$ seqkit fq2fa example.fq
>seq0
ATGCATGCAT
>seq1
ATGCATGCATG
>seq2
ATGCATGCATGC
>seq3
ATGCATGCATGCA
>seq4
ATGCATGCATGCAT
>seq5
ATGCATGCATGCATG
seqkit seq -Q ${quality_value} ${path_fastq} で平均クオリティーによるフィルタリングが可能。
以下の例では、平均クオリティが25以上のリードのみを抽出している。
$ seqkit seq -Q 25 example.fq
@seq1
ATGCATGCATG
+
:::::::::::
@seq2
ATGCATGCATGC
+
;;;;;;;;;;;;
@seq3
ATGCATGCATGCA
+
<<<<<<<<<<<<<
@seq4
ATGCATGCATGCAT
+
==============
@seq5
ATGCATGCATGCATG
+
>>>>>>>>>>>>>>>
想定通り、1本目のリードが除かれ、その他リードは残った。
seqkit seq -m ${min_len} -M ${max_len} ${path_fastq} で長さによるフィルタリングが可能。
11以上、13以下のリード数を取得してみる。
$ seqkit seq -m 11 -M 13 example.fq
@seq1
ATGCATGCATG
+
:::::::::::
@seq2
ATGCATGCATGC
+
;;;;;;;;;;;;
@seq3
ATGCATGCATGCA
+
<<<<<<<<<<<<<
想定通り、1本目,5本目のリードが除かれ、その他リードは残った。
seqkit seq -n ${read_num} ${path_fastq} で本数を指定してサンプリングが可能。
ただし、 seqkitは指定した本数に近づけるだけで、指定した本数通りサンプリングしてくれるわけではない ことに注意。
$ seqkit sample -n 3 -s 1 example.fq
[INFO] sample by number
[INFO] 2 sequences outputted
@seq3
ATGCATGCATGCA
+
<<<<<<<<<<<<<
@seq4
ATGCATGCATGCAT
+
==============
本数を厳密に指定したい場合、shuffleしてheadなどが考えられる。
以下のコマンドで実行可能。
$ seqkit shuffle -s 0 example.fq | seqkit head -n 3 # -sはseed。
[INFO] read sequences ...
[INFO] 6 sequences loaded
[INFO] shuffle ...
[INFO] output ...
@seq4
ATGCATGCATGCAT
+
==============
@seq2
ATGCATGCATGC
+
;;;;;;;;;;;;
@seq3
ATGCATGCATGCA
+
<<<<<<<<<<<<<
seqkit fx2tab -n -g ${path_fastq} でGC含量の計算ができる。
$ seqkit fx2tab -n -g example.fq
seq0 40.00
seq1 45.45
seq2 50.00
seq3 46.15
seq4 42.86
seq5 46.67
Shen, W., Le, S., Li, Y. & Hu, F. SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. PLoS One 11, e0163962 (2016). ↩