モチーフ検索による酵素の論理的探索をしたい

# はじめに
天然には多様な糖質が存在しているため、その合成、分解に関わる<b>糖質関連酵素</b>の種類も多様化している。同じファミリーに分類される酵素であっても基質特異性が多様化していることが多く、アミノ酸配列から基質を推定することが困難なことも多い。

基質特異性の違いは、基質認識に関わるアミノ酸残基の違いによって生じる。
グローバルな配列が異なれば、ローカルな（基質結合部位の）配列も異なり、基質特異性が異なる。
そのため、アミノ酸配列相同性が高い酵素は類似した活性を示す（たまに例外もある）。

では、基質結合部位を形成するアミノ酸残基は完全にランダムなのだろうか？
ということが気になった。

### 基質結合部位に共通点がある酵素の例

同じファミリーに異なる基質特異性を示す酵素が存在する一方で（適応放散？）、異なるファミリーに同基質特異性を持つ酵素が存在していることも多い（収斂進化？）。

例えば、細菌の細胞壁を分解するリゾチーム（ムラミダーゼ）という酵素に関して、異なるファミリーに分類される酵素同士で、活性部位を形成する重要なアミノ酸残基が共通している例が報告されていた[^1]。

このように、同じ基質を認識する酵素で、重要なアミノ酸残基がある程度のパターンに収斂するという現象は、他の酵素にも起こっているのではないかと思う。

そうならば、基質結合部位を形成しうるアミノ酸残基の候補は無限ではなく、合理的に推定可能なのではないか。特定の基質の結合に要求される構造を類推し、その構造をモチーフとして検索できないかと考えた。

これを検証するため、まず、ターゲットとして糖質加水分解酵素ファミリー31（GH31）を選んだ。


### GH31ファミリー

GH31ファミリーは、同一ファミリー内での基質特異性と反応機構の多様化が進んだファミリーの一つである。このファミリーでは、一対の触媒残基（アスパラギン、D）が既報のすべての酵素で保存されている。また、基質は、触媒残基と適切な電子の授受が出来る、ちょうどいい位置に結合する必要があるため、基質の非還元末端側（図中「-1」）が酵素に結合する角度もほぼ同じになっている。

![qiita.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3508057/ca05bae8-2c12-9d76-5f84-e31f4c6b579a.png)


ここに何か法則性を見出せないかとGH31の構造を眺めているうち、構造既知の酵素では、２位の水酸基（2OH）の近くに必ずアルギニン残基（R）があることに気づいた。

一方で、GH31には、2OHを持たない（ガラクトースの 2OH が N-アセチル化されている）糖である N-アセチルガラクトサミン（GalNAc）に作用する α-GalNAcase という酵素が存在している。
α-GalNAcase と GalNAc の複合体の立体構造を見てみると、この R に相当する残基（以下 R <sub>2OH</sub>）が存在していなかった。

![qiita2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3508057/fa621fd9-c1ff-6138-feb1-84157219d6bc.png)


現在、GH31 に登録されている N-アセチル化糖に作用する酵素は α-GalNAcase のみである。もし、R<sub>2OH</sub> を持たない別の配列を見つけられれば、他の〇〇NAc系の基質（GlcNAc、ManNAc 等）に作用する新しい酵素の可能性がある。

さらに構造と配列アラインメントをよく見てみると、 R<sub>2OH</sub> はすべて酸塩基触媒残基の 16 か 17 残基上流にあることが分かった。

なので、これを使ってアミノ酸配列を選別してみる事にした。


# 方法

### データセットの用意

アミノ酸配列相同性が高い配列同士をグループ化したデータセットが必要になる。

