類縁配列検索ソフトウェアHH-Suiteのインストール方法

Johannes Söding Labが開発しているという、タンパク質配列探索プログラムHHsuiteのインストール方法について。オープンソフトウェア。

GitHubはここ https://github.com/soedinglab/hh-suite

## HH-Suiteとは（Google翻訳）
HH-Suiteには、HHsearchとHHblitsが含まれています。HHsearchは、遠隔関連配列の検出およびタンパク質構造予測のための最も一般的な方法であり、Google Scholarによれば現在2000回以上引用されています。HHsearchプログラムとHHblitsプログラムは、クエリとデータベースシーケンスの両方がマルチプルシーケンスアライメント（Multiple Sequence Alignment）によって表されるという事実に、その強みがあります。これらのMSAにおいて、クエリ配列またはデータベース配列は、同種の（関連する）配列と一緒に表に書き込まれ、各列は、同族のアミノ酸残基、すなわち先祖配列中の同じ残基に由来する残基を含みます。そのようなMSAのカラム中のアミノ酸の頻度は、その位置でさらなる相同配列内のアミノ酸を観察する確率として解釈することができます。潜在的な配列のMSA内の配列との関連性のための自動スコア付けを容易にするために、MSAはプロファイル隠れマルコフモデル（HMM）によって簡潔に記述されます。これらは、位置別スコアリングマトリックス（PSSM）の拡張と言えます。

**HHsearch**は、複数の配列アライメントまたはプロファイル隠れマルコフモデル（HMM）を入力として受け取り、同種の（関連する）タンパク質のプロファイルHMMのデータベースを検索します。HHsearchは、ホモロジーモデリング、すなわち、配列のみが既知であるクエリータンパク質の構造のモデルを構築するためにしばしば使用されます。例として、Protein Data Bank (PDB)のような既知の構造を有するタンパク質のデータベースがクエリタンパク質に類似した「鋳型」タンパク質をHHsearch中で検索対象のデータベースとして指定されます。そのような鋳型タンパク質が見出された場合、クエリと鋳型タンパク質配列との対の配列アライメントに基づいて、目的タンパク質の構造を予測することができます。2010年のCASP9タンパク質構造予測コンペティションでは、HHsearchとHHblitsに基づく完全に自動化されたHHpredが、テンプレートベースの構造予測CASP9 TBM / FMで81台のサーバーの中で最も評価されました。

**HHblits**は2011年にHHsuiteに追加されました。単一のクエリシーケンスまたはMSAから高品質の複数シーケンスアライメント（MSA）を構築できます。クエリから、プロファイルHMMを計算することができます。単一配列の代わりにMSAを使用することにより、配列検索の感度および得られた配列アライメントの品質を劇的に改善することができる。
MSAはまた、タンパク質の二次構造および三次構造を予測し、分子機能または細胞経路を予測し、酵素活性に寄与する配列または構造の位置を予測するための方法など、多数の下流計算方法の出発点でもある進化的に保存された残基、疾患を引き起こす変異と中立的な変異、タンパク質の細胞局在化などを予測することができます。これは、最高品質のMSAを製造することの重要性を説明しています。

**HHblitsは、最も一般的な反復配列検索方法であるPSI-BLASTと同様に機能します。** HHblitsは、クエリシーケンスからプロファイルHMMを生成し、HHスイートのuniprot20データベースなどのプロファイルHMMの大きなデータベースを繰り返し検索します。uniprot20データベースには、UniProtデータベースで収集されたすべてのパブリックな高品質タンパク質配列が含まれています。これらの配列はクラスタ化され、複数の配列アラインメントに整列され、そこからuniprot20のプロファイルHMMが生成されます。以前の検索と有意に類似した配列が、次の反復検索にむけてクエリプロファイルHMMに追加されます。**PSI-BLASTおよびHMMERと比較して、HHblitsはより速く、最大2倍の感度を持ち、より正確なアラインメントを生成します。** HHblitsは、HHsearchと同じHMM-HMMアラインメントアルゴリズムを使用しますが、低速HMM-HMM比較を実行するのに使われるデータベースHMMの数を数千万からわずか数千に減らす高速プレフィルタを使用しています。

## HH-Suiteはどんなときに役立つの？（Google翻訳）

