Introduction_to_how_ONT's_medaka_works.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Introduction to how ONT's medaka works",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/epi2me-labs/tutorials/blob/master/Introduction_to_how_ONT's_medaka_works.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y8OphKkrzRXq",
        "colab_type": "text"
      },
      "source": [
        "<h1>How medaka works</h1>\n",
        "\n",
        "The following is a relatively short document describing how Oxford Nanopore Technologies' program for consensus calling of sequencing data, `medaka`, functions internally. We will demonstrate the core functionality required to process alignment data, how it is presented to a recurrent neural network, and how a consensus sequence is formed.\n",
        "\n",
        "The document you are reading is not a static web page, but an interactive environment called a Colab notebook that lets you write and execute code. You can inspect, modify, and run any of the code on this page."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xrD1AW3-0O1w",
        "colab_type": "text"
      },
      "source": [
        "## Installation\n",
        "\n",
        "Before getting started with how `medaka` works, we will install it into the Colab environment. This will enable us to both inspect and run the code. To install the `medaka` run the code cell below by clicking the \"play\" icon to the left or pressing on the cell and typing `<shift>-<enter>`"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PJZWryaduqV5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# setup some prerequisites - medaka model files are stored in git-lfs, which\n",
        "# we need to install\n",
        "%cd /content\n",
        "!wget https://github.com/git-lfs/git-lfs/releases/download/v2.10.0/git-lfs-linux-amd64-v2.10.0.tar.gz\n",
        "!tar -xzvf git-lfs-linux-amd64-v2.10.0.tar.gz\n",
        "!mv git-lfs /bin/\n",
        "!git lfs install\n",
        "\n",
        "# install medaka - we'll do this from source so we can make a couple of\n",
        "# modifications and use some extra bits later on\n",
        "!apt-get install file\n",
        "!rm -rf /content/medaka\n",
        "!git clone https://github.com/nanoporetech/medaka.git\n",
        "%cd /content/medaka\n",
        "!sed -i 's/tensorflow==/tensorflow-gpu==/' requirements.txt\n",
        "!pip install -r requirements.txt\n",
        "!make scripts/mini_align\n",
        "!python setup.py install\n",
        "%cd /content\n",
        "\n",
        "# install minimap2 and pomoxis for later\n",
        "!git clone https://github.com/lh3/minimap2\n",
        "!cd minimap2 && make && cp minimap2 /bin\n",
        "!pip install git+https://github.com/rrwick/Porechop\n",
        "!pip install pomoxis"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HyyI0VZ_zHd8",
        "colab_type": "text"
      },
      "source": [
        "After running the above, it is necessary to select `Runtime>Restart runtime...` from the menu at the top of the page, after which we should be able to import medaka"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LJDS3xRUu9A-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import medaka\n",
        "help(medaka)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DhwqwvbD1e5s",
        "colab_type": "text"
      },
      "source": [
        "## Medaka's input\n",
        "\n",
        "As input the core `medaka` algorithm accepts sequencing reads aligned to an assembly sequence. If you have run the `medaka_consensus` pipeline you will have given as input an assembly sequence and your sequencing data. The pipeline simply runs [`minimap2`](https://https://github.com/lh3/minimap2) to calculate alignments of the reads to the assembly.\n",
        "\n",
        "For the purposes of this demonstration we will download pre-aligned data from an R9.4.1 MinION sequencing run:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xiEyr06dzFmH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!mkdir -p /content/data && cd /content/data/ \\\n",
        "    && wget https://ont-research.s3-eu-west-1.amazonaws.com/datasets/r941_zymo/references.fasta \\\n",
        "    && wget https://ont-research.s3-eu-west-1.amazonaws.com/labs_resources/misc/saureus.bam \\\n",
        "    && wget https://ont-research.s3-eu-west-1.amazonaws.com/labs_resources/misc/saureus.bam.bai \\\n",
        "    && wget https://ont-research.s3-eu-west-1.amazonaws.com/labs_resources/misc/saureus_canu.fasta"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VOcZ6zBD90zz",
        "colab_type": "text"
      },
      "source": [
        "The downloaded `saureus.bam` file contains alignments of sequencing reads to the downloaded `saureus_canu.fasta`. The depth of sequencing has been reduced to around 150-fold coverage of the genome."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gy7qtmCJ-nyH",
        "colab_type": "text"
      },
      "source": [
        "# Diving in: counting bases\n",
        "\n",
        "The first step of `medaka`'s calculation is to parse the alignment data into a base counts table ready for input to the neural network. In this section we explore the functions responsible for doing this, how exactly counting is performed and what the results may represent.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ea9FOGnhBUJO",
        "colab_type": "text"
      },
      "source": [
        "## Pileup interface\n",
        "\n",
        "At the heart of `medaka` resides a straight-forward base-counting procedure. From the alignment data comparing sequencing reads to the reference sequence a pileup is created, much like the display and alignment viewer such as [IGV](https://software.broadinstitute.org/software/igv/) would display.\n",
        "\n",
        "The pileup is summarise by counting the different base types contained within its columns. The function responsible for this counting excercise is called `pileup_counts` in the [`features`](https://github.com/nanoporetech/medaka/blob/d195b9cc1ee7681a121be9fe4fb016a00744ef47/medaka/features.py#L109) module:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RGv0datQAbDH",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 377
        },
        "outputId": "e0b5b0f7-09b2-40e5-dfb7-10208f4fceeb"
      },
      "source": [
        "from medaka.features import pileup_counts\n",
        "help(pileup_counts)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Help on function pileup_counts in module medaka.features:\n",
            "\n",
            "pileup_counts(region, bam, dtype_prefixes=None, region_split=100000, workers=8, tag_name=None, tag_value=None, keep_missing=False, num_qstrat=1, weibull_summation=False)\n",
            "    Create pileup counts feature array for region.\n",
            "    \n",
            "    :param region: `medaka.common.Region` object\n",
            "    :param bam: .bam file with alignments.\n",
            "    :param dtype_prefixes: prefixes for query names which to separate counts.\n",
            "        If `None` (or of length 1), counts are not split.\n",
            "    :param region_split: largest region to process in single thread.\n",
            "    :param workers: worker threads for calculating pileup.\n",
            "    :param tag_name: two letter tag name by which to filter reads.\n",
            "    :param tag_value: integer value of tag for reads to keep.\n",
            "    :param keep_missing: whether to keep reads when tag is missing.\n",
            "    :param num_qstrat: number of layers for qscore stratification.\n",
            "    :param weibull_summation: use a Weibull partial-counts approach,\n",
            "        requires 'WL' and 'WK' float-array tags.\n",
            "    \n",
            "    :returns: pileup counts array, reference positions, insertion postions\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qq3dTUG1Wd-R",
        "colab_type": "text"
      },
      "source": [
        "The `pileup_counts` function above has various arguments, most of which are advanced options and not used within the default operation of medaka. To create a counts matrix we call the function with a Samtools-style region string (`medaka` uses 0-based end exclusive co-ordinates) and a filepath to our alignment file:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pGYE6i8yB85_",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "73d10276-5394-4a88-fb1e-a2e2b676ae4e"
      },
      "source": [
        "from timeit import default_timer as now\n",
        "from medaka.common import Region\n",
        "\n",
        "t0 = now()\n",
        "region = Region.from_string('tig00000061:0-1499707')\n",
        "bam_file = '/content/data/saureus.bam'\n",
        "pileup_data = pileup_counts(region, bam_file)\n",
        "pileup_data = pileup_data[0]  # implementation detail that need not trouble us\n",
        "counts, positions = pileup_data\n",
        "t1 = now()\n",
        "print(\"{:.2f}s to form pileup counts.\".format(t1 - t0))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "10.62s to form pileup counts.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WUaJ09Eq5g7J",
        "colab_type": "text"
      },
      "source": [
        "### The counts matrix\n",
        "\n",
        "The `pileup_counts` function returned two structures. The latter of these is a positions table, this records which pileup columns are reference positions and which are caused by inserted bases in one or more reads:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xlbWIC3iMoPJ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 68
        },
        "outputId": "1371aebb-8d8a-4230-c1b6-da0d772755b8"
      },
      "source": [
        "display(positions)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "array([(      0, 0), (      1, 0), (      2, 0), ..., (1499704, 0),\n",
              "       (1499705, 0), (1499706, 0)],\n",
              "      dtype=[('major', '<i8'), ('minor', '<i8')])"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lwG1D7GkMoYT",
        "colab_type": "text"
      },
      "source": [
        "The field `minor` in the above array indicates reference and insertion columns: it takes a value `0` for a reference position and counts upwards for all following insertion events. The `major` field keeps track of the reference base co-ordinate.\n",
        "\n",
        "The base counts themselves from the alignment pileup are stored separately:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ncxpeIgSDNv4",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 153
        },
        "outputId": "239b4aec-e861-49af-d254-0cf239a38575"
      },
      "source": [
        "display(counts.shape)\n",
        "display(counts)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "(3508694, 10)"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "array([[0, 0, 0, ..., 0, 0, 0],\n",
              "       [0, 0, 0, ..., 0, 0, 0],\n",
              "       [0, 0, 0, ..., 4, 0, 0],\n",
              "       ...,\n",
              "       [0, 0, 0, ..., 2, 0, 0],\n",
              "       [0, 0, 0, ..., 2, 0, 0],\n",
              "       [0, 0, 0, ..., 0, 0, 0]], dtype=uint64)"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cJQvo-n9_3H3",
        "colab_type": "text"
      },
      "source": [
        "The matrix is of shape (# pileup columns, 10), each row of the matrix corresponds to the counts of bases and gaps in the pileup columns (yes, the rows and columns get confusing). There are 10 entries one each for the fours base types and gap, multiplied by two as reads on the forward and reverse strand are counted separately. The ordering of the entries is given by: "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YsIWk-gl3YVZ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "3ab75675-cd9f-4492-e128-25a6ff2c70fb"
      },
      "source": [
        "from medaka.features import libmedaka\n",
        "ffi, lib = libmedaka.ffi, libmedaka.lib\n",
        "plp_bases = lib.plp_bases\n",
        "codes = ffi.string(plp_bases).decode()\n",
        "display(','.join(codes))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "'a,c,g,t,A,C,G,T,d,D'"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SoYRd_TZJiWP",
        "colab_type": "text"
      },
      "source": [
        "in which lower-case letters denote reverse strand counts (upper case, forward) and 'd' and 'D' count deletions. A point of note is that this counting strategy it itself makes a distinction between bases which are deleted in reads with respect to reference sequence and bases which are deleted in reads with respect to other reads (the bases in the other reads being insertions with respect to the reference). Previous versions of `medaka` have performed a symmetrization here: by adding in deletion counts for all read that span a pileup column, whether that pileup column is a reference position (`minor=0`) or an insertion column (`minor>0`).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a2xFr4JDFo5O",
        "colab_type": "text"
      },
      "source": [
        "### An aside on performance\n",
        "\n",
        "Although medaka is mainly written in python, this first calculation is performed in C for speed; `pileup_counts` defers to a `C` implementation of the base counting which makes use of the pileup API from [`htslib`](https://github.com/samtools/htslib). This `C` function is 1-2 orders of magnitude faster than a previous python implementation, using [`pysam`](https://github.com/pysam-developers/pysam). Nevertheless this seemingly trivial step is often the performance bottleneck in `medaka`, particularly when using GPUs to run the neural network calculations. Further it is within the `htslib` function [`resolve_cigar2`](https://github.com/samtools/htslib/blob/9279d76e1186d7155ceea9db9db8c9298f6139bd/sam.c#L3956) that the code is bogged down:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "izsOqORK54Hs",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 527
        },
        "outputId": "95f51c0b-2a2b-43b0-df1d-277d1ef1829e"
      },
      "source": [
        "%cd /content/medaka/\n",
        "print(\"Compiling standalone pileup counts program.\")\n",
        "!rm -rf pileup\n",
        "!sed -i \"s/gcc -pthread/gcc -pg -pthread/\" Makefile\n",
        "!make pileup\n",
        "print(\"\\nRunning pileup...\")\n",
        "!time ./pileup /content/data/saureus.bam tig00000061:1-1499707 2>/dev/null > pileup.txt\n",
        "!gprof pileup gmon.out | head -n 17"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/content/medaka\n",
            "Compiling standalone pileup counts program.\n",
            "gcc -pg -pthread  -g -Wall -fstack-protector-strong -D_FORTIFY_SOURCE=2 -fPIC -std=c99 -msse3 -O3 \\\n",
            "\t-Isrc -Isubmodules/samtools-1.9/htslib-1.9 \\\n",
            "\tsrc/medaka_common.c src/medaka_counts.c src/medaka_bamiter.c libhts.a \\\n",
            "\t-lm -lz -llzma -lbz2 -lpthread -lcurl -lcrypto \\\n",
            "\t-o pileup -std=c99 -msse3 -O3\n",
            "\n",
            "Running pileup...\n",
            "\n",
            "real\t0m28.472s\n",
            "user\t0m26.345s\n",
            "sys\t0m1.940s\n",
            "Flat profile:\n",
            "\n",
            "Each sample counts as 0.01 seconds.\n",
            "  %   cumulative   self              self     total           \n",
            " time   seconds   seconds    calls  Ts/call  Ts/call  name    \n",
            " 49.58      5.85     5.85                             resolve_cigar2\n",
            " 27.97      9.15     3.30                             calculate_pileup\n",
            " 16.02     11.04     1.89                             bam_plp_next\n",
            "  2.88     11.38     0.34                             print_pileup_data\n",
            "  1.10     11.51     0.13                             bam_cigar2rqlens\n",
            "  1.10     11.64     0.13                             bam_cigar2rlen\n",
            "  0.59     11.71     0.07                             bam_mplp_auto\n",
            "  0.25     11.74     0.03                             bam_plp_auto\n",
            "  0.25     11.77     0.03                             kh_init_olap_hash\n",
            "  0.25     11.80     0.03                             mp_free\n",
            "  0.00     11.80     0.00    58568     0.00     0.00  read_bam\n",
            "  0.00     11.80     0.00        7     0.00     0.00  xalloc\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4lP9SAHMHjeP",
        "colab_type": "text"
      },
      "source": [
        "Because of this `medaka` uses multiple worker threads to perform the base counting: the calculation is split into 1 Mbase shards and then further subdivided into 100kbase chunks with the chunks being reassembled for each shard. In conclusion, note the time taken for the simple pure `C` program `pileup` above and that for running the full `medaka consensus` program below:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tNg8ia2eIsxr",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 462
        },
        "outputId": "8954ec9f-65b7-4cc1-e213-37c73eb41df9"
      },
      "source": [
        "!rm -rf /content/data/saureus.hdf\n",
        "!time medaka consensus /content/data/saureus.bam /content/data/saureus.hdf --region tig00000061 --batch_size 50"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[23:06:01 - Predict] Processing region(s): tig00000061:0-1499707\n",
            "[23:06:01 - Predict] Setting tensorflow threads to 1.\n",
            "[23:06:01 - Predict] Found a GPU.\n",
            "[23:06:01 - Predict] If cuDNN errors are observed, try setting the environment variable `TF_FORCE_GPU_ALLOW_GROWTH=true`. To explicitely disable use of cuDNN use the commandline option `--disable_cudnn. If OOM (out of memory) errors are found please reduce batch size.\n",
            "[23:06:01 - Predict] Processing 2 long region(s) with batching.\n",
            "[23:06:01 - Predict] Using model: /usr/local/lib/python3.6/dist-packages/medaka-0.11.5-py3.6-linux-x86_64.egg/medaka/data/r941_min_high_g344_model.hdf5.\n",
            "[23:06:01 - ModelLoad] Building model with cudnn optimization: True\n",
            "[23:06:02 - DLoader] Initializing data loader\n",
            "[23:06:02 - PWorker] Running inference for 1.5M draft bases.\n",
            "[23:06:02 - Sampler] Initializing sampler for consensus of region tig00000061:0-1000000.\n",
            "[23:06:02 - Sampler] Initializing sampler for consensus of region tig00000061:999000-1499707.\n",
            "[23:06:12 - Feature] Processed tig00000061:999000.0-1499706.0 (median depth 150.0)\n",
            "[23:06:12 - Sampler] Took 9.58s to make features.\n",
            "[23:06:12 - PWorker] Samples in cache: 81.\n",
            "[23:06:14 - Feature] Processed tig00000061:0.0-999999.3 (median depth 150.0)\n",
            "[23:06:14 - Sampler] Took 11.76s to make features.\n",
            "[23:06:15 - PWorker] 12.8% Done (0.2/1.5 Mbases) in 12.7s\n",
            "[23:06:17 - PWorker] Samples in cache: 150.\n",
            "[23:06:23 - PWorker] Samples in cache: 0.\n",
            "[23:06:24 - PWorker] All done, 0 remainder regions.\n",
            "[23:06:24 - Predict] Finished processing all regions.\n",
            "\n",
            "real\t0m25.525s\n",
            "user\t0m34.076s\n",
            "sys\t0m3.456s\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IuaZe5SyYAMr",
        "colab_type": "text"
      },
      "source": [
        "On the standard Colab environment with two CPU cores and an NVIDIA P4 GPU, we see that wallclock time for `medaka consensus` is less that the single-threaded `pileup` program. That this speed can be achieved is down to the asynchronous multi-threaded queuing that `medaka` implements for the generation of the counts matrices and the raw power of GPUs to process the data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wGACVcMQVukJ",
        "colab_type": "text"
      },
      "source": [
        "## Normalization\n",
        "\n",
        "After obtained the base-counts matrix produced in the section above `medaka` performs a normalization of the counts. Across the pileup columns, all count vectors with equal corresponding `major` position index are normalized by the total count for the column with `minor=0` (the reference position). This choice of normalization accounts for the lack of symmetry described above, and that whilst consensus insertions are typically rare, isolated insertions may still occur within any one read spanning two reference positions. There are on average up to three pileup columns for every input reference position.\n",
        "\n",
        "Ordinarily this normalization is performed in a [post-processing](https://github.com/nanoporetech/medaka/blob/d195b9cc1ee7681a121be9fe4fb016a00744ef47/medaka/features.py#L372) method of the `CountsFeatureEncoder` class, for the purposes of exposition the operation under normal behaviour is:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oGd1JXR4YemM",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 136
        },
        "outputId": "d91f690c-584c-40a2-d56b-5a4f152501a2"
      },
      "source": [
        "import numpy as np\n",
        "minor_inds = np.where(positions['minor'] > 0)\n",
        "major_pos_at_minor_inds = positions['major'][minor_inds]\n",
        "major_ind_at_minor_inds = np.searchsorted(\n",
        "    positions['major'], major_pos_at_minor_inds, side='left')\n",
        "\n",
        "depth = np.sum(counts, axis=1)\n",
        "depth[minor_inds] = depth[major_ind_at_minor_inds]\n",
        "\n",
        "feature_array = counts / np.maximum(1, depth).reshape((-1, 1))\n",
        "display(feature_array)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "array([[0., 0., 0., ..., 0., 0., 0.],\n",
              "       [0., 0., 0., ..., 0., 0., 0.],\n",
              "       [0., 0., 0., ..., 1., 0., 0.],\n",
              "       ...,\n",
              "       [0., 0., 0., ..., 1., 0., 0.],\n",
              "       [0., 0., 0., ..., 1., 0., 0.],\n",
              "       [0., 0., 0., ..., 0., 0., 0.]])"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4YjHmTTxaEnt",
        "colab_type": "text"
      },
      "source": [
        "The normalization is across all bases, it is not split by strand; splitting the normalization by strand would potentially lose important information with respect to strand bias and relative errors. The plot below visualizes the final input to the neural network used in `medaka`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kUKQRFL6a6mV",
        "colab_type": "code",
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 317
        },
        "outputId": "2f98ba5e-9907-4f9e-9a41-d7c65790255d"
      },
      "source": [
        "#@markdown ***Plot feature array*** *(click to show code)*\n",
        "\n",
        "from bokeh.plotting import figure\n",
        "from bokeh.models import Range1d\n",
        "import bokeh.io as bkio\n",
        "\n",
        "# select just a region to plot\n",
        "reg = slice(10000,10050)\n",
        "pdata = feature_array[reg] * 255\n",
        "ppos = positions[reg]\n",
        "# create RGBA image\n",
        "pdata = np.stack([pdata]*4, axis=-1)\n",
        "pdata[:,:,3] = 0 # alpha channel\n",
        "pdata[:,[x.upper() == 'A' for x in codes],0] *= 0\n",
        "pdata[:,[x.upper() == 'C' for x in codes],1] *= 0\n",
        "pdata[:,[x.upper() == 'G' for x in codes],2] *= 0\n",
        "pdata[:,[x.upper() == 'T' for x in codes],1:3] *= 0\n",
        "pdata[:,[x.upper() == 'D' for x in codes],0:3] *= 10\n",
        "pdata = pdata.astype(dtype=np.uint8)\n",
        "pdata = np.transpose(pdata, axes=[1,0,2])[::-1,:,:]\n",
        "pdata = 255 - pdata\n",
        "\n",
        "# create a figure\n",
        "p = figure(\n",
        "    title=\"Base counts\",\n",
        "    plot_height=300, plot_width=800)\n",
        "\n",
        "p.x_range.range_padding = p.y_range.range_padding = 0\n",
        "p.image_rgba(image=[pdata], x=0, y=0, dw=pdata.shape[1], dh=pdata.shape[0])\n",
        "ylabels = np.arange(0.5,10.5)\n",
        "p.yaxis.ticker = ylabels\n",
        "p.yaxis.major_label_overrides = dict(zip(ylabels, reversed(codes)))\n",
        "p.y_range = Range1d(\n",
        "    start=0, end=10,\n",
        "    bounds=(0, 10))\n",
        "xlabels = np.arange(0.5, pdata.shape[1])\n",
        "p.xaxis.ticker = xlabels\n",
        "p.xaxis.major_label_overrides = dict(zip(\n",
        "    xlabels, ('{}.{}'.format(x['major'], x['minor']) for x in ppos)\n",
        "))\n",
        "p.xaxis.major_label_orientation = 3.14/2\n",
        "bkio.output_notebook(hide_banner=True)\n",
        "bkio.show(p)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "application/javascript": [
              "\n",
              "(function(root) {\n",
              "  function now() {\n",
              "    return new Date();\n",
              "  }\n",
              "\n",
              "  var force = true;\n",
              "\n",
              "  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
              "    root._bokeh_onload_callbacks = [];\n",
              "    root._bokeh_is_loading = undefined;\n",
              "  }\n",
              "\n",
              "  var JS_MIME_TYPE = 'application/javascript';\n",
              "  var HTML_MIME_TYPE = 'text/html';\n",
              "  var EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
              "  var CLASS_NAME = 'output_bokeh rendered_html';\n",
              "\n",
              "  /**\n",
              "   * Render data to the DOM node\n",
              "   */\n",
              "  function render(props, node) {\n",
              "    var script = document.createElement(\"script\");\n",
              "    node.appendChild(script);\n",
              "  }\n",
              "\n",
              "  /**\n",
              "   * Handle when an output is cleared or removed\n",
              "   */\n",
              "  function handleClearOutput(event, handle) {\n",
              "    var cell = handle.cell;\n",
              "\n",
              "    var id = cell.output_area._bokeh_element_id;\n",
              "    var server_id = cell.output_area._bokeh_server_id;\n",
              "    // Clean up Bokeh references\n",
              "    if (id != null && id in Bokeh.index) {\n",
              "      Bokeh.index[id].model.document.clear();\n",
              "      delete Bokeh.index[id];\n",
              "    }\n",
              "\n",
              "    if (server_id !== undefined) {\n",
              "      // Clean up Bokeh references\n",
              "      var cmd = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
              "      cell.notebook.kernel.execute(cmd, {\n",
              "        iopub: {\n",
              "          output: function(msg) {\n",
              "            var id = msg.content.text.trim();\n",
              "            if (id in Bokeh.index) {\n",
              "              Bokeh.index[id].model.document.clear();\n",
              "              delete Bokeh.index[id];\n",
              "            }\n",
              "          }\n",
              "        }\n",
              "      });\n",
              "      // Destroy server and session\n",
              "      var cmd = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
              "      cell.notebook.kernel.execute(cmd);\n",
              "    }\n",
              "  }\n",
              "\n",
              "  /**\n",
              "   * Handle when a new output is added\n",
              "   */\n",
              "  function handleAddOutput(event, handle) {\n",
              "    var output_area = handle.output_area;\n",
              "    var output = handle.output;\n",
              "\n",
              "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
              "    if ((output.output_type != \"display_data\") || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n",
              "      return\n",
              "    }\n",
              "\n",
              "    var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
              "\n",
              "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
              "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
              "      // store reference to embed id on output_area\n",
              "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
              "    }\n",
              "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
              "      var bk_div = document.createElement(\"div\");\n",
              "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
              "      var script_attrs = bk_div.children[0].attributes;\n",
              "      for (var i = 0; i < script_attrs.length; i++) {\n",
              "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
              "      }\n",
              "      // store reference to server id on output_area\n",
              "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
              "    }\n",
              "  }\n",
              "\n",
              "  function register_renderer(events, OutputArea) {\n",
              "\n",
              "    function append_mime(data, metadata, element) {\n",
              "      // create a DOM node to render to\n",
              "      var toinsert = this.create_output_subarea(\n",
              "        metadata,\n",
              "        CLASS_NAME,\n",
              "        EXEC_MIME_TYPE\n",
              "      );\n",
              "      this.keyboard_manager.register_events(toinsert);\n",
              "      // Render to node\n",
              "      var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
              "      render(props, toinsert[toinsert.length - 1]);\n",
              "      element.append(toinsert);\n",
              "      return toinsert\n",
              "    }\n",
              "\n",
              "    /* Handle when an output is cleared or removed */\n",
              "    events.on('clear_output.CodeCell', handleClearOutput);\n",
              "    events.on('delete.Cell', handleClearOutput);\n",
              "\n",
              "    /* Handle when a new output is added */\n",
              "    events.on('output_added.OutputArea', handleAddOutput);\n",
              "\n",
              "    /**\n",
              "     * Register the mime type and append_mime function with output_area\n",
              "     */\n",
              "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
              "      /* Is output safe? */\n",
              "      safe: true,\n",
              "      /* Index of renderer in `output_area.display_order` */\n",
              "      index: 0\n",
              "    });\n",
              "  }\n",
              "\n",
              "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
              "  if (root.Jupyter !== undefined) {\n",
              "    var events = require('base/js/events');\n",
              "    var OutputArea = require('notebook/js/outputarea').OutputArea;\n",
              "\n",
              "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
              "      register_renderer(events, OutputArea);\n",
              "    }\n",
              "  }\n",
              "\n",
              "  \n",
              "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
              "    root._bokeh_timeout = Date.now() + 5000;\n",
              "    root._bokeh_failed_load = false;\n",
              "  }\n",
              "\n",
              "  var NB_LOAD_WARNING = {'data': {'text/html':\n",
              "     \"<div style='background-color: #fdd'>\\n\"+\n",
              "     \"<p>\\n\"+\n",
              "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
              "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
              "     \"</p>\\n\"+\n",
              "     \"<ul>\\n\"+\n",
              "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
              "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
              "     \"</ul>\\n\"+\n",
              "     \"<code>\\n\"+\n",
              "     \"from bokeh.resources import INLINE\\n\"+\n",
              "     \"output_notebook(resources=INLINE)\\n\"+\n",
              "     \"</code>\\n\"+\n",
              "     \"</div>\"}};\n",
              "\n",
              "  function display_loaded() {\n",
              "    var el = document.getElementById(null);\n",
              "    if (el != null) {\n",
              "      el.textContent = \"BokehJS is loading...\";\n",
              "    }\n",
              "    if (root.Bokeh !== undefined) {\n",
              "      if (el != null) {\n",
              "        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n",
              "      }\n",
              "    } else if (Date.now() < root._bokeh_timeout) {\n",
              "      setTimeout(display_loaded, 100)\n",
              "    }\n",
              "  }\n",
              "\n",
              "\n",
              "  function run_callbacks() {\n",
              "    try {\n",
              "      root._bokeh_onload_callbacks.forEach(function(callback) {\n",
              "        if (callback != null)\n",
              "          callback();\n",
              "      });\n",
              "    } finally {\n",
              "      delete root._bokeh_onload_callbacks\n",
              "    }\n",
              "    console.debug(\"Bokeh: all callbacks have finished\");\n",
              "  }\n",
              "\n",
              "  function load_libs(css_urls, js_urls, callback) {\n",
              "    if (css_urls == null) css_urls = [];\n",
              "    if (js_urls == null) js_urls = [];\n",
              "\n",
              "    root._bokeh_onload_callbacks.push(callback);\n",
              "    if (root._bokeh_is_loading > 0) {\n",
              "      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
              "      return null;\n",
              "    }\n",
              "    if (js_urls == null || js_urls.length === 0) {\n",
              "      run_callbacks();\n",
              "      return null;\n",
              "    }\n",
              "    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
              "    root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
              "\n",
              "    function on_load() {\n",
              "      root._bokeh_is_loading--;\n",
              "      if (root._bokeh_is_loading === 0) {\n",
              "        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
              "        run_callbacks()\n",
              "      }\n",
              "    }\n",
              "\n",
              "    function on_error() {\n",
              "      console.error(\"failed to load \" + url);\n",
              "    }\n",
              "\n",
              "    for (var i = 0; i < css_urls.length; i++) {\n",
              "      var url = css_urls[i];\n",
              "      const element = document.createElement(\"link\");\n",
              "      element.onload = on_load;\n",
              "      element.onerror = on_error;\n",
              "      element.rel = \"stylesheet\";\n",
              "      element.type = \"text/css\";\n",
              "      element.href = url;\n",
              "      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
              "      document.body.appendChild(element);\n",
              "    }\n",
              "\n",
              "    for (var i = 0; i < js_urls.length; i++) {\n",
              "      var url = js_urls[i];\n",
              "      var element = document.createElement('script');\n",
              "      element.onload = on_load;\n",
              "      element.onerror = on_error;\n",
              "      element.async = false;\n",
              "      element.src = url;\n",
              "      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
              "      document.head.appendChild(element);\n",
              "    }\n",
              "  };\n",
              "\n",
              "  function inject_raw_css(css) {\n",
              "    const element = document.createElement(\"style\");\n",
              "    element.appendChild(document.createTextNode(css));\n",
              "    document.body.appendChild(element);\n",
              "  }\n",
              "\n",
              "  \n",
              "  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.4.0.min.js\"];\n",
              "  var css_urls = [];\n",
              "  \n",
              "\n",
              "  var inline_js = [\n",
              "    function(Bokeh) {\n",
              "      Bokeh.set_log_level(\"info\");\n",
              "    },\n",
              "    function(Bokeh) {\n",
              "    \n",
              "    \n",
              "    }\n",
              "  ];\n",
              "\n",
              "  function run_inline_js() {\n",
              "    \n",
              "    if (root.Bokeh !== undefined || force === true) {\n",
              "      \n",
              "    for (var i = 0; i < inline_js.length; i++) {\n",
              "      inline_js[i].call(root, root.Bokeh);\n",
              "    }\n",
              "    } else if (Date.now() < root._bokeh_timeout) {\n",
              "      setTimeout(run_inline_js, 100);\n",
              "    } else if (!root._bokeh_failed_load) {\n",
              "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
              "      root._bokeh_failed_load = true;\n",
              "    } else if (force !== true) {\n",
              "      var cell = $(document.getElementById(null)).parents('.cell').data().cell;\n",
              "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
              "    }\n",
              "\n",
              "  }\n",
              "\n",
              "  if (root._bokeh_is_loading === 0) {\n",
              "    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
              "    run_inline_js();\n",
              "  } else {\n",
              "    load_libs(css_urls, js_urls, function() {\n",
              "      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n",
              "      run_inline_js();\n",
              "    });\n",
              "  }\n",
              "}(window));"
            ],
            "application/vnd.bokehjs_load.v0+json": "\n(function(root) {\n  function now() {\n    return new Date();\n  }\n\n  var force = true;\n\n  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n    root._bokeh_onload_callbacks = [];\n    root._bokeh_is_loading = undefined;\n  }\n\n  \n\n  \n  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n    root._bokeh_timeout = Date.now() + 5000;\n    root._bokeh_failed_load = false;\n  }\n\n  var NB_LOAD_WARNING = {'data': {'text/html':\n     \"<div style='background-color: #fdd'>\\n\"+\n     \"<p>\\n\"+\n     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n     \"</p>\\n\"+\n     \"<ul>\\n\"+\n     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n     \"</ul>\\n\"+\n     \"<code>\\n\"+\n     \"from bokeh.resources import INLINE\\n\"+\n     \"output_notebook(resources=INLINE)\\n\"+\n     \"</code>\\n\"+\n     \"</div>\"}};\n\n  function display_loaded() {\n    var el = document.getElementById(null);\n    if (el != null) {\n      el.textContent = \"BokehJS is loading...\";\n    }\n    if (root.Bokeh !== undefined) {\n      if (el != null) {\n        el.textContent = \"BokehJS \" + root.Bokeh.version + \" successfully loaded.\";\n      }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(display_loaded, 100)\n    }\n  }\n\n\n  function run_callbacks() {\n    try {\n      root._bokeh_onload_callbacks.forEach(function(callback) {\n        if (callback != null)\n          callback();\n      });\n    } finally {\n      delete root._bokeh_onload_callbacks\n    }\n    console.debug(\"Bokeh: all callbacks have finished\");\n  }\n\n  function load_libs(css_urls, js_urls, callback) {\n    if (css_urls == null) css_urls = [];\n    if (js_urls == null) js_urls = [];\n\n    root._bokeh_onload_callbacks.push(callback);\n    if (root._bokeh_is_loading > 0) {\n      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n      return null;\n    }\n    if (js_urls == null || js_urls.length === 0) {\n      run_callbacks();\n      return null;\n    }\n    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n    root._bokeh_is_loading = css_urls.length + js_urls.length;\n\n    function on_load() {\n      root._bokeh_is_loading--;\n      if (root._bokeh_is_loading === 0) {\n        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n        run_callbacks()\n      }\n    }\n\n    function on_error() {\n      console.error(\"failed to load \" + url);\n    }\n\n    for (var i = 0; i < css_urls.length; i++) {\n      var url = css_urls[i];\n      const element = document.createElement(\"link\");\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.rel = \"stylesheet\";\n      element.type = \"text/css\";\n      element.href = url;\n      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n      document.body.appendChild(element);\n    }\n\n    for (var i = 0; i < js_urls.length; i++) {\n      var url = js_urls[i];\n      var element = document.createElement('script');\n      element.onload = on_load;\n      element.onerror = on_error;\n      element.async = false;\n      element.src = url;\n      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n      document.head.appendChild(element);\n    }\n  };\n\n  function inject_raw_css(css) {\n    const element = document.createElement(\"style\");\n    element.appendChild(document.createTextNode(css));\n    document.body.appendChild(element);\n  }\n\n  \n  var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-tables-1.4.0.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-gl-1.4.0.min.js\"];\n  var css_urls = [];\n  \n\n  var inline_js = [\n    function(Bokeh) {\n      Bokeh.set_log_level(\"info\");\n    },\n    function(Bokeh) {\n    \n    \n    }\n  ];\n\n  function run_inline_js() {\n    \n    if (root.Bokeh !== undefined || force === true) {\n      \n    for (var i = 0; i < inline_js.length; i++) {\n      inline_js[i].call(root, root.Bokeh);\n    }\n    } else if (Date.now() < root._bokeh_timeout) {\n      setTimeout(run_inline_js, 100);\n    } else if (!root._bokeh_failed_load) {\n      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n      root._bokeh_failed_load = true;\n    } else if (force !== true) {\n      var cell = $(document.getElementById(null)).parents('.cell').data().cell;\n      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n    }\n\n  }\n\n  if (root._bokeh_is_loading === 0) {\n    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n    run_inline_js();\n  } else {\n    load_libs(css_urls, js_urls, function() {\n      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n      run_inline_js();\n    });\n  }\n}(window));"
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "\n",
              "\n",
              "\n",
              "\n",
              "\n",
              "\n",
              "  <div class=\"bk-root\" id=\"7fe1f8f3-5c78-43aa-9344-53f8312344ce\" data-root-id=\"1001\"></div>\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "application/javascript": [
              "(function(root) {\n",
              "  function embed_document(root) {\n",
              "    \n",
              "  var docs_json = {\"d6396d6f-7dfe-4144-80f0-09095cf34615\":{\"roots\":{\"references\":[{\"attributes\":{\"below\":[{\"id\":\"1012\",\"type\":\"LinearAxis\"}],\"center\":[{\"id\":\"1016\",\"type\":\"Grid\"},{\"id\":\"1021\",\"type\":\"Grid\"}],\"left\":[{\"id\":\"1017\",\"type\":\"LinearAxis\"}],\"plot_height\":300,\"plot_width\":800,\"renderers\":[{\"id\":\"1038\",\"type\":\"GlyphRenderer\"}],\"title\":{\"id\":\"1002\",\"type\":\"Title\"},\"toolbar\":{\"id\":\"1028\",\"type\":\"Toolbar\"},\"x_range\":{\"id\":\"1004\",\"type\":\"DataRange1d\"},\"x_scale\":{\"id\":\"1008\",\"type\":\"LinearScale\"},\"y_range\":{\"id\":\"1042\",\"type\":\"Range1d\"},\"y_scale\":{\"id\":\"1010\",\"type\":\"LinearScale\"}},\"id\":\"1001\",\"subtype\":\"Figure\",\"type\":\"Plot\"},{\"attributes\":{\"formatter\":{\"id\":\"1046\",\"type\":\"BasicTickFormatter\"},\"major_label_overrides\":{\"0.5\":\"D\",\"1.5\":\"d\",\"2.5\":\"T\",\"3.5\":\"G\",\"4.5\":\"C\",\"5.5\":\"A\",\"6.5\":\"t\",\"7.5\":\"g\",\"8.5\":\"c\",\"9.5\":\"a\"},\"ticker\":{\"id\":\"1040\",\"type\":\"FixedTicker\"}},\"id\":\"1017\",\"type\":\"LinearAxis\"},{\"attributes\":{},\"id\":\"1018\",\"type\":\"BasicTicker\"},{\"attributes\":{\"dh\":{\"units\":\"data\",\"value\":10},\"dw\":{\"units\":\"data\",\"value\":50},\"image\":{\"field\":\"image\"},\"x\":{\"value\":0},\"y\":{\"value\":0}},\"id\":\"1036\",\"type\":\"ImageRGBA\"},{\"attributes\":{\"dimension\":1,\"ticker\":{\"id\":\"1018\",\"type\":\"BasicTicker\"}},\"id\":\"1021\",\"type\":\"Grid\"},{\"attributes\":{\"active_drag\":\"auto\",\"active_inspect\":\"auto\",\"active_multi\":null,\"active_scroll\":\"auto\",\"active_tap\":\"auto\",\"tools\":[{\"id\":\"1022\",\"type\":\"PanTool\"},{\"id\":\"1023\",\"type\":\"WheelZoomTool\"},{\"id\":\"1024\",\"type\":\"BoxZoomTool\"},{\"id\":\"1025\",\"type\":\"SaveTool\"},{\"id\":\"1026\",\"type\":\"ResetTool\"},{\"id\":\"1027\",\"type\":\"HelpTool\"}]},\"id\":\"1028\",\"type\":\"Toolbar\"},{\"attributes\":{\"ticks\":[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5]},\"id\":\"1040\",\"type\":\"FixedTicker\"},{\"attributes\":{\"source\":{\"id\":\"1035\",\"type\":\"ColumnDataSource\"}},\"id\":\"1039\",\"type\":\"CDSView\"},{\"attributes\":{},\"id\":\"1022\",\"type\":\"PanTool\"},{\"attributes\":{},\"id\":\"1023\",\"type\":\"WheelZoomTool\"},{\"attributes\":{\"dh\":{\"units\":\"data\",\"value\":10},\"dw\":{\"units\":\"data\",\"value\":50},\"image\":{\"field\":\"image\"},\"x\":{\"value\":0},\"y\":{\"value\":0}},\"id\":\"1037\",\"type\":\"ImageRGBA\"},{\"attributes\":{\"overlay\":{\"id\":\"1048\",\"type\":\"BoxAnnotation\"}},\"id\":\"1024\",\"type\":\"BoxZoomTool\"},{\"attributes\":{},\"id\":\"1025\",\"type\":\"SaveTool\"},{\"attributes\":{},\"id\":\"1047\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"1026\",\"type\":\"ResetTool\"},{\"attributes\":{\"callback\":null,\"data\":{\"image\":[{\"__ndarray__\":\"///////////////////////////Ly8v//////////////////////////////////////+7u7v///////////////////////////////////////////7a2tv/////////////////d3d3////////////////////////////////////////////////////////////d3d3/////////////////7u7u/+7u7v//////3d3d/7q6uv/u7u7//////////////////////+7u7v///////////7q6uv/////////////////Ly8v/////////////////7u7u//////+pqan//////+7u7v/u7u7/////////////////Tk5O/////////////////8vLy//u7u7/7u7u/+7u7v/////////////////u7u7//////8vLy////////////////////////////+7u7v+Xl5f/hoaG////////////7u7u/////////////////////////////v////7///////////////7////+/////v////7///////////////////////////////////////////////n///////////////////////////////////////////////////////////////7//////////v///3z//////////////////////////v////f////+////ev////////////////////z/////////fP////7///////////////////9/f/////////////////////////////////////////f3///8/P////////Pz////////g4P///////99ff///////////////////////////////////v7///7+///+/v////////////////////////////////////////z8//99ff///////3x8////////////////////////iIj///7+////////gYH///z8/////////Pz///7+//////////////7//v/+//7//////4H/gf////////////////99/33///////////////////////////////////////7//v+B/4H////////////////////////////////////////////////////////////////////////////+//7//////////////////////////////////v/+////////////f/9///////////////////////////////////////////////7+/////////v7///7+///////////////////+/v///v7///r6///+/v//iIj////////39////v7///7+////////9fX///r6///+/v//m5v///7+/////////////39///96ev//enr//3p6/////////v7///////////////////z8///+/v///v7//39////////////////////z8////v7///z8////////hIT//////////////Pz////////////////////////////8/////v//////////////9/////////////////////7////+/////P///////////////P////////////////////7//////////////////////////v////7//////////v//////////////iP////////////////////7///////////////////+I//////////////////////////////+I///////////////+/v///////4qK///////////////////////////////////+/v//////////////////+vr///r6//+Skv///////4qK//////////////////////////////////////////////7+///+/v////////7+/////////v7//////////////////4uL////////iIj///////////////////////+QkP////////7+//+Ghv/////////////8/P////////7//v/+//7/////////////////kP+Q/////////////////5b/lv///////////////////////////////////////v/+/43/jf////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////+S/5L//////////////////////////////////////////////v7///////////////////////////////////////////////////////+Pj/////////////////////////7+///6+v////////////+zs///////////////////jY3//4uL//+Kiv//ior/////////////////////////////////////////////hob///////////////////////////////////////+IiP///v7////////8/P8=\",\"dtype\":\"uint8\",\"shape\":[10,50,4]}]},\"selected\":{\"id\":\"1050\",\"type\":\"Selection\"},\"selection_policy\":{\"id\":\"1049\",\"type\":\"UnionRenderers\"}},\"id\":\"1035\",\"type\":\"ColumnDataSource\"},{\"attributes\":{},\"id\":\"1027\",\"type\":\"HelpTool\"},{\"attributes\":{},\"id\":\"1049\",\"type\":\"UnionRenderers\"},{\"attributes\":{},\"id\":\"1008\",\"type\":\"LinearScale\"},{\"attributes\":{\"bottom_units\":\"screen\",\"fill_alpha\":{\"value\":0.5},\"fill_color\":{\"value\":\"lightgrey\"},\"left_units\":\"screen\",\"level\":\"overlay\",\"line_alpha\":{\"value\":1.0},\"line_color\":{\"value\":\"black\"},\"line_dash\":[4,4],\"line_width\":{\"value\":2},\"render_mode\":\"css\",\"right_units\":\"screen\",\"top_units\":\"screen\"},\"id\":\"1048\",\"type\":\"BoxAnnotation\"},{\"attributes\":{\"data_source\":{\"id\":\"1035\",\"type\":\"ColumnDataSource\"},\"glyph\":{\"id\":\"1036\",\"type\":\"ImageRGBA\"},\"hover_glyph\":null,\"muted_glyph\":null,\"nonselection_glyph\":{\"id\":\"1037\",\"type\":\"ImageRGBA\"},\"selection_glyph\":null,\"view\":{\"id\":\"1039\",\"type\":\"CDSView\"}},\"id\":\"1038\",\"type\":\"GlyphRenderer\"},{\"attributes\":{},\"id\":\"1010\",\"type\":\"LinearScale\"},{\"attributes\":{},\"id\":\"1046\",\"type\":\"BasicTickFormatter\"},{\"attributes\":{},\"id\":\"1050\",\"type\":\"Selection\"},{\"attributes\":{\"text\":\"Base counts\"},\"id\":\"1002\",\"type\":\"Title\"},{\"attributes\":{\"formatter\":{\"id\":\"1047\",\"type\":\"BasicTickFormatter\"},\"major_label_orientation\":1.57,\"major_label_overrides\":{\"0.5\":\"4918.2\",\"1.5\":\"4918.3\",\"10.5\":\"4921.1\",\"11.5\":\"4921.2\",\"12.5\":\"4921.3\",\"13.5\":\"4922.0\",\"14.5\":\"4922.1\",\"15.5\":\"4923.0\",\"16.5\":\"4923.1\",\"17.5\":\"4924.0\",\"18.5\":\"4925.0\",\"19.5\":\"4925.1\",\"2.5\":\"4919.0\",\"20.5\":\"4925.2\",\"21.5\":\"4925.3\",\"22.5\":\"4926.0\",\"23.5\":\"4926.1\",\"24.5\":\"4926.2\",\"25.5\":\"4926.3\",\"26.5\":\"4927.0\",\"27.5\":\"4928.0\",\"28.5\":\"4929.0\",\"29.5\":\"4930.0\",\"3.5\":\"4919.1\",\"30.5\":\"4930.1\",\"31.5\":\"4930.2\",\"32.5\":\"4930.3\",\"33.5\":\"4931.0\",\"34.5\":\"4931.1\",\"35.5\":\"4932.0\",\"36.5\":\"4932.1\",\"37.5\":\"4933.0\",\"38.5\":\"4934.0\",\"39.5\":\"4934.1\",\"4.5\":\"4919.2\",\"40.5\":\"4934.2\",\"41.5\":\"4935.0\",\"42.5\":\"4936.0\",\"43.5\":\"4937.0\",\"44.5\":\"4937.1\",\"45.5\":\"4938.0\",\"46.5\":\"4939.0\",\"47.5\":\"4940.0\",\"48.5\":\"4940.1\",\"49.5\":\"4940.2\",\"5.5\":\"4920.0\",\"6.5\":\"4920.1\",\"7.5\":\"4920.2\",\"8.5\":\"4920.3\",\"9.5\":\"4921.0\"},\"ticker\":{\"id\":\"1043\",\"type\":\"FixedTicker\"}},\"id\":\"1012\",\"type\":\"LinearAxis\"},{\"attributes\":{\"bounds\":[0,10],\"callback\":null,\"end\":10},\"id\":\"1042\",\"type\":\"Range1d\"},{\"attributes\":{},\"id\":\"1013\",\"type\":\"BasicTicker\"},{\"attributes\":{\"ticker\":{\"id\":\"1013\",\"type\":\"BasicTicker\"}},\"id\":\"1016\",\"type\":\"Grid\"},{\"attributes\":{\"callback\":null,\"range_padding\":0},\"id\":\"1004\",\"type\":\"DataRange1d\"},{\"attributes\":{\"ticks\":[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5,22.5,23.5,24.5,25.5,26.5,27.5,28.5,29.5,30.5,31.5,32.5,33.5,34.5,35.5,36.5,37.5,38.5,39.5,40.5,41.5,42.5,43.5,44.5,45.5,46.5,47.5,48.5,49.5]},\"id\":\"1043\",\"type\":\"FixedTicker\"}],\"root_ids\":[\"1001\"]},\"title\":\"Bokeh Application\",\"version\":\"1.4.0\"}};\n",
              "  var render_items = [{\"docid\":\"d6396d6f-7dfe-4144-80f0-09095cf34615\",\"roots\":{\"1001\":\"7fe1f8f3-5c78-43aa-9344-53f8312344ce\"}}];\n",
              "  root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
              "\n",
              "  }\n",
              "  if (root.Bokeh !== undefined) {\n",
              "    embed_document(root);\n",
              "  } else {\n",
              "    var attempts = 0;\n",
              "    var timer = setInterval(function(root) {\n",
              "      if (root.Bokeh !== undefined) {\n",
              "        clearInterval(timer);\n",
              "        embed_document(root);\n",
              "      } else {\n",
              "        attempts++;\n",
              "        if (attempts > 100) {\n",
              "          clearInterval(timer);\n",
              "          console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
              "        }\n",
              "      }\n",
              "    }, 10, root)\n",
              "  }\n",
              "})(window);"
            ],
            "application/vnd.bokehjs_exec.v0+json": ""
          },
          "metadata": {
            "tags": [],
            "application/vnd.bokehjs_exec.v0+json": {
              "id": "1001"
            }
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mzklY-kWVR8N",
        "colab_type": "text"
      },
      "source": [
        "# The neural network\n",
        "\n",
        "Having counted bases in an alignment pileup `medaka` proceeds to analyse these counts using a [Recurrent Neural Network](https://en.wikipedia.org/wiki/Recurrent_neural_network), (RNN). A full discussion of such algorithms is beyond the scope of this discussion, this section demonstrates their use in calculating a consensus sequence from the base counts array. When `medaka` is used as a variant caller different methods are used.\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sBtKgLVAvpeG",
        "colab_type": "text"
      },
      "source": [
        "## The model\n",
        "\n",
        "In order to construct a consensus sequence `medaka` uses a multi-layer bidirection RNN. This is defined using the [`keras`](https://www.tensorflow.org/guide/keras) API in `tensorflow`. The following code is adapted from the [`models`](https://github.com/nanoporetech/medaka/blob/d195b9cc1ee7681a121be9fe4fb016a00744ef47/medaka/models.py#L31) module of `medaka`, it has been simplified to show only the essential parts:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7pmHJnPhwwfv",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 255
        },
        "outputId": "c1cbf254-742e-4469-a01a-98faa70e0b79"
      },
      "source": [
        "from pkg_resources import resource_filename\n",
        "\n",
        "import tensorflow as tf\n",
        "from tensorflow.keras.models import Sequential\n",
        "from tensorflow.keras.layers import Dense, GRU, CuDNNGRU, Bidirectional\n",
        "\n",
        "from medaka.labels import BaseLabelScheme\n",
        "\n",
        "# parameters of the model\n",
        "gru_size = 128\n",
        "time_steps, feature_len = (1000, counts.shape[1])\n",
        "symbols = BaseLabelScheme.symbols  # [-, A, C, G, T]\n",
        "num_classes = len(symbols)\n",
        "\n",
        "# build the model\n",
        "model = Sequential(name='medaka')\n",
        "input_shape = (time_steps, feature_len)\n",
        "for i in range(2):\n",
        "    gru = CuDNNGRU(gru_size, return_sequences=True, name=\"gru_{}\".format(i))\n",
        "    model.add(Bidirectional(gru, input_shape=input_shape))\n",
        "model.add(Dense(\n",
        "    num_classes, activation='softmax', name='classify',\n",
        "    input_shape=(time_steps, 2 * gru_size)))\n",
        "\n",
        "# add pre-trained model weights\n",
        "weight_file = resource_filename('medaka','data/r941_min_high_g344_model.hdf5')\n",
        "model.load_weights(weight_file)\n",
        "model.summary()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Model: \"medaka\"\n",
            "_________________________________________________________________\n",
            "Layer (type)                 Output Shape              Param #   \n",
            "=================================================================\n",
            "bidirectional_16 (Bidirectio (None, 1000, 256)         107520    \n",
            "_________________________________________________________________\n",
            "bidirectional_17 (Bidirectio (None, 1000, 256)         296448    \n",
            "_________________________________________________________________\n",
            "classify (Dense)             (None, 1000, 5)           1285      \n",
            "=================================================================\n",
            "Total params: 405,253\n",
            "Trainable params: 405,253\n",
            "Non-trainable params: 0\n",
            "_________________________________________________________________\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7-Hojb583Ida",
        "colab_type": "text"
      },
      "source": [
        "The model takes the count matrix as input and outputs for each corresponding column of the counts matrix a set of five scores. The five scores express the possibility that the consensus sequence should contain one of the four bases A, C, G, or T, or a gap '-' character at the pileup column under consideration."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EWt1pOLd9yPc",
        "colab_type": "text"
      },
      "source": [
        "## Making predictions\n",
        "\n",
        "In order to make predictions using the RNN, `medaka` splits the normalized counts array into overlapping chunks before processing by the model. Chunking the array allows for more efficient parallel computation while overlapping is a mitigation against edge-effects at the boundaries of chunks.\n",
        "\n",
        "As mentioned above in the aside on performance, `medaka` has a somewhat elaborate system for managing data chunks. For the purposes of exposition the code below implements a simple chunking and batching of the data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "G6BtUH800yGN",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        },
        "outputId": "f665cb09-f538-486f-def7-c00e98f148fc"
      },
      "source": [
        "from functools import partial\n",
        "from medaka.common import sliding_window, grouper\n",
        "\n",
        "# create a function to perform windowing on an array\n",
        "overlap = 200\n",
        "window = partial(\n",
        "    sliding_window,\n",
        "    window=time_steps, step=time_steps - overlap, axis=0)\n",
        "\n",
        "# run the network on input data\n",
        "def get_predictions(data, batch_size=40):\n",
        "    for batch in grouper(data, batch_size=batch_size):\n",
        "        batch = np.stack(batch)\n",
        "        results = model.predict_on_batch(batch)\n",
        "        yield from results\n",
        "\n",
        "t0 = now()\n",
        "predictions = get_predictions(window(feature_array))\n",
        "seq_chunks = list()\n",
        "for pred predictions:\n",
        "    # remove half the overlapping region of chunks\n",
        "    pred = pred[overlap // 2:-overlap // 2]\n",
        "    # find the most likely base at each position and form the sequence\n",
        "    mp = np.argmax(pred, -1)\n",
        "    seq = ''.join((symbols[x] for x in mp))\n",
        "    seq = seq.replace('*', '')\n",
        "    seq_chunks.append(seq)\n",
        "sequence = ''.join(seq_chunks)\n",
        "t1 = now()\n",
        "print(\"{:.2f} to run predictions\".format(t1 - t0))\n",
        "print(\"Total sequence length: {}.\".format(len(sequence)))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "7.50 to run predictions\n",
            "Total sequence length: 1500668.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M50sfRA70uX3",
        "colab_type": "text"
      },
      "source": [
        "The code performs a simple undoing of the overlapping before stitching the consensus sequence pieces back together. This is sufficient to obtain results here; the full `medaka` implementation also keeps track of the `positions` array to ensure the sequence stitching is performed correctly with respect to the original input reference sequence.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6MWZ7xM0QEbj",
        "colab_type": "text"
      },
      "source": [
        "### Checking our results\n",
        "\n",
        "We can write out the full consensus sequence derived above and compare it to the truth sequence by using the `assess_assembly` program from the [`pomoxis`](https://github.com/nanoporetech/pomoxis) package. By also examining the original draft sequence, we can see the improvement in quality from `medaka`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IzSkq4pd6E4j",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 731
        },
        "outputId": "fc31bf38-80d7-4a2d-f589-edefaba856ce"
      },
      "source": [
        "output = \"/content/data/output.fasta\"\n",
        "with open(output, 'w') as fh:\n",
        "    fh.write(\">seq\\n{}\\n\".format(sequence))\n",
        "for fname in (\"/content/data/saureus_canu.fasta\", output):\n",
        "    print(\"Analysing: {}.\".format(fname))\n",
        "    !assess_assembly -r /content/data/references.fasta -i \"$fname\" 2>/dev/null\n",
        "    print(\"\\n\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Analysing: /content/data/saureus_canu.fasta.\n",
            "Writing list of indels 100 bases and longer to assm_indel_ge100.txt.\n",
            "#  Percentage Errors\n",
            "  name     mean     q10      q50      q90   \n",
            " err_ont  0.133%   0.056%   0.068%   0.290% \n",
            " err_bal  0.133%   0.056%   0.068%   0.290% \n",
            "    iden  0.004%   0.000%   0.001%   0.024% \n",
            "     del  0.117%   0.049%   0.059%   0.166% \n",
            "     ins  0.012%   0.003%   0.006%   0.046% \n",
            "\n",
            "#  Q Scores\n",
            "  name     mean      q10      q50      q90  \n",
            " err_ont  28.76    32.54    31.65    25.38  \n",
            " err_bal  28.76    32.54    31.65    25.38  \n",
            "    iden  44.39      inf    50.00    36.20  \n",
            "     del  29.30    33.10    32.29    27.79  \n",
            "     ins  39.20    45.23    42.22    33.41  \n",
            "\n",
            "All done, output written to assm_stats.txt, assm_summ.txt and assm_indel_ge100.txt\n",
            "\n",
            "\n",
            "Analysing: /content/data/output.fasta.\n",
            "Writing list of indels 100 bases and longer to assm_indel_ge100.txt.\n",
            "#  Percentage Errors\n",
            "  name     mean     q10      q50      q90   \n",
            " err_ont  0.007%   0.003%   0.005%   0.013% \n",
            " err_bal  0.007%   0.003%   0.006%   0.013% \n",
            "    iden  0.001%   0.000%   0.000%   0.003% \n",
            "     del  0.002%   0.000%   0.001%   0.004% \n",
            "     ins  0.005%   0.002%   0.005%   0.008% \n",
            "\n",
            "#  Q Scores\n",
            "  name     mean      q10      q50      q90  \n",
            " err_ont  41.55    45.23    42.60    39.03  \n",
            " err_bal  41.55    45.23    42.60    39.03  \n",
            "    iden  51.35      inf      inf    46.02  \n",
            "     del  48.15      inf    50.00    44.56  \n",
            "     ins  43.25    46.99    43.01    41.25  \n",
            "\n",
            "All done, output written to assm_stats.txt, assm_summ.txt and assm_indel_ge100.txt\n",
            "\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o3FW8O5xQ8GD",
        "colab_type": "text"
      },
      "source": [
        "# Remarks\n",
        "\n",
        "In this short walkthrough we have examined some of the internals of Oxford Nanopore Technologies' `medaka` program performs GPU accelerated consensus calculations from aligned sequencing data. The public `medaka` codebase implements various alternative forms of the algorithms presented here including run length compression and support for multiple datatypes. Hopefully this guide will prove useful to anyone wishing to implement algorithms similar to that implemented in `medaka`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2JZFTwhZbgfm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}