2022/03/14追記
これまで数回に分けてseqkitのコマンドを紹介して来ましたが（リンク）、バージョンアップが続いていて、ありがたいことに新しいコマンドも追加されています（谢谢您）。久しぶりに新機能を確認してみます。
この記事を書いたすぐ後にv2.2が公開されたので、そちらも追記しています。
2022/03/14
seqkit v2.2.0 brings many changes:
— Wei Shen 沈伟 (@shenwei356) March 14, 2022
1. support of reading/writing xz (.xz) and zstd(.zst) files, 2. the new md5-like tool: `seqkit sum`,
3. rewritten `concat`; `split2` with custom out prefix & suffix.
4. other new features proposed by users. https://t.co/n4tgFL8LFG
HP
https://bioinf.shenwei.me/seqkit/
インストール
Githubリリース （armのM1 (apple silicon)のmac向けリリースもあります）
#version2.1を導入（conda）2.2はリリースよりダウンロード（リンク）
mamba install -c bioconda seqkit=2.1 -y
新機能
まず、v0.9.4の時点であったコマンドを復習しておきます。
アルファベット順
- duplicate 配列を複製
- common IDや配列が共通する配列抽出
- concat 同じIDの配列を指定長ずつ合体
- convert fastqのクオリティエンコードの変換 (Solexa形式 ⇔サンガーイルミナ)
- duplicate 配列を指定回、複製
- faidx FASTAのindexを作成、特定の配列を抽出（samttools faidxより高機能）
- fq2fa fastqをfastaに変換
- fx2tab fasta/fastqをtab形式に変換し、さらにGCなどの情報を表示
- genautocomplete シェルオートコンプリートスクリプト
- grep 特定の配列をミスマッチ許容で検索
- head 指定数だけ先頭からリードを抽出
- locate 配列をミスマッチ許容でリプレース
- range 指定間のリードを抽出 (headでもtailでも対応できない時)
- rename duplicateしたIDをユニークな名前に修復
- replace 名前/文字列を正規表現で置き換える
- rmdup duplication除去
- restart 環状ゲノムのスタート位置変更
- sample リードをランダムサンプリング（10%だけランダムに取り出す）
- sana 壊れた一行だけの fastq ファイルをサニタイズ(置き換える)
- seq 配列を変形（逆相補鎖、DNA<=>RNA、指定サイズ以上/以下）
- shuffle ランダムな順番にする
- sliding 配列を指定の単位ずつずらしながら部分重複ありで取り出す
- sort 長さでソート
- split 分割（fastqを均等にn個に分けて出力）
- split2 入力配列を分割
- stats fasta/fastqの要約統計
- subseq gtf/bedの領域に隣接するなど特定の配列だけ抽出
- tab2fx FASTA/FASTQをタブ区切り（TSV）フォーマットに変換
- translate DNA/RNAをアミノ酸配列に翻訳
- watch 配列特徴量のモニタリングとオンラインヒストグラム
その後もバージョンアップが続き、バージョン番号v0.16まであがった後、v2.0になりました (Tags)。2,022年3月現在、バグ修正もされたv2.1.0が最新となっています。v0.9.4以後に追加されたコマンドや当時紹介していなかったコマンドをアルファベット順で見ていきます。
- amplicon - アンプリコン プライマーを介してアンプリコンを取得
フォワードとリバース（逆相補）のプライマー配列を指定する。
echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt
#bedをつけると配列の代わりにポジションが出力される
echo -ne ">seq\nacgcccactgaaatga\n" \
| seqkit amplicon -F ccc -R ttt --bed
出力
出力領域を調整することもできます。詳細はリンク先（赤文字がリンクになっています）で確認して下さい。
- bam - BAMファイルの特徴量のモニタリング
bamファイルを指定する。
#bamの要約統計
seqkit bam -s *.bam
#indexed bamの要約統計（-sよりずっと高速、出力内容は異なる）
seqkit bam -i *.bam
#保存（STDERR）
seqkit bam -i *.bam 2> result
#BAMのindexを使用してリファレンスの染色体ごとににマッピングされたトータルのリード数をカウント（高速）
seqkit bam -C sorted_indexed.bam
#リファレンスにマッピングされたリードのデプス分布をプリント
cat sample.bam | seqkit bam -c counts.tsv -
- -s print BAM satistics of the input files
- -i fast statistics based on the BAM index
cat sample.bam | seqkit bam -c counts.tsv -
seqkit bam -s *.bam
seqkit bam -i *.bam
- fish - ローカルアライメントを利用して、巨大な配列から短い配列のマッチを探す
クエリの配列とターゲットの配列を指定する。
#GGCGGCTGTGACCを探してポジションをプリント(STDOUT)
seqkit fish -F GGCGGCTGTGACC genome.fna
#-gをつけるとアラインメントもプリントされる（STDERR）
seqkit fish -q 4.7 -F GGCGGCTGTGACC -g mouse-p53-cds.fna
#結果をbam形式で保存。-fでクエリの配列を指定することもできる。
seqkit fish -a -q 4.67 -f query.fas -b alignments.bam -g mouse-p53-cds.fna
- -q minimum mapping quality (default 5)
- -F query sequences
- -g print sequence alignments
- -b save aligmnets to this BAM file (memory intensive)
- -f query fasta
出力
出力座標はBEDライクな0ポジションベース、左閉じ、右開き
- mutate - 配列を編集する（点突然変異、挿入、欠失）
クエリの配列とターゲットの配列を指定する。
#-pで点変異；最初（1）をxに置換
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -p 1:x --quiet
#-iで挿入変異；最後(-1)にxxを挿入
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -i -1:xx --quiet
#-dで欠失変異；最初（1）の塩基を欠失
echo -ne ">1\nACTGNactgn\n>2\nactgnACTGN\n" \
| seqkit mutate -d 1:1 --quiet
- -p point mutation: changing base at given position. e.g., -p 2:C for setting 2nd base as C, -p -1:A for change last base as A
- -i insertion mutation: inserting bases behind of given position, e.g., -i 0:ACGT for inserting ACGT at the beginning, -1:* for add * to the end
- -d deletion mutation: deleting subsequence in a range. e.g., -d 1:2 for deleting leading two bases, -d -3:-1 for removing last 3 bases
- pair - 2 つの fastq ファイルからペアエンドのリードをマッチング（同期）させる。
ペアエンドfastqを指定する。出力ディレクトリを指定しなければ現在のパスに書き出される。
seqkit pair -1 reads_1.fq.gz -2 reads_2.fq.gz -O out_dir
#ペアでないリードも保存する
seqkit pair -1 reads_1.fq.gz -2 reads_2.fq.gz -O result -u
- -O output directory
- -u save unpaired reads if there are
reads_1.paired.fq.gzとreads_2.paired.fq.gzが書き出される。
#fastq_dir下の全てのfastqを連結
seqkit scat -j 4 -f fastq_dir > all_records.fq
#ディレクトリを監視し、その下に5秒間書き込みがなくなるまでリアルタイムでFastqレコードをストリーム処理
seqkit scat -j 4 -T "5s" fastq_dir > all_records.fq
#
指定したPIDのプロセスがアライブするまで、ディレクトリを監視してリアルタイムでfastqレコードをストリーム処理
seqkit scat -j 4 -p $PID fastq_dir > all_records.fq
壊れたfastqを指定する。
seqkit sana broken.fq.gz -o rescued.fq.gz
- -i input and output format: fastq or fasta (default "fastq")
壊れたレコードを除くリードがrescued.fq.gzとして出力される。
環状ゲノムのスタート位置だけが違う全く配列かどうか確認するために利用できる。
seqkit sum -c -k 51 virus-*.fasta
- -k k-mer size for processing circular genomes (default 1000)
- -c the file contains a single cicular genome sequence
- -a show all information, including the sequences length and the number of sequences
- -s only consider the positive strand of a circular genome, e.g., ssRNA virus genomes
環状ゲノムのスタート位置だけが違う配列同士（-cを付ける）や、逆補数鎖だけで同じ配列同士（- 一本鎖ゲノムの場合は、-sを使用）は同じ値がプリントされる。 配列のヘッダと品質はスキップされ、配列のみが分析対象になる。
- watch - 配列特徴量の分析と監視
リード長分布のヒストグラム。
#ONTのリードを分析
seqkit watch --fields ReadLen nanopore.fq.gz -O len.png
出力を指定しなければコンソール（STDOUT）に出力される。
seqkit watch --log --fields ReadLen nanopore.fq.gz
出力を指定しなければコンソール（STDOUT）に出力される。デスクトップがないCUI環境で便利。
seqkit watch -p 500 -O qhist.png -f MeanQual nanopore.fq.gz
他にもhead-genomeなどのコマンドがあります。
引用
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
