![[PukiWiki]](image/picc.png "[PukiWiki]")
#author("2025-11-15T14:59:01+00:00","default:sysosa","sysosa") #author("2025-11-18T02:40:50+00:00","default:sysosa","sysosa") 自前でquery配列からMultiple Sequence Alignment(MSA)を作るサーバを作ってみた ***MSA作成ツール [#x51e265e] 複数存在して用途目的にあわせて選ぶ -[[MMseqs2]] https://github.com/soedinglab/MMseqs2 ColabFoldが使用している高速クラスタリング・検索ツール -JackHMMER/HHblits AlphaFold2でも使われるツール データベースと絡めると |BGCOLOR(YELLOW):ツール|BGCOLOR(YELLOW):データベース|BGCOLOR(YELLOW):メモ|BGCOLOR(YELLOW):用途・特徴| |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の導入 [#b55d2635] ここでは ColabFold が採用している [[MMseqs2]] を入れます. [[https://github.com/soedinglab/MMseqs2/releases>+https://github.com/soedinglab/MMseqs2/releases]] に既にコンパイル済みのプログラムが配布されています. コンパイルせず、これをそのまま使いたいと思います. 最新の「MMseqs2 Release 18-8cc5c」で 「mmseqs-linux-avx2.tar.gz」(16.6 MB)を使います. #code(nonumber){{ [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」 #code(nonumber){{ #%Module1.0 set root /apps/mmseqs prepend-path PATH $root/bin }} ***データベース準備 [#j76843e9] UniRef（UniProt Reference Clusters）を用意します. アミノ酸配列が90%以上一致するタンパク質が１つのクラスターに纏められた UniRef90 か 配列が30%以上で似ているタンパク質が１つのクラスターに纏められた UniRef30 のどちらかが採用するみたい. 推奨は UniRef30 っが、uniprot.orgでの提供はない模様で ColabFold で提供されている UniRef30 を使用します 展開先は SSD なドライブが望ましいそうな. まぁー確かに. #code(nonumber){{ [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]# }} これで データベース の準備が出来たので、試しに検索してみる ***検索実行 [#r2fa1a4e] 検索ファイルはfasta形式でこんな感じ #code(nonumber){{ [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 への検索を実施します. #code(nonumber){{ [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 形式へ変換します #code(nonumber){{ [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サーバを作る [#u3c20764] boltzで 「boltz predict --msa_server http://localhost:8000/msa input.fasta」 とされる「--msa_server」の宛先を作ります dockerとかNginX/apache httpを用意せず、必要になったら上げて、使い終わったら落とす. そんな感じに. あと基本的に前段のコマンドを順繰り流すだけなので、ある意味 サーバ に昇格させるほどもでないかなぁ.. 複数のデータベースを使ってマージしたmsaを作るとかなら サーバ にしたほうが便利かも. まあ大規模になりますね. pythonのFastAPIを使います. #code(nonumber){{ [root@rockylinux9 ~]# dnf install python3-pip }} msaサーバを実行するユーザになって pip でパッケージを入れます. 一般ユーザでpipするとインストール先は$HOMEの .local 内になります. #code(nonumber){{ [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」として保存して #code(python,nonumber){{ 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」と実行します #code(nonumber){{ [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」を用意して、下記を実行します #code(nonumber){{ [saber@rockylinux9 ~]$ curl -X POST -F "file=@query.fasta" http://localhost:8000/msa }} これで立ち上げた msa サーバ経由で検索してくれます 結果ですが、 #code(nonumber){{ [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"で繋がった文字列を吐き出してくれます. ***めも [#v9859a9a] こちらはあくまでも１配列から msa を用意しているだけです。こちらのツールは複数配列には対応してないです。 それをするなら colabfold(+配列データベースも用意) で colab_search コマンドで検索ですかね.
