2020年6月9日火曜日

Fastq/Fasta からのリードのリサンプリング




fastqファイルやfastaファイルとしてシークエンスデータを持っているときに、より小さいデータ、リード数の少ないデータにしたい場合がある。このようなリードのダウンサンプリング（あるいは、サブサンプリング）は、データ全体からランダムに一部を取り出すことによって行われる。小さいデータセットを作りたくなる理由は色々あるだろうが、思いつくところでは、以下のようなものが挙げられる。


  • [QCなど] データの性質についてQCを行いたいときに、全部のデータでは負荷がかかるので、一部を取り出して傾向を見たいと言ったとき
  • [条件検討など] 解析の結果がデータ量で変化するかどうか見てみたいとき
  • [ツール側の問題] 用いるプログラムやスクリプト、解析環境が処理できるサイズにデータ量を変更してやらないといけないとき
  • [手軽なtest fileを作りたいとき] ツールの使い方を確かめたり、オプションを試したりしたいときに取り回しの良いファイルを作りたい時がある


Fastqファイルやfastaファイルから、指定した割合やリード数だけリードをサンプリングしたい場合にどうするのかをまとめた。何通りかやり方はあるだろうが、ここではバイオインフォマティクス ツールのseqkitを利用することを想定した。また、同じくよく使われるseqtkのコマンドを併記した。




用いるツール



fastaやfastqファイルを操作するためのツールで、バイナリファイルをダウンロードする他、condaでも導入できるようになっている。



こちらもfastq/fastaを操作するツールで、condaでも導入可能。



Fastq/fasta ファイルの内容確認



リサンプリングの前に、もしそのデータについて概要を知らない場合、リード数などの統計量をあらかじめ知っておきたいかもしれない。その場合は、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



(2) seqtkを使う場合のコマンド：seqtk sample
$ 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



(2) seqtkを使う場合：seqtk sampleで、割合を指定する

 $ seqtk sample -s 100 bud_L001.fastq.gz 0.01 > subseq_bud.fastq




注意点
  • 元ファイルに含まれている以上のリード数を指定した場合、出力ファイルは元のファイルと同一になる。
  • proportionで1以上を指定するとエラーになる。
  • seqtkで読み込んで標準出力すると、ファイルのフォーマットはASCII txtになる (うっかり.gzとか拡張子をつけてもcompressされていない）
  • 多分読み込みにかなりのメモリを食う作業なので、下手なパソコンで始めるとメモリ不足でフリーズする。
ラベル: , , , , , ,