msa_to_hmm.ipynb 5.35 KB
Newer Older
1
2
3
4
5
6
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
7
    "# Build an HMM from a multiple sequence alignment"
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyhmmer\n",
    "pyhmmer.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphabet = pyhmmer.easel.Alphabet.amino()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A new HMM can be built from a single sequence, or from a multiple sequence alignment. Let's load an alignment in digital mode so that we can build our HMM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
49
    "with pyhmmer.easel.MSAFile(\"data/msa/LuxC.sto\", digital=True, alphabet=alphabet) as msa_file:\n",
50
51
52
53
54
55
56
57
58
59
60
    "    msa = next(msa_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Note\n",
    "    \n",
61
    "In this example, we load a multiple sequence alignment from a file, but if your program produces alignment and you wish to make an HMM out of them, you can instantiate a `TextMSA` object yourself, e.g.:\n",
62
63
    "    \n",
    "```python\n",
64
65
66
    "seq1 = pyhmmer.easel.TextSequence(name=b\"seq1\", sequence=\"WVPKQDFT\")\n",
    "seq2 = pyhmmer.easel.TextSequence(name=b\"seq2\", sequence=\"WL--PQGE\")\n",
    "msa  = pyhmmer.easel.TextMSA(name=b\"msa\", sequences=[seq1, seq2])\n",
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    "```\n",
    "\n",
    "Because we need a `DigitalMSA` to build the HMM, you will have to convert it first:\n",
    "    \n",
    "```python\n",
    "msa_d = msa.digitize(alphabet) \n",
    "```\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building an HMM\n",
    "\n",
    "Now that we have a multiple alignment loaded in memory, we can build a pHMM using a `pyhmmer.plan7.Builder`. This also requires a Plan7 background model to compute the transition probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "builder = pyhmmer.plan7.Builder(alphabet)\n",
    "background = pyhmmer.plan7.Background(alphabet)\n",
    "hmm, _, _ = builder.build_msa(msa, background)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a look at the consensus sequence of the HMM with the `consensus` property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hmm.consensus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving the resulting HMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
124
    "Now that we have an HMM, we can save it to a file to avoid having to rebuild it every time. Using the `HMM.write` method lets us write the HMM in ASCII or binary format to an arbitrary file. The resulting file will also be compatible with the `hmmsearch` binary if you wish to use that one instead of pyHMMER."
125
126
127
128
129
130
131
132
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
133
    "with open(\"LuxC.hmm\", \"wb\") as output_file:\n",
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    "    hmm.write(output_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying the HMM to a sequence database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once a pHMM has been obtained, it can be applied to a sequence database with the `pyhmmer.plan7.Pipeline` object. Let's iterate over the protein sequences in a FASTA to see if our new HMM gets any hits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
157
    "with pyhmmer.easel.SequenceFile(\"data/seqs/LuxC.faa\", digital=True, alphabet=alphabet) as seq_file:\n",
158
159
160
161
    "    sequences = list(seq_file)\n",
    "\n",
    "pipeline = pyhmmer.plan7.Pipeline(alphabet, background=background)\n",
    "hits = pipeline.search_hmm(query=hmm, sequences=sequences)"
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then query the `TopHits` object to access the domain hits in the sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ali = hits[0].domains[0].alignment\n",
    "\n",
    "print(\" \"*3, ali.target_name.decode())\n",
    "print(\"{:3}\".format(ali.hmm_from), ali.hmm_sequence[:80] + \"...\")\n",
    "print(\" \"*3, ali.identity_sequence[:80] + \"...\")\n",
    "print(\"{:3}\".format(ali.target_from), ali.target_sequence[:80] + \"...\")\n",
    "print(\" \"*3, ali.hmm_name.decode())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
189
   "display_name": "Python 3 (ipykernel)",
190
191
192
193
194
195
196
197
198
199
200
201
202
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
203
   "version": "3.10.4"
204
205
206
207
208
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}