Introduction

タンパク質の構造情報は、創薬や研究の計画を練る上で現代ではとても重要になってきています。アミノ酸配列が決定されているタンパク質は、2018年12月現在およそ1億くらいあるそうですが（Genbank参照）、それに対してタンパク質構造が決定されているものはおよそ145,000個程度です（重複含む、Protein Data Bank参照）。

いくつかの創薬ターゲットとして有名なタンパク質もまだまだ構造が未決定のものは多いため、これらに対しどうやって阻害剤を探索するかということは1つの困難な課題となっています。

HHPredは、上記のような構造未決定タンパク質のうち、すでに配列や構造上の特徴が似ているということがわかっているタンパク質の場合にとても有効なホモロジーモデリング手法として、とても重宝しています。これはウェブサーバーでも公開されているのですが（参考：http://blog.livedoor.jp/ag_plusplus/archives/68581830.html ←今は画面表示が変わってしまいましたが、やっていることは同じです）、この分野にちょっとだけ専門になってくると、もっと細かな調整を行いながらモデリングをいろいろやってみたいという要望も出てきます。また、このウェブサーバーは世界中からアクセスが多いため、結果が出るまでに少々時間がかかりますので、大量に条件検討してみたい時にはあまり向きません。

そこで、ここでは、HHPredを自分のマシンで完全に動作させるためのセットアップ方法を述べてみます。ここでは、ちょっとMac/Linuxのコマンドライン操作やターミナルの知識が必要となってくるので、こういった操作に慣れていない方には難しいかもしれませんが、一応わかりやすく書いてみようと思います。

環境

マシン環境

以下の環境を想定しています。

macOS High Sierra以降またはLinux OS, CentOS 7.4以降
100GB以上のディスク容量を持つストレージ。SSDであればさらに動作が早い。
類縁配列検索ソフトウェアHH-Suiteのインストールをしてある
Modellerをインストールしてある

HHPredを動作させるためには、これに加えて以下のデータベースのダウンロードが必須です。

データベースのダウンロードとセットアップ

HHPredのために、以下のデータベースをダウンロードしておきます。1つはPDB70のデータベースです。こちらは
ここからpdb70_from_mmcif_latest.tar.gzを選んでダウンロードします。これは1週間に1度更新されるようで、Protein Data Bankのタンパク質構造の更新に合わせてこちらもアップデートしているようです。

もう1つはここにあるuniclust30_2018_08_hhsuite.tar.gzのファイルです。こちらは

Archive containing Uniclust multiple sequence alignments for all clusters in a3m format, generated with Clustal Omega, and additional support files for use with legacy HH-suite version 2 and current version 3.

とあるように、uniclustによって30％で冗長性をカットしたデータファイル、さらにHH-suite（HHPredもこのプログラムスートからの派生）向けにちょっと改変したデータベースになっているそうです。

HHPredプログラムはこの2つのデータベースを参照しながら計算することになるので、必ずこれら2つのファイルのダウンロードが必要です。しかし、この2つはドイツにデータサーバーがある上に、ミラーサーバーが日本にはあまりないようなので、ダウンロードするだけでも丸1日くらい時間がかかります。がんばってダウンロードしましょう。

ダウンロードが終わったら、マシン内の適当な場所にこの2つのデータベースを解凍して置いておきます。ここでは、/Users/Agsmith/database/hhsuiteというディレクトリを作って、ここに解凍したpdb70_latest (最新のファイルなので、この記事の時点ではpdb70_20181125とします)とuniclust30_2018_08というディレクトリを作り、解凍した中身を入れておきます。ダウンロードしてきたtar.gzファイルの解凍のコマンドはtar zxvf uniclust30_2018_08_hhsuite.tar.gzのような形で、pdb70もuniclust30_2018_08の方も行えるはずです。

ディレクトリ/Users/Agsmith/database/hhsuite/uniclust30_2018_08の中身はそれぞれこんな感じ。

