2022/03/06 duplicated IDの配列の修正にseqkit renameを使うように修正
mm2-fastについて紹介しましたが、上手く導入できなかったたため一旦非公開にしました。失礼しました。代わりに簡単な記事を書きます。
メタゲノムのビニングが終わってbin配列を手に入れたら、ファイル名とfastaのヘッダをリネームして、後々の作業で扱いやすくしておくのが重要になります。bin配列があるディレクトリで以下のコマンドを実行すれば、リネームしたfastaを素早く入手できます。
コマンド
seqkitを使ってそれぞれのbin配列(.fna)のmulti-fastaを長さ順にソートし、1000-bp以上のみ抽出、続いてawkとsedを使って、ヘッダーをbin<INT>_にリネーム、bin<INT>.faとして出力する。
#!/bin/bash
count=1
for file in `\find *fna -maxdepth 1 -type f`; do
genome=${file}
seqkit rename $genome |seqkit sort --quiet -r -l -L 1000 | awk '/>/ {$0 = ">bin_" ++n} 1' - > select.fa
sed -e "s/bin/bin$count/g" select.fa > bin${count}_selected.fa
echo $genome bin${count}_selected.fa >> list
count=`expr $count + 1`
done
rm select.fa
元のファイル名とリネーム後のファイル名の関係はlistファイルに保存。出力をbin${count}_selected.faと書きましたが、”_selected”は無い方がスッキリしますね。
テスト
ubuntu18.04でテストした。
1、下の写真は元のファイル。MAG.faの名前が統一されていません。加えて、multi-fastaのそれぞれの配列のヘッダ名は長く、特殊文字も存在している可能性もあります。
2、次は上のコマンドを実行後の写真。元のファイルは残したまま、bin<INT>.faにリネームされたファイルが出力されています。ファイルの中のヘッダも修正しているので、bin1.faファイルなら>bin1_1, >bin1_2, >bin1_3・・・というヘッダになっているはずです。
listにはファイルの対応関係を保存しています。
より実践的には、この結果に分類系統をアサインして種名(無いならそれより上の階級)、にリネームしておくと進めやすくなります。対応関係のリストを作っておく必要があります。
引用
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.