RefSeqから任意の細菌種のfasta配列を網羅的に取得する

背景

以前にも似たような記事を書いたのですが、その際はブラウザを使用していたので、今回はawkとsedを使用してコマンドベースで取得します。

参考：Download Genomic Sequences From NCBI

方法

アセンブリデータ一覧の取得

NCBIのFTPサイトから現在登録されている細菌のアセンブリデータ一覧を取得します。

curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt' > assembly_summary_bacteria.txt

標的細菌ゲノムのFTPパスを抽出

ダウンロードしたassembly_summary_bacteria.txtは行が登録されている細菌を示し、名前やアセンブリレベル、ファイルへパスなど22列で構成されたタブ区切りテキストファイルです。2020年12月24日現在で202,813のアセンブリデータが登録されていました。
ここで8列目のorganism_name、12列目のassembly_level、20列目のftp_pathに着目します。organism_nameはそのまま細菌種名で「属名 種名」が記載されます。assembly_levelはアセンブリがCompleteしているのか、あるいはcontigやscaffoldの状態であるかを示しています。ftp_pathはFTPサイト内のそのアセンブリデータのURLです。
ここではグラム陰性で海洋に広く存在するシェワネラ属（Shewanella）のコンプリートゲノムをawkで抽出してみます。

awk '{FS="\t"} !/^#/ $8 ~ /Shewanella/ && $12 ~ "Complete Genome" {print $20} ' assembly_summary_bacteria.txt > Shewanella_complete_ftp.txt

head -n3 Shewanella_complete_ftp.txt
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/016/585/GCF_000016585.1_ASM1658v1
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/169/215/GCF_000169215.2_ASM16921v2
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/044/255/GCF_003044255.1_ASM304425v1

wc -l Shewanella_complete_ftp.txt
59 Shewanella_complete_ftp.txt

awkコマンドの意味は下記の通りです。
{FS="\t"} FSはField Separator＝区切り文字で、awkはデフォルトでスペースが区切り文字なのでタブ区切りに変更
!/^#/ 先頭が「#」の行=コメント行の1行目と列名の2行目を除外
$8 ~ /Shewanella/ && $12 ~ "Complete Genome"
8列目（organism_name）に「Shewanella」を含み、かつ12列目（assembly_level）が「Complete Genome」の行を抽出
{print $20} 20列目（ftp_path）をプリント

シェワネラ属のコンプリートゲノムはRefSeqに59個登録されているようです。FTPサイトへのURLも問題なく抽出できています。

fasta配列へのURL作成

Shewanella_complete_ftp.txtのURLはアセンブリデータフォルダへのURLです。フォルダ内にはゲノム配列やタンパク配列など複数のデータが存在します。今回はゲノムのfasta配列のみをダウンロードしたいため、URLを修正します。

sed -E 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' Shewanella_complete_ftp.txt > genome_fasta_Shewanella.txt

※-EはBSD sedで拡張正規表現を有効にするオプション。GNU sedの場合は-Eを-rにする。
参考：GNU/BSDでのsedにおける正規表現の扱いの違い

なお以上のコマンドはパイプを用いて下記のようにも書けます。

curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt' | \
awk '{FS="\t"} !/^#/ $8 ~ /Shewanella/ && $12 ~ "Complete Genome" {print $20} ' | \
sed -E 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genome_fasta_Shewanella.txt

標的細菌ゲノムのfasta配列の取得

wgetの--inputオプションでgenome_fastaShewanella.txtを指定し、ダウンロードします。

wget --input genome_fasta_Shewanella.txt

なお、fastaファイルのURL作成まではエクセルのフィルタ機能や関数を用いてもできますので、不安な場合はエクセルでassembly_summary_bacteria.txtを開き、同様の結果になるかを確認するのも良いかと思います（かなりメモリ大きくないと固まりそうですが。。。）。

参照：GenBankファイル（gbff）をダウンロードしたいとき

fastaファイルと同様にassembly_summary_bacteria.txtの20列目をawkとsedで加工していく。ただし、GCFをGCAに、genomicをgbffに変換することに注意。作成したファイルを使ってwgetコマンドでダウンロードする。

awk '{FS="\t"} !/^#/ $8 ~ /Shewanella/ && $12 ~ "Complete Genome" {print $20} ' assembly_summary_bacteria.txt| \
sed -E 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' | sed s/GCF/GCA/g | sed s/fna/gbff/g > gbff_path_Shewanella.txt