【seqkitの使い方】FASTQファイルをダウンサンプリングする

FASTQファイルは一般に数GBを超えるような大容量のファイルになります。

データサイズが大きすぎると、取り回しが大変なので、リード数の少ないデータにして解析に使用したい場合があります。このようなリードのダウンサンプリング(サブサンプリング)は、データ全体からランダムにサンプリングすることで実現できます。

本記事では、seqkitを使用して、FASTQファイルをダウンサンプリングする方法について解説します。

MacBook Air (Apple M1チップ)
macOS Ventura 13.4
SeqKit : v2.3.0

seqkitとは

seqkitは、FASTQファイルやFASTAファイルを操作するためのツールであり、軽量かつ高速に動作します。

公式のドキュメントは以下のリンクを参照ください。

seqkitの使い方

seqkitのインストール方法は簡単で公式ページからバイナリをダウンロードして使用することもできますが、本記事ではdockerイメージを使ってseqkitを利用します。dockerを使える準備ができていればそのままお進みください。dockerがまだ使える状態でない方は下記の記事からdockerをインストールしてから先に進みましょう!

Macの場合

バイオインフォマティクスを始めるときに、初心者が最初にぶつかる壁が環境構築です。解析を行うためには非常にたくさんのソフトウェアをインストールして使わないといけないケースが多くあります。 この記事では、簡単にバイオイ[…]

Windowsの場合

バイオインフォマティクスを始めるときに、初心者が最初にぶつかる壁が環境構築です。解析を行うためには非常にたくさんのソフトウェアをインストールして使わないといけないケースが多くあります。 この記事では、Windows[…]

FASTQファイルの準備

まず、FASTQファイルを準備します。ターミナルで以下のコマンドを実行し、SRR5208824をダウンロードします。

docker container run --rm -v $PWD:/output -w /output pegi3s/sratoolkit:3.0.5 fasterq-dump -e 2 -p SRR5208824

FASTQファイルのサイズが4.3GBあるため、ダウンロードには少し時間がかかりますが、気長に待ちましょう!

もし、スレッド数に余裕があれば、-e オプションのスレッド数を増やすと速くなります。

seqkitのdockerイメージ

dockerイメージは下記のリンクのものを使用します。M1 MacのCPUアーキテクチャであるarm64に対応しているイメージを選択しました。

seqkitのdockerイメージ

seqkitのdockerイメージ

FASTQファイルの入出力と圧縮

seqkitは、pgzipパッケージを内部で使用しており、gzip圧縮されたファイルの入出力に対応しています。pgzipを利用しているため、gzipやpigzといったツールよりも高速にファイルの読み書きができると公式ページに記載されています。

fasterq-dumpでダウンロードされたFASTQファイルは圧縮されていないため、seqkitを使って、FASTQファイルを圧縮します。

ターミナルで以下のコマンドを実行してください。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
seqkit seq SRR5208824.fastq -o SRR5208824.fastq.gz

ここで、バックスラッシュは改行をエスケープするためのコマンドなので、全てを1行で記載して実行しても問題ありません。

seqkitのseqコマンドの-o オプションでgzip圧縮した出力ファイルを指定しています。4.3GBあったFASTQファイルが800MB程度に圧縮できました。

公式にpigzよりも速いと記載されていたので、gzip, pigz, seqkitで実際に圧縮にかかった時間を比較してみました。

gzip 290秒
pigz 56秒
seqkit 35秒

確かにseqkitは高速に圧縮できていました。

先頭のリードを表示する

seqkit headコマンドで、FASTQファイルの先頭のリードを簡単に確認できます。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
seqkit head -n 3 SRR5208824.fastq.gz

以下のように先頭3リードの情報が出力されます。

@SRR5208824.1 1 length=75

AAGTGNTAGCTGAAATGTATGCTTTAGCCNTGCNCNNNNNNNAGGCANAATATCTCACCATAAATTGGCATGCCA

+

/AAAA#EEEEEEEAEEEEEEEEEEEEEEE#EEE#E#######EEAEE#EEEEEEEAAEEEEEEEEEEEEEEEEEA

@SRR5208824.2 2 length=75

ACCAANGTGTGTGTGCTTTAGTCCTTCTTAGAANGNNNNNANAATACTCACAGGAGCAAATATGGAGATAAAGTC

+

/A/AA#EAEEEEEEEEEAEE/EAA//EAE/AAE#E#####E#/EEEEEEEEEEEE6EEEEEE/E/EEEEEEE/A/

@SRR5208824.3 3 length=75

AGGAANTGGGAAAGTTTCAGTGTCAACTGTTGCNANNNNTANCACACAGAATAAGAAAGCACTGGAAATGGAGAG

+

AAAAA#/EEEEEEEEEEEEEEEEEEEEEEEEEE#A####EE#EEEEEEE/EEAEEEEEEEEEEAEEAEEEEEEAE

先頭何リードを表示させるかは -n オプションで指定でき、デフォルトは10となっています。

統計量を確認する

seqkit statコマンドで、FASTQファイルのリード数やリード長などの統計値を確認できます。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
seqkit stat -a SRR5208824.fastq.gz

FASTQファイルが複数ある場合は、以下のようにワイルドカードを使って複数指定して同時に確認することもできます。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
seqkit stat -a *.fastq.gz

3つのFASTQファイルが存在していた場合の出力結果の例を以下に示しています。リード数や全長、平均長、GC含量など、FASTQファイルに関する様々な統計値を確認することができます。-aオプションを付けない場合は出力項目が減ります。

seqkit stat出力結果

ダウンサンプリングする

seqkit sample コマンドでランダムサンプリングすることができます。サンプリングする量は、割合で指定する場合は-pオプション、リード数で指定する場合は-nオプションで指定できます。

容量の大きいFASTQファイルに対して、-nオプションでリード数を指定するコマンドを実行するとメモリ不足でパソコンがフリーズするので、使用しないようにしましょう。

ダウンサンプリングしたい場合は、以下のように実行します。5万リードをランダムサンプリングします。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
sh -c "seqkit sample -s 77 -p 0.1 SRR5208824.fastq.gz | seqkit head -n 50000 -o downsample.fastq.gz"

| (パイプ)を使って一連の処理を実行します。コンテナ内でパイプを使う場合は、上記のように、sh -c “コマンド | コマンド” のように記載することで実行できます。-sオプションで指定している値はランダムサンプリングする時のシード値で、この値を同じにすることで何度やっても同じサンプリング結果を得ることができます。

seqkit statコマンドを使って、ダウンサンプリングできているか確認してみます。

docker container run --rm -v $PWD:/output -w /output bnovak32/seqkit:v2.3.0 \
seqkit stat -a SRR5208824.fastq.gz downsample.fastq.gz

実行結果は下記になります。

seqkitダウンサンプリング結果

num_seqsの列を確認すると、元々、2000万リードあったものから、5万リードにダウンサンプリングできていることがわかります。

おわりに

この記事では、seqkitを使用してFASTQファイルをダウンサンプリングする方法を解説しました。

seqkitはFASTQファイルだけではなく、FASTAファイルも扱うことができ、非常に多機能なツールとなっています。

ぜひ、公式のドキュメントも参考して他の機能も試してみてください。