macでインフォマティクス

macでインフォマティクス

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

fastq / fastaの操作ツール seqkit

fasta/fastqの操作 format conversion 2016 PLoS ONE

2019 4/15 Githubリンク追加 

2019 6/21 seqmit sample コマンド追記

2019 8/7 help追加

2019 8/8 stats追記

2020 3/18 help更新

2021 ツイート追加（対応するバージョンを使っている人は注意）

 

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

 

  

マニュアル

Usage - SeqKit - Ultrafast FASTA/Q kit

チュートリアル

Tutorial - SeqKit - Ultrafast FASTA/Q kit

 

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

f:id:kazumaxneo:20170801220752j:plain

公式マニュアルより。

 

2023/03/22

2023/03/20

2021 8/28

さらにコマンドが増えました！すごい！！

  

インストール 

Github


ダウンロード

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

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

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

chmod u+x seqtk #実行権
ls -l seqtk #確認
mv seqtk /usr/local/bin #パスが通ったディレクトリに移動。ln -sでリンクしても良い。

> seqkit #v0.12.0

$ seqkit --help

SeqKit -- a cross-platform and ultrafast toolkit for FASTA/Q file manipulation

 

Version: 0.12.0

 

Author: Wei Shen <shenwei356@gmail.com>

 

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

Source code: https://github.com/shenwei356/seqkit

Please cite: https://doi.org/10.1371/journal.pone.0163962

 

Usage:

  seqkit [command]

 

Available Commands:

  amplicon        retrieve amplicon (or specific region around it) via primer(s)

  bam             monitoring and online histograms of BAM record features

  common          find common sequences of multiple files by id/name/sequence

  concat          concatenate sequences with same ID from multiple files

  convert         convert FASTQ quality encoding between Sanger, Solexa and Illumina

  duplicate       duplicate sequences N times

  faidx           create FASTA index file and extract subsequence

  fish            look for short sequences in larger sequences using local alignment

  fq2fa           convert FASTQ to FASTA

  fx2tab          convert FASTA/Q to tabular format (with length/GC content/GC skew)

  genautocomplete generate shell autocompletion script

  grep            search sequences by ID/name/sequence/sequence motifs, mismatch allowed

  head            print first N FASTA/Q records

  help            Help about any command

  locate          locate subsequences/motifs, mismatch allowed

  mutate          edit sequence (point mutation, insertion, deletion)

  range           print FASTA/Q records in a range (start:end)

  rename          rename duplicated IDs

  replace         replace name/sequence by regular expression

  restart         reset start position for circular genome

  rmdup           remove duplicated sequences by id/name/sequence

  sample          sample sequences by number or proportion

  sana            sanitize broken single line fastq files

  seq             transform sequences (revserse, complement, extract ID...)

  shuffle         shuffle sequences

  sliding         sliding sequences, circular genome supported

  sort            sort sequences by id/name/sequence/length

  split           split sequences into files by id/seq region/size/parts (mainly for FASTA)

  split2          split sequences into files by size/parts (FASTA, PE/SE FASTQ)

  stats           simple statistics of FASTA/Q files

  subseq          get subsequences by region/gtf/bed, including flanking sequences

  tab2fx          convert tabular format to FASTA/Q format

  translate       translate DNA/RNA to protein sequence (supporting ambiguous bases)

  version         print version information and check for update

  watch           monitoring and online histograms of sequence features

 

Flags:

      --alphabet-guess-seq-length int   length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)

  -h, --help                            help for seqkit

      --id-ncbi                         FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...

      --id-regexp string                regular expression for parsing ID (default "^(\\S+)\\s?")

      --infile-list string              file of input files list (one file per line), if given, they are appended to files from cli arguments

  -w, --line-width int                  line width when outputing FASTA format (0 for no wrap) (default 60)

  -o, --out-file string                 out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

      --quiet                           be quiet and do not show extra information

  -t, --seq-type string                 sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")

  -j, --threads int                     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

 

Use "seqkit [command] --help" for more information about a command.

 

実行方法

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

 

v2になっていくつかコマンドが変化しています（v２リリース）。下に書いてあるのは、v２より前のバージョンなので、新しいバージョンを使っている方は注意して下さい。 

 

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

-Tをつけるとtabular formatで表示。チュートリアルにはcsvtkなどに渡してさらに整形する例も記載されている。

並列化と追加情報つき。csvtkに渡す

seqkit stats *.fa.gz -j 8 -a | csvtk pretty -t

file                                   format   type   num_seqs   sum_len    min_len   avg_len   max_len   Q1    Q2    Q3    sum_gap   N50   Q20(%)   Q30(%)

WS-tuber-1_S276_L001_R1_001.fastq      FASTQ    DNA    28630      8580372    35        299.7     301       299   300   301   0         300   98.00    94.16

WS-tuber-1_S276_L001_R1_001.fastq.gz   FASTQ    DNA    28630      8580372    35        299.7     301       299   300   301   0         300   98.00    94.16

WS-tuber-1_S276_L001_R2_001.fastq      FASTQ    DNA    28630      8598293    35        300.3     301       300   301   301   0         301   93.54    85.00

WS-tuber-1_S276_L001_R2_001.fastq.gz   FASTQ    DNA    28630      8598293    35        300.3     301       300   301   301   0         301   93.54    85.00

WS-tuber-2_S277_L001_R1_001.fastq.gz   FASTQ    DNA    28547      8555935    35        299.7     301       299   300   301   0         300   98.00    94.24

WS-tuber-2_S277_L001_R2_001.fastq.gz   FASTQ    DNA    28547      8575923    35        300.4     301       300   301   301   0         301   92.74    83.63

WS-tuber-3_S278_L001_R1_001.fastq.gz   FASTQ    DNA    35579      10656382   35        299.5     301       299   300   301   0         300   98.13    94.37

WS-tuber-3_S278_L001_R2_001.fastq.gz   FASTQ    DNA    35579      10680068   35        300.2     301       300   301   301   0         301   93.04    84.03

 

EMBOSSのinfoseqも便利（link）

 

 

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　ソート。

長さでソートし、10000-bp以上の配列だけ出力。

seqkit sort --quiet -r -l conitg.fa -L 10000 > 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

#pigzに渡しgzip出力する。
seqkit sample -p 0.1 input.fq.gz |pigz -c - -p 8 > output.fq.gz
  • -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
  • -j     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

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

 

 

split　分割する。

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

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

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

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

 

 

 

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

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

 

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

Usage - SeqKit - Ultrafast FASTA/Q kit

 

追加コマンド

 

  引用

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.

 

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