{ "cells": [ { "cell_type": "markdown", "id": "c9467c38", "metadata": {}, "source": [ "# Fetch marker genes from a genome" ] }, { "cell_type": "markdown", "id": "094d9734", "metadata": {}, "source": [ "This example is adapted from the [fetchMGs](https://github.com/motu-tool/fetchMGs/blob/master/fetchMGs.pl) Perl script used in [mOTU](https://motu-tool.org/) to extract the 40 single-copy universal marker genes from a genome, annotating proteins to find the highest scoring proteins mapping to each of these marker genes.\n", "\n", "In this notebook, we show how to reproduce this kind of analysis, using `pyhmmer` instead of HMMER3 to perform the alignments and extract the bit scores.\n", "\n", "
\n", "\n", "References\n", " \n", "* [Ciccarelli FD, Doerks T, von Mering C, Creevey CJ, Snel B, Bork P. *Toward automatic reconstruction of a highly resolved tree of life.* Science. 2006 Mar 3;311(5765):1283-7. Erratum in: Science. 2006 May 5;312(5774):697.](https://pubmed.ncbi.nlm.nih.gov/16513982/)\n", "* [Sorek R, Zhu Y, Creevey CJ, Francino MP, Bork P, Rubin EM. *Genome-wide experimental determination of barriers to horizontal gene transfer.* Science. 2007 Nov 30;318(5855):1449-52.](https://pubmed.ncbi.nlm.nih.gov/17947550/)\n", " \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "f21f9ece", "metadata": {}, "outputs": [], "source": [ "import pyhmmer\n", "pyhmmer.__version__" ] }, { "cell_type": "markdown", "id": "37db1b70", "metadata": {}, "source": [ "## Getting the cutoffs\n", "\n", "Each HMM has been calibrated and contains custom cutoffs, but they are not in Pfam format, so we need to use them externaly. Let's start by downloading the file with these cutoffs from the GitHub repository of `fetchMG`:" ] }, { "cell_type": "code", "execution_count": null, "id": "77c04216", "metadata": {}, "outputs": [], "source": [ "import csv\n", "import io\n", "import urllib.request\n", "\n", "url = \"https://github.com/motu-tool/fetchMGs/raw/master/lib/MG_BitScoreCutoffs.allhits.txt\"\n", "\n", "cutoffs = {}\n", "with urllib.request.urlopen(url) as f:\n", " for line in csv.reader(io.TextIOWrapper(f), dialect=\"excel-tab\"):\n", " if not line[0].startswith(\"#\"):\n", " cutoffs[line[0]] = float(line[1])" ] }, { "cell_type": "markdown", "id": "cc40f524", "metadata": {}, "source": [ "## Downloading the HMMs\n", "\n", "Since the HMMs for the universal marker genes are also hosted on the `fetchMG` GitHub repository, we can download them from there too. `pyhmmer.plan7.HMMFile` supports reading from a file-handle, so we can parse each HMM as we download it. We also use the occasion to update the bitscore cutoffs specific to each HMM that we just obtained in the previous step.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cb1219de", "metadata": {}, "outputs": [], "source": [ "import urllib.request\n", "import pyhmmer.plan7\n", "\n", "baseurl = \"https://github.com/motu-tool/fetchMGs/raw/master/lib/{}.hmm\"\n", "\n", "hmms = []\n", "for cog in cutoffs:\n", " with urllib.request.urlopen(baseurl.format(cog)) as f:\n", " with pyhmmer.plan7.HMMFile(f) as hmm_file:\n", " hmm = hmm_file.read()\n", " cutoff = cutoffs[hmm.name.decode()]\n", " hmm.cutoffs.trusted = (cutoff, cutoff)\n", " hmms.append(hmm)" ] }, { "cell_type": "markdown", "id": "d4eff896", "metadata": {}, "source": [ "## Loading the sequences\n", "\n", "Now we need protein sequences to annotate. Let's use our set of protein sequences identified in the chromosome of [Anaerococcus provencensis](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=938293). " ] }, { "cell_type": "code", "execution_count": null, "id": "d0e20139", "metadata": {}, "outputs": [], "source": [ "import pyhmmer.easel\n", "with pyhmmer.easel.SequenceFile(\"data/seqs/938293.PRJEB85.HG003687.faa\", digital=True) as seqs_file:\n", " proteins = seqs_file.read_block()" ] }, { "cell_type": "markdown", "id": "f1d021ae", "metadata": {}, "source": [ "## Running the search pipeline\n", "\n", "With the proteins loaded, let's make a `namedtuple` that will contain the data we need to extract for each `hmmsearch` hit: the name of the query gene, the name of the marker gene which produced a hit, and the bitscore for the alignment." ] }, { "cell_type": "code", "execution_count": null, "id": "a5d53e11", "metadata": {}, "outputs": [], "source": [ "import collections\n", "Result = collections.namedtuple(\"Result\", [\"query\", \"cog\", \"bitscore\"])" ] }, { "cell_type": "markdown", "id": "f5baee6c", "metadata": {}, "source": [ "Now we can run the search pipeline: we annotate all the proteins will all the HMMs using the `pyhmmer.hmmsearch` function. Each HMM gives us a `TopHits` instance to process. Note that we run the search pipeline using the *trusted* bitscore cutoffs so that the cutoffs that we manually set on each HMMs are used to filter the results. We only keep hits that fall under the inclusion threshold for each HMM." ] }, { "cell_type": "code", "execution_count": null, "id": "00ae2763", "metadata": {}, "outputs": [], "source": [ "results = []\n", "for hits in pyhmmer.hmmsearch(hmms, proteins, bit_cutoffs=\"trusted\"):\n", " cog = hits.query_name.decode()\n", " for hit in hits:\n", " if hit.included:\n", " results.append(Result(hit.name.decode(), cog, hit.score))" ] }, { "cell_type": "markdown", "id": "814642fa", "metadata": {}, "source": [ "## Filtering results\n", "\n", "Now that we have all the hits that pass the bitscore thresholds, we can create a dictionary that maps each query protein to its highest scoring bitscore alignment, like in the original script. If a protein has alignments to two different marker genes with the same score, that query is ignored." ] }, { "cell_type": "code", "execution_count": null, "id": "d03a1ab3", "metadata": {}, "outputs": [], "source": [ "best_results = {}\n", "keep_query = set()\n", "for result in results:\n", " if result.query in best_results:\n", " previous_bitscore = best_results[result.query].bitscore\n", " if result.bitscore > previous_bitscore:\n", " best_results[result.query] = result\n", " keep_query.add(result.query)\n", " elif result.bitscore == previous_bitscore:\n", " if best_results[result.query].cog != hit.cog:\n", " keep_query.remove(result.query)\n", " else:\n", " best_results[result.query] = result\n", " keep_query.add(result.query)" ] }, { "cell_type": "markdown", "id": "ee538a94", "metadata": {}, "source": [ "Now we can get our final filtered results:" ] }, { "cell_type": "code", "execution_count": null, "id": "85484ed1", "metadata": {}, "outputs": [], "source": [ "filtered_results = [best_results[k] for k in sorted(best_results) if k in keep_query]" ] }, { "cell_type": "markdown", "id": "bd8701f5", "metadata": {}, "source": [ "We can print them to see which gene maps to which marker gene, with the score for each alignment. Let's look at the top 10 results:" ] }, { "cell_type": "code", "execution_count": null, "id": "b4d59cb9", "metadata": {}, "outputs": [], "source": [ "for result in filtered_results[:10]:\n", " print(result.query, \"{:.1f}\".format(result.bitscore), result.cog, sep=\"\\t\")" ] }, { "cell_type": "markdown", "id": "74691d9c", "metadata": {}, "source": [ "All set! You can now use this list of identifiers to filter the initial protein array, to write proteins grouped by marker genes in a FASTA file, compute alignment against the [mOTUs database](https://zenodo.org/record/3364101), or just save it for future use." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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", "version": "3.10.8" } }, "nbformat": 4, "nbformat_minor": 5 }