lrwxrwxrwx 1 root   root          29 Oct 18 18:01 uniclust30_2018_08_hhm_db -> uniclust30_2018_08_hhm.ffdata
lrwxrwxrwx 1 root   root          29 Oct 18 18:01 uniclust30_2018_08_a3m_db -> uniclust30_2018_08_a3m.ffdata
-rw-r--r-- 1 root   root         767 Oct 12 07:13 uniclust30_2018_08_md5sum
-rw-r--r-- 1 root   root     9427452 Oct 12 07:10 uniclust30_2018_08_hhm_db.index
-rw-r--r-- 1 root   root   436453157 Oct 12 07:10 uniclust30_2018_08_a3m_db.index
-rw-r--r-- 1 root   root  4044015407 Oct 12 07:09 uniclust30_2018_08.cs219
-rw-r--r-- 1 root   root          19 Oct 12 07:09 uniclust30_2018_08.cs219.sizes
-rw-r--r-- 1 root   root     8167108 Oct 12 07:06 uniclust30_2018_08_hhm.ffindex
-rw-r--r-- 1 root   root 14171782777 Oct 12 07:06 uniclust30_2018_08_hhm.ffdata
-rw-r--r-- 1 root   root   375805833 Oct 12 07:05 uniclust30_2018_08_a3m.ffindex
-rw-r--r-- 1 root   root 69490273138 Oct 12 07:05 uniclust30_2018_08_a3m.ffdata
-rw-r--r-- 1 root   root   356545673 Oct 12 07:02 uniclust30_2018_08_cs219.ffindex
-rw-r--r-- 1 root   root  3860466139 Oct 12 07:02 uniclust30_2018_08_cs219.ffdata

と/Users/Agsmith/database/hhsuite/pdb70_20181125の方は

lrwxrwxrwx 1 root   root          16 Nov 26 00:39 pdb70_a3m_db -> pdb70_a3m.ffdata
lrwxrwxrwx 1 root   root          16 Nov 26 00:39 pdb70_hhm_db -> pdb70_hhm.ffdata
-rw-r--r-- 1 root   root 15199373148 Nov 22 03:03 pdb70_from_mmcif_latest.tar.gz
-rw-r--r-- 1 root   root         614 Nov 22 02:58 md5sum
-rw-r--r-- 1 root   root     1658050 Nov 22 02:54 pdb70_hhm_db.index
-rw-r--r-- 1 root   root    16623622 Nov 22 02:54 pdb70.cs219
-rw-r--r-- 1 root   root          14 Nov 22 02:54 pdb70.cs219.sizes
-rw-r--r-- 1 root   root     1783338 Nov 22 02:54 pdb70_a3m_db.index
-rw-r--r-- 1 root   root  2503656590 Nov 22 02:54 pdb70_hhm.ffdata
-rw-r--r-- 1 root   root     1413782 Nov 22 02:54 pdb70_hhm.ffindex
-rw-r--r-- 1 root   root     1539070 Nov 22 02:54 pdb70_a3m.ffindex
-rw-r--r-- 1 root   root 39682186995 Nov 22 02:54 pdb70_a3m.ffdata
-rw-r--r-- 1 root   root    16005165 Nov 22 02:40 pdb70_cs219.ffdata
-rw-r--r-- 1 root   root     1168208 Nov 22 02:40 pdb70_cs219.ffindex
-rw-r--r-- 1 root   root     6844342 Nov 21 10:33 pdb70_clu.tsv
-rw-r--r-- 1 root   root    17932168 Nov 21 10:33 pdb_filter.dat

こんな感じになるはず。時刻の部分はともかく、ファイルが過不足なく存在することと、ファイルパーミッションが全部-rw-r--r--とかlrwxrwxrwxとなっていることも確認しておきましょう。

自前マシンでのHHPredを使ったホモロジーモデリング

Step.1 HH-suiteを使ったテンプレートタンパク質構造の検索

