{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Curating taxonomy using tax2tree" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[**tax2tree**](https://github.com/wasade/tax2tree) ([McDonald et al., 2012](https://www.nature.com/articles/ismej2011139)) automatically curates taxonomic annotations of tree tips and internal nodes based on tree topology. This notebook explains how tax2tree is executed and parses its output files to generate more information." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Usage of tax2tree" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Installation" ] }, { "cell_type": "markdown", "metadata": { "code_folding": [] }, "source": [ "```conda create -n tax2tree numpy Click future==0.13.1 scikit-bio==0.2.3 python=2\n", "source activate tax2tree\n", "pip install git+https://github.com/wasade/tax2tree.git\n", "source activate tax2tree```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Input files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A tree in Newick format:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "((((G000830275,G000402355)N8:1.380276217937432,(G000830295,G001871415)N9:0.3229406720234299)N4:0.3513445134276098,((((G001940645,(G001940755,G001940725)N42:0.0861855240487885)N25:0.33699296079495406,(G001940655,(G001563325,(G001563335,G001940705)N67:1.0079577138891038)N43:1.137424498768485)N26:0.08418174824796781)N15:0.06949376998728227,(((G001515215,G001515205)N44:1.0827559306916261,((G001587675,\r\n" ] } ], "source": [ "!cat tree.nwk | cut -c1-400" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A genome ID to lineage str (taxa at the seven ranks) map:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "G000005825\tk__Bacteria; p__Firmicutes; c__Bacilli; o__Bacillales; f__Bacillaceae_D; g__Bacillus_S; s\r\n", "G000006175\tk__Archaea; p__Euryarchaeota; c__Methanococci; o__Methanococcales; f__Methanococcaceae; g\r\n", "G000006605\tk__Bacteria; p__Actinobacteriota; c__Actinobacteria; o__Corynebacteriales; f__Corynebacte\r\n", "G000006725\tk__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Xanthomonadales; f__Xanthomona\r\n", "G000006745\tk__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Enterobacterales; f__Vibrionac\r\n" ] } ], "source": [ "!head -n5 lineages.txt | cut -c1-100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run tax2tree" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Decorate the tree with taxonomy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```t2t decorate -t tree.nwk -m lineages.txt -o decorated.nwk > tax2tree.stdout\n", "mv decorated.nwk-consensus-strings consensus_strings.txt```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assess the consistency between taxdump and topology." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```t2t consistency -t tree.nwk -m lineages.txt -o consistency.tsv >> tax2tree.stdout```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Output files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tree with selected nodes decorated with lineage strings:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "((((G000830275,G000402355)'p__Micrarchaeota_2; c__Iainarchaeia; o__Iainarchaeales':1.38027621794,(G000830295,G001871415)'p__Nanoarchaeota; c__Woesearchaeia':0.322940672023):0.351344513428,((((G001940645,(G001940755,G001940725):0.0861855240488)'c__Heimdallarchaeia':0.336992960795,(G001940655,(G001563325,(G001563335,G001940705):1.00795771389)'f__Thorarchaeaceae':1.13742449877)'c__Lokiarchaeia; o__Th\r\n" ] } ], "source": [ "!cat decorated.nwk | cut -c1-400" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consensus strings (taxon to lineage map):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "G000830275\tk__Archaea; p__Micrarchaeota_2; c__Iainarchaeia; o__Iainarchaeales; f__; g__; s__\r\n", "G000402355\tk__Archaea; p__Micrarchaeota_2; c__Iainarchaeia; o__Iainarchaeales; f__; g__; s__\r\n", "G000830295\tk__Archaea; p__Nanoarchaeota; c__Woesearchaeia; o__; f__; g__; s__\r\n", "G001871415\tk__Archaea; p__Nanoarchaeota; c__Woesearchaeia; o__; f__; g__; s__\r\n", "G001940645\tk__Archaea; p__Asgardarchaeota; c__Heimdallarchaeia; o__; f__; g__; s__\r\n" ] } ], "source": [ "!head -n5 consensus_strings.txt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consistency analysis result:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Taxon\tCount\tConsistency\r\n", "k__Archaea\t645\t1.000\r\n", "k__Bacteria\t9086\t1.000\r\n", "p__Hydrogenedentota\t3\t1.000\r\n", "p__SM23-31\t1\t1.000\r\n" ] } ], "source": [ "!head -n5 consistency.tsv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preparation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dependencies" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from skbio import TreeNode" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Input files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Original tree (with unique internal node IDs)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "tree_ori_fp = 'tree.nwk'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Decorated tree (internal node IDs removed and some nodes appended lineage strings)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "tree_dec_fp = 'decorated.nwk'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Original genome to lineage map." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "lineage_ori_fp = 'lineages.txt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consensus strings (genom to lineage string map)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "cons_str_fp = 'consensus_strings.txt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consistency analysis result." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "consist_fp = 'consistency.tsv'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ranks and codes." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "ranks = ['kingdom', 'phylum', 'class', 'order', 'family', 'genus', 'species']\n", "code2rank = {x[0]: x for x in ranks}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Utilities" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read genome to lineage map into dataframe." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def lineages_to_df(lineage_fp):\n", " g2lin = {}\n", " with open(lineage_fp, 'r') as f:\n", " for line in f:\n", " g, lineage = line.rstrip('\\r\\n').split('\\t')\n", " g2lin[g] = dict(zip(ranks, [x.split('__')[1] for x in lineage.split('; ')]))\n", " df = pd.DataFrame.from_dict(g2lin, orient='index')\n", " return df.loc[sorted(g2lin, key=lambda x: int(x[1:]))]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Export by-rank dataframe as lineage strings." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def df_to_lineages(df, f):\n", " for row in df.itertuples():\n", " lin = []\n", " for i, rank in enumerate(ranks):\n", " lin.append('%s__%s' % (rank[0], row[i + 1]))\n", " f.write('%s\\t%s\\n' % (row.Index, '; '.join(lin)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Curated taxonomy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read consensus strings." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
kingdomphylumclassorderfamilygenusspecies
G000005825BacteriaFirmicutesBacilliBacillales_3Bacillaceae_D
G000006175ArchaeaEuryarchaeota_1MethanococciMethanococcalesMethanococcaceaeMethanococcus
G000006605BacteriaActinobacteriotaActinobacteriaCorynebacterialesCorynebacteriaceaeCorynebacterium
G000006725BacteriaProteobacteriaGammaproteobacteriaXanthomonadalesXanthomonadaceaeXylella
G000006745BacteriaProteobacteriaGammaproteobacteriaEnterobacteralesVibrionaceaeVibrio
\n", "
" ], "text/plain": [ " kingdom phylum class \\\n", "G000005825 Bacteria Firmicutes Bacilli \n", "G000006175 Archaea Euryarchaeota_1 Methanococci \n", "G000006605 Bacteria Actinobacteriota Actinobacteria \n", "G000006725 Bacteria Proteobacteria Gammaproteobacteria \n", "G000006745 Bacteria Proteobacteria Gammaproteobacteria \n", "\n", " order family genus species \n", "G000005825 Bacillales_3 Bacillaceae_D \n", "G000006175 Methanococcales Methanococcaceae Methanococcus \n", "G000006605 Corynebacteriales Corynebacteriaceae Corynebacterium \n", "G000006725 Xanthomonadales Xanthomonadaceae Xylella \n", "G000006745 Enterobacterales Vibrionaceae Vibrio " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfc = lineages_to_df(cons_str_fp)\n", "dfc.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note: In tax2tree's assumption, all lineages should be continuous (i.e., a continuous series of ranks should all contain taxa). However in reality some lineages contain caps. For example:\n", "\n", "`G001458535 k__Bacteria; p__Candidatus Kryptonia; c__; o__; f__; g__Candidatus Kryptonium; s__Candidatus Kryptonium thompsoni`\n", "\n", "In those cases tax2tree may produced gapped consensus strings." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convert to table." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "dfc.to_csv('consensus_ranks.tsv', sep='\\t')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyze polytomies" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Taxon names ending with `__#` are the ones identified as polytomies and split into multiple taxa by tax2tree." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "splits = {}\n", "for rank in ranks:\n", " split_names = {}\n", " for name in set(x for x in dfc[rank] if x):\n", " x = name.split('_')\n", " # convert \"taxonA_1\" into \"taxonA\"\n", " if len(x) > 1 and x[-1].isdigit():\n", " ori_name = '_'.join(x[:-1])\n", " # get all appended indices for each taxon\n", " if ori_name in split_names:\n", " split_names[ori_name].append(int(x[-1]))\n", " else:\n", " split_names[ori_name] = [int(x[-1])]\n", "\n", " splits[rank] = {}\n", " for name in sorted(split_names):\n", " indices = split_names[name]\n", " # check if the indices are incremental integers\n", " if sorted(indices) != list(range(1, max(indices) + 1)):\n", " print('%s: %s is not a split.' % (rank, name))\n", " continue\n", " splits[rank][name] = len(indices)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get original names that appear in the curated table." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "kingdom: 2, phylum: 69, class: 157, order: 347, family: 683, genus: 1443, species: 748\n" ] } ], "source": [ "ori_names = {}\n", "for rank in ranks:\n", " ori_names[rank] = set()\n", " for name in set(x for x in dfc[rank] if x):\n", " x = name.split('_')\n", " if len(x) > 1 and x[-1].isdigit():\n", " name_ = '_'.join(x[:-1])\n", " if name_ in splits[rank]:\n", " ori_names[rank].add(name_)\n", " else:\n", " ori_names[rank].add(name)\n", " else:\n", " ori_names[rank].add(name)\n", "print(', '.join('%s: %d' % (x, len(ori_names[x])) for x in ranks))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Export results." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "phylum: 3 taxa split, max = 2 (Desulfobacterota)\n", "class: 4 taxa split, max = 2 (Elusimicrobia)\n", "order: 2 taxa split, max = 5 (Bacillales)\n", "family: 2 taxa split, max = 2 (Cyanobacteriaceae)\n" ] } ], "source": [ "with open('split_taxa.tsv', 'w') as f:\n", " f.write('rank\\ttaxon\\tsplits\\n')\n", " for rank in ranks:\n", " d = splits[rank]\n", " if len(d) > 0:\n", " max_name = max(d.keys(), key=(lambda k: d[k]))\n", " print('%s: %d taxa split, max = %s (%s)'\n", " % (rank, len(d), d[max_name], max_name))\n", " for name in sorted(d):\n", " f.write('%s\\t%s\\t%d\\n' % (rank, name, d[name]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyze consistency" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
TaxonCountConsistency
0k__Archaea6451.0
1k__Bacteria90861.0
2p__Hydrogenedentota31.0
3p__SM23-3111.0
4p__Elusimicrobiota271.0
\n", "
" ], "text/plain": [ " Taxon Count Consistency\n", "0 k__Archaea 645 1.0\n", "1 k__Bacteria 9086 1.0\n", "2 p__Hydrogenedentota 3 1.0\n", "3 p__SM23-31 1 1.0\n", "4 p__Elusimicrobiota 27 1.0" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfp = pd.read_csv(consist_fp, sep='\\t')\n", "dfp.head()" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CountConsistencyNameRank
Taxon
k__Archaea6451.0Archaeakingdom
k__Bacteria90861.0Bacteriakingdom
p__Hydrogenedentota31.0Hydrogenedentotaphylum
p__SM23-3111.0SM23-31phylum
p__Elusimicrobiota271.0Elusimicrobiotaphylum
\n", "
" ], "text/plain": [ " Count Consistency Name Rank\n", "Taxon \n", "k__Archaea 645 1.0 Archaea kingdom\n", "k__Bacteria 9086 1.0 Bacteria kingdom\n", "p__Hydrogenedentota 3 1.0 Hydrogenedentota phylum\n", "p__SM23-31 1 1.0 SM23-31 phylum\n", "p__Elusimicrobiota 27 1.0 Elusimicrobiota phylum" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfp['Name'] = dfp['Taxon'].str.split('__').str.get(1)\n", "dfp['Rank'] = dfp['Taxon'].str.split('__').str.get(0).map(code2rank)\n", "dfp.set_index('Taxon', inplace=True)\n", "dfp.head()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.8, 1.05)" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGxdJREFUeJzt3Xu8VXWd//HXGxQxREU5U42KqGFKaZrnp6NYYpmp04z30jTTceLh5GXSrJ/+rHQoR/NS/kwz0SHUSsfBpjD5jZoKNuIFSEHAG5LK8ZIXlPqZWcBn/ljfHdvNOee7OGcvzubwfj4e+7HXfX/WOXuv916371ZEYGZm1p0BfV2AmZm1PoeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZVmVhIWmipJclzetivCRdLmmhpLmSPlw3brmkR9JjSlU1mplZOVXuWUwCDuhm/IHAqPQYB1xVN+6tiNglPf6+uhLNzKyMysIiIu4FlnQzycHA9VF4ANhU0nurqsfMzHpuvT587S2AxXX9HWnYi8BgSbOAZcCFEfGzzhYgaRzFXglDhgzZbYcddqi2YjOzfmb27NmvRkRbbrq+DAt1MqzW9siIiHhB0rbA3ZIejYinV5k4YgIwAaC9vT1mzZpVXbVmZv2QpGfLTNeXV0N1AFvV9W8JvAAQEbXnRcA0YNc1XZyZma3Ul2ExBTguXRX1N8DSiHhR0jBJGwBIGg6MARb0YZ1mZuu8yg5DSboRGAsMl9QBnAusDxARPwCmAgcBC4E/ACekWXcErpa0giLMLowIh4WZWR+qLCwi4ujM+ABO7mT4DGCnquoyM7PV5zu4zcwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVmWw8LMzLIcFmZmluWwMDOzLIeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLqiwsJE2U9LKkeV2Ml6TLJS2UNFfSh+vGfV7SU+nx+apqNDOzcqrcs5gEHNDN+AOBUekxDrgKQNJmwLnAHsDuwLmShlVYp5mZZaxX1YIj4l5JI7uZ5GDg+ogI4AFJm0p6LzAWuDMilgBIupMidG6sqlYz67+++tWv8tJLL/Ge97yHiy66qK/LWWtVFhYlbAEsruvvSMO6Gr4KSeMo9koYMWLEKuN3+8r1TSq1uWZffFxfl2BW2tq+sX3ppZd4/vnn1+hrPnb+3Wv09cra8ZyP9XjevgwLdTIsuhm+6sCICcAEgPb29k6nsbX/w762W9v//n2xsT3/2COatqwlLy8tnl96sSnLPedHk3u9jLVRX4ZFB7BVXf+WwAtp+NiG4dPWWFX9UF982JvJG1uzvteXYTEFOEXSTRQns5dGxIuSbgf+te6k9v7A2X1VpPU9b2x75oov39qU5bzx6pt/eW7GMk+59O96vYzVMXjggHc8W89UFhaSbqTYQxguqYPiCqf1ASLiB8BU4CBgIfAH4IQ0bomkbwIz06LG1052m5mtrl03H9rXJfQLVV4NdXRmfAAndzFuIjCxirrMzGz19eVhKLOWNv2j+zRlOW+tNxAk3uroaMoy97l3ehOqKm/IoI3f8WzrJoeFmXVrzHaH9XUJ1gIcFlaZMd8b05TlDHpjEAMYwOI3Fjdlmfedel8TqjJbt/jyADMzy3JYmJlZlg9DmVVs04h3PJutjRwWZhU7dvmKvi7BrNd8GMrMzLK8Z9HCnhu/U1OWs2zJZsB6LFvybFOWOeIbj/a+KDNbq3jPwszMsrxnYS0v3hWsYAXxLp8gNusrDgtreX8e8+e+LsFsnefDUGZmluWwMDOzLIeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZbltqHXA8MErgGXp2cxs9Tks1gFn7vxGX5dgZms5H4YyM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVlWqbCQNLDqQszMrHWV3bNYKOliSaMrrcbMzFpS2bDYGXgSuFbSA5LGSdq4wrrMzKyFlAqLiPh9RFwTEXsBXwXOBV6UdJ2k93U1n6QDJD0haaGkszoZv7WkuyTNlTRN0pZ145ZLeiQ9pvRg3czMrElKNfeRzln8LXACMBK4FPgx8BFgKrB9F/NcCXwC6ABmSpoSEQvqJrsEuD4irpP0MeAC4HNp3FsRsUtPVsrMzJqrbNtQTwH3ABdHxIy64ZMlfbSLeXYHFkbEIgBJNwEHA/VhMRo4PXXfA/ysbOFmZrbmlD5nEREnNgQFABFxWhfzbAEsruvvSMPqzQEOT92HAkMlbZ76B0ualc6RHNLZC6RzJ7MkzXrllVdKroqZma2usmFxpaRNaz2ShkmamJlHnQyLhv4zgX0kPQzsAzwPLEvjRkREO/BZ4DJJ262ysIgJEdEeEe1tbW0lV8XMzFZX2cNQO0fEX9q5jojXJe2amacD2Kquf0vghfoJIuIF4DAASRsBh0fE0rpxRMQiSdOAXYGnS9ZrZmZNVHbPYoCkYbUeSZuRD5qZwChJ20gaBBwFvOOqJknDJdVqOBuYmIYPk7RBbRpgDO8812FmZmtQ2T2LS4EZkian/iOB87ubISKWSToFuB0YCEyMiPmSxgOzImIKMBa4QFIA9wInp9l3BK6WtIIi0C5suIrKzMzWoFJhERHXS5oN7EtxLuKwMhvviJhKcWlt/bBv1HVPBiZ3Mt8MYKcytZmZWfVW52dVHwder80jaUREPFdJVWZm1lLK3pR3KsVd278FllPsXQRFMyBmZtbPld2z+Gfg/RHxWpXFmJlZayp7NdRiYGmVhZiZWesqu2exCJgm6Tbg7drAiPhOJVWZmVlLKRsWz6XHoPQwM7N1SNlLZ/8FQNKQiHiz2pLMzKzVlP1Z1T0lLQAeS/0fkvT9SiszM7OWUfYE92XAJ4HXACJiDtBV0+RmZtbPlA0LImJxw6DlTa7FzMxaVNkT3Isl7QVEahTwNNIhKTMz6//K7lmcRNHI3xYUTY/vAnyxqqLMzKy1lN2zeH9EHFM/QNIY4L7ml2RmZq2m7J7F90oOMzOzfqjbPQtJewJ7AW2SzqgbtTHFb1SYmdk6IHcYahCwUZpuaN3w3wFHVFWUmZm1lm7DIiKmA9MlTYqIZwHSz6BuFBG/WxMFmplZ3yt7zuICSRtLGkLxW9hPSPpKhXWZmVkLKRsWo9OexCEUP5M6AvhcZVWZmVlLKRsW60tanyIsfh4Rf6b4pTwzM1sHlA2Lq4FngCHAvZK2pjjJbWZm64CyTZRfDlxeN+hZSftWU5KZmbWa3H0Wx0bEjxrusajnX8ozM1sH5PYshqTnod1OZWZm/VruPour0/O/rJlyzMysFZX9pbyL0n0W60u6S9Krko6tujgzM2sNZa+G2j/dZ/EpiibKtwd8U56Z2Tqi9H0W6fkg4MaIWFJRPWZm1oLK/p7FrZIeB94CviipDfhjdWWZmVkrKbVnERFnAXsC7enu7TeBg6sszMzMWkfZPQuAHYGRkurnub7J9ZiZWQsqFRaSbgC2Ax4BlqfBgcPCzGydUHbPop2i5Vk3Hmhmtg4qezXUPOA9VRZiZmatq+yexXBggaSHgLdrAyPi7yupyszMWkrZsDivJwuXdADwf4GBwLURcWHD+K2BiUAbsAQ4NiI60rjPA19Lk34rIq7rSQ1mZtZ7ZS+dnQ48TtGg4FDgsTSsS5IGAlcCBwKjgaMljW6Y7BLg+ojYGRgPXJDm3Qw4F9gD2B04V9KwsitlZmbNVbZtqE8DDwFHAp8GHpR0RGa23YGFEbEoIv4E3MSq92aMBu5K3ffUjf8kcGdELImI14E7gQPK1GpmZs1X9jDUOcD/ioiXAdId3L8EJnczzxbA4rr+Doo9hXpzgMMpDlUdCgyVtHkX827R+AKSxgHjAEaMGFFyVczMbHWVvRpqQC0oktdKzKtOhjVeensmsI+kh4F9gOeBZSXnJSImRER7RLS3tbVlyjEzs54qu2fxX5JuB25M/Z8Bpmbm6QC2quvfEnihfoKIeAE4DEDSRsDhEbFUUgcwtmHeaSVrNTOzJut270DS+ySNiYivAFcDOwMfAu4HJmSWPRMYJWkbSYOAo4ApDcsfLqlWw9kUV0YB3A7sL2lYOrG9fxpmZmZ9IHco6TLg9wAR8dOIOCMiTqfYq7isuxkjYhlwCsVG/jHg5oiYL2m8pNr9GWOBJyQ9CbwbOD/NuwT4JkXgzATGu1l0M7O+kzsMNTIi5jYOjIhZkkbmFh4RU2k4XBUR36jrnkwXJ8kjYiIr9zTMzKwP5fYsBnczbsNmFmJmZq0rFxYzJX2hcaCkE4HZ1ZRkZmatJncY6kvAf0o6hpXh0A4MorgvwszM1gHdhkVE/BbYS9K+wAfT4Nsi4u7KKzMzs5ZR6j6LiLiHojkOMzNbB5W9g9vMzNZhDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVmWw8LMzLIcFmZmluWwMDOzLIeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaWVWlYSDpA0hOSFko6q5PxIyTdI+lhSXMlHZSGj5T0lqRH0uMHVdZpZmbdW6+qBUsaCFwJfALoAGZKmhIRC+om+xpwc0RcJWk0MBUYmcY9HRG7VFWfmZmVV+Wexe7AwohYFBF/Am4CDm6YJoCNU/cmwAsV1mNmZj1UZVhsASyu6+9Iw+qdBxwrqYNir+LUunHbpMNT0yV9pLMXkDRO0ixJs1555ZUmlm5mZvWqDAt1Miwa+o8GJkXElsBBwA2SBgAvAiMiYlfgDOAnkjZumJeImBAR7RHR3tbW1uTyzcyspsqw6AC2quvfklUPM50I3AwQEfcDg4HhEfF2RLyWhs8Gnga2r7BWMzPrRpVhMRMYJWkbSYOAo4ApDdM8B3wcQNKOFGHxiqS2dIIcSdsCo4BFFdZqZmbdqOxqqIhYJukU4HZgIDAxIuZLGg/MiogpwJeBaySdTnGI6viICEkfBcZLWgYsB06KiCVV1WpmZt2rLCwAImIqxYnr+mHfqOteAIzpZL5bgFuqrM3MzMrzHdxmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWU5LMzMLMthYWZmWQ4LMzPLcliYmVmWw8LMzLIcFmZmluWwMDOzLIeFmZllOSzMzCzLYWFmZlkOCzMzy3JYmJlZlsPCzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8tyWJiZWZbDwszMshwWZmaW5bAwM7Msh4WZmWVVGhaSDpD0hKSFks7qZPwISfdIeljSXEkH1Y07O833hKRPVlmnmZl1b72qFixpIHAl8AmgA5gpaUpELKib7GvAzRFxlaTRwFRgZOo+CvgA8NfALyVtHxHLq6rXzMy6VuWexe7AwohYFBF/Am4CDm6YJoCNU/cmwAup+2Dgpoh4OyJ+AyxMyzMzsz6giKhmwdIRwAER8Y+p/3PAHhFxSt007wXuAIYBQ4D9ImK2pCuAByLiR2m6fwP+X0RMbniNccC41Pt+4IlKVqYwHHi1wuVXzfX3Ldfft9bm+quufeuIaMtNVNlhKECdDGtMpqOBSRFxqaQ9gRskfbDkvETEBGBCrystQdKsiGhfE69VBdfft1x/31qb62+V2qsMiw5gq7r+LVl5mKnmROAAgIi4X9JgihQtM6+Zma0hVZ6zmAmMkrSNpEEUJ6ynNEzzHPBxAEk7AoOBV9J0R0naQNI2wCjgoQprNTOzblS2ZxERyySdAtwODAQmRsR8SeOBWRExBfgycI2k0ykOMx0fxUmU+ZJuBhYAy4CTW+BKqDVyuKtCrr9vuf6+tTbX3xK1V3aC28zM+g/fwW1mZlkOCzMzy+q3YSFppKR5DcPaJV1e1fLXJEnPSBq+GtOfJ+nMKmvqrbWhxp6SdHy6f6gvazhN0mOSftzL5YyXtF/qniapzy/rXNdImipp0zX5mlVeOttyImIWMKuv67D+TZIozgeu6MUy1ouIZU0sC+CLwIGpVYQei4hvNKke66GIOCg/VXP12z2LepK2TY0VfkXSL9Kw8yRNTN+MFkk6rW76r0t6XNKdkm6sfduVtJukOZLuB06um36wpB9KejS9zr5p+PGSfibpVkm/kXSKpDPSNA9I2qxE7SNTLdelxhYnS3pXGn2qpF+n191B0gBJT0lqS/MOSI0xDm9Y5l++DUoaLumZZtW7OiQdl9ZpjqQbGsZ9QdLMNO6W2jpLOlLSvDT83jTsA5IekvRIWt6oZtbZRe1npDrmSfpS+j89Jun7wK+BrSSdIOlJSdOBMXXztqV1mpkeY9Lw8yRNkHQHcH2T6/0BsC0wRdL/ljQj/V9nSHp/mqbU/1/SJBUtNNQv/0RJ363r/4Kk7zRzHdJyV/lsStpO0n9Jmi3pV5J2qKvz8rSOi2o1SxqrtB1I/VdIOj51XyhpQXofXdLLWodIui29V+dJ+oyKIwLfTu/XhyS9L03b1XtiI63ctsyVdHga/pcjC5KOrXv/Xy1pYHpMSq/7qIorTnsnIvrlAxgJzKNoBuRhYBdgLPCLNP48YAawAcWNgK8B6wPtwCPAhsBQ4CngzDTPXGCf1H0xMC91fxn4YeregeL+kcHA8RTtWg0F2oClwElpuu8CXyq5HgGMSf0TgTOBZ4BT07AvAtem7nNrywX2B26pW9/aekwD2lP3cOCZ1N3relfj//MBiuZZhqf+zRpq3Lxu2m/VreujwBape9P0/D3gmNQ9CNiw4vfWbqmOIcBGwHxgV2AF8Ddpmvem90Fbquk+4Io07ifA3ql7BPBY3f9odlX1p/fMcIr22NZLw/are4+U+v8Dk4Aj6t9L6W/xNLB+Gj4D2KnJ9Xf62QTuAkalafYA7q6r8z8ovhSPpmirDuq2A6n/irTum6X3ZO0q0U17We/hwDV1/Zuk/8E5qf84Vm6PunpPfBu4rG4Zwxr+lzsCt9b93b+flrsbcGfdfL1al4jo94eh2oCfA4dHcY/H2Ibxt0XE28Dbkl4G3g3sDfw8It4CkHRret6E4g8+Pc17A3Bg6t6bYoNFRDwu6Vlg+zTunoj4PfB7SUsp/rFQbGx2LrkeiyPivtT9I6C2F/TT9DwbOCx1T0zrfBnwD8APS75GTTPqLeNjwOSIeBUgIpZI72jl5YOSvgVsSrFBvj0Nvw+YpOI+nNr63w+cI2lL4KcR8VQT6+zM3sB/RsSbAJJ+CnwEeDYiHkjT7AFMi4hX0jT/zsr3xH7A6Lr13VjS0NQ9pfbeq9AmwHVpDywoviTV9Oj/HxFvSrob+JSkxyg2Xo82ue7OPpuDgb2A/6j7e25QN8/PojgcuEDSuzPL/x3wR+BaSbcBv8hMn/MocImkb1OEwq9SjTem8TdShDB0/Z7Yj+KGZgAi4vWG1/g4RTDMTPNuCLxM8X/bVtL3gNso2uDrlf4eFkuBxRSHAOZ3Mv7tuu7lFH+PztqlIg3v6qaUruZpfI0Vdf0rKP/3b3zdWn9tWbXaiYjFkn4r6WMUG6xjOlneMlYeghxcQb1ldPf3hOJb4SERMScdIhgLEBEnSdoD+FvgEUm7RMRPJD2Yht0u6R8j4u4m1tpZ7Z15s6G/q/UbAOzZGArpw964jCp8kyIUDpU0kmLvoKY3//9rgf8DPM7qf0kpo7O/+wDgjYjYpYt56tenNn/9+x/SZyCKG4l3p9gAHwWcQvGlpkci4klJuwEHARekw4vwzvdFrbur90TucyLguog4e5UR0oeAT1IcMv80xZfHHuvv5yz+BBwCHCfpsyXn+W/g71Sch9iIYgNERLwBLJW0d5qufiN8b61f0vYUu5HNbAF3hIqGFqFofPG/M9NfS7EHcnN0fuf7MxTfRgCO6GT8mnAX8GlJmwNo1fMhQ4EXJa1P3d9a0nYR8WAUJ1lfpTg3sC2wKCIup2gqppl7QJ25FzhE0rskDQEOBX7VMM2DwFhJm6d1OLJu3B0UGyIAJHW1oavKJsDzqfv4Zi00Ih6kaNPts6z89txMnX02/wD8RtKRUGxc00ayO89SfIvfIB0xqDU5tBGwSURMBb5Ecei6xyT9NfCHKFrPvgT4cBr1mbrn+1N3V++JxuHDGl7mLuAISX+Vxm8maet0PmNARNwCfL3utXusv4cF6VDBp4DTKT4kuelnUmxw5lAc5phFsYcCcAJwpYoT3PXfAL4PDJT0KPDvFM2W1H+j6a3HgM9LmktxXPWqzPRTKA7ddPXt7hLgnyTNoDjuucZFxHzgfGC6pDlA48nQr1NscO+k+KZac3E6YTePYqM9h+JDN0/SIxTnjJp6criT2n9NsefzUKrxWuD1hmlepDgHcT/wS4qT3jWnAe3phOUC4KQq6+3ERRTfdO+jaIqnmW4G7uvkcEmvdfPZPAY4Mb2P5rPq7+Y0LmdxqnMu8GOKc5pQfEH5RfqcTafYZvTGTsBD6X15DsW5N4AN0p7wP9e9RlfviW8Bw9KJ6jnAvg3rsoDiR+TuSHXfSXG+bAtgWnrtScAqex6ry819dELSRhHx/1VcgXMvMC5tIPqilpEUxzs/uBrztAPfjYiPVFWXWWdUXGX03Yi4q6Llt8xnsydUXHnYXjtXtzbp7+csemqCip92HUxxPHBtejOeBfwTnZ+rMKuEihvEHgLmVBUUyVr72Vzbec/CzMyy+v05CzMz6z2HhZmZZTkszMwsy2Fh1gOSlqe2eOapaEupxy2Ayi232lrAYWHWM29FxC7pkuYl1DUsadYfOSzMeu9+ipugaq2E3qWVrQEfnIbXWqW9RtJ8SXdI2rB+ISpaCb4utYll1lIcFma9IGkgRXMRU9KgPwKHRsSHKe62vVQrW4cbBVwZER8A3qBolbRmPYq7iZ+MiK+tkeLNVoPDwqxnNkxNKbxG0QTLnWm4gH9NTS/8kmKPo9ba6W8i4pHUPZui+fmaqymavD+/6sLNesJhYdYzb6WWTrem+L2K2jmLYyiaxt8tjf8tK1v27ayV45oZwL6SGlsBNmsJDguzXoiIpRSNwJ2ZWpfdBHg5Iv6s4hcTty65qH8DplL8LoOb4bGW47Aw66WIeJiiJdSjKM47tEuaRbGX8Xh38zYs5zsUrdPeIMmfTWspbhvKzMyy/O3FzMyyHBZmZpblsDAzsyyHhZmZZTkszMwsy2FhZmZZDgszM8v6H+vqW93q28iqAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "sns.barplot(x='Rank', y='Consistency', data=dfp, ax=ax)\n", "ax.set_ylim(0.8,)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "fig.tight_layout()\n", "fig.savefig('consistency.pdf', bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fill spaces with unique taxa" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "tax2tree-curated taxa are the ones that occur more than once in the table. However many genomes are the only representatives of some taxa, especially the species level. Therefore, the empty cells will be filled with unique taxa when applicable.\n", "\n", "Note the those added taxa are NOT tax2tree-curated taxa. They are but unique identifiers of certain genomes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read original lineages." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
kingdomphylumclassorderfamilygenusspecies
G000005825BacteriaFirmicutesBacilliBacillalesBacillaceae_DBacillus_SBacillus_S pseudofirmus
G000006175ArchaeaEuryarchaeotaMethanococciMethanococcalesMethanococcaceaeMethanococcusMethanococcus voltae
G000006605BacteriaActinobacteriotaActinobacteriaCorynebacterialesCorynebacteriaceaeCorynebacteriumCorynebacterium jeikeium
G000006725BacteriaProteobacteriaGammaproteobacteriaXanthomonadalesXanthomonadaceaeXylellaXylella fastidiosa
G000006745BacteriaProteobacteriaGammaproteobacteriaEnterobacteralesVibrionaceaeVibrioVibrio cholerae
\n", "
" ], "text/plain": [ " kingdom phylum class \\\n", "G000005825 Bacteria Firmicutes Bacilli \n", "G000006175 Archaea Euryarchaeota Methanococci \n", "G000006605 Bacteria Actinobacteriota Actinobacteria \n", "G000006725 Bacteria Proteobacteria Gammaproteobacteria \n", "G000006745 Bacteria Proteobacteria Gammaproteobacteria \n", "\n", " order family genus \\\n", "G000005825 Bacillales Bacillaceae_D Bacillus_S \n", "G000006175 Methanococcales Methanococcaceae Methanococcus \n", "G000006605 Corynebacteriales Corynebacteriaceae Corynebacterium \n", "G000006725 Xanthomonadales Xanthomonadaceae Xylella \n", "G000006745 Enterobacterales Vibrionaceae Vibrio \n", "\n", " species \n", "G000005825 Bacillus_S pseudofirmus \n", "G000006175 Methanococcus voltae \n", "G000006605 Corynebacterium jeikeium \n", "G000006725 Xylella fastidiosa \n", "G000006745 Vibrio cholerae " ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfo = lineages_to_df(lineage_ori_fp)\n", "dfo.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Identify taxa from the original lineages that can be filled to trailing spaces following consensus strings." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "kingdom: 0, phylum: 20, class: 39, order: 133, family: 377, genus: 2033, species: 4751\n" ] } ], "source": [ "taxa_to_add = {x: {} for x in ranks}\n", "for g in dfc.index:\n", " for rank in ranks:\n", " if dfc[rank][g] == '' and g in dfo.index:\n", " taxon = dfo[rank][g]\n", " if taxon != '' and not taxon in ori_names[rank]:\n", " taxa_to_add[rank][g] = taxon\n", "print(', '.join('%s: %d' % (x, len(taxa_to_add[x])) for x in ranks))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exclude taxa that occur multiple times at each rank." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "kingdom: 0, phylum: 20, class: 39, order: 125, family: 367, genus: 2027, species: 4749\n" ] } ], "source": [ "for rank in ranks:\n", " used, dups = set(), set()\n", " for g, taxon in taxa_to_add[rank].items():\n", " if taxon in used:\n", " dups.add(taxon)\n", " used.add(taxon)\n", " new_dict = {}\n", " for g, taxon in taxa_to_add[rank].items():\n", " if taxon not in dups:\n", " new_dict[g] = taxon\n", " taxa_to_add[rank] = new_dict\n", "print(', '.join('%s: %d' % (x, len(taxa_to_add[x])) for x in ranks))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add those taxa to the lineage table." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
kingdomphylumclassorderfamilygenusspecies
G000005825BacteriaFirmicutesBacilliBacillales_3Bacillaceae_DBacillus_SBacillus_S pseudofirmus
G000006175ArchaeaEuryarchaeota_1MethanococciMethanococcalesMethanococcaceaeMethanococcusMethanococcus voltae
G000006605BacteriaActinobacteriotaActinobacteriaCorynebacterialesCorynebacteriaceaeCorynebacteriumCorynebacterium jeikeium
G000006725BacteriaProteobacteriaGammaproteobacteriaXanthomonadalesXanthomonadaceaeXylellaXylella fastidiosa
G000006745BacteriaProteobacteriaGammaproteobacteriaEnterobacteralesVibrionaceaeVibrioVibrio cholerae
\n", "
" ], "text/plain": [ " kingdom phylum class \\\n", "G000005825 Bacteria Firmicutes Bacilli \n", "G000006175 Archaea Euryarchaeota_1 Methanococci \n", "G000006605 Bacteria Actinobacteriota Actinobacteria \n", "G000006725 Bacteria Proteobacteria Gammaproteobacteria \n", "G000006745 Bacteria Proteobacteria Gammaproteobacteria \n", "\n", " order family genus \\\n", "G000005825 Bacillales_3 Bacillaceae_D Bacillus_S \n", "G000006175 Methanococcales Methanococcaceae Methanococcus \n", "G000006605 Corynebacteriales Corynebacteriaceae Corynebacterium \n", "G000006725 Xanthomonadales Xanthomonadaceae Xylella \n", "G000006745 Enterobacterales Vibrionaceae Vibrio \n", "\n", " species \n", "G000005825 Bacillus_S pseudofirmus \n", "G000006175 Methanococcus voltae \n", "G000006605 Corynebacterium jeikeium \n", "G000006725 Xylella fastidiosa \n", "G000006745 Vibrio cholerae " ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfl = dfc.copy(deep=True)\n", "for i, rank in enumerate(ranks):\n", " for g, taxon in taxa_to_add[rank].items():\n", " dfl[rank][g] = taxon\n", "dfl.head()" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "dfl.to_csv('filled_ranks.tsv', sep='\\t')" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "with open('filled_lineages.txt', 'w') as f:\n", " df_to_lineages(dfl, f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analyze node decoration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Match decorations and node IDs." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10575" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tree_ori = TreeNode.read(tree_ori_fp)\n", "tree_ori.count(tips=True)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10575" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tree_dec = TreeNode.read(tree_dec_fp)\n", "tree_dec.count(tips=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confirm that the decorated tree has the same topology and order of taxa as the original tree." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tree_ori.compare_rfd(tree_dec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read internal node IDs from the original tree." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "nids = []\n", "for node in tree_ori.levelorder():\n", " if not node.is_tip():\n", " nids.append(node.name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read node decorations." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[None, 'k__Archaea', 'k__Bacteria', None, None]" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "decos = []\n", "for node in tree_dec.levelorder():\n", " if not node.is_tip():\n", " decos.append(node.name)\n", "decos[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generate a rank name table." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "nid2label = {}\n", "nid2ranks = {}\n", "for i in range(len(nids)):\n", " nid = nids[i]\n", " deco = decos[i]\n", " nid2ranks[nid] = {x: '' for x in ranks}\n", " if deco is not None and '__' in deco:\n", " deco = deco.split(':')[-1] # remove node support in case there is one\n", " nid2label[nid] = deco\n", " for taxon in deco.split('; '):\n", " code, name = taxon.split('__')\n", " nid2ranks[nid][code2rank[code]] = name" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "with open('decorations.txt', 'w') as f:\n", " for nid in sorted(nid2label, key=lambda x: int(x.lstrip('N'))):\n", " f.write('%s\\t%s\\n' % (nid, nid2label[nid]))" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
kingdomphylumclassorderfamilygenusspecies
N1
N2Archaea
N3Bacteria
N4
N5
\n", "
" ], "text/plain": [ " kingdom phylum class order family genus species\n", "N1 \n", "N2 Archaea \n", "N3 Bacteria \n", "N4 \n", "N5 " ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfd = pd.DataFrame.from_dict(nid2ranks, orient='index')\n", "dfd = dfd.loc[sorted(nid2ranks, key=lambda x: int(x[1:]))]\n", "dfd.head()" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "dfd.to_csv('decorations_by_rank.tsv', sep='\\t')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate lineage for every node" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At a given node, at each rank, check if one taxon can describe all descendants." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "g2ranks = dfl.to_dict('index')" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "nid2ranks, nid2lineage = {}, {}\n", "for nid in nids:\n", " nid2ranks[nid] = {x: '' for x in ranks}\n", " nid2lineage[nid] = []\n", " gs = tree_ori.find(nid).subset()\n", " for rank in ranks:\n", " nid2lineage[nid].append('%s__' % rank[0])\n", " names = set(g2ranks[x][rank] for x in gs)\n", " if len(names) == 1:\n", " name = max(names)\n", " nid2ranks[nid][rank] = name\n", " nid2lineage[nid][-1] += name" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "with open('node_lineages.txt', 'w') as f:\n", " for nid in sorted(nid2lineage, key=lambda x: int(x[1:])):\n", " f.write('%s\\t%s\\n' % (nid, '; '.join(nid2lineage[nid])))" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
kingdomphylumclassorderfamilygenusspecies
N1
N2Archaea
N3Bacteria
N4Archaea
N5Archaea
\n", "
" ], "text/plain": [ " kingdom phylum class order family genus species\n", "N1 \n", "N2 Archaea \n", "N3 Bacteria \n", "N4 Archaea \n", "N5 Archaea " ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfn = pd.DataFrame.from_dict(nid2ranks, orient='index')\n", "dfn = dfn.loc[sorted(nid2ranks, key=lambda x: int(x[1:]))]\n", "dfn.head()" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "dfn.to_csv('node_ranks.tsv', sep='\\t')" ] } ], "metadata": { "hide_input": false, "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.7.2" } }, "nbformat": 4, "nbformat_minor": 2 }