特定の細菌種のゲノムを網羅的に取得する

# 背景
[NCBI](https://www.ncbi.nlm.nih.gov/)のgenomeデータベースを利用すると任意の細菌のゲノム配列やアノテーション情報（GenBankファイル）をブラウザからダウンロードできます。ただ、同じ細菌種の株ごとの違いを検証したいときなど、複数のデータをブラウザから手動でダウンロードするのは大変です。そこで、NCBIから細菌ゲノムの登録情報が記載されたファイルを取得し、ターミナルのコマンドでGenBankファイルを自動・網羅的に取得してみます。

# 手法
#### 任意の細菌種を検索し、登録されている株の登録情報を取得する。
1. NCBIのgenomeデータベースでダウンロードしたい細菌名を検索します。ここでは大腸菌（Escherichia coli）とします。
<img width="1214" alt="National_Center_for_Biotechnology_Information.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/22a8b379-5240-45c4-8d19-1eff8f801603.png">

2. 基準株の情報が画面上部に出てきます。「Browse the list」のlistをクリックします。
<img width="1187" alt="Escherichia_coli__ID_167__-_Genome_-_NCBI.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/1c69b62b-a93e-3783-6e5d-04821b7b1f8d.png">

3. 登録されている大腸菌ゲノムの一覧が表示されます。2020年6月23日現在で20,034の大腸菌ゲノムが登録されているようです。検索ボックスにキーワードを入力すれば、さらに絞り込むこともできます。
4. 「Download」をクリックすると、20,034の大腸菌ゲノムの登録情報が記載されたcsvファイルがダウンロードされます。
<img width="1435" alt="Genome_List_-_Genome_-_NCBI.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/07d83787-a9fc-f2c6-93b9-ce7e87d61681.png">

5. エクセルなどで開きます。株名やゲノムサイズ、CDS数などが一目で分かります。今回使用するのは右から2番目にある「GenBank FTP」列です。
![prokaryotes_csv.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/0a10a5f2-82d9-1d18-2224-851b02c30dd2.png)
![prokaryotes_csv-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/d80ff2ac-c518-6894-27fe-7045f1cfb3a3.png)

6. 試しに1番上の「Escherichia coli str. K-12 substr. MG1655」の「ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/005/845/GCA_000005845.2_ASM584v2 」をブラウザに入力してみます。

7. 「Escherichia coli str. K-12 substr. MG1655」に関するファイルが格納されたディレクトリにアクセスできます。「GCA_000005845.2_ASM584v2_genomic.fna.gz」をクリックすればゲノム配列（fasta）、「GCA_000005845.2_ASM584v2_genomic.gbff.gz」をクリックすればGenBankファイルがダウンロードできます。
<img width="1381" alt="_genomes_all_GCA_000_005_845_GCA_000005845_2_ASM584v2_のインデックス.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/984b4e6c-c26a-7dab-c231-63740f91bdde.png">

8. 以上から、csvファイルに記載のあるFTPサーバー情報から入手したいファイルのURLを作成することで、ある細菌種のデータファイルを取得できることが分かります。

#### 取得するファイルのURLを作成する。
1. ここではGenBankファイル（gbff）を取得することにします。csvファイルに記載されているFTPサーバー情報は「ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/005/845/GCA_000005845.2_ASM584v2
 」、実際のGenBnakファイルのアドレスは（ブラウザ上で右クリック→「リンクのアドレスをコピー」で取得できます）「ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/005/845/GCA_000005845.2_ASM584v2/GCA_000005845.2_ASM584v2_genomic.gbff.gz
」です。
2. 見比べるとcsvファイルのアドレスに「/GCA_000005845.2_ASM584v2_genomic.gbff.gz」を追記すれば良いことが分かります。このうち「GCA_000005845.2_ASM584v2」は本ファイルが格納されているディレクトリ、すなわち最下層のディレクトリ名です。そのため、作成したいURLは「csvファイルのFTPサーバーアドレス + "/" + 最下層のディレクトリ名 + "_genomic.gbff.gz"」となります。
3. エクセルの関数を使用しても作成できますが、ここではRを用いて作成してみます。

```r
#環境
sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

#データの読み込み
dt <- read.csv("ダウンロードしたcsvファイル", stringsAsFactors = F)

#最も下層のディレクトリ名の抽出
dirVec <- sapply(dt$GenBank.FTP, FUN = function(x) {
        split <- strsplit(x, "/")[[1]]
        return(split[length(split)])
})

#取得するファイルのURLを作成
gbffVec <- paste(dt$GenBank.FTP, "/", dirVec, 
                 "_genomic.gbff.gz", sep = "") 

#書き出し
write.table(gbffVec, "gbff_table.txt", quote = F, row.names = F, 
            col.names = F)
```

gbff_table.txtの中身は下記の通りです（20,034は多いので一部間引いています）。
![gbff_table_txt.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/e1a372ca-e340-0a08-c83a-287600084af2.png)


#### ファイルのダウンロード
1. R上でダウンロードすることも可能ですが、ここではターミナルのwgetコマンドを用いてダウンロードしてみます。
2. wgetコマンドは-iオプションの後ろにファイル名を記載することで、ファイルに記載のあるURLにアクセスしてくれます。ここで先ほど作成したgbff_table.txtを指定します。また、-Pオプションでダウンロード先のディレクトリを指定することができます。

    ```sh
wget -i gbff_table.txt -P gbff/
    ```
<img width="1461" alt="1__bash.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/363491/201bea4e-c40e-e0e8-fd4c-601c72751adc.png">

3. 後はgunzipコマンドでファイルを解凍すれば目的のGenBankファイルを利用することができます。ファイル数が多い場合はgunzipを並列に行えるpigzコマンドがお勧めです。

[【pigz】Linux環境でマルチコアをフル活用して圧縮・解凍する方法](http://www.lesstep.jp/step_on_board/linux/572/#pigz)


# 補足
もっと良い方法あれば教えて下さい。