AlphaFoldを真似た論文ProSPrの距離行列予測の一部解析

論文 ProSPr: Democratized Implementation of Alphafold Protein Distance Prediction Network についての解析。論文URLは
<https://biorxiv.org/content/10.1101/830273v2.full>

## 概要

実装が公開されていないAlphaFold（正確にはAlphaFoldの距離分布を考慮したコンタクトマップの生成部分のみ）を再実装してみました論文。GitHubで公開され、Dockerも置いてある。
https://github.com/dellacortelab/prospr

Proteinsに掲載された[AlphaFoldの論文](https://doi.org/10.1002/prot.25834)はこのProSPrの論文とほぼ同時期に発表されたので、実際にはAlphaFold論文を読んで真似したものではなく、それまでのDeepMind社のプレスリリースなどの図から想像して再実装したものであることに注意。

この2ヶ月後、Natureの方にAlphaFoldの論文が改めて掲載された。これを読むと、まぁまぁ実装が似ているような似ていないような……。 <https://www.nature.com/articles/s41586-019-1923-7>

このProSPr論文は距離行列生成部分に注力しており、他の様々なソフトウェアを使う構造予測パートについては触れられていないので注意。

## MATLABとplmDCAのインストール

ProSPrは処理の途中で[plmDCAというプログラム](https://github.com/magnusekeberg/plmDCA)を利用することになる（ちなみにAlphaFoldではplmDCAと目的が似ている[CCMPred](https://qiita.com/Ag_smith/items/cd426f9582b4d9c79835)を使い、かつ同様にフロベニウスノルム出力のためのコード改変を加えているようだ）。このプログラムはMATLABで書かれており、実行するためにはMATLAB環境が必要である。ただ、コンパイルが終了すれば無料版のMATLAB runtimeで実行することができるようになる（？） 以下、MATLABをインストールしてplmDCAが使えるようになるまでの記録を残す。

### MATLABのインストール

Linux OSの場合について説明する。まず、なんとかしてMATLABライセンスをゲットする（UTokyoは全学生が使えるMATLABライセンスがあったので助かりました）。
その後、引き続いて、インストーラーを選ぶことができる。ここではprosprが動作確認しているバージョン**R2018a**を選択する。インストーラー名は`matlab_R2018a_glnxa64.zip`。ただし、重要なアップデートがあるため、そのファイルもダウンロードする(`R2018a_Update_6_glnxa64.sh`)。

MATLABのファイル展開は

```bash
mkdir -p ~/Downloads/temp ; cd ~/Downloads/temp
mv /path/to/{matlab_R2018a_glnxa64.zip,R2018a_Update_6_glnxa64.sh} .
unzip matlab_R2018a_glnxa64.zip
```

これを展開したところに現れる`install`というファイルがインストーラーだが、**ライセンス認証にはどうしてもGUI環境が必要っぽいので**GUIを用意した上でターミナルから`./install`を打つとインストーラーが起動する。

インストールの手順として

1.  ログインの種類は**Log in to your Mathworks
    Account**を選んで大学包括ライセンスで登録したメールアドレスとパスワードを入力する。
2.  **Choose installation Folder**のところでは、インストール先を選ぶことができる。私の場合は`/home/apps/MATLAB/R2018a`にした。
3.  **Select products to install**のところでは、様々な追加パッケージのインストールオプションが存在するが、**MATLAB,
    Simulink, Bioinformatics Toolbox, MATLAB Compiler, Parallel
    computing toolbox, Statics and Machine Learning Toolbox, Symbolic
    Math Toolbox**を選択する。

こののち、インストールを開始する。インストールが終わると続けてアクティベーションが問われる。このライセンスは1台のマシンにつき1アカウントしか利用することができないので、インストールしたいマシンの中のユーザーアカウント名を入力してネットワーク越しのライセンス認証を行う。これで完了。

環境設定はこんな感じ

``` bash
export PATH="/path/to/MATLAB/R2018a/bin:$PATH"
export LD_LIBRARY_PATH="/path/to/MATLAB/R2018a/bin/glnxa64:/path/to/MATLAB/R2018a/sys/os/glnxa64:/path/to/MATLAB/R2018a/bin/glnxa64:$PATH"
```

### 改変したplmDCAのインストール

オリジナルは[plmDCA](https://github.com/magnusekeberg/plmDCA)だが、ProSPrにはこれをかなり改変したものを用いているため、いくつかの操作を行う。

まずはplmDCAオリジナルをGitでクローンしてくる。

``` bash
git clone https://github.com/magnusekeberg/plmDCA.git
# またはzipで落としてきて展開
```

次に私がオリジナルprosprから修正したprosprをGitでクローンしてくる。

``` bash
git clone https://github.com/YoshitakaMo/prospr_mod.git
```

この`prospr_mod`の中に`potts/plmDCA_asymmetric`というディレクトリがあるので、**このディレクトリの中にmagnusekeberg氏のplmDCAの`plmDCA_asymmetric_v2`の中身をすべてコピーする。ただし、`plmDCA_asymmetric.m`ファイルだけは`YoshitakaMo/prospr_mod`のものに入れ替える。**

次に、`plmDCA_asymmetric/3rd_party_code/minFunc/autoDiff/autoTensor.m`コードの一部を以下のように書き換える。

```patch
- [~ ~ diff(:,:,j)] = funObj(x + mu*e_j,varargin{:});
+ [~, ~, diff(:,:,j)] = funObj(x + mu*e_j,varargin{:});
```

これで、`prospr_mod/potts/plmDCA_asymmetric`ディレクトリはこのようにファイル構成になっているはず

```shell:prospr_mod/potts/plmDCA_asymmetric
-rw-rw-r--. 1 moriwaki moriwaki      8665 Jan  8 16:00 plmDCA_asymmetric.m
drwxrwxr-x. 2 moriwaki moriwaki       104 Dec 14 02:09 functions
drwxrwxr-x. 3 moriwaki moriwaki        50 Oct 18  2016 3rd_party_code
-rw-rw-r--. 1 moriwaki moriwaki      3591 Oct 18  2016 readme.txt
-rw-rw-r--. 1 moriwaki moriwaki       178 Oct 18  2016 mexAll.m
drwxrwxr-x. 2 moriwaki moriwaki       143 Oct 18  2016 test_examples
```

この状態でMATLABのコンパイラを使って、`plmDCA_asymmetric.m`をコンパイルする。

```bash
mcc -m plmDCA_asymmetric.m -a functions -a 3rd_party_code
```
ちょっと時間がかかるが、エラーなく終了すれば成功。すると、`plmDCA_asymmetric`と`run_plmDCA_asymmetric.sh`という実行可能ファイルが生成されるはず。前者はスタンドアロンで動作するMATLABのコンパイル済みバイナリであり、他の環境に移すことも可能らしく、MATLAB Runtimeさえ入っていればMATLABは必要なく動作するらしい……。後者はターミナルのシェル上でこのコンパイルされた`plmDCA_asymmetric`を動作させるためのシェルスクリプトとなっており、以下のコマンドでplmDCAを動作させることができる。

```bash
./run_plmDCA_asymmetric.sh /path/to/MATLAB/R2018a ./test_examples/example1_alignment.txt test_score.mat 0.2 8
```
第1引数`/path/to/MATLAB/R2018a`はMATLABのバイナリへのPATHを指定する。MATLABを持っていない場合はMATLAB Runtimeをダウンロードしてきてインストールし、`/path/to/MATLAB_Runtime/v94`のように指定する。第2引数`./test_examples/example1_alignment.txt`はインプットとなる`a2m`形式の配列アライメントファイルを指定する。これは後述の**hhblitsを利用したMSAの取得**部分を参照されたし。第3引数はアウトプットファイルのファイル名を指定する。第4引数`0.2`はしきい値、第5引数`8`は計算に使用するCPUコア数を指定する。

これは以下のplmDCA公式のREADME （https://github.com/magnusekeberg/plmDCA/tree/master/plmDCA_asymmetric_v2） でのMATLAB用実行コマンドと等価である。

```MATLAB
plmDCA_asymmetric('./test_examples/example1_alignment.txt','test_score.mat',0.2,8)
```
よって、上記の`run_plmDCA_asymmetric.sh`コマンドで正しく`test_score.mat`ファイルが生成されれば、ProSPr用改変済plmDCAのテストは正常に終了である。


## ProSPrを使ってみる

ProSPrオリジナルは https://github.com/dellacortelab/prospr.git においてあるが、このコードだけでは色々動作しないようになっているので、以下に記す実装解析を参考にしてなんとか動かすことを考える。

あらかじめオリジナルprosprから修正したprosprをGitでクローンしておき、上の改変済みplmDCAがターミナルから動作することを確認しておく。

``` bash
git clone https://github.com/YoshitakaMo/prospr_mod.git
```

## 実装の解析

### ProSPrのbuild処理の全体の流れ
以下では例としてファイル名（本来は[CATH domain list S35](http://download.cathdb.info/cath/releases/latest-release/cath-classification-data/cath-domain-list-S35.txt)に登録されているドメイン区分名）を`2o6pA00`とする。

`python3.6 prospr.py build 2o6pA00`でやっている中身は以下の通り。

1.  アミノ酸配列情報が含まれているFASTA形式の`2o6pA00.fasta`を用意する。**アミノ酸配列部分に改行が含まれていると正しくpklファイルが作られないので注意**。<br>
<pre>
>2O6P:A|PDBID|CHAIN|SEQUENCE
GSDSGTLNYEVYKYNTNDTSIANDYFNKPAKYIKKNGKLYVQITVNHSHWITGMSIEGHKENIISKNTAKDERTSEFEVSKLNGKIDGKIDVYIDEKVNGKPFKYDHHYNITYKFNGPTDVAG
</pre>
1.  `hhblits.py`の内部で、`hhblits`を利用したMSAの取得を行う。出てきた`2o6pA00.a3m`ファイルから`reformat.pl`, `hhmake`コマンドを使ってそれぞれ`2o6pA00.a2m`, `2o6pA00.hhm`を得る。
1.  （並行して）`psiblast`を行い、pssmファイルの生成を行う。`2o6pA00.pssm`を得る。
1.  `2o6pA00.a2m`の結果を用いて、改変したplmDCAを使って`2o6pA00.mat`を生成する。
1.  `2o6pA00.fasta`, `2o6pA00.pssm`, `2o6pA00.mat`, `2o6pA00.hhm`などを用いて、Pythonのpickleを使って特徴量のpickleファイル`2o6pA00.pkl`を生成する。

これでbuild処理は終了。要はこの手順さえ踏めば目的の`2o6pA00.pkl`が得られるので、下の**hhblitsを利用したMSAの取得**, **psiblastを利用したPSSMの取得**, **plmDCAを利用したmatファイルの取得**, **pklファイルの生成**を順に達成すれば良い。

### hhblitsを利用したMSAの取得

実際にプログラムを動かすコマンド。例として、`/mnt/HDD/database/2o6pA00/`においた`2o6pA00.fasta`を利用する。
動作には`uniclust30_2018_08_hhsuite.tar.gz`ファイルが必要。`hhblits`の使い方の詳細や準備方法は[類縁配列検索ソフトウェアHH-Suiteのインストール方法](https://qiita.com/Ag_smith/items/401291c22fca19d61af9)を読んでください。

``` bash
# 1. hhblitsでMSAを取得
hhblits -i /mnt/HDD/database/2o6pA00/2o6pA00.fasta \
        -oa3m /mnt/HDD/database/2o6pA00/2o6pA00.a3m \
        -d /mnt/HDD/database/hhblits/uniclust30_2018_08/uniclust30_2018_08 \
        -o /mnt/HDD/database/2o6pA00/2o6pA00.hhr -e 0.01 -n 3 -B 100000 -Z 100000 \
        -maxmem 4.0 -v 0 -cpu 8
```

ここで`-maxmem`はメモリ使用量を表している。例えば`less /proc/meminfo`で使っているマシンのメモリを確認することができるが、64GBメモリを積んでいるマシンであれば`-maxmem 60.0`としても大丈夫だろう（少し下げておいたほうが良い）。同様に、`-cpu 6`も使っているのCPUのコア数の最大値を設定するとよい。こちらは`less /proc/cpuinfo`で確認することができる。

しばらくするとアウトプットとして`2o6pA00.a3m`ファイルが出てくる。これは

```bash
{HH-suite}/scripts/reformat.pl /mnt/HDD/database/2o6pA00/2o6pA00.a3m /mnt/HDD/database/2o6pA00/2o6pA00.a2m
```

とすることでa2m形式に変換できる。これを後のplmDCA解析のインプットとして与えるMSA情報になる。

また、

```bash
hhmake -i 2o6pA00.a3m -o 2o6pA00.hhm
```

で出てくる`2o6pA00.hhm`ファイルも後で必要なので、このコマンドを実行しておく。

### psiblastを利用したPSSMの取得

処理としては`psiblast`コマンドで単純にPSSMファイルを取得している。以下はそのコマンド。

``` bash
DOMAINNAME="2o6pA00"
BASEDIR="/mnt/HDD/database/${DOMAINNAME}"
DBDIR="/mnt/HDD/database/blast"
psiblast \
  -query ${BASEDIR}/${DOMAINNAME}.fasta \
  -db ${DBDIR}/nr/nr \
  -out_ascii_pssm ${BASEDIR}/${DOMAINNAME}.pssm \
  -out ${BASEDIR}/${DOMAINNAME}.out \
  -evalue 0.001 -num_iterations 3 \
  -save_pssm_after_last_round -num_threads 8
```

### plmDCAを利用したmatファイルの取得
先述の改変済みplmDCAを使って以下のコマンドでplmDCAを起動する。

```
./run_plmDCA_asymmetric.sh /path/to/MATLAB/R2018a /mnt/HDD/database/2o6pA00/2o6pA00.a2m /mnt/HDD/database/2o6pA00/2o6pA00.mat 0.2 8
```

これによって`2o6pA00.mat`ファイルが生成される。このファイルにはフロベニウスノルムの情報が入っている。

### pklファイルの生成
最後に、`2o6pA00.fasta`, `2o6pA00.pssm`, `2o6pA00.mat`の結果ファイルを使って`2o6pA00.pkl`ファイルを生成する。この処理は`./prospr/io.py`の中の以下の`build`関数にある。

```python
class Sequence(object):
    def __init__(s, name, **kwargs):
        """Accepts name as input and expects a corresponding .pssm, hhm, mat, and fastfa txt file (dashes represent gaps), and outputs a pkl file that can be used by ProSPr to predict contacts"""
        s.name = name
        s.base = pconf.basedir + s.name + "/" + s.name
        s.data = dict()
        s.outfile = s.base + ".pkl"

    def build(s, args):
        """Builds the pkl file for the protein sequence and saves it.  `args` is an argparse object and it will look, specifically, for a boolean args.potts"""
        s.seq_name()
        s.pssm()
        s.hh()
        s.potts()
        filename = s.base + ".pkl"
        save(s.data, filename)
```


### ProSPrのrun処理の流れ
buildされた`2o6pA00.pkl`ファイルが存在する状態で、`python3.6 prospr.py run 2o6pA00`とすると、学習済みのモデル（これは本家ProSPrのリンクをたどっていけばダウンロードできます）を用いて、2o6pA00についての各残基の距離確率分布の生成を行う。この可視化スクリプトが`prediction_visualization.ipynb`である。これは最も確率が高かった距離の区分（2〜22 Åで0.32Åごとのbinのどこが最も確率が高いか）を示している。

## 正解の距離行列との比較

正解の距離行列を取得するPython3スクリプト

``` python
import numpy as np
from scipy.spatial import distance_matrix
from Bio.PDB.mmtf import *

class DistanceMatrix():
    # 変数 struc に入る情報を持っておいてinit
    # 実体はまだ
    def __init__(self, resname):
        self.startres = -1
        self.seqlen = 0
        self.tempcoords = None
        self.resname = resname

    # タンパク残基っぽいかどうか。False Positiveは例が思い浮かばない。
    def is_sane(self, res):
        return "C" in res and "N" in res and "O" in res and "CA" in res

    # 一定残基 targetlen ごとに
    def save_matrix(self):
        '''
        与えられたnumpy形式の座標データに対して距離行列を生成する関数
        scipy.spatialのdistance_matrixを使用
        変数matに保存
        保存ファイル名は2Z6FAという入力に対して"npys/2Z6FA.npy"みたいな形になる
        '''
        templen = len(self.tempcoords)
        mat = distance_matrix(self.tempcoords, self.tempcoords)
        savename = "npys/" + self.resname[:5]
        np.save(savename, mat)

    # run関数で計算開始
    def run(self):
        mp = MMTFParser()
        try:
            # 入力としてとってきたresnameの最初の4文字をPDB IDとしてインターネットから取得する
            # 例： 2Z6FAならば、self.resname[:4]は2Z6Fとなる
            struc = mp.get_structure_from_url(self.resname[:4])
        except:
            return
        for res in struc[0][ self.resname[4:] ]:
            '''
            struc[0][ self.resname[4:] ]はstruc[0]["A"]のような感じ
            resは以下のような構造体になっている
            <Residue LEU het=  resseq=544 icode= >
            <Residue THR het=  resseq=545 icode= >
            res.get_idはこんな感じ
            (' ', 544, ' ')
            (' ', 545, ' ')
            '''
            resid = res.get_id()[1]
            if not self.is_sane(res): continue
            '''
            startresの初期値は-1なので最初の1回だけstartresを読み込むファイルに合わせて
            更新する。
            res["CA"].get_coord()とすることで
            [-11.45   27.282   5.276]
            [-9.645 30.579  5.574]
            [-7.091 32.439  7.742]
            のようにCα原子の座標を取得できる。
            '''
            if self.startres < 0:
                self.startres = resid
                self.tempcoords = np.array( [res["CA"].get_coord()] )
            else:
                self.tempcoords = np.concatenate(
                    (self.tempcoords, np.array( [res["CA"].get_coord()] )
                    )
                )

            '''
            self.tempcoordsは
            [[-11.45   27.282   5.276]
            [-9.645 30.579  5.574]
            [-7.091 32.439  7.742]
            ...
            [a b c]]
            みたいな配列
            '''
            self.seqlen += 1
        self.save_matrix()
        # np.saveで"ans_coords_2Z6FA.npy"みたいなファイル名で座標情報を保存する
        np.save( "ans_coords_" + self.resname[:5], self.tempcoords )
        print(self.resname)

# ここにIDとChainを5文字で入れる。例：2Z6FA
#DM = DistanceMatrix("2Z6FA")
DM = DistanceMatrix("2O6PA")
DM.run()
```

これで`./npys/2O6PA.npy`というファイルが作成される。次に、これをmatplotlibで22 Åカットオフで表示させる。

``` python
import numpy as np
import matplotlib.pyplot as plt
import pickle as pkl  
import pandas as pd

# 入力
file = "./npys/2O6PA.npy"

truemap = np.load(file)
fig = plt.figure()
im = plt.imshow(truemap, cmap='viridis_r', interpolation='none', vmin='0.0', vmax='22.0')
fig.colorbar(im)
```

正解の距離行列アウトプットはこんな感じ。

![5e6a5ef750f43c340c8b6472.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/fc647d31-a47f-b622-ff94-2e8090c5c13b.png)

ProSPrによる予測結果はこんな感じ。

![5e6a5efb50f43c340c8b6473.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/73b0cb27-9de5-4d3e-d371-e2afde47d10a.png)

大体良さげな感じ？