BLAST、FASTA、またはPSI-BLASTなどの配列検索方法は、配列検索で同定された相同タンパク質または遺伝子からタンパク質または遺伝子の機能情報が推測できる理由で、生物学的研究にとって最も重要です。しかし、既知の機能を有するタンパク質との重要な関係性を確立することはできないということがしばしばあります。これはオーソログがまだ研究されていない、最も興味深いタンパク質のグループの場合にまさにその通りと言えるでしょう。従来の配列検索方法が失敗する場合、HHblitsおよびHHsearchは、たいていの場合において、より遠隔の相同的関係からの推論を可能にします。
HHblitsはより良好なMSAを構築し、HHsearchまたはHHblitsを使用して（PDBまたはPfamなどのドメインデータベースを検索したりして）より遠隔な相同体を見つけることができます。関係が非常に遠く離れていて共通の機能が仮定できない場合でも、可能性のあるメカニズム、活性部位の位置および残基、または基質結合のクラスについての仮説を引き出せることがしばしばあります。約30億年前に共通の祖先を共有していたタンパク質ドメインでさえほとんど同じような性質を持つため、既知の構造を持つ相同タンパク質が同定できれば、その構造は目的のタンパク質の3D構造をモデル化するためのテンプレートとして使用できます。3Dモデルは、実験を導くための仮説を生成するのに役立ちます。

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

  - macOS High Sierra / MojaveまたはCentOS 7.4, Ubuntu 18.04以降
  - **100GB以上の空き容量を持つハードディスク**

HH-Suiteを動作させるためには後述のHHsuite用データベースをインストールしておくだけの十分な空き容量を持つマシンが必要です。ちなみに**SSD上にそのデータベースを置き、十分大きなメモリがあるマシンならば、HH-suiteは快適に動作します**。