HH-suiteのインストールとModellerのインストールは、上の環境のところでリンクを張ってインストール方法を紹介していますので、ここでは省略します。

以上のセットアップが揃ったら、これでHHPredの実行をマシン上で行うことができるようになっています。早速やってみましょう。

例として以下のようなquerytest.faのFASTA形式ファイルを用意してみます。

>hogehoge
PRGSHMASMTGGMIEVPPFWCPLPIAIHPAADQAEKDARAWAERYGVRLRIAD
QVQPGRLGAYWAPHGTYEGMLAVGCWNFWAFAFDDHLDEPLPLDVPVTTSLVQQAVDIPSPPITDDPW
AAGAQAVFNMFRDLATPTQVRYCADNHRRWLHGACWRHSNHVNRRLPPLAEYIPLRMQDAAAQATCLI
AVLIGSDISVPEQEMDSPRVRALLETASWTATIDSDLHSFQLEDTQRPVSQHIVSVLMHERGIGVDEA
LRQSVALRDRFMTRFLHLQQECARTGSSELARFAHTLGYVISGYLQWAVDTSRYGQTEATFSFTDTPR
DDTPEPPPGIPSVEWLWTL

次に以下のコマンドを入力します。

hhblits -cpu 12 -i querytest.fa \
        -d /Users/Agsmith/database/hhsuite/uniclust30_2018_08/uniclust30_2018_08 \
        -oa3m query.a3m \
        -n 3 -e 0.001 -loc -norealign \
        -cov 20 -qid 0

-iにはインプットとなるクエリファイル、-dには参照するデータベース、-oa3mには結果のマルチプルシーケンスアライメントをアウトプットさせるa3mフォーマット形式ファイル名(write result MSA with significant matches in a3m format)、-cpuは計算させるコア数、-nにはiteration回数(1-8)の指定です。最高のパフォーマンスのためには-cpuにそのマシンの最大CPU数を与えていることが良いでしょう。

続いて、hhsearchコマンドを用いて、上述のアウトプットのa3mファイルをインプットとして、PDB70データベースからそれに合致するタンパク質構造を取ってきます。

export HHLIB="/Users/Agsmith/apps/hh-suite/3.0b3"
# add Secondary structure information into query.a3m
${HHLIB}/scripts/addss.pl query.a3m
hhsearch -cpu 12 -i query.a3m -o results.hhr -oa3m results.a3m \
         -d /Users/Agsmith/database/hhsuite/pdb70_20181125/pdb70 \
         -p 20 -Z 250 -loc -z 1 -b 1 -B 250 -ssm 2 -sc 1 \
         -seq 1 -dbstrlen 10000 -norealign -maxres 32000 \
         -contxt ${HHLIB}/data/context_data.crf \
         -cov 20 -qid 0

${HHLIB}の部分はHH-suiteがインストールされているディレクトリになります。これにより得られたresults.hhrに、すでに構造が解かれているタンパク質構造と類縁関係にあると思われるPDB
IDのリストが表示されます。

得られたresults.hhr結果は以下の通り(2018年7月30日時点のもの)。

