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

## 背景
以前にも[似たような記事](https://qiita.com/otsukimi200624/items/e9e1164f48d6f1728ea7)を書いたのですが、その際はブラウザを使用していたので、今回はawkとsedを使用してコマンドベースで取得します。

参考：[Download Genomic Sequences From NCBI](https://ftp.ncbi.nlm.nih.gov/pub/factsheets/HowTo_Downloading_Genomic_Data.pdf)

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

```sh
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で抽出してみます。![assembly_summary_bacteria_1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/40f752a6-e4ca-d1da-44a4-5ce162074a0d.png)

```sh
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を修正します。

```sh
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における正規表現の扱いの違い](https://rcmdnk.com/blog/2015/07/19/computer-gnu-bsd/)

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

```sh
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を指定し、ダウンロードします。

```sh
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コマンドでダウンロードする。

```sh

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
```