fasta, fastq ファイルを扱うときに便利なseqkit

fasta, fastqファイルに含まれている配列数、各配列の長さ、N50 の基本情報がほしい
全部の配列名を取得したい、名前を変えたい、特定の領域の配列だけ取り出したい

そんなときに激烈便利なのがseqkitです。
condaでインストール可能　https://anaconda.org/bioconda/seqkit
seqkitの使い方（公式）​https://bioinf.shenwei.me/seqkit/usage/

使用例

fastaファイルの基本情報が知りたい（配列数、平均長、N50など）

$ seqkit stat hoge.fastafile
format  type  num_seqs  sum_len  min_len  avg_len  max_lentest
FASTA   DNA    3      160       49     53.3       58​

$ seqkit stat -a hoge.fastafile
format  type  num_seqs  sum_len  min_len  avg_len  max_len  Q1  Q2  Q3  sum_gap  N50  Q20(%)  Q30(%)
FASTA   DNA          3      160       49     53.3       58  49  53  58        0   53       0       0
Bash

option
-a N50とかクオリティ（fastaqを扱うと値がでる）とか出してくれる

fastaファイルから各配列の情報のリストが欲しい（配列ごとの長さ、GC含量など）

$ seqkit fx2tab -nl hoge.fasta > hoge.name.length.list
Bash

option
  -B　ある塩基が含まれる量を出力。-B ATと指定すれば配列あたりのATの数を出力。
  -g　GC含量を出力
  -G　GC-Skew （つまり、各配列の(C-G)/(C+G)）を出力
  -H　print header line
  -l　配列長を出力
  -n　配列名だけ出力。このオプションを入れないと、配列とそのクオリティも出力される
  -i　print ID instead of full head

ある配列だけほしい

$ seqkit grep -np contig01 hoge.fasta > hoge.tig01.fasta
Bash

option
-n 正規表現を利用する。（完全一致のみを抜き出したいときは不要）

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig02
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

→

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG

ある配列のある領域だけほしい (部分的に抽出)

$ seqkit subseq --chr contig01 -r 15:30 hoge.fasta > hoge.chr1_15_30.fasta
Bash

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig02
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

→

​>contig01_15-30 contig01​
GTCATCGATCGTCAG

一覧を使って複数の配列を抜き出したい

$ seqkit grep -f list hoge.fasta > hoge.tig0103.fasta
Bash

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig02
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

> list

contig01
contig02

→

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

全配列の名前を名前一覧リストをつかってrenameしたい

$ seqkit replace -p '^(\S+)(.+?)$' -r '{kv}' -k rename.tsv  hoge.fasta > hoge.renamed.fasta

名前*だけ書き換えたいとき（名前の末尾にスペースやタブが入ってるときも有効）　
$ seqkit replace -p '^(\S+)' -r '{kv}' -k rename.tsv  hoge.fasta > hoge.renamed.fasta

説明部分*を書き換えたいとき
$ seqkit replace -p ' (.+)$ -r ' {kv}' -k rename.tsv  hoge.fasta > hoge.renamed.fasta
Bash

*　>seq01 hogehoge
のとき、
 　seq01が名前、hogehogeが説明部分に対応

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig02
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

> rename.tsv (tab区切り)

contig01 seqA
contig02 seqB
​contig03 seqC

→

​>seqA
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>seqB
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>seqC
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

配列名をHoge0001の形に一括で変換したい

$ seqkit replace -p '^(\S+)(.+?)$' -r "Hoge{nr}" --nr-width 4 test.fa
Bash

option
–nr-width 桁数

​>contig01
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>contig02
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>contig03
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

→

​>Hoge0001
AACTGCTAATGCTTGTCATCGATCGTCAGTCGATCGATCGCTAGCATCGATCGATCAG
>Hoge0002
ACTGCTTTGCTTTTCCCCCCGCTCGTTGCTAAATCGTGCTGATGTGCTC
>Hoge0003
TGCTGGCTTCCGCTGTACTCGAGCACGATCGTCGTACGTGGTCGTATCGGATC

コメント