CAZyデータベースからの配列の取得とFASTAファイルの取り扱い

:::note info
大学院生の備忘録です
:::

# はじめに

[Carbohydrate-Active enZYmes Database](http://www.cazy.org/ "CAZy") （CAZy Database）は、糖質関連酵素の分類を行っているデータベースである。CAZy の分類は、生化学的あるいは構造生物学的実験結果に紐づけられているため信頼性が高く、世界中の研究者に参照されている。

しかし、
* CAZy はタンパク質配列データを一括ダウンロードする方法を提供していない。
* タンパク質のファミリー分類はアミノ酸配列相同性に基づいて行われているが、CAZy はそのアルゴリズムを公開していない。

という問題点も存在する[^1]。

前者に関しては、第三者機関によって CAZy からデータを取得するプログラムが開発されている。


<b>SACCHARIS</b>[^2] は、おそらく最初に作られた CAZy データ取得プログラムで、CAZy のページをウェブスクレイピングし、配列取得、系統解析を行える。CAZy のサイトにタンパク質ID のリストをダウンロードできる機能がまだない時に作られ、非常に便利だったが、色々なサードパーティーソフトウェアが必要で、インストールが煩雑だった。
[GitHub](https://github.com/DallasThomas/SACCHARIS "SACCHARIS")

私はまだ試していないが、<b>cazy_webscraper</b>[^3] というものもあるらしい。
[GitHub](https://github.com/HobnobMancer/cazy_webscraper "cazy_webscraper")


2022年頃に、各ファミリーのページから NCBI の accession 番号のリストをダウンロードできる機能が追加され、GenBank ID の取得にウェブスクレイピングをする必要がなくなった。そのため、私は現在は Python で配列を取得する方法に落ち着いている。

後輩（たまに先輩）によくCAZy からの配列の取得方法を訊かれるのだが、コードの画面を見せると😓こういう顔をされたりする。

でもやっぱり Python できると便利なので、コードに触れた事のない後輩にも使ってもらえることを目指してこの記事を作成した。

また、単一のファミリーに数千、数万の配列が登録されている場合も多く、配列の取得と取捨選択がセットで必要となることがほとんどであるため、CD-HIT を用いて配列を絞る方法にも触れる。


# Python で NCBI ID のリストから配列を取得する

CAZy の各ファミリーのページで `Download GH○○` のような表示をクリックすると .txt ファイルがダウンロードできる。

このリストから配列を取得するためのコードをGoogle Colab で作成した。
コードが何をやっているか分かるようにノートブック中にできるだけコメントを残した。
[GitHub](https://github.com/lomon978410/fetch_cazy/tree/main)

### 使い方
#### 準備
1)  Google drive に解析用のディレクトリを作る。
2) [GitHub](https://github.com/lomon978410/fetch_cazy/tree/main) > Code > Download zip  からファイルをダウンロードし、解凍した fetch_fasta.ipynb ファイルを自分の Google drive に保存する。
3) 保存したファイルを右クリック、アプリで開く から Google Colaboratory を選択して開く。選択肢にない場合はアプリを接続する。

#### 実行
1) ダウンロードした ID の .txt ファイルを解析用のディレクトリにアップロードする。
2) コード中の `Entrez.email = "xxxx@yyyy.zz"` を自分のアドレスに書き換える。
3) `ID_list="/content/drive/MyDrive/xxx/yyy/zzz.txt"` このパスを ID の .txt ファイルのものに書き換える。
4) ランタイム > すべてのセルを実行 をクリックして実行。または、三角のマークをクリックしてセルを一つずつ実行。

:::note warn
注意
配列はNCBIのサーバーから取得している。一万を超えるような大量の配列を取得する場合、非常に時間がかかるため、途中でランタイムの接続が切れる可能性がある。その場合、ローカル環境で行う、事前にファイルから行をランダム抽出する、後述の 「CUPP」 を利用するといった対応が必要になる。
:::


# CD-HIT を使って冗長性を排除する

