msa

自前でquery配列からMultiple Sequence Alignment(MSA)を作るサーバを作ってみた

MSA作成ツール†

複数存在して用途目的にあわせて選ぶ

MMseqs2 https://github.com/soedinglab/MMseqs2
ColabFoldが使用している高速クラスタリング・検索ツール
JackHMMER/HHblits
AlphaFold2でも使われるツール

データベースと絡めると

ツール	データベース	メモ	用途・特徴
MMseqs2	UniRef30		高速・省リソース(Boltz/ColabFold用途)
JackHMMER	UniRef90 + MGnify	UniRef90(ホモログ探索)向けとMGnify(追加探索)	伝統的・精度重視
HHblits	BFD / UniClust30	どちらか一方. alphafold2では２つとも検索対象	超大規模メタゲノムDBを検索(BFD[1.4TB])

alphafold2では JackHMMER でUniRef90とMGnifyを検索して、HHblitsでBFDとUniClust30を検索している. っで出てきたMSAを全てマージしている.
LocalColabFoldをローカルで実行する場合は、UniRef30、ColabFold_envdb(メタゲノム系配列)、PDB70を検索するみたい.

ここでは UniRef30 を対象に MMseqs2 で MSA を作る事にした

マルチマーな蛋白質に対しては素直に localcolabfold を導入して、その中にある colabfold_search コマンドを使った方がいいかな

↑

MMseqs2の導入†

ここでは ColabFold が採用している MMseqs2 を入れます.
https://github.com/soedinglab/MMseqs2/releases に既にコンパイル済みのプログラムが配布されています.
コンパイルせず、これをそのまま使いたいと思います.

最新の「MMseqs2 Release 18-8cc5c」で「mmseqs-linux-avx2.tar.gz」(16.6 MB)を使います.

[root@rockylinux9 ~]# cd /apps/src/
[root@rockylinux9 src]# wget https://github.com/soedinglab/MMseqs2/releases/download/18-8cc5c/mmseqs-linux-avx2.tar.gz
[root@rockylinux9 src]# cd ..
[root@rockylinux9 apps]# tar xf src/mmseqs-linux-avx2.tar.gz
 
(オプション)
[root@rockylinux9 apps]# cp mmseqs/util/bash-completion.sh /usr/share/bash-completion/completions/mmseqs

EnvironmentModulesは「/apps/modulefiles/mmseqs」

#%Module1.0
set root /apps/mmseqs
prepend-path PATH $root/bin

↑

データベース準備†

UniRef（UniProt Reference Clusters）を用意します.
アミノ酸配列が90%以上一致するタンパク質が１つのクラスターに纏められた UniRef90 か
配列が30%以上で似ているタンパク質が１つのクラスターに纏められた UniRef30 のどちらかが採用するみたい.

推奨は UniRef30
っが、uniprot.orgでの提供はない模様で ColabFold で提供されている UniRef30 を使用します

展開先は SSD なドライブが望ましいそうな. まぁー確かに.

[root@rockylinux9 ~]# mkdir -p /db/{uniref30,src} 
[root@rockylinux9 ~]# cd /db/src 
 
[root@rockylinux9 src]# wget https://wwwuser.gwdg.de/~compbiol/colabfold/uniref30_2302.tar.gz
 
[root@rockylinux9 src]# cd ../uniref30
 
[root@rockylinux9 uniref30]# tar xf /db/src/uniref30_2302.tar.gz
 
[root@rockylinux9 uniref30]# ls -l
total 224591556
-rw-------. 1 528745 9100  30961144274 May 16  2023 uniref30_2302_aln.tsv
-rw-------. 1 528745 9100   5797891705 May 22  2023 uniref30_2302_db_mapping
-rw-------. 1 528745 9100    667957493 May 22  2023 uniref30_2302_db_taxonomy
-rw-------. 1 528745 9100  46247602628 May 16  2023 uniref30_2302_h.tsv
-rw-------. 1 528745 9100          337 May 22  2023 uniref30_2302.md5sum
-rw-------. 1 528745 9100 137235400133 May 16  2023 uniref30_2302_seq.tsv
-rw-------. 1 528745 9100   9071701972 May 16  2023 uniref30_2302.tsv
 
[root@rockylinux9 uniref30]# du -hs .
215G    .
 
[root@rockylinux9 uniref30]#
 
(MMseqs2 用プロファイルDB変換)
[root@rockylinux9 uniref30]# /apps/mmseqs/bin/mmseqs  tsv2exprofiledb  uniref30_2302  uniref30_2302_db
 
