アミノ酸配列のみからタンパク質構造のコンタクトマップを予想するCCMPredプログラムの利用法

## Introduction
**タンパク質のコンタクトマップ** とは、タンパク質が折り畳まれていったときに、どのアミノ酸残基とどのアミノ酸とが近くなったかという距離の情報を2次元マップで表したものです。例として、157残基長のPDB: 5DIRのChain Aのコンタクトマップは次の白黒画像のようになります。

<img src="https://qiita-image-store.s3.amazonaws.com/0/224327/b8390178-f42a-d59e-6865-25602945bf3c.png" width="400px"><img src="https://qiita-image-store.s3.amazonaws.com/0/224327/d3f12a96-ea90-cb3e-d466-8e2bee220c1c.png" width="400px">

この記事では、なんと、**アミノ酸配列だけから、このコンタクトマップの予想を高精度で返してくれる** プログラムである **GREMLIN/CCMPred** の使い方を紹介します。
GREMLINは、David Baker研究室から2013年に発表された、**共進化情報を利用した強力なコンタクトマップ予測法** であり（DOI: 10.1073/pnas.1314045110）、CCMPredはGREMLINを並列化して大幅に高速化したプログラムとなっています（DOI: 10.1103/PhysRevE.87.012707 , DOI: 10.1002/prot.22934）。
GREMLIN/CCMPredは**Direct coupling analysis (DCA)**の一種として分類されることもあり、類似のウェブサービスは[EVcouplings](http://evfold.org/evfold-web/evfold.do)やライス大学の[DCA Webservice](http://dca.rice.edu/portal/dca/)などがあります。

はじめに、ウェブサーバー版GREMLINを使った結果を紹介しながら、それを自前のパソコン・サーバー上で動かす方法を紹介していきます。
## 何に使えるの？
**正確なコンタクトマップが得られるのならば、配列からのタンパク質構造予測に強力に役立ちます**。

## 計算環境準備
### マシン環境
-   macOS 10.13.6 (High Sierra), 2.8GHz Intel Core i7, メモリ16GB またはLinux OS
-   macOSの場合[ターミナルからHomebrewをインストール](http://brew.sh/index_ja.html)してある
-   **100GB以上のディスク容量を持つストレージ**。SSDであればさらに動作が早い。
-   [類縁配列検索ソフトウェアHH-Suiteのインストール](https://qiita.com/Ag_smith/items/401291c22fca19d61af9)をしてある

### 配列検索のソフトウェアHMMERのインストール
#### macOSの場合
Homebrewで簡単にインストールできる。 `brew install hmmer` で終了
#### Linux OSの場合
MPIを用いた並列計算をしたい場合は./configure時に `--enable-mpi` を付けるべき。 `--prefix` でインストール先を指定する。今回の例では`/usr/local/package/hmmer/3.2.1` にしておきます。

    wget http://eddylab.org/software/hmmer/hmmer-3.2.1.tar.gz
    tar zxvf hmmer-3.2.1 && cd hmmer-3.2.1
    ./configure --prefix=/usr/local/package/hmmer/3.2.1 --enable-sse --enable-threads
    make -j8
    make -j8 install
    # miniappsの追加インストール
    cd easel
    make install

### map_alignのインストール
`map_align`は2つのコンタクトマップが得られたときに、それらの中でオーバーラップしている領域を最大化し、ギャップを最小化するための簡単なプログラムです。あとでちょっと必要になるので、こいつもインストールしておきます。GitHubは https://github.com/sokrypton/map_align

> map_align takes two contact maps and returns an alignment that attempts to maximize the number of overlapping contacts while minimizing the number of gaps.

インストール方法は

    mkdir -p /path/to/apps/map_align #インストール先ディレクトリを適当に作る
    cd /path/to/apps
    git clone https://github.com/sokrypton/map_align.git
    cd map_align
    g++ -O3 -std=c++0x -o map_align main.cpp # map_alignプログラムのコンパイル
    # binディレクトリを作り、その中にプログラムを入れておく
    mkdir bin
    cp mk_map/*.pl ./bin
    cp map_align ./bin
    # 実行権限を付与(permission denied対策)
    chmod +x ./bin/*

上記の処理後、`export PATH="/path/to/apps/map_align/bin:$PATH"`として、ここにPATHを通しておきます。

### CCMPredのインストール
GREMLINの並列＆GPU計算実装版のCCMpred ver.0.3.2をインストールします（現行は0.3.3 alpha）。NVIDIA GPU/CUDAによる高速化に対応しており、インストール時にGPU対応オプションを指定しておくと5倍速くらいになりますが、CPUだけでも数分の処理時間なので、まあなくてもどうにかなります。

macOS, LinuxOSの場合も、いずれもCMakeをインストールしておく必要があります。**CMakeのバージョンは3.14以上の最新版を使っている方が無難です**。Ubuntu 18.04に`apt`でインストールできるCMakeはver.3.10.2ですが、`nvcc fatal   : Unknown option 'fopenmp'`と出てしまいインストールに失敗します。macOSならば `brew install cmake` でインストールできます。

インストール方法は

    # gitでダウンロード
    git clone --recursive https://github.com/soedinglab/CCMpred.git
    cd CCMpred
    mkdir build ; cd build
    # GPU/CUDAを使いたい場合は指定……しなくても勝手に検知してくれることがある。
    # -DCUDA_TOOLKIT_ROOT_DIR= でCUDAドライバの場所を指定。
    # CUDAを使わない場合は-DWITH_CUDA=OFFをcmake時のオプションに追加する
    # ここでインストール先は/usr/local/package/CCMpred/0.3.2とする（任意）
    cmake .. -DCMAKE_C_COMPILER=/bin/gcc -DCMAKE_CXX_COMPILER=/bin/g++ -DCMAKE_INSTALL_PREFIX=/usr/local/package/CCMpred/0.3.2 -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda
    make -j4
    make install

これでインストール完了。簡単な使い方は次のような感じ

    $ ccmpred -h
    Usage: ccmpred [options] input.aln output.mat

    Options:
            -d DEVICE       Calculate on CUDA device number DEVICE (set to -1 to use CPU) [default: 0]
            -t THREADS      Calculate using THREADS threads on the CPU (automatically disables CUDA if available) [default: 1]
            -n NUMITER      Compute a maximum of NUMITER operations [default: 50]
            -e EPSILON      Set convergence criterion for minimum decrease in the last K iterations to EPSILON [default: 0.01]
            -k LASTK        Set K parameter for convergence criterion to LASTK [default: 5]

            -i INIFILE      Read initial weights from INIFILE
            -r RAWFILE      Store raw prediction matrix in RAWFILE

            -w IDTHRES      Set sequence reweighting identity threshold to IDTHRES [default: 0.8]
            -l LFACTOR      Set pairwise regularization coefficients to LFACTOR * (L-1) [default: 0.2]
            -A              Disable average product correction (APC)
            -R              Re-normalize output matrix to [0,1]
            -h              Display help

ヘルプを読んだら実行。`-t`にはCPU数、`-n`にMax Iteration回数、`-d`でGPUのナンバー（通常は0）を指定する。その後インプット・アウトプットファイル名の指定。

### HH-suiteのインストール
上の [類縁配列検索ソフトウェアHH-Suiteのインストール](https://qiita.com/Ag_smith/items/401291c22fca19d61af9)を参考に、プログラムをインストールしておきます。

### データベースの構築
Uniclust30データベースと、Uniref100データベースの2つのセットアップが必要です。
#### uniclust30_2018_08データベースの用意
[類縁配列検索ソフトウェアHH-Suiteのインストール](https://qiita.com/Ag_smith/items/401291c22fca19d61af9)の中で紹介してありますので、そちらを参照してください。
#### Uniref100データベースの用意
ftp://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref100 からダウンロードし、BLAST+のプログラム`makeblastdb`でデータベース構築をします。BLAST+のインストール方法はここでは割愛します。

    mkdir -p /path/to/database/uniref100 && cd /path/to/database/uniref100
    # uniref100は34.6 GB以上のファイルサイズ
    wget ftp://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref100/uniref100.fasta.gz
    gunzip uniref100.fasta.gz
    makeblastdb -title uniref100 -dbtype prot -in uniref100.fasta -max_file_sz 2GB

## GREMLIN webサーバー版を使ってコンタクトマップ予測をしてみる

Baker研究室がGREMLINのwebサーバー版を公開してくれています。ここに投げたクエリ配列とその結果をもとに、自前マシンでの再現実験のときの指標にします。

<http://gremlin.bakerlab.org/submit.php> に向けて配列を投げてみます。テストケースとして、膜タンパク質PDB: 5DIRのChain Aについてやってみます。これは、当時**Pfamで新規ドメインとして認識され、類似構造が存在していないタンパク質構造とされていたものです**。なので、この構造が報告されるまではHHPredなどのホモロジーモデリング手法ではうまくいかなかったそうです。

    >querytest(5DIR_A)
    LVPRGSHMPDVDRFGRLPWLWITVLVFVLDQVSKAFFQAELSMYQQIVVIPDLFSWTLAYNTGAAFSFLADSSGWQRWLF
    ALIAIVVSASLVVWLKRLKKGETWLAIALALVLGGALGNLYDRMVLGHVVDFILVHWQNRWYFPAFNLADSAITVGAVML
    ALDMF

（見る人が見たら、最初に`LVPRGS`っていうthrombin認識配列がついていることがわかります）これをGREMLINのフォームに入力しますが、オプション設定は以下の画像のように設定しました。  
![600px-Gremlin1.png](https://qiita-image-store.s3.amazonaws.com/0/224327/27e552fc-0be0-d9dd-20fe-511f005b821d.png)
 
Generate MSA using:`HHblits`にし、E-value:は`1E-20`、Iteration:は`8`にした。Filter
MSAについてはCoverage: `75`、Remove gaps:`50`とした。

10日ほど待つと（このときはジョブ待ちが非常に長かった）、結果はこんな感じで返ってきました。  
![600px-Gremlin3.png](https://qiita-image-store.s3.amazonaws.com/0/224327/68691fd0-3b74-f4b1-6d65-f899f328bd77.png)

右上のDownload alignmentにはフィルタリングされる前のアライン結果が表示されます。Seq/Lenというのは、`hhblits`で得られたMSAの配列数がSeq, クエリの配列長をLenとして割った値です。()内の数字はおそらく実効の値で、Lenでは最初の15アミノ酸くらいでMSA中で50％以上のギャップが存在していたため、その部分の配列は実効値にカウントされなかったと考えられる。Seqの実効値は上のCoverage: 75%フィルターで減ったのかな？ちょっと自信がないですが、多分そういうことだと思います（適当）。

![600px-Gremlin4.png](https://qiita-image-store.s3.amazonaws.com/0/224327/a49195a5-08e9-f8b1-7e74-4be42cdfcd6a.png)
左にGREMLINが出したContact Mapが示されています。最初に示した白黒画像でのコンタクトマップとは違って、上下反転しているので気をつけて比較してください。

![12.png](https://qiita-image-store.s3.amazonaws.com/0/224327/0cf6dcfd-98af-8386-b571-f10bfde254d5.png)
並べてみると、だいたいマッチしていることがわかりますね。

その下に共進化情報に基づいてコンタクトすると考えられる残基ペアのリストが書かれています。

下のNew Features！のところで**Import to GREMLIN BETA**を押すと、以下の画面が表示されました。

![600px-Gremlin2.png](https://qiita-image-store.s3.amazonaws.com/0/224327/4cf8a8e4-d4b8-cd70-d55c-cffa9c940cb9.png)
情報に大きく違いはないように見えますが……？

さて、まあこうした内容を、*in house*で出力させるにはどうしたらいいかを考えてみます。

## 自前マシンでGREMLIN/CCMPredによるコンタクトマップ生成を行う
上のWebサーバーになげたときのクエリ配列を使って再現実験をしてみます。
まずはHH-suiteとHMMER, map_alignのコマンドが使えることを確認しておきます。

    export PATH="/path/to/apps/map_align/bin:$PATH"
    which map_align
    which hhblits
    which hmmbuild
    echo $HHLIB

としてみて、PATHが通っていることと、環境変数`HHLIB`が存在することを確認します。上のセットアップ通りにやっていれば大丈夫なはずです。

まず`hhblits`をクエリ配列が入っているファイル`query.fa`について、以下のオプションとともに行います。

    hhblits -cpu 12 -i query.fa \
            -d /Users/Agsmith/database/hhsuite/uniclust30_2018_08/uniclust30_2018_08 \
            -oa3m query.a3m -n 8 -e 1E-20 -neffmax 20 -nodiff -maxfilt 999999 \
            -realign_max 999999

ここで`-cpu`はCPU数を指定します。使っているマシンの最大CPU数を指定すると良いでしょう。結果として得られる`query.a3m`を上のGREMLINウェブサーバーの結果と比較してみたのですが、この時点ですでに10683本の類似配列が釣れてきました（ウェブサーバーは2724本）。`-neffmax 20 -nodiff -maxfilt 999999 -realign_max 999999`このへんのオプションが効いているんでしょうかね？

続いてこれを`hhfilter`でフィルタリングします。

    hhfilter -i query.a3m -o query_filt.a3m -id 90 -cov 75

得られたこのファイルについて、`reformat.pl`で`a2m`ファイル形式に変換します。これは続く`hmmbuild`でHMMの生成を行うためです。

    ${HHLIB}/scripts/reformat.pl query_filt.a3m query_filt.a2m
    hmmbuild --cpu 12 hmmout query_filt.a2m

これにより、作成されたHMMが`hmmout`ファイルに出力されます。

続く`hmmsearch`のステップでは、この`hmmout`と対象データベース（ここではUniref100）をもとに、`hmmsearch [options]  <hmmfile> <seqdb>`の形式で指定します。

    hmmsearch --cpu 28 -A hmmsearch.sto --tblout tblout.txt --incT 27 \
    hmmout /Users/Agsmith/database/blast/db/FASTA/uniref100/uniref100 > hmmsearch.out

`--incT 27`というのがPfam境界定義にも使われるという、bit-scoreのしきい値27の指定です。また、得られるMSAを`-A hmmsearch.sto`で出力させます。この出力形式はstockholmフォーマットになっていますため、続くフィルタリングでは`reformat.pl`でa2m形式に変換してから処理を行います。

得られた`hmmsearch.sto`に対し、再びフィルタリングをかけます。

    ${HHLIB}/scripts/reformat.pl hmmsearch.sto hmmsearch.a2m
    hhfilter -i hmmsearch.a2m -o hmmsearch_filt.a3m -id 90 -cov 75

さて、これからCCMpredにかけてcontact mapを得たいわけなのですが、この`hmmsearch_filt.a3m`は現時点では

    >UniRef100_A0A1W9WY37/34-184
    SNANALIWLWLSLLVVVLDQVTKVWVSATLELHQSIAVLPFFDLRLVHNEGAAWSFLADASGWQRWFLSGLAILISGVIVVWLSRLERQQRWLACALALILGGALGNVMDRAMYGYVVDFIDIYYQQMHWPVFNIADSAISVGAVMLLIDA
    >UniRef100_A0A148NB86/5-151
    ----HVKWLGLSVAVIVLDQLTKWWVVSTLRLHESIPVWPFFNLTYVRNTGAAFSFLAGAGGWQRWFFLGLAVVASMAILVWMWRLKDDEVLTAAGLSLVLGGALGNFIDRAAYGYVIDFLDFYIQNWHWPAFNVADSAITIGVILLLWDG
    >UniRef100_UPI0009EEDF12/24-173
    -VSRYRLLAIVTAVVLVLDQATKLYVDSHFRLHESVTVLEnFFHFTYVRNKGAAFGIFADS-AIRIPFFITVSLVAAVGILWYLRRLKDNQKLVCFSLALIFSGALGNLIDRVRFGEVIDFIDVHWYQYHWPAFNVADSAISVGVTLLLVDL
    >UniRef100_A0A0A6P993/4-156
    kSSEKALVWFWLSLLVVILDQVTKVWISANLALYQSIVVLPFFNLTLLHNEGAAWSILATAGGWQRWFLSGLAIVISGVLMIWLSRLKRQQRWLACALALILGGALGNVIDRVIYGYVVDFIDVYYNhEWHWPAFNIADSAISVGAVMLLIDA
    >UniRef100_UPI00041D792D/6-151
    -----RLLLIISAVVLVLDQATKLYIDSRFALHQSLTVVEnFFHITYVRNKGAAFGIFADS-AVRIPFFITVSIVAAIGIVWYLHRLNEQQRLLHVALSLIFAGAVGNLIDRIRLGEVIDFIDVHWYQYHWPAFNIADSAITVGVGILLLDL

こんな感じとなっており、FASTA形式でかつ配列長が揃っていません。CCMpredの計算処理をさせるためには、**PSICOV式フォーマット**になっていることと、**配列長が揃っていること** が条件だそうです（参考：https://github.com/soedinglab/CCMpred/wiki/FAQ ）。

CCMpredのディレクトリの中にも`scripts/`ディレクトリの中にその問題を解決してくれるスクリプトがあるのですが、それよりも便利なものが**map\_align**の中の`mk_map`ディレクトリに同梱されていました。ここの中にある`a3m2alncut.pl, mk_map.pl, pdb2map.pl`はこのcontact map作成に非常に便利なものです。

まずは <https://github.com/sokrypton/map_align/blob/master/mk_map/cmd> を参考にして、`hmmsearch_filt.a3m`からPSICOVインプットファイルを生成します。

    input=hmmsearch_filt ; a3m2alncut.pl ${input}.a3m ${input}.aln ${input}.cut

この`a3m2alncut.pl`内部で50％以上のギャップ位置の削除と75％未満のcoverage削除を行ってくれています。こうして、CCMpredのインプットとなる`hmmsearch_filt.aln`（下のMSA）と`hmmsearch_filt.cut`が生成されます。

    ALIWLWLSLLVVVLDQVTKVWVSATLELHQSIAVLPFFDLRLVHNEGAAWSFLADASGWQRWFLSGLAILISGVIVVWLSRLERQQRWLACALALILGGALGNVMDRAMYGYVVDFIDIYYQQMHWPVFNIADSAISVGAVMLLIDA
    HVKWLGLSVAVIVLDQLTKWWVVSTLRLHESIPVWPFFNLTYVRNTGAAFSFLAGAGGWQRWFFLGLAVVASMAILVWMWRLKDDEVLTAAGLSLVLGGALGNFIDRAAYGYVIDFLDFYIQNWHWPAFNVADSAITIGVILLLWDG
    YRLLAIVTAVVLVLDQATKLYVDSHFRLHESVTVLEFFHFTYVRNKGAAFGIFADS-AIRIPFFITVSLVAAVGILWYLRRLKDNQKLVCFSLALIFSGALGNLIDRVRFGEVIDFIDVHWYQYHWPAFNVADSAISVGVTLLLVDL
    ALVWFWLSLLVVILDQVTKVWISANLALYQSIVVLPFFNLTLLHNEGAAWSILATAGGWQRWFLSGLAIVISGVLMIWLSRLKRQQRWLACALALILGGALGNVIDRVIYGYVVDFIDVYYNEWHWPAFNIADSAISVGAVMLLIDA
    -RLLLIISAVVLVLDQATKLYIDSRFALHQSLTVVEFFHITYVRNKGAAFGIFADS-AVRIPFFITVSIVAAIGIVWYLHRLNEQQRLLHVALSLIFAGAVGNLIDRIRLGEVIDFIDVHWYQYHWPAFNIADSAITVGVGILLLDL
    GLRWLWISFIVLLLDQATKWWASSALELYQPVPVLPYLNFTLLHNPGAAFSFLADAGGWQRWFFTGLAVVISAVIVMWLARLSPQQRMLGIALALVLGGAVGNVIDRFIYGYVIDFIDVYYQTWHWPAFNIADSAITIGAVLLLIDV

`hmmsearch_filt.cut`の方はあとでcontact mapを生成する時に使います。

さて、いよいよCCMpredを使った計算を行います。`-t 12`というのはCPU数を指定します。収束を確実にするために`-n 100`を付けておきます。

    $ ccmpred -t 12 -n 100 hmmsearch_filt.aln hmmsearch_filt.mtx
     _____ _____ _____               _
    |     |     |     |___ ___ ___ _| |
    |   --|   --| | | | . |  _| -_| . |
    |_____|_____|_|_|_|  _|_| |___|___|
                      |_|

    using CPU (12 thread(s))

    Reweighted 12998 sequences with threshold 0.8 to Beff=9296.36 weight mean=0.715215, min=0.0144928, max=1

    Will optimize 9532509 32-bit variables

    iter    eval    f(x)            ║x║             ║g║             step
    1       1       2.41303e+06     25531.6         8.7357614e+09   9.75e-06
    2       1       2.36412e+06     25531           7.3350098e+09   6.19e-06
    3       1       2.31537e+06     25531           6.0129638e+09   5.16e-06
    4       1       2.26682e+06     25533.3         4.817707e+09    4.99e-06
    5       1       2.21858e+06     25538.5         3.8217572e+09   5.31e-06
    6       1       2.17071e+06     25549.4         3.0241239e+09   6.01e-06
    7       1       2.1231e+06      25568           2.3991114e+09   7.08e-06

    .....

    55      1       1.33715e+06     41937.3         55085000        1.95e-05
    56      1       1.3343e+06      41937.2         50249676        2.1e-05
    57      1       1.33154e+06     41931.3         51348760        2.46e-05
    58      1       1.32892e+06     41920.5         48426468        2.7e-05
    59      1       1.32657e+06     41904.9         56314588        3.43e-05

    Done with status code 0 - Success!

    Final fx = 1325345.875000

    xnorm = 93.2845
    Output can be found in hmmsearch_filt.mtx

こんな感じで、2\~5分くらいで返ってきました。GPU acceleratedにするとさらに5倍くらいは早くなりますが、このくらいの計算時間でしたらCPU onlyでも問題にならないくらいです（ただし**CPU並列数指定オプションを忘れると悲惨**）。`hmmsearch_filt.mtx`は例によって対称行列のようになっているはずです。ちなみに、`ccmpred`の`-d`オプションのところでNf計算（number of sequenses \>80% identity）の80%しきい値を他の値に変えられるようです。デフォルトは`-d 0.8`となっているので、ここでは明示的に指定しませんでした。

この`hmmsearch_filt.mtx`を`mk_map.pl`でmap形式に変換します。

    input=hmmsearch_filt ; mk_map.pl -aln ${input}.aln -chk ${input}.chk -cut ${input}.cut -mtx ${input}.mtx -map ${input}.map -do_apc

得られた`hmmsearch_filt.map`はこんな感じ

    LEN 151
    CON 14 138 1.000
    CON 21 119 1.000
    CON 47 115 1.000
    CON 121 128 1.000
    CON 14 101 1.000
    CON 26 117 1.000
    CON 15 108 1.000
    CON 24 28 1.000
    CON 41 120 1.000
    CON 42 121 1.000
    ……

`LEN 151`の151は実効的な配列長です。また、**この残基番号の数字は0-origin**になっていることに気をつけてください。query配列が165アミノ酸であり、50%ギャップの位置情報は失われているので、14アミノ酸分が消し飛んでつながっていることになります。つまり、これだけでは本来のクエリ配列長と照らし合わせることができません。しかし、最初の方で作成した`query_filt.a3m`にはqueryの配列情報がまだ残っています。これを使って上の流れと同様に`query_filt.mtx`と`query_filt.map`を作成します。この`query_filt.map`を見てみると

    LEN 165
    CON 32 131 1.000
    CON 25 151 1.000
    CON 133 141 1.000
    CON 26 120 1.000
    CON 59 127 1.000
    CON 53 132 1.000
    CON 25 113 1.000
    CON 35 39 1.000
    CON 37 59 1.000
    CON 40 46 1.000
    CON 36 40 1.000
    ……

こんな感じでした。**これはGREMLIN ウェブサーバの結果と（0-originのために残基番号が1ずれていることを考慮すれば）同一となっています**。

というわけで、2つのContact map：`hmmsearch_filt.map`と`query_filt.map`が得られました。ここで、`query_filt.map`には**query配列の情報が完全に残っている**ので、この2つをContact mapをうまいことマッチングさせるのに`map_align`を使います。

    map_align -a query_filt.map -b hmmsearch_filt.map > mapalignout.txt

`-a`と`-b`には2つのmapファイルを指定します（順不同）。結果は標準出力で返ってくるので、`mapalignout.txt`に入れておきます。これを見てみると、

    MAX 0_1_0       query_filt.map        hmmsearch_filt.map      119.165 -1      118.165 141     20:9    21:10   22:11   23:12   24:13   25:14   26:15
       27:16   28:17   29:18   30:19   31:20   32:21   33:22   34:23   35:24   36:25   37:26   38:27   39:28   40:29   41:30   42:31   43:32   44:33   45:34   46:35   47:36   48:37   50:38   51:39   52:40   53:41   54:42   55:43   56:44   57:45   58:46   59:47   60:48   61:49   62:50   63:51   64:52   65:53   66:54   67:55   68:56   69:57   70:58   71:59   72:60   73:61   74:62   75:63   76:64   77:65   78:66   79:67   80:68   81:69   82:70   83:71
       84:72   85:73   86:74   87:75   88:76   89:77   90:78   91:79   92:80   93:81   94:82   95:83   96:84   97:85   98:86   99:87   100:88  101:89  102:90  103:91  104:92  105:93  106:94  107:95  108:96  109:97  110:98  111:99  112:100 113:101 114:102 115:103 116:104 117:105 118:106 119:107 120:108 121:109 122:110 123:111 124:112 125:113 126:114 127:115 128:116 129:117 130:118 131:119 132:120 133:121 134:122 135:123 136:124 138:125 139:126 140:127 141:128 142:129 143:130 144:131 145:132 146:133 147:134 148:135 149:136 150:137 151:138 152:139 153:140 154:141 155:142 156:143 157:144 158:145 159:146 160:147 161:148 162:149

こんな感じで出力されます。最後のこの`20:9 21:10`というのが、要はこの2つのcontact mapの対応表となるわけです。これに準じて`hmmsearch_filt.map`の番号を振り直すと……query側の`CON 32 131 1.000`がhmmsearch側の`CON 21 119 1.000`といった具合に、ほぼすべてのquery側の情報がhmmsearch側に含まれていることがわかりました！！

こうして、Uniref100データベースを用いたContact mapをクエリ配列に対応させることができるようになりました。

#### weight scoreの値について

GREMLINウェブサーバーのときに、Import to GREMLIN BETAのところで登場した

    33_S 132_F 4.630 1.00

という行があるのですが、ここのScoreの`4.630`に相当する値は最後の`mk_map.pl`を使ったときに表示されるProbの計算に使うので（GREMLINのWebページの解説にも出てる）、そのアウトプットファイル`query_filt.map`や`hmmsearch_filt.map`を表示させる処理中で存在しているはずです。ということで`mk_map.pl`を調べてみたら、score変数が隠されていました。それを再び表示するためには以下の修正をする必要があります。

`mk_map.pl`を開き、48行目あたりにある次の箇所を修正します。

``` perl
while(exists $CST[$c]){
    my ($i,$j,$sc) = @{$CST[$c]};
    my $prob = sprintf("%.3f",$N_v_H/(1+exp(-$N_v_S * ($sc - $N_v_B))) + $N_v_Y);   
    print MAP "CON $i $j $prob\n";
    $c++;
}
```

これをこう

``` perl
while(exists $CST[$c]){
        my ($i,$j,$sc) = @{$CST[$c]};
        my $prob = sprintf("%.3f",$N_v_H/(1+exp(-$N_v_S * ($sc - $N_v_B))) + $N_v_Y);
        #print MAP "CON $i $j $prob\n";
        # add $newsc to print score for Rosetta distance restraints
        my $newsc = sprintf("%.3f", $sc);
        print MAP "CON $i $j $newsc $prob\n";
        $c++;
}
```

して、改めて`mk_map.pl`を使ったcontactmapの計算を実行すると、例えば上のクエリ例の`query_filt.map`では次のように表示されるようになります。

    LEN 165
    CON 32 131 4.774 1.000
    CON 25 151 3.894 1.000
    CON 133 141 3.490 1.000
    CON 26 120 3.145 1.000
    CON 59 127 3.017 1.000
    CON 53 132 2.924 1.000
    CON 25 113 2.892 1.000

まあGREMLINの結果とちょっとだけ値が少し違うんですけれど、これはデータベースで新しいものを使っているせいだと思いますたぶん。

これでscoreの値を得ることができました。**必ずこの修正をしていたほうが良いです**。

## 一括処理のシェルスクリプト

以下のファイルを一括でやりたい場合は、これらを単に連結するだけでいけるはずです。この処理をまとめたシェルスクリプトを`contactmap.sh`という名前で作成します。そしてコンタクトマップを予測したい配列の入ったFASTA形式のファイルを`query.fa`として、 `./contactmap.sh query.fa` とすれば実行されます。

順調に行けば1ファイルあたり40分〜80分程度で終了します。ファイル名が決めうちされている部分があるので、複数のfastaファイルを同一ディレクトリで実行するのはご遠慮ください。

後の一手のために、Process 9とProcess 10でPythonスクリプトによるファイル変換を行っています。これは https://github.com/BILAB/scripts にてスクリプトを公開していますので、ダウンロードしてお使いください。`pymol_mapalign.py` は次の項のPyMOL表示用プラグインを動作させる上で必要です。

``` bash
#!/bin/sh

PROGNAME=$(basename $0)
query=$1

# set the path to uniclustDB and uniref100DB
uniclustDB="/Users/Agsmith/database/hhsuite/uniclust30_2018_08/uniclust30_2018_08"
uniref100DB="/Users/Agsmith/database/blast/db/FASTA/uniref100/uniref100"

# CPU number. Default is 12.
test ${PBS_NP} || PBS_NP=12

# main process
echo "Process 1 : hhblits to search for sequences"
hhblits -cpu ${PBS_NP} -i $query \
    -d ${uniclustDB} \
    -oa3m query.a3m \
    -n 8 \
    -e 1E-20 -neffmax 20 -nodiff -maxfilt 999999 -realign_max 999999 || exit 1

echo "Process 2 : filtering the sequences obtained from Proc 1"
hhfilter -i query.a3m -o query_filt.a3m -id 90 -cov 75 || exit 1

echo "Process 3 : reformat to a2m format to build hmm model using hmmbuild"
${HHLIB}/scripts/reformat.pl query_filt.a3m query_filt.a2m
hmmbuild --cpu ${PBS_NP} hmmout query_filt.a2m || exit 1

echo "Process 4: search for all similar sequences from uniref100 database and write in hmmsearch.a2m."
hmmsearch --cpu ${PBS_NP} -A hmmsearch.sto \
          --incT 27 hmmout ${uniref100DB} > hmmsearch.out

echo "Process 5: filtering hmmsearch.sto"
${HHLIB}/scripts/reformat.pl hmmsearch.sto hmmsearch.a2m
hhfilter -i hmmsearch.a2m -o hmmsearch_filt.a3m -id 90 -cov 75 || exit 1

echo "Process 6-1: reformat hmmsearch_filt for CCMPred"
input=hmmsearch_filt
a3m2alncut.pl ${input}.a3m ${input}.aln ${input}.cut || exit 1

echo "Process 6-2: execute CCMPred to generate contact map for hmmsearch_filt"
ccmpred -t ${PBS_NP} -n 100 ${input}.aln ${input}.mtx || exit 1
input=hmmsearch_filt
mk_map.pl -aln ${input}.aln -chk ${input}.chk -cut ${input}.cut \
          -mtx ${input}.mtx -map ${input}.map -do_apc

echo "Process 7-1: reformat query_filt for CCMPred"
input=query_filt
a3m2alncut.pl ${input}.a3m ${input}.aln ${input}.cut || exit 1

echo "Process 7-2: execute CCMPred to generate contact map for query_filt"
ccmpred -t ${PBS_NP} -n 100 ${input}.aln ${input}.mtx
input=query_filt
mk_map.pl -aln ${input}.aln -chk ${input}.chk -cut ${input}.cut \
          -mtx ${input}.mtx -map ${input}.map -do_apc || exit 1

echo "Process 8: do map_align to merge the alignments"
map_align -a query_filt.map -b hmmsearch_filt.map > mapalignout.txt

echo "Process 9: renumbering the hmmsearch_filt.map using mapalignout.txt"
# pymol_mapalign.py is called by 'module load map_align'
pymol_mapalign.py hmmsearch_filt.map mapalignout.txt 1> hmmrenumbered.map 2>/dev/null

echo "Process 10: create distance restraint files for Rosetta AbinitioRelax"
# makecst.py is called by 'module load map_align'
makecst.py -f $query
```

## PyMOL表示用プラグイン

この結果をPyMOLで表示するためのプラグインを作成しました。 <https://github.com/BILAB/psico>に便利なPyMOLプラグインを入れております。この中の`contactmap.py`が表示プラグインです。
PyMOLプラグインが入っているpsicoディレクトリの使い方は、まず
    
    git clone https://github.com/BILAB/psico.git

コマンドで中の`psico`ディレクトリを適当な場所にダウンロードしたあと、これについてのPATHをPyMOL側から通しておくことです。例えば`~/apps`にpsicoディレクトリを生成した場合は、`~/.pymolrc`を作成して

    # ~/appsはpsicoディレクトリの場所に合わせる
    sys.path.append(os.path.expanduser('~/apps'))
    import psico.fullinit

としておきます。この状態でPyMOLを立ち上げれば、`contactmap`を始めとしたpsicoの中のプラグインに由来する追加コマンドが使えるようになっているはずです。

PyMOLの窓に`help contactmap`とすると、使い方が表示されます。

    PyMOL>help contactmap

    DESCRIPTION

        visualize the results from map_align algorithm in PyMOL.
        https://github.com/sokrypton/map_align

    ARGUMENTS

        selection = string: atom selection {default: all}

        chain = string: chain ID you want to map {default: ''}

        align = string: name of map_align result txt file {default: 'hmmrenumbered.map'}

        pcut = float: probability cutoff to display {default: 0.999}

        offset = int: slide the sequence number to fit its original number manually {default: 0}

    EXAMPLE

        contactmap 5dir, A, /path/to/hmmrenumbered.map, offset=-7

` contactmap <object name>, <chain ID>`とし、上記Process 9で作成した`hmmrenumbered.map`のパスを適切に入れてやると動きます。`pcut`は表示させるときのProbabilityのカットオフ値です（Probabilityの低いものまで全部表示させていると大変なので）。`offset`は、PDBの番号とFASTAのナンバリングが合わない時にのみ、適当な整数値を与えてやってください。以下の図はPDB: 5DIRでやったときの図です。PDBのナンバリングととquery.faの残基数は7ずれていたため、`offset=-7`にしてあります。

![600px-Gremlin5.png](https://qiita-image-store.s3.amazonaws.com/0/224327/b177d33e-2396-78eb-7f7c-174e1f37a216.png)


こうしてみると、だいたいのCβ-Cβ距離は6\~8 Å以内にあることがわかります。

## 特性
GREMLIN/CCMPredは以下のような特性があります。

- **配列空間が大きければ大きいほど、精度が上昇する**：つまり、タンパク質をコードするアミノ酸配列データベースが増えれば増えるほど、精度が上がっていきます。現在も配列空間はどんどん更新されており、Uniref100に集約されていますので、新しいUniref100データベースを利用すれば精度が増す仕様となっています。
- **クエリ配列に対する類縁配列が多いほど精度が上昇する**: `hhblits`などでマルチプルシーケンスアラインメントが取れた本数だけ精度が上がる仕組みになっています。類縁配列が乏しいタンパク質だとよい精度が出ません。