Query         hogehoge
Match_columns 344
No_of_seqs    155 out of 1362
Neff          10.4876
Searched_HMMs 59613
Date          Mon Jul 30 21:54:30 2018
Command       hhsearch -cpu 12 -i query.a3m -o results.hhr -oa3m results.a3m -d /Users/Agsmith/database/hhsuite/pdb70_latest/pdb70 -p 20 -Z 250 -loc -z 1 -b 1 -B 250 -ssm 2 -sc 1 -seq 1 -dbstrlen 10000 -norealign -maxres 32000 -contxt /Users/Agsmith/apps/hh-suite/3.0b3/data/context_data.crf -cov 20 -qid 0

 No Hit                             Prob E-value P-value  Score    SS Cols Query HMM  Template HMM
  1 4OKZ_C Terpene synthase metal- 100.0 5.5E-49 9.2E-54  354.0  36.2  332   11-344     2-347 (365)
  2 4OKZ_A Terpene synthase metal- 100.0 5.5E-49 9.2E-54  354.0  36.2  332   11-344     2-347 (365)
  3 4MC3_A Putative sesquiterpene  100.0 2.8E-44 4.8E-49  321.4  33.2  317   12-331     2-336 (346)
  4 4MC8_A Putative sesquiterpene  100.0 2.8E-44 4.8E-49  321.4  33.2  317   12-331     2-336 (346)
  5 5I1U_B Germacradien-4-ol synth 100.0   3E-44   5E-49  319.9  30.9  308    4-314    11-335 (335)
  6 5I1U_A Germacradien-4-ol synth 100.0   3E-44   5E-49  319.9  30.9  308    4-314    11-335 (335)
  7 5DZ2_B Germacradienol/geosmin  100.0 8.3E-44 1.4E-48  317.5  30.9  306    9-317     1-332 (338)
（以下省略）

ここにはhhblitsやhhsearchコマンドによって、querytest.faに入力していたタンパク質に近いと判定された構造既知のタンパク質のリストが現れます。4OKZ_Cとか4MC3_Aとかの名前が現れていますが、これはPDB IDと、チェインIDが続いた表記になっています。つまり、PDB: 4OKZのCチェインということです。

ここで、続くmodellerを使ったホモロジーモデリングに向けて、PDB_chain名前のモデル番号を控えておきます。ここでは例として1, 3, 5を使用するとします。上の例のように、E-valueのとても小さいものがヒットしている場合は1つだけでもOKです。

Step.2 Modellerを使ったホモロジーモデリング

HHPredによって構造既知の類縁タンパク質を推定した後は、Modellerの機能を使って、そのアミノ酸配列を使ってアライメントを取り、モデリングを行います。

startmodeller.shは自作のModeller用前処理＆モデリング実行の汎用スクリプトです。下記のシェルスクリプトをコピーしてこの名前のスクリプトを作成した後、chmod +x startmodeller.shを実行しておきます（一度やれば以降は不要）。

#!/bin/sh

# 各自のHH-SuiteへのPATHを指定する
export HHLIB="/Users/Agsmith/apps/hh-suite/3.0b3"

PROGNAME=$(basename $0)
VERSION="1.0"
HELP_MSG="'$PROGNAME -h'と指定することでヘルプを見ることができます"

red=31
green=32
yellow=33
blue=34

function cecho {
    color=$1
    shift
    echo -e "\033[${color}m$@\033[m"
}

# ヘルプメッセージ
usage() {
  echo "Usage: $PROGNAME -i input.hhr -p output.pir [-d arg] param"
  echo
  echo "Options:"
  echo "  -h, --help"
  echo "      --version"
  echo "  -i, --inputhhr <ARG>     <必須>"
  echo "  -p, --pir   <ARG>     <必須>"
  echo "  -d, --direcotry <string>  [任意] *オプション指定時は引数を設定"
  echo "  -m, --model   テンプレートにするモデル番号の選択 e.g.\"-m 1 3\" <必須>"
  echo "  -n, --num     作成するモデルの個数。default is 1."
  echo
  exit 1
}

ARG_D="template"
ARG_N=1

