docs/023.peroba_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 023. Peroba pipeline\n",
    "The `peroba` suite consists of 3 programs:\n",
    "1. `peroba_database`: collects information from several sources and generates a set of `perobaDB` files\n",
    "2. `peroba_backbone`: selects a set of \"global\" sequences from `perobaDB` to be analysed together with the local ones (NORW). It finds local sequences within the database, but the user should also include other local sequences. \n",
    "3. `peroba_report`: once the user finishes the analysis (i.e. has a phylogenetic tree using suggestions from `peroba_backbone`), this script will estimate ancestral states and generate a PDF report. Importantly it will estimate (impute) the lineages, uk_lineages and phylotypes for \"all\" local sequences that were not estimated by COGUK (except those failing QC, like excess of `N`s...)\n",
    "\n",
    "### how to run the pipeline\n",
    "There are two main usages for the pipeline. If we do not have new local sequences (i.e. that were not incorporated into the COGUK phylogenetics pipeline), then:\n",
    "1. run `peroba_database` whenever the databases upstream (on COGUK) are updated. This usually happens by the beginning of the week (Tue or Wed)\n",
    "2. run `peroba_report` using `perobaDB` global tree.\n",
    "This will generate the report but will not impute lineages (since we are not giving extra information from our tree)\n",
    "\n",
    "However if we have new sequences, we don't need to wait until they are analysed by COGUK:\n",
    "1. Assume we have already ran `peroba_database` so that we have a `perobaDB` for the most recent data available (a few days old, from last Wednesday let's say). \n",
    "2. run `peroba_backbone` adding the local sequences. It will generate alignments and backbone trees to be used in the phylogenetic inference. It generates a NJ tree which can be used as a starting point but it's expected that the user then runs iqtree or raxml by hand. \n",
    "3. run `peroba_report` using the tree estimated in step 2 above.\n",
    "\n",
    "The latter strategy has the advantage that `peroba` is more permissive than COGUK &mdash; that is, NORW sequences that were rejected by the phylogenetic pipeline at COGUK may still be included by `peroba_backbone` and thus have their lineages inferred. \n",
    "\n",
    "### caveats\n",
    "\n",
    "* Please keep in mind that some sequences will never have a UK_lineage (or other classification) because they cannot be included in the phylogenetic analysis (e.g. too many `N`s confuse MAFFT and make the sample become a \"rogue taxon\"). \n",
    "* For more information about te COGUK phylogenetic pipeline, see further below on the notebook.\n",
    "* Both the software and this tutorial are a work-in-progress. The command options and output may not be as shown here.\n",
    "* this is a `bash` notebook (not `python`). Paths etc. won't be the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/notebooks01/0618\n"
     ]
    }
   ],
   "source": [
    "# mkdir 0618 # please do this by hand outside the notebook; (`mkdir`+`cd` in loop is a bad combo)\n",
    "cd 0618\n",
    "pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## peroba_database\n",
    "\n",
    "<div>\n",
    "<img src=\"023_peroba_database.png\" width=\"800\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: peroba_database [options] \n",
      "\n",
      "peroba_database is the script that generates from scratch the integrated data\n",
      "from COGUK and GISAID\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -d, --debug           Print debugging statements\n",
      "  -v, --verbose         Add verbosity\n",
      "  -m csv [csv ...], --metadata csv [csv ...]\n",
      "                        csv/tsv files with metadata information (from COGUK or\n",
      "                        GISAID)\n",
      "  -t treefile, --tree treefile\n",
      "                        single treefile in newick format\n",
      "  -s fas [fas ...], --fasta fas [fas ...]\n",
      "                        fasta files with unaligned genomes from COGUK, GISAID\n",
      "  -a aln [aln ...], --alignment aln [aln ...]\n",
      "                        aligned sequences from previous iteration (speeds up\n",
      "                        calculations)\n",
      "  -i INPUT, --input INPUT\n",
      "                        Directory where input files are. Default: working\n",
      "                        directory\n",
      "  -o OUTPUT, --output OUTPUT\n",
      "                        Output database directory. Default: working directory\n",
      "  -f, --force-pruning   Prune tree, removing leaves absent from metadata and\n",
      "                        sequences. Slow, few benefits?\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_database.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/alignments\u001b[0m\u001b[K/\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/civet\u001b[0m\u001b[K/\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/microreact\u001b[0m\u001b[K/\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/public\u001b[0m\u001b[K/\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/reports\u001b[0m\u001b[K/\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/trees\u001b[0m\u001b[K/\n",
      "\u001b[0m\u001b[01;31m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/0324.nextstrain.aligned.fasta.bz2\u001b[0m\u001b[K\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/0324.nextstrain.tree_raw.nwk\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/0324.nextstrain.tree_raw.nwk.log\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gadm36_GBR_1_sp.rds\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gadm36_GBR_2_sp.rds\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gadm36_GBR_3_sp.rds\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gadm36_GBR.gpkg\n",
      "\u001b[01;31m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gadm36_GBR_shp.zip\u001b[0m\u001b[K\n",
      "\u001b[01;31m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.metadata_2020-06-03_06-11.tsv.gz\u001b[0m\u001b[K\n",
      "\u001b[01;31m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.sequences_2020-06-03_06-11.fasta.xz\u001b[0m\u001b[K\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/readme.txt\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/simplemaps_ukcities.csv\n",
      "/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/simplemaps_worldcities.csv\n",
      "\u001b[01;34m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/UK-postcode-boundaries\u001b[0m\u001b[K/\n",
      "\u001b[01;31m/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/UK-postcode-boundaries.txz\u001b[0m\u001b[K\n"
     ]
    }
   ],
   "source": [
    "# brief look at input files, downloaded from COGUK and GISAID\n",
    "workfolder=\"${HOME}/Academic/Quadram/021.ncov\"\n",
    "\n",
    "ls -d ${workfolder}/01.COGUK/phylogenetics/*  # downloaded from CLIMB server\n",
    "ls -d ${workfolder}/00.data/*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "`peroba_database` can receive several tables (metadata) from both GISAID and COGUK and it will merge them (based on colun called `sequence_name`). The tree must come from COGUK, and the sequences from GISAID and COGUK. \n",
    "\n",
    "The COGUK sequences must be unaligned and from the phylogenetics pipeline (the raw sequences have very non-standard names)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_tree.newick\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/alignments/cog_2020-06-12_all_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/civet/cog_global_2020-06-12_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.metadata_2020-06-03_06-11.tsv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.sequences_2020-06-03_06-11.fasta.xz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/public/cog_2020-06-12.fasta.xz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Found file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200519.peroba/perobaDB.0608.sequences.aln.xz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Starting database '/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618'\n",
      "peroba 2020-06-18 14:33 [INFO] Reading tree file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_tree.newick (only first tree will be used) and checking for duplicate names\n",
      "peroba 2020-06-18 14:33 [INFO] 39904 leaves in treefile id=0\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading metadata file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/alignments/cog_2020-06-12_all_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading metadata file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/civet/cog_global_2020-06-12_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading metadata file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading metadata file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.metadata_2020-06-03_06-11.tsv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Finished reading metadata files. Merged (raw) dataframe has shape (47821, 94) (rows x cols)\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading fasta sequence files (if not set individually below)\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading sequence file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.sequences_2020-06-03_06-11.fasta.xz\n",
      "peroba 2020-06-18 14:33 [INFO] Read 37155 sequences from file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/00.data/gisaid.sequences_2020-06-03_06-11.fasta.xz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Reading sequence file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/public/cog_2020-06-12.fasta.xz\n",
      "peroba 2020-06-18 14:33 [INFO] Read 26243 sequences from file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/phylogenetics/public/cog_2020-06-12.fasta.xz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Database now has 47517 valid sequences\n",
      "perobaDB 2020-06-18 14:33 [INFO] DataSeq:: New dataframe size, after removing missing sequences: (47517, 94) (r x c)\n",
      "perobaDB 2020-06-18 14:33 [INFO] Full merge between sequences, metadata, and tree before saving database\n",
      "perobaDB 2020-06-18 14:33 [INFO] DataSeq:: New dataframe size, after removing missing sequences: (47517, 96) (r x c)\n",
      "perobaDB 2020-06-18 14:33 [INFO] Saving metadata to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.metadata.csv.gz\n",
      "perobaDB 2020-06-18 14:33 [INFO] Generating an HTML report of metadata variables in /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.html\n",
      "Summarize dataset: 100%|███████████| 109/109 [00:38<00:00,  2.86it/s, Completed]\n",
      "Generate report structure: 100%|██████████████████| 1/1 [00:19<00:00, 19.51s/it]\n",
      "Render HTML: 100%|████████████████████████████████| 1/1 [00:05<00:00,  5.12s/it]\n",
      "Export report to file: 100%|██████████████████████| 1/1 [00:00<00:00, 73.66it/s]\n",
      "perobaDB 2020-06-18 14:34 [INFO] Saving raw (metadata) table to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.raw.csv.gz\n",
      "perobaDB 2020-06-18 14:34 [INFO] Saving tree to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.tree.nhx\n",
      "perobaDB 2020-06-18 14:34 [INFO] Saving unaligned sequences to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.sequences.fasta.xz\n",
      "perobaDB 2020-06-18 14:37 [INFO] Reading Alignment file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200519.peroba/perobaDB.0608.sequences.aln.xz\n",
      "peroba 2020-06-18 14:37 [INFO] Read 45115 sequences from file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200519.peroba/perobaDB.0608.sequences.aln.xz\n",
      "perobaDB 2020-06-18 14:38 [INFO] From 47517 sequences, 45102 were found in alignment (originally with 45115 sequences)\n",
      "perobaDB 2020-06-18 14:38 [INFO] Aligning remaining sequences with mafft\n",
      "OS = linux\n",
      "The number of physical cores =  6\n",
      "nadd =  1000\n",
      "npair =  1000\n",
      "nseq =  1001\n",
      "nlen =  29906\n",
      "use ktuples, size=6!\n",
      "nadd = 1000\n",
      "ppenalty_ex = -10\n",
      "nthread = 6\n",
      "blosum 62 / kimura 200\n",
      "sueff_global = 0.100000\n",
      "norg = 1\n",
      "njobc = 2\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 1029173 ambiguous characters\n",
      "    1 / 1 (thread    0)\n",
      "done.\n",
      "\n",
      "fTEP 499 / 1000 (thread 4)                    \n",
      "STEP 500 / 1000 (thread 3)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 600 / 1000 (thread 5)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 700 / 1000 (thread 3)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 800 / 1000 (thread 0)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 900 / 1000 (thread 5)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "\n",
      "Combining ..\n",
      "   done.                      \n",
      "\n",
      "   done.                      \n",
      "\n",
      "addsingle (nuc) Version 7.467\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "\n",
      "To keep the alignment length, 55 letters were DELETED.\n",
      "To know the positions of deleted letters, rerun the same command with the --mapout option.\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-fragment (Not tested.)\n",
      " ?\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "peroba 2020-06-18 14:38 [INFO] First 1000 sequences aligned\n",
      "OS = linux\n",
      "The number of physical cores =  6\n",
      "nadd =  1000\n",
      "npair =  1000\n",
      "nseq =  1001\n",
      "nlen =  29905\n",
      "use ktuples, size=6!\n",
      "nadd = 1000\n",
      "ppenalty_ex = -10\n",
      "nthread = 6\n",
      "blosum 62 / kimura 200\n",
      "sueff_global = 0.100000\n",
      "norg = 1\n",
      "njobc = 2\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 1022373 ambiguous characters\n",
      "    1 / 1 (thread    0)\n",
      "done.\n",
      "\n",
      "fTEP 499 / 1000 (thread 2)                    \n",
      "STEP 500 / 1000 (thread 1)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 600 / 1000 (thread 5)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 700 / 1000 (thread 0)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 800 / 1000 (thread 5)                    \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "STEP 900 / 1000 (thread 0)                    \n",
      "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\n",
      "\n",
      "Combining ..\n",
      "   done.                      \n",
      "\n",
      "   done.                      \n",
      "\n",
      "addsingle (nuc) Version 7.467\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "\n",
      "To keep the alignment length, 45 letters were DELETED.\n",
      "To know the positions of deleted letters, rerun the same command with the --mapout option.\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-fragment (Not tested.)\n",
      " ?\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "peroba 2020-06-18 14:38 [INFO] First 2000 sequences aligned\n",
      "OS = linux\n",
      "The number of physical cores =  6\n",
      "nadd =  415\n",
      "npair =  415\n",
      "nseq =  416\n",
      "nlen =  29906\n",
      "use ktuples, size=6!\n",
      "nadd = 415\n",
      "ppenalty_ex = -10\n",
      "nthread = 6\n",
      "blosum 62 / kimura 200\n",
      "sueff_global = 0.100000\n",
      "norg = 1\n",
      "njobc = 2\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 443180 ambiguous characters\n",
      "    1 / 1 (thread    0)\n",
      "done.\n",
      "\n",
      "fTEP 414 / 415 (thread 1)                    \n",
      "\n",
      "Combining ..\n",
      "   done.                      \n",
      "\n",
      "   done.                      \n",
      "\n",
      "addsingle (nuc) Version 7.467\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "\n",
      "To keep the alignment length, 23 letters were DELETED.\n",
      "To know the positions of deleted letters, rerun the same command with the --mapout option.\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-fragment (Not tested.)\n",
      " ?\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "peroba 2020-06-18 14:38 [INFO] First 2415 sequences aligned\n",
      "perobaDB 2020-06-18 14:38 [INFO] Saving alignment to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.sequences.aln.xz\n",
      "perobaDB 2020-06-18 14:41 [INFO] Finished saving alignment\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_database.py -v \\\n",
    " -m ${workfolder}/01.COGUK/phylogenetics/alignments/cog_2020-06-12_all_metadata.csv.gz \\\n",
    "    ${workfolder}/01.COGUK/phylogenetics/civet/cog_global_2020-06-12_metadata.csv.gz  \\\n",
    "    ${workfolder}/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_metadata.csv.gz \\\n",
    "    ${workfolder}/00.data/gisaid.metadata_2020-06-03_06-11.tsv.gz \\\n",
    " -t ${workfolder}/01.COGUK/phylogenetics/trees/cog_global_2020-06-12_tree.newick \\\n",
    " -s ${workfolder}/00.data/gisaid.sequences_2020-06-03_06-11.fasta.xz \\\n",
    "    ${workfolder}/01.COGUK/phylogenetics/public/cog_2020-06-12.fasta.xz \\\n",
    " -a ${workfolder}/200519.peroba/perobaDB.0608.sequences.aln.xz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "The output of `peroba_database` is what I call `perobaDB`, a set of files with a day timestamp:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "perobaDB.0618.html             \u001b[0m\u001b[01;31mperobaDB.0618.sequences.aln.xz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.metadata.csv.gz\u001b[0m  \u001b[01;31mperobaDB.0618.sequences.fasta.xz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.raw.csv.gz\u001b[0m       perobaDB.0618.tree.nhx\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## peroba_backbone\n",
    "\n",
    "A common pitfall here is that the local sequences must be named `NORW-Exxxx` (where `Exxxx` is the sample hex code), which is not their original name when sequenced. \n",
    "\n",
    "<div>\n",
    "<img src=\"023_peroba_backbone.png\" width=\"800\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: peroba_backbone <perobaDB> [options]\n",
      "\n",
      "peroba_backbone is the script that generates a global backbone data set\n",
      "(COGUK+GISAID) given a local one (NORW). It depends on the prefix for a\n",
      "perobaDB set of files (from `peroba_database`), like \"perobaDB.0519\". It's\n",
      "recommended that you also local sequences, even without CSV metadata. You can\n",
      "furthermore add a newick file with extra trees (the tree from previous run is\n",
      "a good choice).\n",
      "\n",
      "positional arguments:\n",
      "  perobaDB\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -d, --debug           Print debugging statements\n",
      "  -v, --verbose         Add verbosity\n",
      "  -i INPUT, --input INPUT\n",
      "                        Directory where perobaDB files are. Default: working\n",
      "                        directory\n",
      "  -c csv, --csv csv     csv table with metadata from NORW\n",
      "  -s fasta.bz2, --sequences fasta.bz2\n",
      "                        extra sequences from NORW\n",
      "  -t , --trees          file with (user-defined) trees in newick format to\n",
      "                        help produce backbone\n",
      "  -o OUTPUT, --output OUTPUT\n",
      "                        Output database directory. Default: working directory\n",
      "  -r, --replace         replace database sequence with local version\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_backbone.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading metadata, sequences, and tree from peroba_database\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading database metadata from '/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.metadata.csv.gz'\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading database tree from '/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.tree.nhx'\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading database sequences from '/usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.sequences.aln.xz'\n",
      "peroba 2020-06-18 15:02 [INFO] Read 47517 sequences from file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/perobaDB.0618.sequences.aln.xz\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Finished loading the database; dataframe has dimensions (47517, 96) and it's assumed we have the same number of sequences; the tree may be smaller\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading fasta file with sequences from NORW\n",
      "peroba 2020-06-18 15:02 [INFO] Read 115 sequences from file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/01.COGUK/NORW/NORW-20200610.fa.xz\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Reading file with current trees and checking for duplicate names\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] 6155 leaves in treefile 0\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] 5173 leaves in treefile 1\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Imported 47517 rows from database\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Updating sequence names if they are on global database\n",
      "peroba_backbone 2020-06-18 15:02 [WARNING] Sequences not found in global database will be added by hand:\n",
      "NORW-EB7F7   NORW-EB1BD   NORW-EB733   NORW-EB6DC   NORW-EB092   NORW-ECC60   NORW-ECDE5   NORW-EBF68   NORW-ED713   NORW-EBF77   NORW-ED740   NORW-EB1AE   NORW-EBFB3   NORW-EC028   NORW-ECF58   NORW-EA5B2   NORW-EB126   NORW-EB153   NORW-EBFE0   NORW-EB180   NORW-ED625   NORW-EBF95   NORW-ED28B   NORW-EBF2C   NORW-ED494   NORW-ECE97   NORW-EBF3B   NORW-EC00A   NORW-EB6EB   NORW-EB117   NORW-EA6FB   NORW-EBFA4   NORW-EB6FA   NORW-EB715   NORW-EB7E8   NORW-ED7AA\n",
      "\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] 0 long-named sequences and 79 short-named sequences found on database. 36 new sequences (not in database).\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Now I'll check for obvious duplicate sequences, with identical names (only first will be used)\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] 0 long-named sequences and 79 short-named sequences found on database. 36 new sequences (not in database).\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Checking for non-obvious duplicate sequences, where both short- and long-named versions appear\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Will update frequency info and align 36 sequences\n",
      "OS = linux\n",
      "The number of physical cores =  6\n",
      "nadd =  36\n",
      "npair =  36\n",
      "nseq =  37\n",
      "nlen =  29903\n",
      "use ktuples, size=6!\n",
      "nadd = 36\n",
      "ppenalty_ex = -10\n",
      "nthread = 6\n",
      "blosum 62 / kimura 200\n",
      "sueff_global = 0.100000\n",
      "norg = 1\n",
      "njobc = 2\n",
      "generating a scoring matrix for nucleotide (dist=200) ... done\n",
      "\n",
      "\n",
      "Making a distance matrix ..\n",
      "\n",
      "There are 59215 ambiguous characters\n",
      "    1 / 1 (thread    0)\n",
      "done.\n",
      "\n",
      "fTEP 35 / 36 (thread 4)                    \n",
      "\n",
      "Combining ..\n",
      "   done.                      \n",
      "\n",
      "   done.                      \n",
      "\n",
      "addsingle (nuc) Version 7.467\n",
      "alg=A, model=DNA200 (2), 1.53 (4.59), -0.00 (-0.00), noshift, amax=0.0\n",
      "6 thread(s)\n",
      "\n",
      "\n",
      "To keep the alignment length, 1 letters were DELETED.\n",
      "To know the positions of deleted letters, rerun the same command with the --mapout option.\n",
      "\n",
      "Strategy:\n",
      " FFT-NS-fragment (Not tested.)\n",
      " ?\n",
      "\n",
      "If unsure which option to use, try 'mafft --auto input > output'.\n",
      "For more information, see 'mafft --help', 'mafft --man' and the mafft page.\n",
      "\n",
      "The default gap scoring scheme has been changed in version 7.110 (2013 Oct).\n",
      "It tends to insert more gaps into gap-rich regions than previous versions.\n",
      "To disable this change, add the --leavegappyregion option.\n",
      "\n",
      "peroba 2020-06-18 15:02 [INFO] First 36 sequences aligned\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Finished aligning 36 local sequences\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Trimming genomes from site 265 to site 29675\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Ordering metadata after categorisation\n",
      "peroba_backbone 2020-06-18 15:02 [INFO] Finding SNPs to speed up calculations\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] In total 23014 SNPs were found (i,e, alignment size)\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] Splitting data into global and local (NORW)\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] split data into 46523 global and 1030 local sequences\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] Storing 2 trees with leaves mapped to sequence names\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] In tree id 0 the number of leaves to be pruned is 1\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] In tree id 1 the number of leaves to be pruned is 1\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] Removing duplicates (identical sequences)\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] Number of unique sequences: 37425\n",
      "peroba_backbone 2020-06-18 15:03 [INFO] removing global sequences after global metadata update\n",
      "peroba_backbone 2020-06-18 15:04 [WARNING] In tree id 0, 8448 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:04 [WARNING] In tree id 1, 22 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:04 [WARNING] In tree id 2, 20 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:04 [INFO] Remove global sequences with proportion of ACGT less than  0.75 or proportion of N higher than 0.1\n",
      "peroba_backbone 2020-06-18 15:04 [INFO] removing global sequences after global metadata update\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] Also removing local sequences with proportion of ACGT less than 0.5 or proportion of N higher than 0.3\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] removing local sequences after metadata update\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] After removal, global data has 34723 samples and local data has 979.\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] Subsampling redundant global samples (i.e. those fom same lineage etc.)\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] After subsampling, global metadata has 18858 samples\n",
      "peroba_backbone 2020-06-18 15:05 [INFO] removing global sequences after global metadata update\n",
      "peroba_backbone 2020-06-18 15:06 [WARNING] In tree id 0, 12322 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:06 [WARNING] In tree id 1, 547 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:06 [WARNING] In tree id 2, 482 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:06 [INFO] Finding neighbours to local sequences, using a distance of 4 to 2000 segments\n",
      "peroba 2020-06-18 15:06 [INFO] Creating a hashed genome with blocks of 11 bases\n",
      "peroba 2020-06-18 15:06 [INFO] And finding neighbours with distance smaller than 0.0020499999999999997\n",
      "peroba_backbone 2020-06-18 15:07 [INFO] Found 3544 neighbours; now will find their 25 closest neighbours on 4000 segments\n",
      "peroba 2020-06-18 15:07 [INFO] Creating a hashed genome with blocks of 5 bases\n",
      "peroba 2020-06-18 15:08 [INFO] And finding 25 closest neighbours\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Found 4728 neighbours\n",
      "peroba_backbone 2020-06-18 15:11 [WARNING] In tree id 0, 13828 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:11 [WARNING] In tree id 1, 386 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:11 [WARNING] In tree id 2, 241 leaves were pruned (absent from metadata or list)\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Saving sequences to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.coguk.aln.xz\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Finished saving alignment\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Finished saving global tree(s) to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.coguk.trees.nhx\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Saving user-defined trees with added sequences\n",
      "peroba_backbone 2020-06-18 15:11 [INFO] Saving sequences to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.user.aln.xz\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peroba_backbone 2020-06-18 15:12 [INFO] Finished saving alignment\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] In total, 6662 sequences are part of the user-defined data set (user trees plus sequences found here)\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Estimating NJ tree with all user-defined sequences (will be the first tree in file)\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Finished saving user-defined trees to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.user.trees.nhx\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Saving data set with all local and selected global sequences\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Saving sequences to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.norw-coguk.aln.xz\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Finished saving alignment\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Total of 5707 sequences will form the local+global data set\n",
      "peroba_backbone 2020-06-18 15:12 [INFO] Estimating NJ tree with local+global sequences\n",
      "peroba_backbone 2020-06-18 15:13 [INFO] Finished saving local+global tree to /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.norw-coguk.trees.nhx\n",
      "peroba_backbone 2020-06-18 15:13 [INFO] Saving sequences to file /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/peroba_backbone.0618_1502.norw.aln.xz\n",
      "peroba_backbone 2020-06-18 15:13 [INFO] Finished saving alignment\n",
      "Finished. The output files are described below, where 'global' means COGUK and GISAID data which were \n",
      "\u001b[0;0;31;1mnot\u001b[0m generated in NORW. Those, together with the extra sequences are being called'local'.\n",
      "Files produced:\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.coguk.aln.xz\u001b[0m\n",
      " Aligned sequences selected from global data\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.coguk.trees.nhx\u001b[0m\n",
      " Pruned trees, with leaves only from global data (please refer to 'user' trees below for phylo inference)\n",
      " First tree is full tree from COGUK, remaining trees are provided from user, if any\n",
      "\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.user.aln.xz\u001b[0m\n",
      " alignment with sequences from all user-defined trees (that I have access to), as well as all \n",
      " global and local sequences.\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.user.trees.nhx\u001b[0m\n",
      " User-defined trees with subset of leaves for which I've found sequences\n",
      " where the first tree is a NJ estimate from all 'user' sequences above\n",
      "\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.norw-coguk.aln.xz\u001b[0m\n",
      " alignment with (all) local and (selected) global sequences. This is main output of program\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.norw-coguk.trees.nhx\u001b[0m\n",
      " NJ estimate of all global and local sequences (useful for phylo inference)\n",
      "\n",
      "\u001b[0;0;33;1mperoba_backbone.0618_1502.norw.aln.xz\u001b[0m\n",
      " alignment with all local sequences only\n",
      "\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_backbone.py -v perobaDB.0618 \\\n",
    "-s ${workfolder}/01.COGUK/NORW/NORW-20200610.fa.xz \\\n",
    "-t ${workfolder}/200519.peroba/peroba_backbone.0612_1119.user.trees.nhx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mperoba_backbone.0618_1502.coguk.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.coguk.trees.nhx\n",
      "\u001b[01;31mperoba_backbone.0618_1502.norw.aln.xz\u001b[0m\n",
      "\u001b[01;31mperoba_backbone.0618_1502.norw-coguk.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.norw-coguk.trees.nhx\n",
      "\u001b[01;31mperoba_backbone.0618_1502.user.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.user.trees.nhx\n",
      "perobaDB.0618.html\n",
      "\u001b[01;31mperobaDB.0618.metadata.csv.gz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.raw.csv.gz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.sequences.aln.xz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.sequences.fasta.xz\u001b[0m\n",
      "perobaDB.0618.tree.nhx\n"
     ]
    }
   ],
   "source": [
    "ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## peroba_report\n",
    "\n",
    "<div>\n",
    "<img src=\"023_peroba_report.png\" width=\"800\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: peroba_report [options]\n",
      "\n",
      "peroba_report is the script that generates a PDF report given a tree and\n",
      "metadata for new genomes\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help            show this help message and exit\n",
      "  -d, --debug           Print debugging statements\n",
      "  -v, --verbose         Add verbosity\n",
      "  -m csv.gz, --metadata csv.gz\n",
      "                        metadata file formated by peroba\n",
      "  -c csv, --csv csv     csv table with metadata from NORW\n",
      "  -t treefile, --tree treefile\n",
      "                        single treefile in newick format\n",
      "  -p PREFIX, --prefix PREFIX\n",
      "                        Date to be added to report\n",
      "  -i INPUT, --input INPUT\n",
      "                        Input directory with pandoc templates, if not using\n",
      "                        default\n",
      "  -o OUTPUT, --output OUTPUT\n",
      "                        Output database directory. Default: working directory\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_report.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peroba_report 2020-06-18 15:15 [INFO] Reading metadata (previously prepared by peroba)\n",
      "peroba_report 2020-06-18 15:15 [INFO] Reading CSV file with metadata from NORW\n",
      "peroba_report 2020-06-18 15:15 [INFO] Reading tree file and checking if there are duplicate names\n",
      "peroba_report 2020-06-18 15:15 [INFO] 5707 leaves in treefile and metadata with shape (47821, 96)\n",
      "peroba 2020-06-18 15:15 [INFO] Merging global metadata (COGUK+GISAID usually) with local one (NORW)\n",
      "peroba_report 2020-06-18 15:15 [WARNING] 513\tsamples in local csv metadata were not found on tree (excluded from analysis due to low quality?)\n",
      "peroba_report 2020-06-18 15:15 [INFO] 972\tsamples from local csv metadata were found on tree and will form the basis of the analysis\n",
      "peroba_report 2020-06-18 15:15 [WARNING] Removing 6 unmapped leaves from tree dictionary\n",
      "peroba_report 2020-06-18 15:15 [INFO] Will now prune these leaves from tree.\n",
      "peroba_report 2020-06-18 15:16 [INFO] Comparing tree leaves to check for duplicates\n",
      "peroba_report 2020-06-18 15:16 [INFO] Merging metadata files by sequence name, after renaming local sequences when possible\n",
      "peroba 2020-06-18 15:16 [INFO] Start estimating ancestral states by DOWNPASS for locality and ACCTRAN for others\n",
      "peroba 2020-06-18 15:16 [INFO] Finished estimating ancestral state for 'locality', which defines clusters\n",
      "peroba 2020-06-18 15:16 [INFO] Follow-up columns found in CSV: uk_lineage lineage phylotype\n",
      "peroba 2020-06-18 15:16 [INFO] Will now estimate ancestral states for uk_lineage lineage phylotype\n",
      "peroba_report 2020-06-18 15:17 [INFO] Decorating tree before plotting all clusters\n",
      "peroba_report 2020-06-18 15:17 [INFO] Entering loop over clusters found; sit tight!\n",
      "peroba_report 2020-06-18 15:22 [INFO] running command:: cd /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618 && pandoc /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618/_pdf_2020-06-18.md -o /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618/report_2020-06-18.pdf --from markdown --template /usr/users/QIB_fr005/deolivl/Academic/Quadram/022.peroba/peroba/data/report/eisvogel --listings\n",
      "peroba_report 2020-06-18 15:22 [INFO] running command:: cd /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618 && pandoc -s /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618/_htm_2020-06-18.md -o /usr/users/QIB_fr005/deolivl/Academic/Quadram/021.ncov/200404.UK/0618/report.0618/report_2020-06-18.html --toc\n"
     ]
    }
   ],
   "source": [
    "python ~/Academic/Quadram/022.peroba/peroba/peroba_report.py -v \\\n",
    "-c ${workfolder}/01.COGUK/NORW/20200617.SARCOV2-Metadata.csv \\\n",
    "-t peroba_backbone.0618_1502.norw-coguk.trees.nhx \\\n",
    "-m perobaDB.0618.raw.csv.gz -o report.0618"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[01;31mperoba_backbone.0618_1502.coguk.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.coguk.trees.nhx\n",
      "\u001b[01;31mperoba_backbone.0618_1502.norw.aln.xz\u001b[0m\n",
      "\u001b[01;31mperoba_backbone.0618_1502.norw-coguk.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.norw-coguk.trees.nhx\n",
      "\u001b[01;31mperoba_backbone.0618_1502.user.aln.xz\u001b[0m\n",
      "peroba_backbone.0618_1502.user.trees.nhx\n",
      "perobaDB.0618.html\n",
      "\u001b[01;31mperobaDB.0618.metadata.csv.gz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.raw.csv.gz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.sequences.aln.xz\u001b[0m\n",
      "\u001b[01;31mperobaDB.0618.sequences.fasta.xz\u001b[0m\n",
      "perobaDB.0618.tree.nhx\n",
      "\n",
      "report.0618:\n",
      "\u001b[01;31mcsv_2020-06-18.csv.gz\u001b[0m  \u001b[01;31mmetadata_2020-06-18.csv.gz\u001b[0m  report_2020-06-18.html\n",
      "\u001b[01;34mfigures\u001b[0m/               pandoc.css                  report_2020-06-18.pdf\n",
      "_htm_2020-06-18.md     _pdf_2020-06-18.md\n"
     ]
    }
   ],
   "source": [
    "ls *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The COGUK pipeline\n",
    "Every Fri at noon the COGUK collects all sequences (with metadata) submitted that week and start processing them &mdash; uploading them to GISAID and feeding it to the phylogenetic pipeline. Every week around Tue/Wed we then have a new files on the CLIMB server. Here is a description of the available files: \n",
    "\n",
    "```bash\n",
    "Phylogenetics pipeline output published to `/cephfs/covid/bham/artifacts/published/latest/phylogenetics/`\n",
    "\n",
    "Reports published to `reports/`\n",
    "\n",
    "Unaligned (deduplicated, clean headers) COG sequences published to `alignments/cog_2020-06-19_all.fasta`\n",
    "\n",
    "Aligned (deduplicated, clean headers) COG sequences published to `alignments/cog_2020-06-19_all_alignment.fasta`\n",
    "Matching metadata published to `alignments/cog_2020-06-19_all_metadata.csv`\n",
    "\n",
    "Filtered, aligned COG sequences published to `alignments/cog_2020-06-19_alignment.fasta`\n",
    "Matching metadata with lineage information published to `alignments/cog_2020-06-19_metadata.csv`\n",
    "\n",
    "Full, annotated tree published to `trees/cog_2020-06-19_tree.nexus`\n",
    "Matching metadata published to `trees/cog_global_2020-06-19_metadata.csv`\n",
    "UK lineage subtrees published in `trees/uk_lineages/`\n",
    "UK lineage timetrees published in `trees/uk_lineages/`\n",
    "\n",
    "Public tree published to `public/cog_global_2020-06-19_tree.newick`\n",
    "Associated unaligned sequences published to `alignments/cog_2020-06-19_all.fasta`\n",
    "Matching metadata with public fields only published to `public/cog_2020-06-19_metadata.csv`\n",
    "\n",
    "Public tree for microreact published to `microreact/cog_global_2020-06-19_tree_public.newick`\n",
    "Public metadata for microreact published to `microreact/cog_global_2020-06-19_metadata_public.csv`\n",
    "Private metadata for microreact published to `microreact/cog_global_2020-06-19_metadata_private.csv`\n",
    "\n",
    "Data for Civet published to `/cephfs/covid/bham/civet-cat/`\n",
    "\n",
    "```\n",
    "The particular files, and sometimes the table columns, change from week to week. As seen above, we do not use them all but I downloaded them anyway. using an LZMA algorithm (like `xz` or `tar -J`) which makes the files a gazillion times smaller. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}