tax_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# `sourmash tax` subcommands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**for integrating taxonomic information**\n",
    "\n",
    "The sourmash `tax` or `taxonomy` commands integrate taxonomic information into the results of `sourmash gather`. `tax` commands  require a properly formatted `taxonomy` csv file that corresponds to the database used for `gather`. For supported databases (e.g. GTDB, NCBI), we provide these files, but they can also be generated for user-generated databases. For more information, see [databases](databases.md).\n",
    "\n",
    "These commands rely upon the fact that `gather` results are non-overlapping: the fraction match for gather on each query will be between 0 (no database matches) and 1 (100% of query matched). We use this property to aggregate gather matches at the desired taxonomic rank. For example, if the gather results for a metagenome include results for 30 different strains of a given species, we can sum the fraction match to each strain to obtain the fraction match to this species.\n",
    "\n",
    "As with all reference-based analysis, results can be affected by the completeness of the reference database. However, summarizing taxonomic results from `gather` minimizes issues associated with increasing size and redundancy of reference databases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All `tax` commands require a lineage file for the sourmash database used with `gather`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### download and look at gtdb-rs202 lineage file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100   483  100   483    0     0    133      0  0:00:03  0:00:03 --:--:--   133\n",
      "100 35.2M  100 35.2M    0     0  3018k      0  0:00:11  0:00:11 --:--:-- 11.3M\n"
     ]
    }
   ],
   "source": [
    "!curl -L https://osf.io/p6z3w/download -o gtdb-rs202.taxonomy.v2.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ident,superkingdom,phylum,class,order,family,genus,species\n",
      "GCF_014075335.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_002310555.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_900013275.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_000168095.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_002459845.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_001614695.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_000356585.2,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_014528595.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n",
      "GCF_900196445.1,d__Bacteria,p__Proteobacteria,c__Gammaproteobacteria,o__Enterobacterales,f__Enterobacteriaceae,g__Escherichia,s__Escherichia flexneri\n"
     ]
    }
   ],
   "source": [
    "!head gtdb-rs202.taxonomy.v2.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All `tax` commands work on the output files from `sourmash gather`\n",
    "\n",
    "example gather output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intersect_bp,f_orig_query,f_match,f_unique_to_query,f_unique_weighted,average_abund,median_abund,std_abund,name,filename,md5,f_match_orig,unique_intersect_bp,gather_result_rank,remaining_bp,query_filename,query_name,query_md5,query_bp\n",
      "442000,0.08815317112086159,0.08438335242458954,0.08815317112086159,0.05815279361459521,1.6153846153846154,1.0,1.1059438185997785,\"GCF_001881345.1 Escherichia coli strain=SF-596, ASM188134v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,683df1ec13872b4b98d59e98b355b52c,0.042779713511420826,442000,0,4572000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n",
      "390000,0.07778220981252493,0.10416666666666667,0.07778220981252493,0.050496823586903404,1.5897435897435896,1.0,0.8804995294906566,\"GCF_009494285.1 Prevotella copri strain=iAK1218, ASM949428v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,1266c86141e3a5603da61f57dd863ed0,0.052236806857755155,390000,1,4182000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n",
      "206000,0.041084962106102914,0.007403148134837921,0.041084962106102914,0.2215344518651246,13.20388349514563,3.0,69.69466823965065,\"GCA_002754635.1 Plasmodium vivax strain=CMB-1, CMB-1_v2\",/home/irber/sourmash_databases/outputs/sbt/genbank-protozoa-x1e6-k31.sbt.zip,8125e7913e0d0b88deb63c9ad28f827c,0.0037419167332703625,206000,2,3976000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n",
      "138000,0.027522935779816515,0.024722321748477247,0.027522935779816515,0.015637726014008795,1.391304347826087,1.0,0.5702120455914782,\"GCF_013368705.1 Bacteroides vulgatus strain=B33, ASM1336870v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,7d5f4ba1d01c8c3f7a520d19faded7cb,0.012648945921173235,138000,3,3838000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n",
      "338000,0.06741124850418827,0.013789581205311542,0.010769844435580374,0.006515719172503665,1.4814814814814814,1.0,0.738886568268889,\"GCF_003471795.1 Prevotella copri strain=AM16-54, ASM347179v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,0ebd36ff45fc2810808789667f4aad84,0.04337782340862423,54000,4,3784000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n",
      "110000,0.021938571998404467,0.000842978957948319,0.010370961308336658,0.023293696041700604,5.5,2.5,7.417494911978758,\"GCA_000256725.2 Toxoplasma gondii TgCatPRC2 strain=TgCatPRC2, TGCATPRC2 v2\",/home/irber/sourmash_databases/outputs/sbt/genbank-protozoa-x1e6-k31.sbt.zip,2a3b1804cf5ea5fe75dde3e153294548,0.0008909768346023004,52000,5,3732000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000\n"
     ]
    }
   ],
   "source": [
    "!head outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `sourmash tax summarize` (for summarizing metagenomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " **for each gather query, summarize gather results by taxonomic lineage.**\n",
    "\n",
    "There are three possible output formats, `summary`, `lineage_summary`, and `krona`.\n",
    "\n",
    "- `summary` is the default output format. This outputs a `csv` with lineage summarization for each taxonomic rank. This output currently consists of four columns, `query_name,rank,fraction,lineage`, where `fraction` is the  fraction of the query matched to the reported rank and lineage.\n",
    "- `krona` format is a tab-separated list of these results at a specific rank. The first column, `fraction` is the fraction of the query matched to the reported rank and lineage. The remaining columns are `superkingdom`, `phylum`, .. etc down to the rank used for summarization. This output can be used directly for summary visualization.\n",
    "- `lineage_summary` - the lineage summary format is most useful when comparing across metagenomes. each row is a lineage at the desired reporting rank. The columns are each query used for gather, with the fraction match reported for each lineage. This format is commonly used as input for many external multi-sample visualization tools.\n",
    "\n",
    "example `lineage_summary`:\n",
    "\n",
    "    lineage    sample1  sample2 sample3\n",
    "    lin_a     0.4    0.17     0.6\n",
    "    lin_b     0.0    0.0      0.1\n",
    "    lin_c     0.3    0.4      0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "usage:  summarize [-h] [--from-file FILE] [-q] [-o OUTPUT_BASE] -t FILE\n",
      "                  [FILE ...] [--keep-full-identifiers]\n",
      "                  [--keep-identifier-versions] [--fail-on-missing-taxonomy]\n",
      "                  [--output-format {summary,krona,lineage_summary} [{summary,krona,lineage_summary} ...]]\n",
      "                  [-r {species,genus,family,order,class,phylum,superkingdom}]\n",
      "                  [-f]\n",
      "                  [gather_results ...]\n",
      "\n",
      "positional arguments:\n",
      "  gather_results\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  --from-file FILE      input many gather results as a text file, with one\n",
      "                        gather csv per line\n",
      "  -q, --quiet           suppress non-error output\n",
      "  -o OUTPUT_BASE, --output-base OUTPUT_BASE\n",
      "                        base filepath for output file(s) (default stdout)\n",
      "  -t FILE [FILE ...], --taxonomy-csv FILE [FILE ...]\n",
      "                        database lineages csv\n",
      "  --keep-full-identifiers\n",
      "                        do not split identifiers on whitespace\n",
      "  --keep-identifier-versions\n",
      "                        after splitting identifiers, do not remove accession\n",
      "                        versions\n",
      "  --fail-on-missing-taxonomy\n",
      "                        fail quickly if taxonomy is not available for an\n",
      "                        identifier\n",
      "  --output-format {summary,krona,lineage_summary} [{summary,krona,lineage_summary} ...]\n",
      "                        choose output format(s)\n",
      "  -r {species,genus,family,order,class,phylum,superkingdom}, --rank {species,genus,family,order,class,phylum,superkingdom}\n",
      "                        For non-default output formats: Summarize genome\n",
      "                        taxonomy at this rank and above\n",
      "  -f, --force           continue past errors in taxonomy database loading\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax summarize -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `summarize` to aggregate taxonomic information for gather results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 2 lineage assignments.\n",
      "\u001b[KThe following are missing from the taxonomy information: GCA_002754635,GCA_000256725\n",
      "query_name,rank,fraction,lineage\n",
      "HSMA33MX,superkingdom,0.131,d__Bacteria\n",
      "HSMA33MX,phylum,0.073,d__Bacteria;p__Bacteroidota\n",
      "HSMA33MX,phylum,0.058,d__Bacteria;p__Proteobacteria\n",
      "HSMA33MX,class,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia\n",
      "HSMA33MX,class,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria\n",
      "HSMA33MX,order,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales\n",
      "HSMA33MX,order,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales\n",
      "HSMA33MX,family,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae\n",
      "HSMA33MX,family,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae\n",
      "HSMA33MX,genus,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia\n",
      "HSMA33MX,genus,0.057,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella\n",
      "HSMA33MX,genus,0.016,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Phocaeicola\n",
      "HSMA33MX,species,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia;s__Escherichia coli\n",
      "HSMA33MX,species,0.057,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella;s__Prevotella copri\n",
      "HSMA33MX,species,0.016,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Phocaeicola;s__Phocaeicola vulgatus\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax summarize outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv -t gtdb-rs202.taxonomy.v2.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the default `summarize` output current reports the summarization at all ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: `of 6, missed 2 lineage assignments.`\n",
    "\n",
    "This gather used additional databases -- we need to add additional taxonomies to properly summarize!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100   459  100   459    0     0    802      0 --:--:-- --:--:-- --:--:--     0 --:--:-- --:--:--   803\n",
      "100 83.1M  100 83.1M    0     0  13.0M      0  0:00:06  0:00:06 --:--:-- 18.0M  0:00:05  0:00:02 12.9M\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100   459  100   459    0     0   1242      0 --:--:-- --:--:-- --:--:--  1240\n",
      "100 24.2M  100 24.2M    0     0  2858k      0  0:00:08  0:00:08 --:--:-- 5074k\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100   459  100   459    0     0   1249      0 --:--:-- --:--:-- --:--:--  1247\n",
      "100  107k  100  107k    0     0  45242      0  0:00:02  0:00:02 --:--:--  125k\n"
     ]
    }
   ],
   "source": [
    "!curl -L https://osf.io/cbhgd/download -o bacteria_genbank_lineages.csv\n",
    "!curl -L https://osf.io/q8j6d/download -o bacteria_refseq_lineages.csv\n",
    "!curl -L https://osf.io/urtfx/download -o protozoa_genbank_lineages.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accession,taxid,superkingdom,phylum,class,order,family,genus,species,strain\n",
      "GCA_004431415,2762,Eukaryota,,Glaucocystophyceae,,Cyanophoraceae,Cyanophora,Cyanophora paradoxa,\n",
      "GCA_000150955,556484,Eukaryota,Bacillariophyta,Bacillariophyceae,Naviculales,Phaeodactylaceae,Phaeodactylum,Phaeodactylum tricornutum,Phaeodactylum tricornutum CCAP 1055/1\n",
      "GCA_000310025,2880,Eukaryota,,Phaeophyceae,Ectocarpales,Ectocarpaceae,Ectocarpus,Ectocarpus siliculosus,\n",
      "GCA_000194455,2898,Eukaryota,,Cryptophyceae,Cryptomonadales,Cryptomonadaceae,Cryptomonas,Cryptomonas paramecium,\n",
      "GCA_000372725,280463,Eukaryota,Haptista,Haptophyta,Isochrysidales,Noelaerhabdaceae,Emiliania,Emiliania huxleyi,Emiliania huxleyi CCMP1516\n",
      "GCA_001939145,2951,Eukaryota,,Dinophyceae,Suessiales,Symbiodiniaceae,Symbiodinium,Symbiodinium microadriaticum,\n",
      "GCA_900617105,2996,Eukaryota,,Chrysophyceae,Hydrurales,Hydruraceae,Hydrurus,Hydrurus foetidus,\n",
      "GCA_001638955,158060,Eukaryota,Euglenozoa,Euglenida,Euglenales,Euglenaceae,Euglena,Euglena gracilis,\n",
      "GCA_900893395,3039,Eukaryota,Euglenozoa,Euglenida,Euglenales,Euglenaceae,Euglena,Euglena gracilis,\n",
      "accession,taxid,superkingdom,phylum,class,order,family,genus,species,strain\n",
      "GCF_000010525,438753,Bacteria,Proteobacteria,Alphaproteobacteria,Rhizobiales,Xanthobacteraceae,Azorhizobium,Azorhizobium caulinodans,Azorhizobium caulinodans ORS 571\n",
      "GCF_000007365,198804,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. Sg (Schizaphis graminum)\n",
      "GCF_000007725,224915,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. Bp (Baizongia pistaciae)\n",
      "GCF_000009605,107806,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. APS (Acyrthosiphon pisum)\n",
      "GCF_000021065,561501,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. Tuc7 (Acyrthosiphon pisum)\n",
      "GCF_000021085,563178,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. 5A (Acyrthosiphon pisum)\n",
      "GCF_000090965,372461,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola BCc\n",
      "GCF_000174075,593275,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. LSR1 (Acyrthosiphon pisum)\n",
      "GCF_000183285,713601,Bacteria,Proteobacteria,Gammaproteobacteria,Enterobacterales,Erwiniaceae,Buchnera,Buchnera aphidicola,Buchnera aphidicola str. JF99 (Acyrthosiphon pisum)\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head protozoa_genbank_lineages.csv\n",
    "head bacteria_refseq_lineages.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and now run `summarize` with multiple taxonomies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 0 lineage assignments.\n",
      "query_name,rank,fraction,lineage\n",
      "HSMA33MX,superkingdom,0.245,Eukaryota\n",
      "HSMA33MX,superkingdom,0.131,d__Bacteria\n",
      "HSMA33MX,phylum,0.245,Eukaryota;Apicomplexa\n",
      "HSMA33MX,phylum,0.073,d__Bacteria;p__Bacteroidota\n",
      "HSMA33MX,phylum,0.058,d__Bacteria;p__Proteobacteria\n",
      "HSMA33MX,class,0.222,Eukaryota;Apicomplexa;Aconoidasida\n",
      "HSMA33MX,class,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia\n",
      "HSMA33MX,class,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria\n",
      "HSMA33MX,class,0.023,Eukaryota;Apicomplexa;Conoidasida\n",
      "HSMA33MX,order,0.222,Eukaryota;Apicomplexa;Aconoidasida;Haemosporida\n",
      "HSMA33MX,order,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales\n",
      "HSMA33MX,order,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales\n",
      "HSMA33MX,order,0.023,Eukaryota;Apicomplexa;Conoidasida;Eucoccidiorida\n",
      "HSMA33MX,family,0.222,Eukaryota;Apicomplexa;Aconoidasida;Haemosporida;Plasmodiidae\n",
      "HSMA33MX,family,0.073,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae\n",
      "HSMA33MX,family,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae\n",
      "HSMA33MX,family,0.023,Eukaryota;Apicomplexa;Conoidasida;Eucoccidiorida;Sarcocystidae\n",
      "HSMA33MX,genus,0.222,Eukaryota;Apicomplexa;Aconoidasida;Haemosporida;Plasmodiidae;Plasmodium\n",
      "HSMA33MX,genus,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia\n",
      "HSMA33MX,genus,0.057,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella\n",
      "HSMA33MX,genus,0.023,Eukaryota;Apicomplexa;Conoidasida;Eucoccidiorida;Sarcocystidae;Toxoplasma\n",
      "HSMA33MX,genus,0.016,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Phocaeicola\n",
      "HSMA33MX,species,0.222,Eukaryota;Apicomplexa;Aconoidasida;Haemosporida;Plasmodiidae;Plasmodium;Plasmodium vivax\n",
      "HSMA33MX,species,0.058,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia;s__Escherichia coli\n",
      "HSMA33MX,species,0.057,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella;s__Prevotella copri\n",
      "HSMA33MX,species,0.023,Eukaryota;Apicomplexa;Conoidasida;Eucoccidiorida;Sarcocystidae;Toxoplasma;Toxoplasma gondii\n",
      "HSMA33MX,species,0.016,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Phocaeicola;s__Phocaeicola vulgatus\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax summarize outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv --taxonomy-csv bacteria_refseq_lineages.csv protozoa_genbank_lineages.csv gtdb-rs202.taxonomy.v2.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Produce Krona-compatible summary at species-level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 2 lineage assignments.\n",
      "\u001b[KThe following are missing from the taxonomy information: GCA_002754635,GCA_000256725\n"
     ]
    }
   ],
   "source": [
    "!mkdir -p outputs/tax\n",
    "!sourmash tax summarize -r species outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv -t gtdb-rs202.taxonomy.v2.csv --output-format krona -o outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing text.krona.html...\n"
     ]
    }
   ],
   "source": [
    "!ktImportText outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.krona.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `sourmash tax classify` (for classifying genomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sourmash tax classify` - for each gather query, report likely classification based on `gather` matches. By default, classification requires at least 10% of the query to be matched. Thus, if 10% of the query was matched to a species, the species-level classification can be reported. However, if 7% of the query was matched to one species, and an additional 5% matched to a different species in the same genus, the genus-level classification will be reported.\n",
    "\n",
    "Optionally, `classify` can instead report classifications at a desired `rank`, regardless of match threshold. \n",
    "\n",
    "Note that these thresholds and strategies are under active testing.\n",
    "\n",
    "There are two possible output formats, `summary` and `krona`.\n",
    "\n",
    "- `summary` is the default output format. This outputs a `csv` with lineage summarization for each taxonomic rank. This output currently consists of four columns, `query_name,rank,fraction,lineage`, where `fraction` is the  fraction of the query matched to the reported rank and lineage.\n",
    "- `krona` format is a tab-separated list of these results at a specific rank. The first column, `fraction` is the fraction of the query matched to the reported rank and lineage. The remaining columns are `superkingdom`, `phylum`, .. etc down to the rank used for summarization. This output can be used directly for summary visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage:  classify [-h] [-q] -t FILE [FILE ...] [--from-file FILE]\n",
      "                 [-o OUTPUT_BASE]\n",
      "                 [-r {species,genus,family,order,class,phylum,superkingdom}]\n",
      "                 [--keep-full-identifiers] [--keep-identifier-versions]\n",
      "                 [--fail-on-missing-taxonomy]\n",
      "                 [--output-format {summary,krona} [{summary,krona} ...]] [-f]\n",
      "                 [--containment-threshold CONTAINMENT_THRESHOLD]\n",
      "                 [gather_results ...]\n",
      "\n",
      "positional arguments:\n",
      "  gather_results\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -q, --quiet           suppress non-error output\n",
      "  -t FILE [FILE ...], --taxonomy-csv FILE [FILE ...]\n",
      "                        database lineages csv\n",
      "  --from-file FILE      input many gather results as a text file, with one\n",
      "                        gather csv per line\n",
      "  -o OUTPUT_BASE, --output-base OUTPUT_BASE\n",
      "                        base filepath for output file(s) (default stdout)\n",
      "  -r {species,genus,family,order,class,phylum,superkingdom}, --rank {species,genus,family,order,class,phylum,superkingdom}\n",
      "                        Summarize genome taxonomy at this rank and above\n",
      "  --keep-full-identifiers\n",
      "                        do not split identifiers on whitespace\n",
      "  --keep-identifier-versions\n",
      "                        after splitting identifiers, do not remove accession\n",
      "                        versions\n",
      "  --fail-on-missing-taxonomy\n",
      "                        fail quickly if taxonomy is not available for an\n",
      "                        identifier\n",
      "  --output-format {summary,krona} [{summary,krona} ...]\n",
      "                        choose output format(s)\n",
      "  -f, --force           continue past survivable errors in loading taxonomy\n",
      "                        database or gather results\n",
      "  --containment-threshold CONTAINMENT_THRESHOLD\n",
      "                        minimum containment threshold for classification;\n",
      "                        default=0.1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "sourmash tax classify -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 2 lineage assignments.\n",
      "\u001b[KThe following are missing from the taxonomy information: GCA_000256725,GCA_002754635\n",
      "\u001b[Kloaded 0 gather files for classification.\n",
      "query_name,status,rank,fraction,lineage\n",
      "HSMA33MX,match,superkingdom,0.131,d__Bacteria\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax classify outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv -t gtdb-rs202.taxonomy.v2.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 2 lineage assignments.\n",
      "\u001b[KThe following are missing from the taxonomy information: GCA_002754635,GCA_000256725\n",
      "\u001b[KWARNING: classifying at desired rank species does not meet containment threshold 0.1\n",
      "\u001b[Kloaded 0 gather files for classification.\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax classify --rank species outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv -t gtdb-rs202.taxonomy.v2.csv --output-format krona -o outputs/tax.HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.classify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `sourmash tax label` (for labeling gather results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sourmash tax label` - for any gather results, add a column with taxonomic lineage information for each database match. Do not summarize or classify. Note that this is not required for either `summarize` or `classify`.\n",
    "\n",
    "By default, `label` uses the name of each input gather csv to write an updated version with lineages information. For example, labeling `sample1.gather.csv` would produce `sample1.gather.with-lineages.csv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage:  label [-h] [-q] [--from-file FILE] -t FILE [FILE ...] [-o OUTPUT_DIR]\n",
      "              [--keep-full-identifiers] [--keep-identifier-versions]\n",
      "              [--fail-on-missing-taxonomy] [-f]\n",
      "              [gather_results ...]\n",
      "\n",
      "positional arguments:\n",
      "  gather_results\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -q, --quiet           suppress non-error output\n",
      "  --from-file FILE      input many gather results as a text file, with one\n",
      "                        gather csv per line\n",
      "  -t FILE [FILE ...], --taxonomy-csv FILE [FILE ...]\n",
      "                        database lineages csv\n",
      "  -o OUTPUT_DIR, --output-dir OUTPUT_DIR\n",
      "                        directory for output csv(s)\n",
      "  --keep-full-identifiers\n",
      "                        do not split identifiers on whitespace\n",
      "  --keep-identifier-versions\n",
      "                        after splitting identifiers, do not remove accession\n",
      "                        versions\n",
      "  --fail-on-missing-taxonomy\n",
      "                        fail quickly if taxonomy is not available for an\n",
      "                        identifier\n",
      "  -f, --force           continue past errors in file and taxonomy loading\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "sourmash tax label -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[K\n",
      "== This is sourmash version 4.1.1.dev84+g7044ae5. ==\n",
      "\u001b[K== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==\n",
      "\n",
      "\u001b[Kloaded 6 gather results.\n",
      "\u001b[Kof 6, missed 2 lineage assignments.\n",
      "\u001b[KThe following are missing from the taxonomy information: GCA_002754635,GCA_000256725\n"
     ]
    }
   ],
   "source": [
    "!sourmash tax label outputs/gather/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.csv -t gtdb-rs202.taxonomy.v2.csv --output-dir outputs/tax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.krona.tsv\n",
      "outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.krona.txt\n",
      "outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.with-lineages.csv\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "ls outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intersect_bp,f_orig_query,f_match,f_unique_to_query,f_unique_weighted,average_abund,median_abund,std_abund,name,filename,md5,f_match_orig,unique_intersect_bp,gather_result_rank,remaining_bp,query_filename,query_name,query_md5,query_bp,lineage\n",
      "442000,0.08815317112086159,0.08438335242458954,0.08815317112086159,0.05815279361459521,1.6153846153846154,1.0,1.1059438185997785,\"GCF_001881345.1 Escherichia coli strain=SF-596, ASM188134v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,683df1ec13872b4b98d59e98b355b52c,0.042779713511420826,442000,0,4572000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia;s__Escherichia coli\n",
      "390000,0.07778220981252493,0.10416666666666667,0.07778220981252493,0.050496823586903404,1.5897435897435896,1.0,0.8804995294906566,\"GCF_009494285.1 Prevotella copri strain=iAK1218, ASM949428v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,1266c86141e3a5603da61f57dd863ed0,0.052236806857755155,390000,1,4182000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella;s__Prevotella copri\n",
      "206000,0.041084962106102914,0.007403148134837921,0.041084962106102914,0.2215344518651246,13.20388349514563,3.0,69.69466823965065,\"GCA_002754635.1 Plasmodium vivax strain=CMB-1, CMB-1_v2\",/home/irber/sourmash_databases/outputs/sbt/genbank-protozoa-x1e6-k31.sbt.zip,8125e7913e0d0b88deb63c9ad28f827c,0.0037419167332703625,206000,2,3976000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,\n",
      "138000,0.027522935779816515,0.024722321748477247,0.027522935779816515,0.015637726014008795,1.391304347826087,1.0,0.5702120455914782,\"GCF_013368705.1 Bacteroides vulgatus strain=B33, ASM1336870v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,7d5f4ba1d01c8c3f7a520d19faded7cb,0.012648945921173235,138000,3,3838000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Phocaeicola;s__Phocaeicola vulgatus\n",
      "338000,0.06741124850418827,0.013789581205311542,0.010769844435580374,0.006515719172503665,1.4814814814814814,1.0,0.738886568268889,\"GCF_003471795.1 Prevotella copri strain=AM16-54, ASM347179v1\",/group/ctbrowngrp/gtdb/databases/ctb/gtdb-rs202.genomic.k31.sbt.zip,0ebd36ff45fc2810808789667f4aad84,0.04337782340862423,54000,4,3784000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,d__Bacteria;p__Bacteroidota;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Prevotella;s__Prevotella copri\n",
      "110000,0.021938571998404467,0.000842978957948319,0.010370961308336658,0.023293696041700604,5.5,2.5,7.417494911978758,\"GCA_000256725.2 Toxoplasma gondii TgCatPRC2 strain=TgCatPRC2, TGCATPRC2 v2\",/home/irber/sourmash_databases/outputs/sbt/genbank-protozoa-x1e6-k31.sbt.zip,2a3b1804cf5ea5fe75dde3e153294548,0.0008909768346023004,52000,5,3732000,outputs/abundtrim/HSMA33MX.abundtrim.fq.gz,HSMA33MX,9687eeed,5014000,\n"
     ]
    }
   ],
   "source": [
    "!head outputs/tax/HSMA33MX_gather_x_gtdbrs202_genbank_euks_k31.with-lineages.csv"
   ]
  }
 ],
 "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}