macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

fastq / fastaの操作ツール seqkit

 

2016年に発表されたfastqの操作ツール。競合ツールより多機能とされる。seqtkと同様、動作は非常に早い。メモリ使用量はseqtkより少ないとされる。

 

  

マニュアル

Usage - SeqKit - Ultrafast FASTA/Q kit

チュートリアル

Tutorial - SeqKit - Ultrafast FASTA/Q kit

 

公式マニュアルの機能比較。競合のツールより多機能に見える。

f:id:kazumaxneo:20170801220752j:plain

公式マニュアルより。

 

  

インストール 

 

ダウンロード

http://bioinf.shenwei.me/seqkit/download/

解凍するとbinaryファイルできる。

実行権をつけ、それを/usr/local/bin/などに移動して使う。

chmod u+x seqtk #実行権
ls -l seqtk #確認
mv seqtk /usr/local/bin #パスがとったところに移動。ln -sしても同じ。

 

ラン

入力はfasta、fastq、.gzなど。

 

seq 配列を変形

reverse complementary (相補鎖にしてさらにリバースにする)

seqkit seq -pr input.fa > output.fa
  • -p      complement sequence (blank for Protein sequence)
  • -r                      reverse sequence
  • --dna2rna         DNA to RNA
  • --rna2dna         RNA to DNA
  • -l                  print sequences in lower case
  • -i                    print ID instead of full head
  • -u           print sequences in upper case
  • -n           only print names
  • -m                    only print sequences longer than the minimum length (-1 for no limit) (default -1)
  • -M,             only print sequences shorter than the maximum length (-1 for no limit) (default -1)

ヘッダーだけ取り出すことでcontigの数をカウントする。

seqkit seq -n contig.fa |wc -l

fastqからヘッダーのIDだけ取り出す。

seqkit seq -ni fastq |less

--id-regexpを使うと、標準正規表現で仕切りを自由にカスタム可能(リンク)。

ギャップを除く。("-"が除かれる)。

seqkit seq -g scaffolds.fa > output.fasta

1000bp以上の配列だけ取り出す。

seqkit seq -m 1000 scaffolds.fa > output.fasta

1000bp以下の配列だけ取り出し、statsコマンドに渡して分析。

seqkit seq -M 1000 scaffolds.fa |seqkit stats

他にもいくつかのオプションがある。

 

 

subseq 特定のリードや配列だけ抽出

最初の100-bpを取り出す。

seqkit subseq -r 1:100 contigs.fa > left.fa

 最後の100-bpを取り出す(マイナスで後ろからカウント)。 

seqkit subseq -r -100:-1 contig.fa > right.fa

 最初と最後の100-bpを除く。

seqkit subseq -r 100:-100 contig.fa > center.fa

 gtfで指定した領域の配列だけ取り出す。

seqkit subseq --gtf cDNA.gtf genome.fa > cdna.fa

 gtfで指定した領域の配列とその100bp上流までを取り出す。

seqkit subseq --gtf cDNA.gtf -u 100 genome.fa > cdna_and100bp_upstream.fa

-fもつけると、上流配列だけ取り出される。

bedファイルの指定領域かつchr1の配列のみ取り出す。

seqkit subseq --bed Homo_sapiens.GRCh38.bed.gz --chr 1 input.fa >  chr1_bed.fa

bedtoolsのツールと似ているが、動作が高速でメモリ使用量も少ないので、こちらの方が使いやすい。

 

 

sliding 配列を指定の単位ずつずらしながら取り出す。

100-bp windowでゲノム配列を切り出す(GC含量カウントに使ったりする)。

seqkit sliding -s 1 -W 100 genome.fa > 100bp_windows.fa

1bpずつずらしながら、100bpの配列が取り出される。

環状ゲノムなら-Cをつける。

GC含量をカウント

seqkit fx2tab genome.fa| head -n 1 | seqkit tab2fx | seqkit sliding -s 1 -W 100 | seqkit fx2tab -n -g > GC.txt

ウィンドウサイズ100でのGC率が出力される。平均GC率などを求めるなら、続けて以下のように打つ。

cut -f 4 GC.txt > GC #GC列を抽出
st --mean GC

 stコマンドはbrewで導入できます。

 

 

 

stats fastaの簡単な分析を行う。

カレントディレクトリの全fa.gzファイルを分析。

seqkit stats *.fa.gz

user$ seqkit stats R1.fastq 

file               format  type  num_seqs     sum_len  min_len  avg_len  max_len

T1second_R1.fastq  FASTQ   DNA    261,774  77,786,418       35    297.2      301

-aをつけるとより詳細な情報を表示。

user$ seqkit stats -a R1.fastq 

file               format  type  num_seqs     sum_len  min_len  avg_len  max_len  sum_gap  N50

T1second_R1.fastq  FASTQ   DNA    261,774  77,786,418       35    297.2      301        0  301

 

 

fq2fa fastqをfastaに変換

seqkit fq2fa input.fastq > output.fa

gz圧縮ファイルにも対応。

seqkit fq2fa input.fastq.gz > output.fa.gz

 

 

 

fx2tab fasta/fastqをtab形式に変換し、さらにGCなどの情報を表示。