(ケーキとお茶)
 
[timeを入れたら下記になった. 実時間で143minの待ち. 4coreなら37/4でまぁrealは10minくらい. それが143minなので14倍も遅い. まぁメモリー不足のswap発生でしょうね]
real    143m4.638s    <--  (各coreが動いていた総時間 user+sys) + (I/O待ち + スワップ待ち + 他プロセス待ち + 休止時間 など)
user    37m45.828s
sys     6m14.018s
 
[root@rockylinux9 uniref30]#
[root@rockylinux9 uniref30]# ls -l
total 381182412
-rw-------. 1 528745 9100  30961144274 May 16  2023 uniref30_2302_aln.tsv
-rw-r--r--. 1 root   root   5787495369 Nov 14 21:23 uniref30_2302_db
-rw-r--r--. 1 root   root   8709887243 Nov 14 21:40 uniref30_2302_db_aln
-rw-r--r--. 1 root   root            4 Nov 14 21:40 uniref30_2302_db_aln.dbtype
-rw-r--r--. 1 root   root    867494002 Nov 14 21:40 uniref30_2302_db_aln.index
-rw-r--r--. 1 root   root            4 Nov 14 21:24 uniref30_2302_db.dbtype
lrwxrwxrwx. 1 root   root           22 Nov 14 21:24 uniref30_2302_db_h -> uniref30_2302_db_seq_h
lrwxrwxrwx. 1 root   root           29 Nov 14 21:24 uniref30_2302_db_h.dbtype -> uniref30_2302_db_seq_h.dbtype
lrwxrwxrwx. 1 root   root           28 Nov 14 21:24 uniref30_2302_db_h.index -> uniref30_2302_db_seq_h.index
-rw-r--r--. 1 root   root    879290728 Nov 14 21:24 uniref30_2302_db.index
-rw-------. 1 528745 9100   5797891705 May 22  2023 uniref30_2302_db_mapping
-rw-r--r--. 1 root   root  83036144795 Nov 14 20:56 uniref30_2302_db_seq
-rw-r--r--. 1 root   root            4 Nov 14 21:00 uniref30_2302_db_seq.dbtype
-rw-r--r--. 1 root   root  43200163261 Nov 14 21:14 uniref30_2302_db_seq_h
-rw-r--r--. 1 root   root            4 Nov 14 21:17 uniref30_2302_db_seq_h.dbtype
-rw-r--r--. 1 root   root   8910693488 Nov 14 21:17 uniref30_2302_db_seq_h.index
-rw-r--r--. 1 root   root   8957791292 Nov 14 21:00 uniref30_2302_db_seq.index
-rw-------. 1 528745 9100    667957493 May 22  2023 uniref30_2302_db_taxonomy
-rw-------. 1 528745 9100  46247602628 May 16  2023 uniref30_2302_h.tsv
-rw-------. 1 528745 9100          337 May 22  2023 uniref30_2302.md5sum
-rw-------. 1 528745 9100 137235400133 May 16  2023 uniref30_2302_seq.tsv
-rw-------. 1 528745 9100   9071701972 May 16  2023 uniref30_2302.tsv
[root@rockylinux9 uniref30]#
 
 
 
(インデックス作成)
[root@rockylinux9 uniref30]# /apps/mmseqs/bin/mmseqs  createindex      uniref30_2302_db    tmp
 
(ケーキとお茶 2nd)
 
real    80m36.952s
user    59m5.088s
sys     2m13.538s
 
[root@rockylinux9 uniref30]#
 