# オプション解析
for OPT in "$@"
do
  case "$OPT" in
    # ヘルプメッセージ
    '-h'|'--help' )
      usage
      exit 1
      ;;
    # バージョンメッセージ
    '--version' )
      echo $VERSION
      exit 1
      ;;
    # オプション-i、--inputhhr, HHRファイルの指定
    '-i'|'--inputhhr' )
      FLG_I=1
      # オプションに引数がなかった場合（必須）
      if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
        cecho $red "$PROGNAME:「$1」オプションにはhhrファイルの引数が必要です。" 1>&2
        exit 1
      fi
      ARG_I="$2"
      shift 2
    ;;
    # オプション-p, --pir, PIRファイルの指定
    '-p'|'--pir' )
      FLG_P=1
      # オプションに引数がなかった場合（必須）
      if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
        cecho $red "$PROGNAME:「$1」オプションにはPIRファイルの引数が必要です。" 1>&2
        exit 1
      fi
      ARG_P="$2"
      shift 2
    ;;
    # オプション-d, --directory, direcotryの指定
    '-d'|'--directory' )
      FLG_D=1
      # オプションに引数がなかった場合
      if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
    shift 1
      else
        # オプションの引数設定
        ARG_D="$2"
        shift 2
      fi
    ;;
    # オプション-m
    '-m'|'--model' )
      # 配列型を引数にとる(-m 1 3 のような感じ)。
      FLG_M=1
      # オプションに引数がなかった場合（必須）
      if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
        cecho $red "$PROGNAME:「$1」オプションにはhhrファイルの引数が必要です。" 1>&2
        exit 1
      fi
      ARG_M="$2"
      shift 2
      ;;
    # オプション-n, --num, modellerのモデル作成個数の指定
    '-n'|'--num' )
      FLG_N=1
      # オプションに引数がなかった場合
      if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
    shift 1
      else
        # オプションの引数設定
        ARG_N="$2"
        shift 2
      fi
    ;;
    '--'|'-' )
      # 「-」か「--」だけ打った時
      shift 1
      param+=( "$@" )
      break
      ;;
    -*)
      cecho $red "$PROGNAME: 「$(echo $1 | sed 's/^-*//')」オプションは存在しません。'$PROGNAME -h'で確認してください" 1>&2
      exit 1
      ;;
    *)
      # コマンド引数（オプション以外のパラメータ）
      if [[ ! -z "$1" ]] && [[ ! "$1" =~ ^-+ ]]; then
        param+=( "$1" )
        shift 1
      fi
      ;;
  esac
done

# 「-i」のオプション指定がない場合
if [ -z $FLG_I ]; then
  cecho $red "$PROGNAME:「-i」オプションは必須です。正しいオプションを指定してください" 1>&2
  echo $HELP_MSG 1>&2
  exit 1
fi
# 「-p」のオプション指定がない場合
if [ -z $FLG_P ]; then
  cecho $red "$PROGNAME:「-p」オプションは必須です。正しいオプションを指定してください" 1>&2
  echo $HELP_MSG 1>&2
  exit 1
fi
# 「-m」のオプション指定がない場合
if [ -z $FLG_M ]; then
  cecho $red "$PROGNAME:「-m」オプションは必須です。正しいオプションを指定してください" 1>&2
  echo $HELP_MSG 1>&2
  exit 1
fi

# copy template mmCIF files
preparepdb() {
    mkdir -p ${ARG_D}
    modelarray=(${ARG_M})
    PDBID_dict=()
    for i in ${modelarray[@]};do
    uprange=`grep "No Hit" ${ARG_I} -n | sed -e "s/:.*//g"`
    line=`sed -n $((${uprange}+${i}))p ${ARG_I}`
    echo "Using: ${line}"
    PDBID=`echo $line | awk '{print $2}' | cut -c 1-4 | tr '[A-Z]' '[a-z]'`
    divided=`echo $PDBID | cut -c 2-3`
    curl --silent ftp://ftp.pdbj.org/pub/pdb/data/structures/divided/mmCIF/${divided}/${PDBID}.cif.gz ./${ARG_D}/
    gzip -df ./${ARG_D}/${PDBID}.cif.gz
    PDBID=`echo $PDBID | tr '[a-z]' '[A-Z]'`
    PDBID_dict+=(${PDBID})
    done
}