[CUPP](https://cupp.info/) は、糖質関連酵素を k-mer 法という手法で分類したデータセットを提供している。系統解析を行い、両者の結果がほぼ一致している事を確認したので、酵素の分類には CUPP のサブファミリーを用いる事にした。

データの前処理として、CUPP から取得したテーブルを一度 fasta 形式のファイルにし、CD-HIT で配列相同性 100% の配列を除いた。

必要な行のみ取り出したところ。
![cupp_table.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3508057/a9c51a61-767e-d1ad-463b-2a7487e4ae65.png)



テーブルの形は、 (8653, 5)。 

`df2 = df1[df1['DomainLength'] > 100]` で不完全長と思われる配列を除くと、

(8557, 5)

CUPP サブファミリーの個数は全部で 70 個だった。


### 酸塩基触媒残基モチーフ

![prosite.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3508057/4f406770-479b-f77a-7002-211452e1ab9a.png)


[prosite](https://prosite.expasy.org/PDOC00120) によると、GH31ファミリーの酸塩基触媒残基は

[GFY]-[LIVMF]-W-x-<font color=red>D</font>-M-[NSA]-E

で表されるらしい。

テーブルと正規表現パターンを引数として、マッチする配列を持つ行をデータフレームとして返す関数を書いた。

```python
def match_yes (regex, df):
    dfnew = df.iloc[0:0]
    regex = re.compile(regex)
    for i in range(len(df)):
        m = re.search(regex, df.iat[i, 3])
        if m:
            dfnew = pd.concat([dfnew, df[i:i+1]])
    return dfnew
```

これを使って本当にすべての GH31 にこのモチーフが含まれているか調べると、ヒットしたのは 2151 個だった。
2151/8557 ……???

詳しく見てみると、このモチーフは動物の α-グルコシダーゼにしか適用できないことが分かった。

全部の GH31 酵素に適応できる触媒残基モチーフはないものかと各サブファミリーのアラインメントを作成して眺めまわした結果、

側鎖に芳香環をもつアミノ酸　-　嵩の小さいアミノ酸1 or 2個　-　<font color=red>D</font>

が、すべての（当時）構造既知の酵素に共通していることが分かった。従って、

[WFY]-x{1,2}-<font color=red>D</font>

がより包括的な酸塩基触媒残基のモチーフと言えそうだ。
あまりに短いので触媒残基と関係のない部分にもマッチしそうだが、これと16、17残基上流に R<sub>2OH</sub> があるという条件を組み合わせれば、何とかなりそうな気がする。


# 結果

### モチーフ検索

逆に、モチーフを持たない配列を調べるため、上の関数の `if m:` を `if not m:` とすることで、マッチしない配列を返す関数（`match_no`）を作成した。

これを使って調べると、[WFY]-x{1,2}-D のモチーフを持たない配列はたった 25 個だった。


続いて本命の、R-D 間の距離 16 or 17 × [WFY]-x{1,2}-D に当てはまらないものをふるいにかける。 

```python
print(len(df3))
dfA = match_no('[R].{13}[WFY][GAVLISTDENQCM]D', df3)#R-D => 16, WxD
print('-> ', len(dfA))
dfB = match_no('[R].{12}[WFY][GAVLISTDENQCM]{2}D', dfA)#R-D => 16, WxxD
print('-> ', len(dfB))
dfC = match_no('[R].{14}[WFY][GAVLISTDENQCM]D', dfB)#R-D => 17, WxD
print('-> ', len(dfC))
dfD = match_no('[R].{13}[WFY][GAVLISTDENQCM]{2}D', dfC)#R-D => 17, WxxD
print('-> ', len(dfD))
```
【結果】

8557
->  6608
->  1816
->  1704
->  708

十分の一以下まで絞れた。
[WFY]-x{1,2}-D というモチーフは触媒残基以外にも存在するが、R<sub>2OH</sub>と セットにするとふるいとして機能してくれた。


# 選別された配列の評価

残った708 配列を各 CUPP サブファミリーごとに集計し、R<sub>2OH</sub> が保存されていない割合順にソートした。

![results_table.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3508057/abbe3f25-c297-6927-5958-d34b7a40f1be.png)



系統樹と照らし合わせて、内訳を見てみる。

既報の GalNAcase が属している CUPPサブファミリー GH31:25.1 、それと近縁なGH31:29.1、GH31:50.1 などは 90 ％ 以上がふるいに残り、「N-アセチル基が結合するならば R<sub>2OH</sub> を持たない」という仮説と合致した。


残りのサブファミリーを見てみると……
あれ、機能解析済みの酵素がある。
既報の α-ガラクトシダーゼとその近縁配列がヒットしていた。

AlphaFold2 による予想構造を見てみると、 R<sub>2OH</sub> は存在している。
なぜだ、と思いよく見ると、R ではなく、触媒残基のモチーフのほうが保存されていなかった……。

結局、他の候補も R<sub>2OH</sub> と 触媒残基の間にループの挿入があるだけだったりと、GalNAcase 以外に R<sub>2OH</sub> を持たないタンパク質を見つけることはできなかった。


# 総括

残念なことに、GH31から新たな〇〇NAcase を発見することはできなかった。

しかし、活性スクリーニングプロジェクトによって発見されたGH31 GalNAcase は、実はアミノ酸配列情報からも発見できた可能性が示された。

ついでに、 GH31 の触媒残基以外の部分が高度に多様化していることを改めて確認できた。

この方法は他のファミリーの酵素に適用できる可能性が高い。必ずしも触媒残基の近くに重要なアミノ酸が存在しているわけではないので、今後はモチーフ検索に依らない手法を使う必要がある。

その手法の候補としては、ペアワイズアラインメントを考えている。（というかコード書いて放置している）

[^1]:Cernooka, E., Rumnieks, J., Zrelovs, N. et al. Diversity of the lysozyme fold: structure of the catalytic domain from an unusual endolysin encoded by phage Enc34. Sci Rep 12, 5005 (2022). https://doi.org/10.1038/s41598-022-08765-1