[root@rockylinux9 uniref30]# ls -l
total 621506080
drwxr-xr-x. 3 root   root         4096 Nov 14 22:12 tmp
-rw-------. 1 528745 9100  30961144274 May 16  2023 uniref30_2302_aln.tsv
-rw-r--r--. 1 root   root   5787495369 Nov 14 21:23 uniref30_2302_db
-rw-r--r--. 1 root   root   8709887243 Nov 14 21:40 uniref30_2302_db_aln
-rw-r--r--. 1 root   root            4 Nov 14 21:40 uniref30_2302_db_aln.dbtype
-rw-r--r--. 1 root   root    867494002 Nov 14 21:40 uniref30_2302_db_aln.index
-rw-r--r--. 1 root   root            4 Nov 14 21:24 uniref30_2302_db.dbtype
lrwxrwxrwx. 1 root   root           22 Nov 14 21:24 uniref30_2302_db_h -> uniref30_2302_db_seq_h
lrwxrwxrwx. 1 root   root           29 Nov 14 21:24 uniref30_2302_db_h.dbtype -> uniref30_2302_db_seq_h.dbtype
lrwxrwxrwx. 1 root   root           28 Nov 14 21:24 uniref30_2302_db_h.index -> uniref30_2302_db_seq_h.index
-rw-r--r--. 1 root   root    388177920 Nov 14 23:32 uniref30_2302_db.idx.0
-rw-r--r--. 1 root   root 159321395200 Nov 14 23:32 uniref30_2302_db.idx.1
-rw-r--r--. 1 root   root  28806189056 Nov 14 23:32 uniref30_2302_db.idx.2
-rw-r--r--. 1 root   root  28789460992 Nov 14 23:32 uniref30_2302_db.idx.3
-rw-r--r--. 1 root   root  28786122752 Nov 14 23:33 uniref30_2302_db.idx.4
-rw-r--r--. 1 root   root            4 Nov 14 23:33 uniref30_2302_db.idx.dbtype
-rw-r--r--. 1 root   root          848 Nov 14 23:33 uniref30_2302_db.idx.index
-rw-r--r--. 1 root   root    879290728 Nov 14 21:24 uniref30_2302_db.index
-rw-------. 1 528745 9100   5797891705 May 22  2023 uniref30_2302_db_mapping
-rw-r--r--. 1 root   root  83036144795 Nov 14 20:56 uniref30_2302_db_seq
-rw-r--r--. 1 root   root            4 Nov 14 21:00 uniref30_2302_db_seq.dbtype
-rw-r--r--. 1 root   root  43200163261 Nov 14 21:14 uniref30_2302_db_seq_h
-rw-r--r--. 1 root   root            4 Nov 14 21:17 uniref30_2302_db_seq_h.dbtype
-rw-r--r--. 1 root   root   8910693488 Nov 14 21:17 uniref30_2302_db_seq_h.index
-rw-r--r--. 1 root   root   8957791292 Nov 14 21:00 uniref30_2302_db_seq.index
-rw-------. 1 528745 9100    667957493 May 22  2023 uniref30_2302_db_taxonomy
-rw-------. 1 528745 9100  46247602628 May 16  2023 uniref30_2302_h.tsv
-rw-------. 1 528745 9100          337 May 22  2023 uniref30_2302.md5sum
-rw-------. 1 528745 9100 137235400133 May 16  2023 uniref30_2302_seq.tsv
-rw-------. 1 528745 9100   9071701972 May 16  2023 uniref30_2302.tsv
[root@rockylinux9 uniref30]#
 
[root@rockylinux9 uniref30]# du -hs .
593G    .
 
[root@rockylinux9 uniref30]#

これでデータベースの準備が出来たので、試しに検索してみる

↑

検索実行†

検索ファイルはfasta形式でこんな感じ

[saber@rockylinux9 ~]$ mkdir msa
[saber@rockylinux9 ~]$ cd msa/
[saber@rockylinux9 msa]$
 
[saber@rockylinux9 msa]$ vi query.fasta
>sample
PIAQIHILEGRSDEQKETLIREVSEAISRSLDAPLTSVRVIITEMAKGHFGIGGELASK
 
[saber@rockylinux9 msa]$ 
[saber@rockylinux9 msa]$ ls -l
total 4
-rw-r--r--. 1 saber saber 68 Nov  5 03:27 query.fasta
 
[saber@rockylinux9 msa]$

これを検索してみます
「mmseq createdb」で検索ファイルのデータベース化を行って
「mmseq search」で uniref30 への検索を実施します.

[saber@rockylinux9 msa]$ module use /apps/modulefiles/
[saber@rockylinux9 msa]$ module load mmseqs
 
[saber@rockylinux9 msa]$ work=$(mktemp -d -p $(pwd))
[saber@rockylinux9 msa]$ mmseqs createdb  query.fasta  $work/query.qdb
createdb query.fasta query.qdb
 :
 :
[saber@rockylinux9 msa]$
[saber@rockylinux9 msa]$ ls -l $work
total 32
-rw-r--r--. 1 saber saber 61 Nov 15 00:31 query.qdb
-rw-r--r--. 1 saber saber  4 Nov 15 00:31 query.qdb.dbtype
-rw-r--r--. 1 saber saber  8 Nov 15 00:31 query.qdb_h
-rw-r--r--. 1 saber saber  4 Nov 15 00:31 query.qdb_h.dbtype
-rw-r--r--. 1 saber saber  6 Nov 15 00:31 query.qdb_h.index
-rw-r--r--. 1 saber saber  7 Nov 15 00:31 query.qdb.index
-rw-r--r--. 1 saber saber 11 Nov 15 00:31 query.qdb.lookup
-rw-r--r--. 1 saber saber 14 Nov 15 00:31 query.qdb.source
 
