notung_pipeline_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "48935755-9202-47dd-8e99-53a3bbab2710",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio import SearchIO, Entrez, SeqIO\n",
    "from ete3 import NCBITaxa, Tree\n",
    "# upload database when necessary\n",
    "#ncbi.update_taxonomy_database()\n",
    "\n",
    "ncbi = NCBITaxa()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fb62139d-e752-4336-8882-439d4dbf6e47",
   "metadata": {},
   "outputs": [],
   "source": [
    "domtblout_file = \"data/OG_343.domtblout\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a467b9ad-a795-4c3a-bd63-db86d7c8e829",
   "metadata": {},
   "source": [
    "### Parse hmmsearch results and taxonomy\n",
    "Parse the file. Get the taxonomy from the last part of the description, ie. `[Siphoviridae sp.]`. Filter out hits that are not Bacteria or Virus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "034044ac-ce8c-4066-8c3a-5b75f229e94c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['DAW89363.1', 'DAW86324.1', 'DAW94066.1', 'DAU94425.1', 'DAF01075.1', 'DAR20966.1', 'DAW20005.1', 'DAP51448.1', 'DAR70384.1', 'DAX21802.1', 'DAR45680.1', 'DAQ45717.1', 'DAW07241.1', 'DAL80299.1', 'DAE28469.1', 'DAX32588.1', 'DAP34063.1', 'DAJ75044.1', 'DAH06617.1', 'DAP92696.1', 'DAJ18143.1', 'DAO43059.1', 'DAH51560.1', 'DAX07875.1', 'WP_153163798.1', 'MBS6796132.1', 'DAV68353.1', 'WP_153170287.1', 'DAL58409.1', 'DAX15220.1', 'DAO36476.1', 'MBR5795345.1', 'MBQ7819940.1', 'MBS7122782.1', 'DAL29697.1', 'DAQ34250.1', 'CCX55664.1', 'DAJ12619.1', 'DAF73599.1', 'NBK99050.1', 'DAH77719.1', 'DAS74269.1', 'DAL04585.1', 'DAS78946.1', 'NBK97982.1']\n"
     ]
    }
   ],
   "source": [
    "target_proteins = list()\n",
    "target_taxids   = list()\n",
    "target_proteins_taxids = dict()\n",
    "target_proteins_lineages = list()\n",
    "\n",
    "records = SearchIO.parse(domtblout_file, \"hmmsearch3-domtab\")\n",
    "for record in records:\n",
    "    og_len = record.seq_len\n",
    "    # iterate the protein hits\n",
    "    for hit in record.hits:\n",
    "        hit_len = hit.seq_len\n",
    "        if hit.evalue < 0.01:\n",
    "            # set up the lists for the coverage calculation\n",
    "            hit_starts = list()\n",
    "            hit_ends   = list()\n",
    "            og_starts = list()\n",
    "            og_ends   = list()\n",
    "            for hsp in hit.hsps:\n",
    "                if hsp.evalue < 0.01:\n",
    "                    hit_starts.append(hsp.env_start) \n",
    "                    hit_ends.append(hsp.env_end)\n",
    "                    og_starts.append(hsp.query_start)\n",
    "                    og_ends.append(hsp.query_end)\n",
    "            #print(hit.id, hit_starts, hit_ends, og_starts, og_ends)\n",
    "                \n",
    "            # check there were hits passing the i-value thrshold\n",
    "            if hit_starts:\n",
    "                # calculate coverage of hit (protein seq) and query (og)\n",
    "                intervals_hit = [[s,e] for s,e in zip(hit_starts, hit_ends)]\n",
    "                intervals_og = [[s,e] for s,e in zip(og_starts, og_ends)]\n",
    "               \n",
    "                intervals_hit.sort(key=lambda interval: interval[0])\n",
    "                intervals_og.sort(key=lambda interval: interval[0])\n",
    "\n",
    "                merged_hit = [intervals_hit[0]]\n",
    "                for current in intervals_hit:\n",
    "                    previous = merged_hit[-1]\n",
    "                    if current[0] <= previous[1]:\n",
    "                        previous[1] = max(previous[1], current[1])\n",
    "                    else:\n",
    "                        merged_hit.append(current)\n",
    "\n",
    "                covered_length_hit = float(0)\n",
    "                for interval in merged_hit:\n",
    "                    covered_length_hit += (interval[1] - interval[0])\n",
    "\n",
    "                coverage_hit = float(covered_length_hit/hit_len)\n",
    "\n",
    "\n",
    "                merged_og = [intervals_og[0]]\n",
    "                for current in intervals_og:\n",
    "                    previous = merged_og[-1]\n",
    "                    if current[0] <= previous[1]:\n",
    "                        previous[1] = max(previous[1], current[1])\n",
    "                    else:\n",
    "                        merged_og.append(current)\n",
    "\n",
    "                covered_length_og = float(0)\n",
    "                for interval in merged_og:\n",
    "                    covered_length_og += (interval[1] - interval[0])\n",
    "\n",
    "                coverage_og = float(covered_length_og/og_len)\n",
    "\n",
    "                    \n",
    "                # check that coverages pass the cutoff\n",
    "                if coverage_og >= 0.5 and coverage_hit >= 0.5:\n",
    "                    # check which is the taxonomy\n",
    "                    taxa_description = hit.description.split(\" [\")[-1][:-1]\n",
    "                    # line below returns a dict, k=taxa_description  v=[tax_id]\n",
    "                    tax_id_dict = ncbi.get_name_translator([taxa_description])\n",
    "                    tax_id = tax_id_dict.get(taxa_description)[0]\n",
    "                    # get the lineage of the tax id\n",
    "                    lineage = ncbi.get_lineage(tax_id)\n",
    "                    # check if it is bacteria or virus\n",
    "                    if 2 in lineage or 10239 in lineage:\n",
    "                        # EXCLUDE HITS FROM CRASS-LIKE PHAGES\n",
    "                        if 1978007 not in lineage:\n",
    "                            target_proteins.append(hit.id)\n",
    "                            target_taxids.append(tax_id)\n",
    "                            target_proteins_taxids[hit.id] = tax_id\n",
    "                            target_proteins_lineages += lineage\n",
    "\n",
    "print(target_proteins)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10b45f5f-1872-414a-965b-bfa44ccca925",
   "metadata": {},
   "source": [
    "### Download proteins from NCBI\n",
    "First download the GenBank records, then parse the files to fasta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8c36dcbc-9b92-40f9-a9da-5d688a73d1c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download hits in GenBank format\n",
    "target_proteins = ['DAW89363.1', 'DAW86324.1', 'DAW94066.1', 'DAU94425.1', 'DAF01075.1', 'DAR20966.1', 'DAW20005.1', 'DAP51448.1', 'DAR70384.1', 'DAX21802.1', 'DAR45680.1', 'DAQ45717.1', 'DAW07241.1', 'DAL80299.1', 'DAE28469.1', 'DAX32588.1', 'DAP34063.1', 'DAJ75044.1', 'DAH06617.1', 'DAP92696.1', 'DAJ18143.1', 'DAO43059.1', 'DAH51560.1', 'DAX07875.1', 'WP_153163798.1', 'MBS6796132.1', 'DAV68353.1', 'WP_153170287.1', 'DAL58409.1', 'DAX15220.1', 'DAO36476.1', 'MBR5795345.1', 'MBQ7819940.1', 'MBS7122782.1', 'DAL29697.1', 'DAQ34250.1', 'CCX55664.1', 'DAJ12619.1', 'DAF73599.1', 'NBK99050.1', 'DAH77719.1', 'DAS74269.1', 'DAL04585.1', 'DAS78946.1', 'NBK97982.1']\n",
    "Entrez.email   = 'dcarrillo.bioinf@gmail.com'\n",
    "Entrez.api_key = 'b6db7fece605d37fcabd4b93749d2e46aa09'\n",
    "\n",
    "with open(\"borrar.gb\", \"w\") as fout:\n",
    "    handle = Entrez.efetch(db=\"protein\", id=target_proteins, rettype=\"gb\", retmode=\"text\")\n",
    "    batch_data = handle.read()\n",
    "    fout.write(batch_data)\n",
    "    handle.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1345b966-45d0-4b61-bead-3472fe4693f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# parse the GenBank file and get FASTA sequences\n",
    "\n",
    "to_faa = list()\n",
    "\n",
    "records = SeqIO.parse(\"data/OG_343_NCBI.gb\", \"gb\")\n",
    "\n",
    "for record in records:\n",
    "    # put the id in this form: <TAXID>|<SEQID>\n",
    "    record.id = f\"{target_proteins_taxids[record.id]}|{record.id}\"\n",
    "    record.description = \"\"\n",
    "    to_faa.append(record)\n",
    "    \n",
    "with open(\"data/OG_343_NCBI.faa\", \"w\") as fout:\n",
    "    SeqIO.write(to_faa, fout, \"fasta\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f5920a4-86d7-47d2-b8f8-0df5d531610a",
   "metadata": {},
   "source": [
    "### Merge OG and NCBI, align, trim and fasttree\n",
    "\n",
    "```\n",
    "(mafft_env) danielc@encode:~/projects/Bas_phages/crassphage_proteins/notung_pipeline/data:cat ../../../3_Broccoli/4_easy_OGs/0_easy_OGs_faa/OG_343.faa OG_343_NCBI.faa > gene_tree.faa\n",
    "\n",
    "(mafft_env) danielc@encode:~/projects/Bas_phages/crassphage_proteins/notung_pipeline/data:mafft-einsi --thread 10 gene_tree.faa > gene_tree.mafft\n",
    "\n",
    "(mafft_env) danielc@encode:~/projects/Bas_phages/crassphage_proteins/notung_pipeline/data:trimal -in gene_tree.mafft -out gene_tree_trim.mafft -gt 0.9\n",
    "\n",
    "(phylogenies) danielc@encode:~/projects/Bas_phages/crassphage_proteins/notung_pipeline/data:fasttree gene_tree_trim.mafft > gene_tree_trim.nwk\n",
    "```\n",
    "\n",
    "Now parse the tree and change leaves names to only the genome_id for the reference crassphages, or the taxid for the NCBI hits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "4deff857-203f-4ca6-aa44-3e47aad0b612",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read TerL tree\n",
    "gene_tree = Tree(\"data/gene_tree_trim.nwk\", format=1)\n",
    "# changes leaves' names\n",
    "for leaf in gene_tree.iter_leaves():\n",
    "    # check if the leaf comes from the reference set\n",
    "    label = leaf.name.split(\"|\")[0]\n",
    "    leaf.name = label\n",
    "# write the final tree\n",
    "gene_tree.write(format=1, outfile=\"data/gene_tree_final.nwk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf755b98-00b7-4e6e-8268-6ca3ad3a4a63",
   "metadata": {},
   "source": [
    "### Create species tree\n",
    "I have stored in a list all the tax_ids for the proteins that passed the cutoffs. Create a tree with them and attach it to the crAssphage TerL tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ca1f42bb-60c1-4d5d-9ece-129e7acf0ccf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                                         /-Bacteriophage sp.\n",
      "           /unclassified bacterial viruses\n",
      "          |                              \\-virus sp. ct9pU4\n",
      "          |\n",
      "          |                                    /-Podoviridae sp. ctY3D12\n",
      "    /Viruses                                  |\n",
      "   |      |            /unclassified Podoviridae-Podoviridae sp.\n",
      "   |      |           |                       |\n",
      "   |      |           |                        \\-crAss-like viruses\n",
      "   |       \\Caudovirales\n",
      "   |                  |--Myoviridae sp.\n",
      "   |                  |\n",
      "   |                   \\-Siphoviridae sp.\n",
      "-root\n",
      "   |                                                    /-Erysipelotrichaceae bacterium\n",
      "   |                                  /Erysipelotrichales\n",
      "   |                  /Erysipelotrichia                 \\-Coprobacillus sp.\n",
      "   |                 |               |\n",
      "   |        /Firmicutes               \\-Erysipelotrichia bacterium\n",
      "   |       |         |\n",
      "   |       |          \\-Veillonella sp. CAG:933\n",
      "    \\Bacteria\n",
      "           |--Thauera sp. 2A1\n",
      "           |\n",
      "            \\-Bacteroidales bacterium\n",
      "\n",
      "         /-38018\n",
      "      /-|\n",
      "     |   \\-2828248\n",
      "     |\n",
      "     |      /-2828134\n",
      "   /-|     |\n",
      "  |  |   /-|--2202567\n",
      "  |  |  |  |\n",
      "  |  |  |   \\-1978007\n",
      "  |   \\-|\n",
      "  |     |--2202564\n",
      "  |     |\n",
      "  |      \\-2170413\n",
      "--|\n",
      "  |            /-2049044\n",
      "  |         /-|\n",
      "  |      /-|   \\-2137881\n",
      "  |     |  |\n",
      "  |   /-|   \\-2184014\n",
      "  |  |  |\n",
      "  |  |   \\-1262980\n",
      "   \\-|\n",
      "     |--2570191\n",
      "     |\n",
      "      \\-2030927\n"
     ]
    }
   ],
   "source": [
    "# make tax_ids unique\n",
    "target_taxids = list(set(target_taxids))\n",
    "\n",
    "# add crAss-like taxid \n",
    "target_taxids.append(1978007)\n",
    "\n",
    "# get the ncbi tree for these tax_ids\n",
    "ncbi_topology = ncbi.get_topology(target_taxids)\n",
    "print(ncbi_topology.get_ascii(attributes=[\"sci_name\"]))\n",
    "# root the tree manually\n",
    "ncbi_topology_nwk = ncbi_topology.write()\n",
    "ncbi_topology_nwk = ncbi_topology_nwk[:-1] + \"root;\"\n",
    "ncbi_tree = Tree(ncbi_topology_nwk, format=1)\n",
    "print(ncbi_tree)\n",
    "ncbi_tree.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "94d289d6-6c28-432b-9b8b-6060ef558e0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read taxonomy file and save it in a dict\n",
    "crass_taxonomy = dict()\n",
    "\n",
    "tax_file = \"data/crass_taxonomy.txt\"\n",
    "lines = [line.strip().split(\"\\t\") for line in open(tax_file).readlines()]\n",
    "for line in lines:\n",
    "    crass_taxonomy[line[0]] = line[2]\n",
    "\n",
    "families = list(set(crass_taxonomy.values()))\n",
    "families.remove(\"outgroup\")\n",
    "\n",
    "# read TerL tree\n",
    "terl_tree = Tree(\"data/species_tree.nwk\", format=1)\n",
    "# assign taxonomy\n",
    "for leaf in terl_tree.iter_leaves():\n",
    "    # check if the leaf comes from the reference set\n",
    "    genome = leaf.name.split(\"|\")[0]   \n",
    "    leaf.add_features(family=crass_taxonomy[genome], genome=genome)\n",
    "\n",
    "\n",
    "\n",
    "outgs_leaves = terl_tree.search_nodes(family=\"outgroup\")\n",
    "outgs_lca = terl_tree.get_common_ancestor(outgs_leaves)\n",
    "# reroot the tree \n",
    "terl_tree.set_outgroup(outgs_lca)\n",
    "\n",
    "# change leaves names to genomes_ids only\n",
    "for leaf in terl_tree.iter_leaves():\n",
    "    leaf.name = leaf.genome\n",
    "    \n",
    "## collapse nodes, just for ease of visualization\n",
    "#for family in families:\n",
    "#    fam_leaves = terl_tree.search_nodes(family=family)\n",
    "#    fam_lca = terl_tree.get_common_ancestor(fam_leaves)\n",
    "#    fam_lca.name = family\n",
    "#    for child in fam_lca.get_children():\n",
    "#        child.detach()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "b15af8f8-d035-4c7c-8a89-c18c0f1aab86",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(terl_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "a543b039-d474-422a-b34d-261583ad679d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Tree node '' (0x7f79f53f5d3)"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# merge both trees in the crAss-like clade of the ncbi tree\n",
    "ncbi_crass_node = ncbi_tree.search_nodes(name=\"1978007\")[0]\n",
    "ncbi_crass_node.add_child(terl_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "959145ba-127e-414f-8d12-12930aec7155",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(ncbi_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "29a09d7b-27dc-4a9a-a6bf-54718b0f73de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write to final species tree file\n",
    "ncbi_tree.write(format=1, outfile=\"data/species_tree_final.nwk\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "094e8a0e-f5a4-4134-b75a-80c03eda3d4e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "878eb31f-7be5-4374-ba9d-e5660102bdf4",
   "metadata": {},
   "source": [
    "## TESTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3314bfd6-d250-47b7-9c66-44e507f7a759",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = list(set(target_proteins_lineages))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "034beac3-43d1-4a70-9c76-ab7ae54687c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "            /-Proteobacteria\n",
      "           |\n",
      "           |--Acidobacteria\n",
      "    /Bacteria\n",
      "   |       |--Firmicutes\n",
      "-root      |\n",
      "   |        \\-Bacteroidetes\n",
      "   |\n",
      "    \\-crAss-like viruses\n"
     ]
    }
   ],
   "source": [
    "\n",
    "target_taxids = [57723, 976, 1239, 1224]\n",
    "\n",
    "# add crAss-like taxid \n",
    "target_taxids.append(1978007)\n",
    "\n",
    "# get the ncbi tree for these tax_ids\n",
    "ncbi_topology = ncbi.get_topology(target_taxids)\n",
    "print(ncbi_topology.get_ascii(attributes=[\"sci_name\"]))\n",
    "# root the tree manually\n",
    "ncbi_topology_nwk = ncbi_topology.write()\n",
    "ncbi_topology_nwk = ncbi_topology_nwk[:-1] + \"root;\"\n",
    "ncbi_tree = Tree(ncbi_topology_nwk, format=1)\n",
    "#print(ncbi_tree)\n",
    "#ncbi_tree.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a300ccf1-bdcb-4b66-9257-013934132bf5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}