fastqファイルやfastaファイルとしてシークエンスデータを持っているときに、より小さいデータ、リード数の少ないデータにしたい場合がある。このようなリードのダウンサンプリング（あるいは、サブサンプリング）は、データ全体からランダムに一部を取り出すことによって行われる。小さいデータセットを作りたくなる理由は色々あるだろうが、思いつくところでは、以下のようなものが挙げられる。
- [QCなど] データの性質についてQCを行いたいときに、全部のデータでは負荷がかかるので、一部を取り出して傾向を見たいと言ったとき
- [条件検討など] 解析の結果がデータ量で変化するかどうか見てみたいとき
- [ツール側の問題] 用いるプログラムやスクリプト、解析環境が処理できるサイズにデータ量を変更してやらないといけないとき
- [手軽なtest fileを作りたいとき] ツールの使い方を確かめたり、オプションを試したりしたいときに取り回しの良いファイルを作りたい時がある
Fastqファイルやfastaファイルから、指定した割合やリード数だけリードをサンプリングしたい場合にどうするのかをまとめた。何通りかやり方はあるだろうが、ここではバイオインフォマティクス ツールのseqkitを利用することを想定した。また、同じくよく使われるseqtkのコマンドを併記した。
fastaやfastqファイルを操作するためのツールで、バイナリファイルをダウンロードする他、condaでも導入できるようになっている。
seqtk: https://github.com/lh3/seqtk
こちらもfastq/fastaを操作するツールで、condaでも導入可能。
リサンプリングの前に、もしそのデータについて概要を知らない場合、リード数などの統計量をあらかじめ知っておきたいかもしれない。その場合は、seqkit statsを使って統計量を出すことができる。
$ seqkit stats input.fastq.gz
複数のファイルがある場合は、以下のように指定してもいいかもしれない。
$ seqkit stats *.fastq
fastaとfastqが混じっていてもいいらしい。
$ seqkit stats *.fa*
結果は以下のように示される。
$ seqkit stats bud_L001.fastq.gz
file format type num_seqs sum_len min_len avg_len max_len
bud_L001.fastq.gz FASTQ DNA 29,342,622 2,195,350,462 50 74.8 75
(1) seqkitを使う場合のコマンド: seqkit sample
bud_L001.fastq.gzから、100リードをサンプルして、bud_L001.n100.fastq.gzとして出力する場合、
$ seqkit sample -n 100 -s 100 bud_L001.fastq.gz -o bud_L001.n100.fastq.gz
実際に確かめると、100リードのファイルができている。
$ seqkit stats bud_L001.n100.fastq.gz
ここで-sで指定されているのは、乱数を発生させるときのrandom seedで、デフォルトは11。同じrandom seedを使えば同じ配列が得られる。
single readではなく、PEのデータの場合は、R1とR2で同じrandom seedを使う必要がある。
$ seqkit sample -n 100 -s 100 seq_R1.fastq.gz -o seq_R1.fastq.gz
$ seqkit sample -n 100 -s 100 seq_R2.fastq.gz -o seq_R2.fastq.gz
使えるスレッドに余裕がある場合は、-jオプションでスレッド数を増やすことができる
$ seqkit sample -100 -s 100 -j 8 bud_L001.fastq.gz -o bud_L001.n100.fastq.gz
$ seqtk sample -s 100 bud_L001.fastq.gz 1000 > subseq_bud.fastq
PEの場合、やはり-sを同一に指定する。
$ seqtk sample -s 100 bud_R1.fastq.gz 1000 > subseq_bud_R1.fastq
$ seqtk sample -s 100 bud_R2.fastq.gz 1000 > subseq_bud_R2.fastq
(1) seqkitを使う場合：seqkit sampleで、オプション-p (0,1)を使う
$ seqkit sample -p 0.01 -s 100 bud_L001.fastq.gz -o bud_L001.n100.fastq.gz
$ seqtk sample -s 100 bud_L001.fastq.gz 0.01 > subseq_bud.fastq
- 元ファイルに含まれている以上のリード数を指定した場合、出力ファイルは元のファイルと同一になる。
- proportionで1以上を指定するとエラーになる。
- seqtkで読み込んで標準出力すると、ファイルのフォーマットはASCII txtになる (うっかり.gzとか拡張子をつけてもcompressされていない）
- 多分読み込みにかなりのメモリを食う作業なので、下手なパソコンで始めるとメモリ不足でフリーズする。
