Yang Zhang研究室で開発されているResPREというコンタクトマップ予想プログラムの紹介です。
論文はBioinformatics誌上で2019年5月9日に発表されたこれ。GitHubはここ。Webサーバーはここ。
コンタクトマップについてのIntroduction
近年のタンパク質の構造予測手法のトレンドとして、2010年代から再注目され始めたタンパク質のコンタクトマップを使った構造予測手法があります。コンタクトマップはタンパク質の3次元構造を2次元で表現する方法の1つです。例として、PDB: 5DIRのA chainはこんな感じの画像になります。
上の右図が左図に対応するコンタクトマップで、タンパク質を構成するアミノ酸のCα原子間の距離（単位：nm）情報を2次元的にヒートマップで示しています。性質上、同一番号の交差点（図における対角線上）の距離は0となり、対角線を挟んで線対称に必ずなります。
コンタクトマップ予想についての過去の研究
コンタクトマップに含まれる情報は、タンパク質構造モデリングを行う上でのガイドとしての距離拘束として使うことができることが2000年から報告されています。そして、タンパク質構造を配列のみからモデリングするde novoモデリングというグランドチャレンジにも使われています。
コンタクトマップの予測アルゴリズムはこの10年で大きく進歩しており、現在の手法の主な流れは3つに分けられます。
1. 共進化情報ベースの手法
2. 機械学習ベースの手法
3. メタベースの手法
です。とはいえ、最近はこれらの境界は曖昧になりつつあり、一例を上げるなら、最近は多くのディープラーニング手法が共進化情報ベースの手法を採用しつつあります。
各手法を紹介していきます
1. 共進化情報ベースの手法
これは「共進化の原理」に着想を持ったもので、もしタンパク質構造のある1ヶ所の残基の変異が入った場合、構造上で近位のアミノ酸残基がその変異に対応するように同時に変異すれば、結果としてそのタンパク質の構造自体が保たれるだろう、という原理に基づいたものです。
現在配列データベースGenbankには2億以上のタンパク質をコードする配列が登録されているので、任意のあるアミノ酸配列に対し、そのデータベースの中からそれに近い配列をpsiblastやhhblitsを使うことでマルチプルシーケンスアライメント（MSA）を作成することができます。この手法の代表例としてpseudo-likelihood maximizationアルゴリズムを用いたplmDCAや、GREMLINとそのGPU/並列実装版のCCMpred（DOI: 10.1093/bioinformatics/btu500）、平均場近似 (mean-field approximation)を使うEVfoldやmfDCA、sparse inverse of covariance matrix（スパースな精度行列、グラフィカルLASSO的なアレ？）を用いるPSICOVなどがあります。
2. 機械学習ベースの手法
コンタクトマップ予測をパターン認識の問題として見て機械学習手法で解くものです。2007〜2008年頃はSupport Vector Machineを用いた手法がこのカテゴリに属していましたが、近年ではディープラーニングの台頭がめざましく、その予測精度を大幅に上昇させることができたことを報告しています。plmConv (2016)やDeepCov (2018)、RaptorX-contact (2017)、DNCON2 (2018)、SPOT-contact (2018), DeepContact (2018), DeepConPred (2018)などが雨後の筍のように報告されています。
3. メタベースの手法
異なる予測手法をうまく補完するように組み合わせる手法のことで、LRcon (2011), R2C (2016), MetaPSICOV (2015), PconsC2 (2014), NeBcon (2017)などが存在する。メタ予測器では元となる予測器よりも良い精度が達成されたことが報告されていますが、その効果は限定的であるような示唆が、それぞれの論文から伺えています。
現在、ほとんどすべてのコンタクトマップ予測手法は十分な数の冗長性を排除したホモログ配列のMSAを取得可能であることに依存している。特に、共進化情報ベースの場合は顕著です。
ResPREのWebサーバー版の利用
ResPREはこのコンタクトマップの予想手法やタンパク質構造予測を研究するのバイオインフォマティクス分野で有名なZhang教授から発表されました。Webサーバー版があるのでさっそくこれを使ってみましょう。 https://zhanglab.ccmb.med.umich.edu/ResPRE/ にアクセスします
やり方は簡単、"Cut and paste your sequence"と書かれているところに予測してもらいたいアミノ酸配列をコピー＆ペーストで貼り付けるだけです。Emailは必須で、結果がここに指定したメールアドレスに送られてきます。だいたい3〜5時間くらいです。IDの欄は任意です。
終わったら、Submitボタンを押して結果を待つだけです。
結果はこんな感じのフォーマットで送られてきます。
PFRMAT RR
TARGET PDB:5go5
AUTHOR ResPRE
REMARK The full prediction can be viewed at https://zhanglab.ccmb.med.umich.edu/ResPRE/output/RPCXX
METHOD contact prediction by deep residual neural networks
MODEL 1
YTYGKPFAVMYIPRLGFTWNKPVLEGTGTEVLKKGLGHYANTARLGQKGN
FAVAGHRRTYGDPFKDFPKLRHGDEVVLTDGTTWFTYVIDTGPYKTVPTD
VEVIDPVPRKSGYEREGRYLTLTTCEPEWGHSHRLIVWAHLDSTQPVEAG
KPEALRR
i j d1 d2 p
87 139 0 8 0.9987920000
89 139 0 8 0.9979270000
76 87 0 8 0.9974700000
75 88 0 8 0.9971590000
90 138 0 8 0.9970550000
87 141 0 8 0.9967210000
……
i j d1 d2 pより下の行で、左2つのi jがコンタクトすると予想されるアミノ酸の番号（入力したアミノ酸配列に対応）、d1, d2部分は0, 8で固定されているようであまり意味はなく、pの列がprobabilityを表しています。基本的にはこの数字が大きいほど自信ありげなようです。
また、一応コンタクトマップ予想の結果画像も作ってくれますが、これはprobabilityが0.5以上のものについて黒、それ以下は白で表してくれるようです。
これだとちょっとわかりにくいので、probabilityの大小に応じて色を付けて正解の距離マップと比較してみます。
ここで左の正解画像はGromacsに付属しているgmx mdmatとgmx xpm2psコマンドで生成することができます。Gromacs 2019だとなぜかこのコマンドがバグっていて使えないので2018シリーズを使いました。
# 正解の距離マップをPDBファイルから生成。
# 先に5go5のA chainだけのPDBファイルを5go5A.pdbとして用意しておく。
# gmx mdmatで距離行列を生成する。-tオプションで距離のカットオフ（1.5 nm = 15 Å）
# コマンド入力後、Select group for analysisを追加で尋ねられるので、'C-alpha'または'3'を選択して終了する
gmx mdmat -f 5go5A.pdb -s 5go5A.pdb -mean 5go5A.xpm -t 1.5 << EOF
C-alpha
EOF
# gmx xpm2psでは-noframeをつけると目盛りを消せます
gmx xpm2ps -f 5go5A.xpm -o 5go5A.eps -rainbow red -bx 3 -by 3 -noframe
色合いのことを無視していただければ、だいたいのところでコンタクトすべき残基ペアがうまく予想できており、模様がだいたいそっくりなんじゃないかな、と言えると思います。
ちなみにPDB: 5go5はこんな感じのβシートリッチな構造です。
PyMOLなどでコンタクトマップと3次元構造を比較してみて、どんな模様がどんな構造に対応しているかを理解するのも良いでしょう。
ResPREをローカルマシンにインストールして利用してみる
この記事の本題です。上記のResPREは割と精度良くコンタクトマップを予想してくれるのですが、なかなか時間がかかってしまうので、大量に予想したいという人には向いていませんし、「アミノ酸配列の情報自体を外に出すのも躊躇われる……」という企業研究者の方もいらっしゃるかもしれません。そこで、作者が公開してくださっているResPREのプログラムを、自前のマシン環境にインストールして動かしてみる方法を紹介します。
マシン環境
OSはLinuxまたはWindows 10のWSL上にインストールされたUbuntuでならば動作します。CentOS 7.4, Ubuntu 18.04で動作を確認しました。macOSの場合はOSの仕様上の関係で、途中で使うことになるプログラムbin/calNf_lyが動作せず、結果としておかしな値を返します。これについてはnoah というパッケージを利用すればmacOSでもbin/calNf_lyをきちんと動作させられるようになるので回避できます（かんたんな使い方はいつか書きます）。
さらにgit, Python3, Pytorchのバージョン0.3.1を入れておきます。Pytorchのバージョンが1.1.0では動きませんでした。PyTorchのインストール方法については https://pytorch.org/get-started/previous-versions/ を読んでください。pip3で入れる場合はこんな感じ。
pip3 install torch==0.3.1 numba
書いてありませんが、 numba パッケージも動作に必要です。
ResPRE自体は上記の環境で動作しますが、任意のアミノ酸配列でコンタクトマップ予測を行いたい場合は、自前でMSAを作るためにHH-suite（に含まれるコマンドhhblits）のインストールが必要になります。これは類縁配列検索ソフトウェアHH-Suiteのインストール方法に書いています（結構時間がかかります）。
ResPREの予測結果の精度はMSAの質に大きく依存します。また、MSAは形式自体が整っていれば究極的にはどうやって生成しても良いので、hhblitsを使ってMSAを生成する以外の方法もあると思います。
インストール
pip3でPyTorchの古いバージョン（推奨は0.3.1付近）とnumba, numpy, scipyをインストールしておきます。一応、pip3 listコマンドでバージョンを確認しておきます。
$ pip3 list
（略）
numba 0.44.0
（略）
torch 0.3.1
（略）
このバージョンなら動作します。
その後、以下のコマンドでResPRE本体のインストールを行います。以下の例では~/appsにResPREをインストールします。
# インストールディレクトリに移動
mkdir -p ~/apps ; cd ~/apps
# gitでResPREをダウンロード
git clone https://github.com/leeyang/ResPRE.git
cd ResPRE
# マニュアルに書いていないですがbinディレクトリにあるプログラムを実行可能にしないと動作してくれない
chmod +x bin/calNf_ly
これで完了です。テスト動作コマンドは、このResPREディレクトリ内で
python3 respre.py test/T0759-D1.aln output.out
を実行するだけです。しばらくすると、output.outが生成されるはずです。この結果とtest/T0759-D1.outが一致していれば成功です。
hhblitsを使ってResPREに与えるMSAを生成する方法
ダウンロードしてきたResPREディレクトリの中にあるtest/T0859-D1.alnはこんな内容になっています。
VVIHPDPGRELSPEEAHRAGLIDWNMFVKLRSQE
VVIHPDTGRELSPEEAHRVGLIDWNMFVKLRSQE
VVIHPDTGRELSPEEAHRAGLIDWNMFVKLRSQE
VVIDPDTGRELSPEEAHRAGLIDWKMFVKLRSQE
S----------STPTRHGAGLIDRAMFDRLRGQE
---------LWSSTLTHRVGLIDWNMFVKLRSQE
VIVDPETNKEMSVQEAYKKGLIDYETFKELCEQE
VIVDPDTGKEMTVREAYHKELIDYDTFLDLSEQE
VVVHKDTGKEMTLREAYHSELIDYNTFLELSEQE
VIVDPETNKEMSVREAYHKDLIDYETFMELSEQE
VVIDPDTGRELSPEEAHRAGLIDWNMFVKLRSQE
VVIDPDTGREMSPEEARRTGLIDWSMFVKLQNQE
VVIDPETGKEMRPEEAYKLGLIDWKMFVNLQSQE
LVIDPGTSKEMTPEEAYSQGLINLAMFVHLQSQE
VVIDPDTGHELSPEEARRAGLIDWSMFVKLRSQE
VVIDPDSGRELSPEEAHRAGLIDWNMFVKLRSQE
（以下略）
こんな感じで長さが揃っているMSAをインプットとしてResPREに与えてやれば、それだけでコンタクトマップを生成してくれるわけです。ちなみに確か一番上はクエリの配列になっていると思います。
任意の配列に対してコンタクトマップ予測をやってみたい場合はこのようなMSAを自前で生成する必要があるのですが、ここではhhblitsを使った方法を紹介します。
準備は類縁配列検索ソフトウェアHH-Suiteのインストール方法に書いています。準備に必要なデータベース（UniclustとかUniref100とか）のダウンロードが丸1日くらいかかるかもしれませんが根気よく待ってください。
下にローカル版ResPREを使ってコンタクトマップを自前で作るスクリプトhhblits.shの一例を紹介しておきます。（データベースのPATHを設定してあげる必要があります。）Step 6のところはMSAの長さを切りそろえておく処理です。
#!/bin/sh
PROGNAME=$(basename $0)
# 解析したい配列が入ったFASTAファイルを指定する。例：query=5go5.fa
query=$1
# querynameにはqueryの名前から拡張子を外した文字列が入る 例：query=5go5.fa -> queryname=5go5
queryname=${query%.*}
# PATHの設定。
uniclustDB="/path/to/database/hhsuite/uniclust30_2018_08/uniclust30_2018_08"
uniref100DB="/path/to/database/FASTA/uniref100/uniref100"
respredir="/path/to/ResPRE" # ResPREをインストールしたディレクトリの場所
# 使用するCPUコア数。環境に合わせて変更する。
PBS_NP=4
# main process
echo "Process 1 : hhblits to search for sequences"
hhblits -cpu ${PBS_NP} -i $query \
-d ${uniclustDB} \
-oa3m ${queryname}.a3m \
-n 3 \
-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 ${queryname}.a3m -o ${queryname}_filt.a3m -id 90 -cov 75 || exit 1
echo "Process 3 : reformat to a2m format to build hmm model using hmmbuild"
reformat.pl ${queryname}_filt.a3m ${queryname}_filt.a2m
hmmbuild --cpu ${PBS_NP} ${queryname}hmmout ${queryname}_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 ${queryname}hmmsearch.sto \
--incT 27 ${queryname}hmmout ${uniref100DB} > ${queryname}hmmsearch.out
echo "Process 5: filtering hmmsearch.sto"
reformat.pl ${queryname}hmmsearch.sto ${queryname}hmmsearch.a2m
hhfilter -i ${queryname}hmmsearch.a2m -o ${queryname}hmmsearch_filt.a3m -id 90 -cov 75 || exit 1
echo "Process 6: generate clean MSA with the same sequence length"
# query_filt.hhbaln will keep the query sequence.
grep -v '^>' ${queryname}_filt.a3m | sed 's/[a-z]//g' > ${queryname}_filt.hhbaln
grep -v '^>' ${queryname}hmmsearch_filt.a3m | sed 's/[a-z]//g' > ${queryname}hmmsearch_filt.hhbaln
echo "Process 7: ResPRE"
python3 ${respredir}/respre.py ${queryname}_filt.hhbaln ${queryname}_pre.map
python3 ${respredir}/respre.py ${queryname}hmmsearch_filt.hhbaln ${queryname}hmmsearch_pre.map
queryseq=`grep -v '^>' ${query} | perl -pe 's/\n//g' | fold`
cat << EOS > ${queryname}.map
PFRMAT RR
AUTHOR ResPRE
METHOD contact prediction by deep residual neural networks
MODEL 1
$queryseq
i j d1 d2 p
EOS
cat ${queryname}_pre.map >> ${queryname}.map
cat << EOS > ${queryname}hmmsearch.map
PFRMAT RR
AUTHOR ResPRE
METHOD contact prediction by deep residual neural networks
MODEL 1
$queryseq
i j d1 d2 p
EOS
cat ${queryname}hmmsearch_pre.map >> ${queryname}hmmsearch.map
これと任意のFASTA形式ファイル（例として5go5.fa）を同じディレクトリに用意します。
> 5go5
YTYGKPFAVMYIPRLGFTWNKPVLEGTGTEVLKKGLGHYANTARLGQKGN
FAVAGHRRTYGDPFKDFPKLRHGDEVVLTDGTTWFTYVIDTGPYKTVPTD
VEVIDPVPRKSGYEREGRYLTLTTCEPEWGHSHRLIVWAHLDSTQPVEAG
KPEALRR
この2つを用意して nohup ./hhblits.sh 5go5.fa > nohup.out & とターミナルから打てば動作します。十分なパワーとメモリを積んだマシンならば30分くらいで結果が帰ってきます。あとはそいつを画像に変換するなりすればOKでしょう。