makealignpy() {
    knowns=""
    for i in ${PDBID_dict[@]};do
    knowns+="'${i}', "
    done
    cat <<EOF > align.py || exit 1
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from modeller import *
from modeller.automodel import *    # Load the automodel class

log.verbose()
env = environ()

# directories for input atom files
env.io.atom_files_directory = ['.', './template']

a = automodel(env,
              alnfile  = '${ARG_P}', # alignment filename
              knowns   = (${knowns}),     # codes of the templates
              sequence = 'UKNP')               # code of the target

a.starting_model= 1                 # index of the first model
a.ending_model  = ${ARG_N}          # index of the last model
a.make()                            # do the actual comparative modeling

EOF
}


preparepdb
${HHLIB}/scripts/hhmakemodel.py ${ARG_I} ${ARG_D}/ ${ARG_P} ./ -m ${ARG_M}
makealignpy
mod9.19 align.py

このstartmodeller.shの処理の流れとしては、

preparepdbの関数定義で、PDBファイルをmmcif形式でダウンロードし、templateディレクトリに解凍する
${HHLIB}/scripts/hhmakemodel.pyを使ってuknp.pirファイルに配列アライメントを取る。環境変数${HHLIB}が正しく設定されていることを確認
makealignpyの関数定義で、Modeller用実行スクリプトをalign.pyファイルに書き出す
Modeller実行コマンドmod9.19で、上記align.pyを読み込んでモデリングを行う。mod9.19の部分はバージョンによって名前が異なるので、自身が使っているコマンド名に書き換えてください。

となっています。
とりあえずここでは、align.pyで作り出す実行スクリプトはお決まりのものになっていますが、Modellerの使い方に熟知していれば、ここでモデリングをする時に2点の原子間に距離拘束を入れた条件を付与できたり、モデルをたくさん作り出したりすることができるようになります。HHPredのウェブサーバー版ではこの機能にはついておらず、またモデル数も1に限定されているはずです。これが、自前でHHPredを用いるメリットの1つと言えるでしょう。

uknp.pirは配列アライメントのファイルとなっています。Modellerはこれを読んで、中にあるそれぞれのテンプレート構造を利用してモデルを作成しようと試みます。

こうして作ったstartmodeller.shとresults.hhrと同じディレクトリに置いたら、ホモロジーモデリングをするためのスクリプトmodeller.shを実行します。同様に下記のスクリプトをコピペして作成後、chmod +x modeller.shを実行しておきます（一度やれば以降は不要）。

#!/bin/sh

test $hhr || hhr="results.hhr"
test $pir || pir="uknp.pir"
test $num || num=1

# startmodeller.sh -i results.hhr -p uknp.pir -m "1 3 5" -n 2
bash startmodeller.sh -i ${hhr} -p ${pir} -m "${model}" -n ${num}

これを作成し、model="1 3 5" num=3 ./modeller.shとして実行します。num=の値を変えると生成されるモデルの数が変化します。

Modellerで動かした後のログファイルalign.logの最後の方には、作り出したモデル構造のスコア値が記載されています。

>> Summary of successfully produced models:
Filename                          molpdf
----------------------------------------
UKNP.B99990001.pdb            9450.61133
UKNP.B99990002.pdb            9513.01562
UKNP.B99990003.pdb            9432.22168

このmolpdfの値が小さいほど、よいモデリングであるとされているようです。この値はモデルのテンプレートの選び方で大きく変わってきますし、距離拘束などをうまく取り入れることで変化します。また、Modellerには少し構造を緩和させるような処理なども追加で加えることもできるので（この辺はModellerのマニュアルやチュートリアルなり他のウェブサイトなりを見てください）、それらをうまく使ってより良い値になるようなモデルが引けるまでいろいろ試してみると良いでしょう。

おわりに

自前マシンでHHPredを使ったホモロジーモデリングの方法を紹介してみました。
これは実際ウェブサーバー版だけで作り出したモデル構造でも（とても近い類縁配列がある場合には）かなり強力で、これひとつで十分に考察することが可能です。
しかし、ホモロジーが存在しないものについては結構無力なこともあるので、ケースバイケースです。
次回の記事では、そういったホモロジー構造がないものでも、強力にタンパク質モデリングをしてくれる方法を紹介したいと思います。