取得した配列を様々な解析に用いる際、配列数が多いと計算に時間が掛かる。また、データベースに登録されている配列は、大腸菌など何度もシークエンス解析されている配列の割合が多くなっており、冗長性がある。そのため、私は [CD-HIT](https://sites.google.com/view/cd-hit "CD-HIT") を使って相同性が高い配列を除いている。
＊CD-HIT Web サーバー版があったのだが、2022年10月以降使えなくなった……[^1]。

### CD-HIT のインストール

* OS が Ubuntu の場合や、 Windows で wsl 等の Ubuntu の仮想環境を使っている場合は、 [こちら](https://www.howtoinstall.me/ubuntu/18-04/cd-hit/) のコマンドで簡単にインストールできる。

* [Anaconda でも簡単にインストール](https://anaconda.org/bioconda/cd-hit) できる。

どちらも使っていない人は、これを機にどちらか始めてみるといいかもしれない。

### 使い方


ターミナルを開いて解析用のディレクトリに移動し、以下のコマンドを実行する。
```bash:
cd-hit -i cazy_all.fasta -o cazy_all_90.fasta -c 0.9
```
オプションは、
`-i` 入力ファイル名 
`-o` 出力ファイル名
`-c` 配列相同性の閾値（デフォルト値は 0.9 = 90%）

これで`-c` オプションで指定した以上の相同性を持つ配列を除ける。
`-c 0.5` など低い閾値を指定するとエラーがでる（どこまでいけるかはデータセットによって異なる）。

他のオプションについては、
```bash:
cd-hit -h
```
でヘルプを確認。

# CAZy から  Characterized 配列を取得する

解析にあたって、機能既知配列のホモログとその他を区別したいケースが多いと思う。しかし、 CAZy データベースは特定のページ（ Characterized など）を指定してのID ファイルのダウンロードには対応していない[^1]。

そこで、CAZy のサイトからウェブスクレイピングを行って GenBank ID を取得し、配列をダウンロードするコードも作成した。

:::note warn
ウェブスクレイピングしている都合上、CAZy のページ構成が変わるとエラーがでる可能性がある。
:::
### 使い方

事前準備は同じ。
[fetch_cazy.ipynb](https://github.com/lomon978410/fetch_cazy) を使う。

#### 実行

1) コード中の `Entrez.email = "xxxx@yyyy.zz"` を自分のアドレスに書き換える。
2) `outdir=` の後ろのパスをファイルを出力したいディレクトリに書き換える。
3) `cazy_url=` の後ろのurl を取得したいCAZy のページのものに書き換える。配列が 100 を越える場合は複数のタブがあるので注意。
4) ランタイム > すべてのセルを実行 をクリックして実行。または、三角のマークをクリックしてセルを一つずつ実行。
:::note warn
表の中に GenBank ID がない場合があるため、CAZy に登録されている配列数と実際にダウンロードできた配列数が一致しているか確認した方が良い。
:::
:::note warn
「Structure」 には対応していない。
:::



# おわりに

登録されている配列の数はファミリーによって異なり、数万に及ぶ配列が登録されているファミリーでは NCBI から配列をダウンロードするのに時間が掛かる。

その場合、私は Uniprot や [CUPP](https://cupp.info/ "CUPP")[^4] を利用することが多い。

CUPP は k-mer 法を用いた手法によって、ファミリーのサブファミリーへの分類を行っており、CUPPプログラムで分類済みの配列の tsv ファイルが簡単にダウンロードできる。
（ダウンロードは　[CUPP](https://cupp.info/ "CUPP") のページで BROUSE CAZY > 見たい Class / Family / SubFamily に✔ > BROUSE > Download Annotations:　の ALL TSV をクリック ）

ただし、Uniprot や CUPP は、各データベース独自の手法（Uniprot は隠れマルコフモデル、CUPP は CUPP プログラム）によって配列を取得しており、登録されている配列は完全に同一ではない。そのため、あるファミリーに分類されるかされないか微妙な境界領域の配列は、アノテーションがCAZy 公式とは異なっている可能性がある。
このような懸念はあるが、CUPP の表形式のデータは扱いが楽だし、サブファミリー分類も素晴らしいと思っているので私はよく利用している。

（そもそも CAZy に前述の問題があるのが原因なので、CAZy 公式も早くこういう形でデータをダウンロードできるようにしてほしい……。）




[^1]:2023年7月7日現在。今後のアップデートに期待。
[^2]:Jones DR, Thomas DK, Alger N, Ghavidel A, Inglis GD, Abbott DW. SACCHARIS: An automated pipeline to inform discover of new carbohydrate active enzyme activities within polyspecific families and de novo sequence datasets. Biotechnology for Biofuels, 11 (2018), p. 27, DOI: 10.1186/s13068-018-1027-x.
[^3]:Hobbs, E. E. M., Gloster, T. M., and Pritchard, L. (2022) 'cazy_webscraper: local compilation and interrogation of comprehensive CAZyme datasets', bioRxiv, https://doi.org/10.1101/2022.12.02.518825
[^4]:Barrett, K., Hunt, C. J., Lange, L., & Meyer, A. S. (2020). Conserved unique peptide patterns (CUPP) online platform: peptide-based functional annotation of carbohydrate active enzymes. Nucleic Acids Research, 48(W1), W110–W115. https://doi.org/10.1093/nar/gkaa375