HMMERの使い方

## HMMERとは
HMMERは、タンパク質やDNA配列のドメインファミリーを確率論的にモデル化するツールを提供するソフトウェアパッケージです。この記事では、HMMERで使用できるプログラムの使い方を紹介します。
## 動作環境
- Windows11
- Ubuntu 22.04.3 LTS
## インストール方法
HMMERはanacondaでインストールできます。(他にも方法はあります)
anacondaの環境構築方法は[他の記事](https://qiita.com/116ryusei/items/adb893a79235bf84fd7b)でまとめています。
```bash:Ubuntu
$ conda install bioconda::hmmer
```
## HMMER プログラム一覧(抜粋)
HMMERで使用できるプログラムを抜粋しました。他のプログラムは、[ユーザーガイド](http://eddylab.org/software/hmmer/Userguide.pdf)から使用方法を確認できます。
**1. hmmstat**
**2. hmmpress**
**3. hmmscan**
**4. hmmfetch**
**5. hmmsearch**


これらプログラムの使用方法の紹介にあたり、HMM profiles, 2種のアミノ酸配列,　大腸菌のゲノム配列を用意しました(以下データセット)。 
## データセット
#### HMM profiles
HMM profilesは[InterPro(Pfam)](https://www.ebi.ac.uk/interpro/download/Pfam/)にあるものをダウンロードしました(Pfam-A.hmm.gz)。ファイルは圧縮されているので`$ gunzip Pfam-A.hmm.gz`で解凍する必要があります(Pfam-A.hmm)。
![interpro.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602355/19218d9b-a7e7-6e68-5a3b-6345134e7f12.png)

### アミノ酸配列
アミノ酸配列を二つ用意しました。これらのアミノ酸は、ドメインを一つまたは複数持っています。
```bash:sample.fa
>singledomain
MSKSDVFHLGLTKNDLQGATLAIVPGDPDRVEKIAALMDKPVKLASHREFTTWRAELDGKPVIVCSTGIGGPSTSIAVEELAQLGIRTFLRIGTTGAIQPHINVGDVLVTTASVRLDGASLHFAPLEFPAVADFECTTALVEAAKSIGATTHVGVTASSDTFYPGQERYDTYSGRVVRHFKGSMEEWQAMGVMNYEMESATLLTMCASQGLRAGMVAGVIVNRTQQEIPNAETMKQTESHAVKIVVEAARRLL
>multipledomain
MADSPLRHLSTFVRYLRGERTAVALLALLVPAGMGLQLVAPYLLRGFIDGALSGDSRKTLLDLAAWSLAAAVGTLVVTAGTEALSSRVAWRSTNRLRADLVEHCLSRPPGFYRKHPPGELVERMDGDVTRLAAVMSTLLLELLAQALLIVGILVALFRLEWRLALVFAPFAAGTLLLLRTLVGRAMPFVTARQRVAADLQGFLEERLAAAEDLRVNGASRYTLRELGDRQDDLYRKARDAARASVRWPATVQGLSAVSVVLALAVSAWLHARGQLSTGTAFASLSYAMLLRRPLLAVTTRFRELEDAAASAQRLRDLLGHGTAAPRTGRGTLPAGLPGVRFDGVSFGYEPDEPVLRDVSFTLRPGERLGVVGRTGSGKSTVVRLLFGLHHPGAGSVSAGGLDLTEIDPRALRSRVALVTQEVHVFYASLRDNLTFFDRSVPDDRLRAALGEAGLGPWLRTLPDGLDTPLGAGARGMSAGEEQQLALARVFLRDPGLVLMDEPTARLDPYSERLLMPALERLLEGRTAVVVEHRPHLLRNVDRILVLEEGKVAEEGERRVLAADPGSRFHALLRTAGATR
```
### ゲノム配列
ゲノム配列は、NCBIから取得しました。今回使用するのは、[Escherichia coli str. K-12 substr. MG1655](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000005845.2/)の配列です。以下スクショのように指定してダウンロードしました。
GCA_000005845.2のフォルダにあるprotein.faaというファイル名を、GCA_000005845.2.faaに変更して使用しました。

![refseq.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602355/8c0ae47d-5534-8f2b-18cd-ef245cfb73f0.png)





## 実施例
### 1. hmmstat
hmmstatは、今回使用するPfam-A.hmm(hmmfile)の要約統計量を出力します。これにより、HMM profilesに含まれているプロファイルが簡単に分かります。
```bash:Ubuntu
$ hmmstat Pfam-A.hmm > pfam_hmmstat.txt
$ head -n 15 pfam_hmmstat.txt
# hmmstat :: display summary statistics for a profile file
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#
# idx  name                 accession        nseq eff_nseq      M relent   info p relE compKL
# ---- -------------------- ------------ -------- -------- ------ ------ ------ ------ ------
1      1-cysPrx_C           PF10417.14         30    15.00     41   1.34   1.28   1.26   0.12
2      10_blade             PF21734.2          14     1.62    345   0.59   0.55   0.54   0.07
3      117-like_vir         PF21578.2           1     0.38    101   0.59   0.54   0.51   0.08
4      120_Rick_ant         PF12574.13          2     0.45    240   0.59   0.60   0.52   0.03
5      12TM_1               PF09847.14          6     2.00    448   0.59   0.62   0.51   0.12
```

### 2. hmmpress
hmmpressは標準的なhmmfileをバイナリ圧縮します。これは、hmmscanを実施する前に必要なステップです。出力として、hmmfile.h3m, hmmfile.h3m, hmmfile.h3f, hmmfile.h3pが得られます。実際にこのコマンドを使用すると以下の出力が得られます。

```bash:Ubuntu
$ hmmpress Pfam-A.hmm
Working...    done.
Pressed and indexed 21979 HMMs (21979 names and 21979 accessions).
Models pressed into binary file:   Pfam-A.hmm.h3m
SSI index for binary model file:   Pfam-A.hmm.h3i
Profiles (MSV part) pressed into:  Pfam-A.hmm.h3f
Profiles (remainder) pressed into: Pfam-A.hmm.h3p
```

### 3. hmmscan
hmmscanは、クエリー配列のファイル(今回はsample.fa)にある各配列をもとに、hmmdb(hmmpressしたhmmfile)を検索し、配列に最もマッチしたプロファイルのリストを出力します。以下出力結果から、単一ドメインの配列には1つのプロファイルがヒットし、それはPNP_UDP_1であることがわかりました。複数ドメインの配列には複数のプロファイルがヒットしました。
```bash:Ubuntu
$ hmmscan --domtblout sample_hmmscan.txt Pfam-A.hmm sample.fa
$ head sample_hmmscan.txt #出力確認
#                                                                            --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
# target name        accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
#------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
PNP_UDP_1            PF01048.25   233 singledomain         -            253     4e-40  138.0   1.5   1   1     2e-44   4.4e-40  137.8   1.5     1   223    20   244    20   253 0.85 Phosphorylase superfamily
ABC_membrane         PF00664.28   274 multipledomain       -            579   1.1e-29  104.3   0.0   1   1   5.3e-33   1.3e-29  104.0   0.0    13   274    35   294    23   294 0.96 ABC transporter transmembrane region
ABC_tran             PF00005.32   137 multipledomain       -            579   3.4e-27   96.0   0.0   1   1   2.3e-30   5.6e-27   95.3   0.0     1   137   355   504   355   504 0.95 ABC transporter
AAA_21               PF13304.11   304 multipledomain       -            579   5.1e-06   27.0   0.5   1   2    0.0037       8.9    6.5   0.0     3    22   369   388   368   414 0.78 AAA domain, putative AbiEii toxin, Type IV TA system
AAA_21               PF13304.11   304 multipledomain       -            579   5.1e-06   27.0   0.5   2   2   7.4e-07    0.0018   18.7   0.1   236   295   474   530   421   537 0.88 AAA domain, putative AbiEii toxin, Type IV TA system
DUF87                PF01935.22   220 multipledomain       -            579    0.0074   16.7   0.1   1   1   6.4e-06     0.016   15.6   0.1    25    44   367   386   361   388 0.92 Helicase HerA, central domain
SMC_N                PF02463.24   220 multipledomain       -            579     0.012   15.3   0.0   1   2      0.01        26    4.4   0.0    28    43   369   384   357   395 0.84 RecF/RecN/SMC N terminal domain
```



### 4. hmmfetch
hmmfetchは、複数のHMMプロファイルがまとまったhmmfileから、目的とするHMMプロファイルを抽出します。ここでは、単一ドメインの配列にヒットした、PNP_UDP_1のプロファイルを取り出します。
```bash:Ubuntu
$ hmmfetch -o pnp_udp_1.hmm Pfam-A.hmm PNP_UDP_1
$ head -n 5 pnp_udp_1.hmm #出力確認
HMMER3/f [3.1b2 | February 2015]
NAME  PNP_UDP_1
ACC   PF01048.25
DESC  Phosphorylase superfamily
LENG  233
$ tail -n 5 pnp_udp_1.hmm #出力確認
          0.02583  6.16670  3.75488  0.61958  0.77255  0.48576  0.95510
    233   1.95585  6.02370  2.96175  2.00713  5.37000  3.63762  3.97967  3.70400  2.10203  3.54825  5.06404  2.88035  3.47748  3.13431  2.06920  2.33703  2.85980  3.92648  6.45565  5.04395    491 a - - H
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00215  6.14302        *  0.61958  0.77255  0.00000        *
//

```
### 5. hmmsearch
hmmsearchは、配列データベースに対して、目的のプロファイルにマッチする配列があるかを検索することができます。hmmfetchで取り出したPNP_UDP_1のプロファイルにマッチする配列が、Escherichia coli str. K-12 substr. MG1655に存在するか調べました。その結果、複数の配列にこのドメインが含まれていることがわかりました。
```bash:Ubuntu
$ hmmsearch --domtblout GCA_000005845.2.txt pnp_udp_1.hmm GCA_000005845.2.faa
$ head -n 8 GCA_000005845.2.txt
#                                                                            --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
# target name        accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
#------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
AAC73270.1           -            232 PNP_UDP_1            PF01048.25   233   8.2e-57  190.2   2.2   1   1   1.3e-59   9.1e-57  190.1   2.2     1   232     2   226     2   227 0.96 5'-methylthioadenosine/S-adenosylhomocysteine nucleosidase [Escherichia coli str. K-12 substr. MG1655]
AAC75460.1           -            277 PNP_UDP_1            PF01048.25   233     7e-45  151.2   0.6   1   1   1.1e-47   7.9e-45  151.0   0.6     1   233    27   274    27   274 0.95 xanthosine phosphorylase [Escherichia coli str. K-12 substr. MG1655]
AAC77337.1           -            239 PNP_UDP_1            PF01048.25   233   4.8e-41  138.6   1.3   1   1   7.8e-44   5.6e-41  138.4   1.3     2   218    16   223    15   234 0.90 purine nucleoside phosphorylase [Escherichia coli str. K-12 substr. MG1655]
AAC76834.1           -            253 PNP_UDP_1            PF01048.25   233   7.7e-41  138.0   1.5   1   1   1.2e-43   8.7e-41  137.8   1.5     1   223    20   244    20   253 0.85 uridine phosphorylase [Escherichia coli str. K-12 substr. MG1655]
AAC75046.1           -            484 PNP_UDP_1            PF01048.25   233   3.5e-22   76.9   0.0   1   1   8.1e-25   5.8e-22   76.2   0.0    37   201   266   431   246   452 0.86 AMP nucleosidase [Escherichia coli str. K-12 substr. MG1655]
```

# まとめ
HMMERの使用方法を簡単にまとめました。以降の記事で、詳細なオプションの紹介や、出力ファイルの説明等していきたいです。