seqkit fx2tab input.fa > output.fa
  • -g print GC content
  • -G print GC-Skew
  • -H print header line (先頭行のコメント出力)
  • -l print sequence length
  • -n only print names (no sequences and qualities)
  • -i print ID instead of full head

出力はタブ区切り形式になっている。 fastqならクオリティも後ろについている。

f:id:kazumaxneo:20170808221310j:plain

GC-skew、GC含量、長さも出力。配列は出力しない。

seqkit fx2tab -Ggln input.fa > output.fa

出力。名前の後ろの列に長さ、GC率、GC-skewがプリントされている。

f:id:kazumaxneo:20170808221724j:plain

タブからfastqなどに戻すtab2fxというコマンドもある。

 

 

head 指定数だけ先頭からリードを取り出す。

seqkit head -10 input.fastq > 10read.fastq

先頭から10リード(つまり40行)出力される。

 

 

rename duplicateしたIDをユニークな名前に修復する。

seqkit rename input.fastq > renamed.fastq
  •  -n check duplicated by full name instead of just id.

duplicateしていると判断されると"_2"とか名前につく。

 

 

restart 環状ゲノムのスタート位置を変える。

seqkit restart -i 100 genome.fa > genome_restart.fa
  • -i new start position (1-base, supporting negative value counting from the end) (default 1).

この例だと100番目の塩基がスタートになる(元々の99塩基は最後につく)。

 

 

shuffle 順番をバラバラにする。

seqkit shuffle input.fastq > shuffled.fastq
  • -2 two-pass mode read files twice to lower memory usage. (only for FASTA format).
  • -s rand seed for shuffle (default 23)

一度メモリーに読み込むので、シーケンスデータが大きすぎるとランできない恐れあり。fasta限定だが、-2をつけるとtwo-passにしてメモリ使用量を削減可能。ペアリードをシャッフルするなら、R1とR2それぞれ-sの値を同じ値にすること。

 

 

sort ソート。

長さでソート。

seqkit sort --quiet -l conitg.fa > sorted.fa

 

  • -l by sequence length
  • -n by full name instead of just id
  • -s by sequence
  • -r reverse the result
  • -L length of sequence prefix on which seqkit sorts by sequences (0 for whole sequence) (default 10000)
  • -2 two-pass mode read files twice to lower memory usage. (only for FASTA format)
  • -i ignore case

配列でソート。

seqkit sort --quiet -s conitg.fa > sorted.fa

長さでソートして、名前と長さ、ヘッダー行だけタブ出力。

seqkit sort --quiet -l conitg.fa | seqkit fx2tab -l -n -i -H > length.txt

 

 

rmdup duplicationを除去。

フルネームが重複したリードを除去。

seqkit rmdup -n input.fq 
  • -n by full name instead of just id
  • -s by seq
  • -D file to save number and list of duplicated seqs
  • -d string file to save duplicated seqs
  • -i ignore case

簡単な動作テストはseqkitのdupコマンドでfastqを複製して行えば良い。

seqkit dup -n 2 input.fastq > duplicated.fastq #全リードを1回ずつコピー。ファイルサイズが倍はなる。
seqkit rmdup -n duplicated.fastq > dup_removed.fastq

user$ seqkit dup -n 2 input.fastq > dup.fq 

user$ seqkit rmdup -n dup.fq > dup_removed.fastq 

[INFO] 261774 duplicated records removed

 

 

common 複数ファイルの共通する配列や名前を分析

カレントディレクトリのfastaファイルから、IDが共通する配列を抽出。

seqkit common *fa > common.fa
  • -n match by full name instead of just id
  • -s match by sequence
  • -i ignore case

カレントディレクトリのfastqファイルから、配列が一致する配列を抽出。

seqkit common -s *fq > common_seq.fq

 

 

sample リードをランダムサンプリング

リードを10%ランダムサンプリングする。

seqkit sample -p 0.1 input.fastq > output.fastq
  • -n int sample by number (result may not exactly match)
  • -p float sample by proportion
  • -s  int rand seed (default 11)
  • -2 2-pass mode read files twice to lower memory usage. Not allowed when reading from stdin

-sでランダムさを制御しているので、ペアリードは-sを同じ値にして抽出する。

 

 

split 分割する。

fastqを均等に10ファイルに分けて出力。

seqkit split -p 10 input.fastq
  • -p split squences into N parts

input.fastq.split/ができ、その中に10ファイルにfastqが分割され出力される。

 巨大なデータを分割して渡したり、サーバーからダウンロードする時に使えるかもしれない。

 

 

 

他にもいくつかのコマンドがあるが、動作は1リードを1行として扱う以外unixの同名コマンドと同じ感じで使える。

  • locate 配列をリプレースする。
  • grep   特定の配列を検索する。
  • dup    配列を指定回、複製させる。

 

紹介したコマンドや引数は全てではありません。使えそうなコマンドがあれば、公式マニュアルで確認してみてください。

Usage - SeqKit - Ultrafast FASTA/Q kit

 

 

似たコマンドにseqtkがある。seqtkも以前紹介しています。

 

引用

SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation

Wei Shen, Shuai Le, Yan Li , Fuquan Hu 

PLOS ONE. 2016. doi:10.1371/journal.pone.0163962.