[saber@rockylinux9 msa]$ mmseqs   search   $work/query.qdb   /db/uniref30/uniref30_2302_db   $work/query.res    $work/tmp

検索が終わって、
得られた検索結果ファイルを「mmseqs result2msa」で a3m 形式へ変換します

[saber@rockylinux9 msa]$ ls -l $work
total 48
-rw-r--r--. 1 saber saber   61 Nov 15 00:31 query.qdb
-rw-r--r--. 1 saber saber    4 Nov 15 00:31 query.qdb.dbtype
-rw-r--r--. 1 saber saber    8 Nov 15 00:31 query.qdb_h
-rw-r--r--. 1 saber saber    4 Nov 15 00:31 query.qdb_h.dbtype
-rw-r--r--. 1 saber saber    6 Nov 15 00:31 query.qdb_h.index
-rw-r--r--. 1 saber saber    7 Nov 15 00:31 query.qdb.index
-rw-r--r--. 1 saber saber   11 Nov 15 00:31 query.qdb.lookup
-rw-r--r--. 1 saber saber   14 Nov 15 00:31 query.qdb.source
-rw-r--r--. 1 saber saber 2237 Nov 15 00:50 query.res          <--- 検索結果ファイル
-rw-r--r--. 1 saber saber    4 Nov 15 00:50 query.res.dbtype   <---
-rw-r--r--. 1 saber saber    9 Nov 15 00:50 query.res.index    <---
drwxr-xr-x. 3 saber saber   48 Nov 15 00:35 tmp
 
[saber@rockylinux9 msa]$
[saber@rockylinux9 msa]$ head $work/query.res
151993339       104     0.833   5.264E-25       0       58      59      0       58      129
76067200        86      0.745   1.157E-18       0       54      59      1       55      88
284902848       73      0.610   3.026E-14       0       58      59      1       59      63
21494430        72      0.605   5.716E-14       0       58      59      1       59      74
43193666        70      0.620   5.302E-13       0       54      59      23      77      87
306434389       66      0.598   1.279E-11       0       53      59      1       54      85
98039153        64      0.574   6.288E-11       0       54      59      1       55      86
208056876       63      0.526   1.635E-10       0       58      59      1       60      68
317866183       61      0.557   8.045E-10       0       53      59      1       54      59
43586680        61      0.519   8.045E-10       0       58      59      1       59      65
[saber@rockylinux9 msa]$ wc -l query.res
49 query.res
[saber@rockylinux9 msa]$
 
[saber@rockylinux9 msa]$ mmseqs result2msa  $work/query.qdb   /db/uniref30/uniref30_2302_db   $work/query.res   $work/query.a3m
 
[saber@rockylinux9 msa]$
[saber@rockylinux9 msa]$ ls -l $work
total 68
-rw-r--r--. 1 saber saber 9046 Nov 15 01:07 query.a3m          <--- 変換された a3m ファイル
-rw-r--r--. 1 saber saber    4 Nov 15 01:07 query.a3m.dbtype   <---
-rw-r--r--. 1 saber saber    9 Nov 15 01:07 query.a3m.index    <---
-rw-r--r--. 1 saber saber   61 Nov 15 00:31 query.qdb
-rw-r--r--. 1 saber saber    4 Nov 15 00:31 query.qdb.dbtype
-rw-r--r--. 1 saber saber    8 Nov 15 00:31 query.qdb_h
-rw-r--r--. 1 saber saber    4 Nov 15 00:31 query.qdb_h.dbtype
-rw-r--r--. 1 saber saber    6 Nov 15 00:31 query.qdb_h.index
-rw-r--r--. 1 saber saber    7 Nov 15 00:31 query.qdb.index
-rw-r--r--. 1 saber saber   11 Nov 15 00:31 query.qdb.lookup
-rw-r--r--. 1 saber saber   14 Nov 15 00:31 query.qdb.source
-rw-r--r--. 1 saber saber 2237 Nov 15 00:50 query.res
-rw-r--r--. 1 saber saber    4 Nov 15 00:50 query.res.dbtype
-rw-r--r--. 1 saber saber    9 Nov 15 00:50 query.res.index
drwxr-xr-x. 3 saber saber   48 Nov 15 00:35 tmp
[saber@rockylinux9 msa]$