## インストール方法
インストール要件に、

  - CMake 2.8.12以降
  - [pdbx python library](https://github.com/soedinglab/pdbx)

と書かれており、後者をまずインストールする必要があります。

    #インストールディレクトリに移動。~/appsとか。
    mkdir -p ~/apps ; cd ~/apps
    # Gitで直接hh-suiteディレクトリのダウンロード
    git clone https://github.com/soedinglab/pdbx.git
    cd pdbx
    mkdir build ; cd build
    cmake .. -DCMAKE_CXX_COMPILER=g++ -DCMAKE_C_COMPILER=gcc (-DUserInstallOption=ON)
    make -j4 && make install

root権限がない場合は`-DUserInstallOption=ON`を入れておくと良いでしょう。root権限がある場合はシステム全体にpdbxライブラリがインストールされます。全体にインストールしても大きな害はないと思われるのでしておいても悪くないと思います。

次に、HHSuiteをGitから直接ダウンロードしてきます。現時点ではバージョンは3.3となっています。これを`/Users/hogehoge/apps/hh-suite/3.3`ディレクトリにインストールすることにします。コンパイラは`gcc`,`g++`を使用(GNU Compiler 9.3.0)。MPIコンパイラはこれに合わせてgcc9.3.0でビルドしたOpenMPIを指定したものを用意します（参考： [OpenMPIのインストール方法](https://qiita.com/Ag_smith/items/0e73d6d10e0b35c3f6d5#%E5%8F%82%E8%80%83openmpi%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB%E6%96%B9%E6%B3%95)）。**macOSの場合、gcc-5, 6系ではコンパイルが問題なく行くのですが、実行するとsegmentation faultを起こして失敗しました**。そのためgcc-9系でビルドしました。（Linux OS系ではgcc-5でも問題ないようです。）

```bash
#インストールディレクトリに移動。下の例ではホームディレクトリ直下の~/appsにします。
mkdir -p ~/apps
cd ~/apps
brew install wget
# Gitで直接hh-suiteディレクトリのダウンロード
wget https://github.com/soedinglab/hh-suite/archive/v3.3.0.tar.gz
# 解凍してディレクトリ内部に入る
tar zxvf v3.3.0.tar.gz ; cd hh-suite-3.3.0
# cmake用buildディレクトリの作成
mkdir build ; cd build
# cmakeを行う。${HOME}は自動的にホームディレクトリに設定されているはずの環境変数。
cmake .. -DCMAKE_INSTALL_PREFIX=${HOME}/apps/hh-suite/3.3 \
-DCMAKE_C_COMPILER=gcc-9 \
-DCMAKE_CXX_COMPILER=g++-9 \
-DMPI_CXX_COMPILER=${HOME}/apps/openmpi/4.0.5_gcc9.3.0/bin/mpicxx \
-DMPI_C_COMPILER=${HOME}/apps/openmpi/4.0.5_gcc9.3.0/bin/mpicc
# インストール開始
make -j4 && make install
```

公式ページによれば`-DCMAKE_BUILD_TYPE=RelWithDebInfo`を付けてインストールしているようですが、他にはDebug, Releaseなどが指定できます。ふつう安定版であればReleaseをつけることが多いのですが、今回はbeta版だからか、`RelWithDebInfo`らしいです。これをつけるとCMAKE\_C\_FLAGS\_RELWITHDEBINFOまたはCMAKE\_CXX\_FLAGS\_RELWITHDEBINFOがつくらしい（参考： https://gitlab.kitware.com/cmake/community/wikis/doc/cmake/Useful-Variables）

インストール後、使用準備のために`$HHLIB/scripts/HHPaths.pm`のディレクトリ中のファイルを設定します。ここで`$HHLIB`変数は上の`-DCMAKE_INSTALL_PREFIX`と同じものにしておきます。今回の例では`HHLIB=${HOME}/apps/hh-suite/3.3`としておきました。

    ##############################################################################################
    # PLEASE COMPLETE THE PATHS ... TO PSIPRED AND OLD-STYLE BLAST (NOT BLAST+) (NEEDED FOR PSIPRED)
    our $execdir = "/usr/local/package/psipred/3.5/bin";  # path to PSIPRED V2 binaries
    our $datadir = "/usr/local/package/psipred/3.5/data"; # path to PSIPRED V2 data files
    our $ncbidir = "/raid1/share/tools/blast/executables/blast/blast-2.2.26/bin";    # path to NCBI binaries (for PSIPRED in addss.pl)
    
    ##############################################################################################
    # PLEASE COMPLETE THE PATHS ... TO YOUR LOCAL PDB FILES, DSSP FILES ETC.
    our $pdbdir  =  "/raid1/share/database/PDB/data/structures/divided"; # where are the pdb files? (pdb/divided directory will also work)
    our $dsspdir =  "/raid1/share/database/DSSPrsync"; # where are the dssp files? Used in addss.pl
    our $dssp    =  "/bin/mkdssp";  # where is the dssp binary? Used in addss.pl
    ##############################################################################################

ここの記述をそれぞれ設定しておきます。注意点として、`$ncbidir`に設定すべきlegacy blastは**ver 2.2.26**をしておいたほうがいいです（2.2.23では`addss.pl`が動作しなかった）。blast-2.2.26は ftp://ftp.ncbi.nlm.nih.gov/blast/executables/legacy.NOTSUPPORTED/2.2.26/ で配布されていますのでゲットしておきましょう。`mkdssp`はCentOS7系では`yum install dssp`を使ってインストールすることが可能。`$pdbdir`と`$dsspdir`はそれぞれの公式からrsyncなどでとってきたものをおいておきます。

## PDBミラーデータベースの構築（？）
タンパク質構造情報のデータベースをサーバー内に作っておく必要があるらしいです（？）。本当に必要なんだっけ……？ まあタンパク構造を扱うバイオインフォ系ラボでしたら入れておいて損はないです。
やり方は[PDBjのウェブサイト](https://pdbj.org/info/archive)に書いてありますので、例えば

    # /raid1/share/database/PDB/data/structuresにPDBデータを置きたい場合の例
    mkdir -p /raid1/share/database/PDB/data/structures/divided
    cd /raid1/share/database/PDB/data/structures/divided
    rsync -avz --delete ftp.pdbj.org::ftp_data/structures/divided/pdb .

とすれば全PDBデータ（pdbフォーマット）をダウンロードすることができます（割と高速にダウンロードできます）。ただしこれは`pdb1000d.ent.gz`のようにgzip圧縮されていますので、利用するためには

    # オプション-kをつけるとgzファイルを残す。つけなければ残らない。
    gunzip -k pdb1000d.ent.gz
    mv pdb1000d.ent 1000d.pdb

としておく必要があるでしょう。
`rsync`コマンドでgzファイルを取得してきているので、元のent.gzファイルを残しておかないといつも全データをダウンロードしてくる羽目になります。更新分だけgunzipしてmvするくらいのうまい設定をしておいたほうがいいかもしれません。

## DSSPデータベースについて
PDB構造情報に対するDSSPの計算結果は定期的にデータベース上にアップロードされています。これは<https://swift.cmbi.umcn.nl/gv/dssp/DSSP_5.html>に書いてあるように、以下のコマンドで簡単に取ってくることができます。

    rsync -avz --delete rsync://rsync.cmbi.umcn.nl/dssp-from-mmcif/ /raid1/share/database/DSSPrsync/

これを定期的にやっておけばOK。

## HHsuite用データベースについて

さて、配列検索をするときにはHH-suite用の以下のデータベースが必要となります。欲しいのは`UniRef30_2020_06_hhsuite.tar.gz`です。最新版はhttp://wwwuser.gwdg.de/~compbiol/uniclust/2020_06/ からダウンロードできます。今は1年に1回ペースで新しいデータベースを公開してくれているようです。
いずれもファイルサイズが数十GBくらいある上に、ドイツにファイル本体があるようなので、ダウンロードには12時間くらいかかると思います。ゆっくり待ちましょう。

ダウンロードが完了して解凍されると中には`UniRef30_2020_06_hhm.ffdata`や`UniRef30_2020_06_a3m.ffdata`などのデータが入っています。まとめて`UniRef30_2020_06`という名前のディレクトリを作成して中に入れておきましょう。忘れがちなのですが、解凍したファイルはパーミッションが`600`になっているので、他のユーザーからも利用したい場合や共用機に置く場合は`chmod 644 *`をして利用可能にしておく必要があります。

## 使い方の例

### Generating a multiple sequence alignment using HHblits

インプットとして指定した`query.seq`に含まれる単独の配列または最初のMSAについてのMSAを生成するには、検索対象のデータベースが配列空間全体をカバーしている必要があります。ここで、データベースの指定`-d`に上で解凍した`UniRef30_2020_06`ディレクトリを指定しましょう。オプション`-oa3m`<msa_file>は、有意なヒットから出力MSAを生成するようにHHblitsに指示します。`-n 1`は、1回だけ検索のiterationをするように指定します。

```
$ hhblits -cpu 4 -i data/query.seq -d dbs/UniRef30_2020_06/UniRef30_2020_06 -oa3m query.a3m -n 1
```

> To generate an MSA for a sequence or initial MSA in query.a3m, the database to be searched should cover the entire sequence space, such as uniprot20 or nr20. The option `-oa3m `<msa_file> tells HHblits to generate an output MSA from the significant hits, -n 1 specifies a single search iteration.

検索の最後に、HHblitsは、MSAを含むpackedデータベースファイルから、閾値以下のE値を有するHMMに属するシーケンスを読み込みます。MSAに含めるE値のしきい値は、`-e `<E-value>オプションを使用して指定できます。検索後、`query.a3m`にはA3M形式のMSAが含まれます。

> At the end of the search, HHblits reads from the packed database file containing the MSAs the sequences belonging to HMMs with E-value below the threshold. The E-value threshold for inclusion into the MSA can be specified using the `-e `<E-value> option. After the search, query.a3m will contain the MSA in A3M format.

以前の検索で得られたMSAからスタートすることで、より多くの配列を追加するために2回目の反復検索を行うことも可能です。以前の検索後に生成されたMSAには`query.seq`に含まれる単一の配列よりも多くの情報が含まれているため、このMSAで検索すると、より多くの相同性のあるデータベースマッチングが得られる可能性があります。

> We could do a second search iteration, starting with the MSA from the previous search, to add more sequences. Since the MSA generated after the previous search contains more information than the single sequence in query.seq, searching with this MSA will probably result in many more homologous database
matches.

```
$ hhblits -cpu 4 -i query.a3m -d dbs/UniRef30_2020_06/UniRef30_2020_06 -oa3m query.a3m -n 1
```

以上のやり方は、以下のようにコマンドを変えれば、2回コマンドを打つことなく、最初の`query.seq`から直接的に反復検索を2回行うことができます。

> Instead, we could directly perform two search iterations starting from query.seq:

```
$ hhblits -cpu 4 -i data/query.seq -d dbs/UniRef30_2020_06/UniRef30_2020_06 -oa3m query.a3m -n 2
```

この2つの結果が微妙に異なることについては公式マニュアルのAbout unexpected HHsuite behavior and troubleshootingの項を読んでください。

実用上は、「信頼性と特異性 VS 遠隔同族配列の感度」間のトレードオフに応じて、**MSAを構築するのに1〜4回ほど反復検索をするのが推奨されます**。 より多くの検索反復が行われるほど、MSAに入る非相同配列または配列セグメントのリスクが高くなり、続くその後の反復でそういったものをより多く集めてくることになります。これは、短い反復を含む配列、アミノ酸組成バイアスを有する領域、およびそれほど劇的ではないが複数のドメインを含む配列を用いて検索する場合に特に問題になります。幸いにも、**この問題は、PSI-BLASTと比較すると、hhblitsの反復回数はより少なく、より堅牢なMaximum Accuracyアラインメントアルゴリズム、およびそのHMM-HMMアライメントのより高い精度のために、hhblitsではあまり顕著になることはありません。**