ビニングして得たfastaのファイル名とヘッダ名を一括リネームする

2022/03/06 duplicated IDの配列の修正にseqkit renameを使うように修正

2023/08/10追記

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でテストした。

１、下の写真は元のファイル。MAG.faの名前が統一されていません。加えて、multi-fastaのそれぞれの配列のヘッダ名は長く、特殊文字も存在している可能性もあります。

f:id:kazumaxneo:20210728115809p:plain

２、次は上のコマンドを実行後の写真。元のファイルは残したまま、bin<INT>.faにリネームされたファイルが出力されています。ファイルの中のヘッダも修正しているので、bin1.faファイルなら>bin1_1, >bin1_2, >bin1_3・・・というヘッダになっているはずです。

f:id:kazumaxneo:20210728115936p:plain

listにはファイルの対応関係を保存しています。

f:id:kazumaxneo:20210728120250p:plain

より実践的には、この結果に分類系統をアサインして種名（無いならそれより上の階級）、にリネームしておくと進めやすくなります。対応関係のリストを作っておく必要があります。

2023/08/10追記

上記ではbin${count}_selected.faのようなファイルにリネームしました。各MAGについて、bin5ならヘッダ名を

> bin5_1

> bin5_3

> bin5_2

のようにリネームしましたが、大規模なMAGセットでは、"binN_<INT++>"だけだとほかのMAGのヘッダ名と重複する可能性があります。ファイル名中のSRA/ENA/DRA名もヘッダに含めるように修正します。カレントのパスに既に上のリネーム処理をしている"bin"をファイル名中に含む.faファイルセットがあるとして、以下のようにリネームを実行。

#!/bin/bash
for file in `\find *bin*.fa -maxdepth 1 -type f`; do
 genome=$file
 name=${genome%.fa}_
 output=${genome%.fa}_renamed.fa
awk -v genome_var="$name" '
  /^>/ {
    count++;
    printf ">" genome_var count "\n";
  }
  !/^>/ {
    print;
  }
' $file > $output

例；このようなファイルがあるとして（"bin"を途中に含む.faの拡張子を持つfasta）、