↑

msaサーバを作る†

boltzで
「boltz predict --msa_server http://localhost:8000/msa input.fasta」
とされる「--msa_server」の宛先を作ります

dockerとかNginX/apache httpを用意せず、必要になったら上げて、使い終わったら落とす. そんな感じに.
あと基本的に前段のコマンドを順繰り流すだけなので、ある意味サーバに昇格させるほどもでないかなぁ.. 複数のデータベースを使ってマージしたmsaを作るとかならサーバにしたほうが便利かも. まあ大規模になりますね.

pythonのFastAPIを使います.

[root@rockylinux9 ~]# dnf install python3-pip

msaサーバを実行するユーザになって pip でパッケージを入れます. 一般ユーザでpipするとインストール先は$HOMEの .local 内になります.

[saber@rockylinux9 ~]$ pip install fastapi uvicorn python-multipart
 
[saber@rockylinux9 ~]$ ls -l .local/bin/
total 8
-rwxr-xr-x. 1 saber saber 210 Nov 15 21:39 fastapi
-rwxr-xr-x. 1 saber saber 211 Nov 15 21:39 uvicorn
[saber@rockylinux9 ~]$

っでmsaサーバのスクリプトですが、下記を使います。これを「msa_server.py」として保存して

from fastapi import FastAPI, UploadFile
import subprocess
import uuid
import os
 
app = FastAPI()
 
DB_PATH = "/db/uniref30/uniref30_2302_db"
TMP_DIR = "/tmp"
 
@app.post("/msa")
async def run_msa(file: UploadFile):
    # 一時ファイル
    job_id = str(uuid.uuid4())
    query_fasta = f"{TMP_DIR}/{job_id}.fasta"
    query_db = f"{TMP_DIR}/{job_id}_qdb"
    result_db = f"{TMP_DIR}/{job_id}_res"
    msa_out = f"{TMP_DIR}/{job_id}.a3m"
 
    # 保存
    with open(query_fasta, "wb") as f:
        f.write(await file.read())
 
    # MMseqs2 実行
    subprocess.run(["mmseqs", "createdb", query_fasta, query_db])
    subprocess.run(["mmseqs", "search", query_db, DB_PATH, result_db, TMP_DIR, "--threads", "16", "-s", "7.5"])
    subprocess.run(["mmseqs", "result2msa", query_db, DB_PATH, result_db, msa_out])
 
    # 結果返却
    with open(msa_out, "r") as f:
        msa = f.read()
 
    # cleanup は省略 (必要なら cron で)
    return {"msa": msa}

PATH環境変数に $HOME/.local/bin がある事を確認して、「msa_server.py」がある場所とかで「uvicorn msa_server:app --host 0.0.0.0 --port 8000」と実行します

[saber@rockylinux9 ~]$ echo $PATH
/home/saber/.local/bin:/home/saber/bin:/usr/share/Modules/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin
 
[saber@rockylinux9 ~]$ ls -l $HOME/msa_server.py
-rw-r--r--. 1 saber saber 960 Nov 15 21:44 /home/saber/msa_server.py
 
[saber@rockylinux9 ~]$ uvicorn msa_server:app --host 0.0.0.0 --port 8000
INFO:     Started server process [43478]
INFO:     Waiting for application startup.
INFO:     Application startup complete.
INFO:     Uvicorn running on http://0.0.0.0:8000 (Press CTRL+C to quit)
 :

msaサーバを停止させるには、CTRL+C (コントロールキーを押しながらCを押下) で終わります

っでこの段階でテストしてみます. 検索用ファイル「query.fasta」を用意して、下記を実行します

[saber@rockylinux9 ~]$ curl -X POST -F "file=@query.fasta" http://localhost:8000/msa

これで立ち上げた msa サーバ経由で検索してくれます

結果ですが、

[saber@rockylinux9 ~]$ curl -X POST -F "file=@query.fasta" http://localhost:8000/msa
{"msa":">sample\nPIAQIHILEGRSDEQKETLIREVSEAISRSLDAPLTSVRVIITEMAKGHFGIGGELASK\n>UniRef100_UPI0 ......AAHWAVGNQ----\n\u0000"}[saber@rockylinux9 msa]$

とjson形式で改行が"\n"で繋がった文字列を吐き出してくれます.

↑