fiddle.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "fiddle.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [
        "eYT-RNR_gGGL",
        "tGklhEqLdkl2",
        "PyyE4IMtP4iM",
        "qARBUY0Sn8Xh",
        "ldOTIObHi_A9",
        "xzzTYPYJjP3e",
        "88e0wDN4XSTw",
        "x3-4W64keKrx",
        "TtAM-2bCq_lS",
        "-LSZGQqvf5jE",
        "MMOoeTdH33oQ",
        "Agxyr5Cx_ZGg",
        "uGwEmm16_a_H",
        "X7aQLpfcE-sE",
        "xIEhmaYK_gDM",
        "vBYh84bRhQRs",
        "91pmO_OOpMou",
        "dk4oJHoEEvTA",
        "NPcWRSFcxJzW",
        "jUdA9dZ8eL4e",
        "911AoZz24Ai1",
        "XM95zaCzuO6v",
        "e-9uvmZzCHkn",
        "6foHpV5LzY1E",
        "WENdY0K8xwUy",
        "t_wnvJc4Yovp",
        "CkOUInXY-A2Z",
        "tar6re3PzUwf",
        "TGqN3e7wx80Z",
        "ZvNrKCzjDWm8",
        "e_v4KYJIBox6",
        "Ki9lBEBjbA3C",
        "7y9xKGTflGui",
        "4vSGLPNmleBT"
      ],
      "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/churchmanlab/FIDDLE/blob/master/fiddle.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eII45S2yk1hQ",
        "colab_type": "text"
      },
      "source": [
        "# SET UP - open via CoLab"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "86Eg7Zu-vjWX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# check present working directory\n",
        "!pwd"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZNaa_e6E5Kh4",
        "colab_type": "text"
      },
      "source": [
        "**_user specific variables_**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "tdveRrr3rgk0",
        "colab": {}
      },
      "source": [
        "# head directory\n",
        "_head_dir = \"< your_head_directory >\"\n",
        "# compute cluster specificities\n",
        "_slurm = {\n",
        "  \"email\": \"< your_compute_cluster_email >\",\n",
        "  \"virtual_env\": \"< your_python3.6_virtual_environment >\"\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m366Aq6f5e_n",
        "colab_type": "text"
      },
      "source": [
        "**_load the following_**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "cellView": "form",
        "id": "5xm1ayEGrgkt",
        "colab": {}
      },
      "source": [
        "#@title <-- imports / settings\n",
        "\n",
        "# native libs\n",
        "import os, re, subprocess, sys\n",
        "# external libs\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "from scipy import signal as sig\n",
        "from scipy import special\n",
        "import scipy.stats as stats\n",
        "# style\n",
        "# sns.set(style=\"white\", rc={\"axes.facecolor\": (0, 0, 0, 0)})\n",
        "sns.set_style(\"ticks\")\n",
        "\n",
        "print(\"imports / settings are loaded\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "RPMqhjtYrgkk",
        "colab": {}
      },
      "source": [
        "#@title <-- backend variables { run: \"auto\", display-mode: \"form\" }\n",
        "\n",
        "# arabic <---> roman conversions\n",
        "_a2r = {\n",
        "  \"chrI\": 1, 1: \"chrI\",\n",
        "  \"chrII\": 2, 2: \"chrII\",\n",
        "  \"chrIII\": 3, 3: \"chrIII\",\n",
        "  \"chrIV\": 4, 4: \"chrIV\",\n",
        "  \"chrV\": 5, 5: \"chrV\",\n",
        "  \"chrVI\": 6, 6: \"chrVI\",\n",
        "  \"chrVII\": 7, 7: \"chrVII\",\n",
        "  \"chrVIII\": 8, 8: \"chrVIII\",\n",
        "  \"chrIX\": 9, 9: \"chrIX\",\n",
        "  \"chrX\": 10, 10: \"chrX\",\n",
        "  \"chrXI\": 11, 11: \"chrXI\",\n",
        "  \"chrXII\": 12, 12: \"chrXII\",\n",
        "  \"chrXIII\": 13, 13: \"chrXIII\",\n",
        "  \"chrXIV\": 14, 14: \"chrXIV\",\n",
        "  \"chrXV\": 15, 15: \"chrXV\",\n",
        "  \"chrXVI\": 16, 16: \"chrXVI\"\n",
        "}\n",
        "\n",
        "# ENSEMBL Fungi genome annotations &\n",
        "# Winston (HMS Genetics) lab's TSS assay-based annotations\n",
        "_annotations = {\n",
        "  \"sc\": {\n",
        "    \"gff\": \"Saccharomyces_cerevisiae.R64-1-1.40.gff3\",\n",
        "    \"tss\": \"Scer_polIItranscripts-adjustedTSS.bed\"\n",
        "  },\n",
        "  \"kl\": {\n",
        "    \"gff\": \"GCF_000002515.2_ASM251v1_genomic.gff\",\n",
        "    \"tss\": {\n",
        "      1: \"Klac1.ctss\",\n",
        "      2: \"Klac2.ctss\"\n",
        "    }\n",
        "  }\n",
        "}\n",
        "\n",
        "# bedgraph files\n",
        "_bedgraphs = {\n",
        "  \"ms\": {\n",
        "    \"sc\": {\n",
        "      \"wt\": {\n",
        "        1: \"WT-30C-2_mnase-midpoint-counts.bedgraph\",\n",
        "        2: \"WT-30C-1_mnase-midpoint-counts.bedgraph\",\n",
        "        3: \"WT-30C-3_mnase-midpoint-counts.bedgraph\",\n",
        "        4: \"GSM2561057_Dyads_H3_CC_rep_1.bw.bedGraph\"\n",
        "      },\n",
        "      \"ko\": {\n",
        "        \"spt6-1004\": \"spt6-1004-37C-2_mnase-midpoint-counts.bedgraph\",\n",
        "        \"spt6-YW\": \"spt6-YW-37C-1_mnase-midpoint-counts.bedgraph\"\n",
        "      }\n",
        "    },\n",
        "    \"kl\": {\n",
        "      1: \"GSM552917_Klac_CM_Dec0908_nucCount_100_edited.wig\",\n",
        "      2: \"GSM552917_Klac_CM_Dec0908_nucCount_100_edited.wig\"\n",
        "    }\n",
        "  },\n",
        "  \"ns\": {\n",
        "    \"sc\": {\n",
        "      \"wt\": {\n",
        "        \"pos\": \"WT-30C-1_netseq-5end-counts-plus.bedgraph\",\n",
        "        \"neg\": \"WT-30C-1_netseq-5end-counts-minus.bedgraph\"\n",
        "      },\n",
        "      \"ko\": {\n",
        "        \"chd1\": {\n",
        "          \"pos\": \"chd1-4.pos.bed\",\n",
        "          \"neg\": \"chd1-4.neg.bed\"\n",
        "        },\n",
        "        \"ino80\": {\n",
        "          \"pos\": \"ino80-2.pos.bed\",\n",
        "          \"neg\": \"ino80-2.neg.bed\"\n",
        "        },\n",
        "        \"isw1\": {\n",
        "          \"pos\": \"isw1-2.pos.bed\",\n",
        "          \"neg\": \"isw1-2.neg.bed\"\n",
        "        },\n",
        "        \"isw2\": {\n",
        "          \"pos\": \"isw2-4.pos.bed\",\n",
        "          \"neg\": \"isw2-4.neg.bed\"\n",
        "        },\n",
        "        \"rsc\": { # rsc30\n",
        "          \"pos\": \"rsc30-2.pos.bed\",\n",
        "          \"neg\": \"rsc30-2.neg.bed\"\n",
        "        },\n",
        "        \"spt6-1004\": {\n",
        "          \"pos\": \"spt6-1004-37C-1_netseq-5end-counts-plus.bedgraph\",\n",
        "          \"neg\": \"spt6-1004-37C-1_netseq-5end-counts-minus.bedgraph\"\n",
        "        },\n",
        "        \"spt6-YW\": {\n",
        "          \"pos\": \"spt6-YW-37C-2_netseq-5end-counts-plus.bedgraph\",\n",
        "          \"neg\": \"spt6-YW-37C-2_netseq-5end-counts-minus.bedgraph\"\n",
        "        }\n",
        "      }\n",
        "    },\n",
        "    \"kl\": \"GSM2339539_YJ160_1.bedgraph\"\n",
        "  }\n",
        "}\n",
        "\n",
        "# chromosome homogenization (gff / gff3 / bedgraph)\n",
        "_chrom_convert = {\n",
        "  # S. cerevisiae\n",
        "  \"I\": \"chrI\",\n",
        "  \"II\": \"chrII\",\n",
        "  \"III\": \"chrIII\",\n",
        "  \"IV\": \"chrIV\",\n",
        "  \"V\": \"chrV\",\n",
        "  \"VI\": \"chrVI\",\n",
        "  \"VII\": \"chrVII\",\n",
        "  \"VIII\": \"chrVIII\",\n",
        "  \"IX\": \"chrIX\",\n",
        "  \"X\": \"chrX\",\n",
        "  \"XI\": \"chrXI\",\n",
        "  \"XII\": \"chrXII\",\n",
        "  \"XIII\": \"chrXIII\",\n",
        "  \"XIV\": \"chrXIV\",\n",
        "  \"XV\": \"chrXV\",\n",
        "  \"XVI\": \"chrXVI\",\n",
        "  # K. lactis\n",
        "  \"klac_Chr_1\": \"chrI\",\n",
        "  \"klac_Chr_2\": \"chrII\",\n",
        "  \"klac_Chr_3\": \"chrIII\",\n",
        "  \"klac_Chr_4\": \"chrIV\",\n",
        "  \"klac_Chr_5\": \"chrV\",\n",
        "  \"klac_Chr_6\": \"chrVI\",\n",
        "  \"NC_006037.1\": \"chrI\",\n",
        "  \"NC_006038.1\": \"chrII\",\n",
        "  \"NC_006039.1\": \"chrIII\",\n",
        "  \"NC_006040.1\": \"chrIV\",\n",
        "  \"NC_006041.1\": \"chrV\",\n",
        "  \"NC_006042.1\": \"chrVI\",\n",
        "  # remove mitochondria and other\n",
        "  \"Mito\": \"remove\",\n",
        "  \"MT\": \"remove\",\n",
        "  \"MTR\": \"remove\",\n",
        "  \"NC_006077.1\": \"remove\"\n",
        "}\n",
        "\n",
        "# species, chromosome lengths\n",
        "_chrom_lengths = {\n",
        "  \"sc\": {\n",
        "    \"chrI\": 230218,\n",
        "    \"chrII\": 813184,\n",
        "    \"chrIII\": 316620,\n",
        "    \"chrIV\": 1531933,\n",
        "    \"chrV\": 576874,\n",
        "    \"chrVI\": 270161,\n",
        "    \"chrVII\": 1090940,\n",
        "    \"chrVIII\": 562643,\n",
        "    \"chrIX\": 439888,\n",
        "    \"chrX\": 745751,\n",
        "    \"chrXI\": 666816,\n",
        "    \"chrXII\": 1078177,\n",
        "    \"chrXIII\": 924431,\n",
        "    \"chrXIV\": 784333,\n",
        "    \"chrXV\": 1091291,\n",
        "    \"chrXVI\": 948066\n",
        "  },\n",
        "  \"kl\": {\n",
        "    \"chrI\": 1062590,\n",
        "    \"chrII\": 1320834,\n",
        "    \"chrIII\": 1753957,\n",
        "    \"chrIV\": 1715506,\n",
        "    \"chrV\": 2234072,\n",
        "    \"chrVI\": 2602197\n",
        "  }\n",
        "}\n",
        "\n",
        "# directories\n",
        "_dirs = {\n",
        "  \"raw\": _head_dir + \"DATA/RAW/\",\n",
        "  \"collated\": _head_dir + \"DATA/COLLATED/\",\n",
        "  \"figs_stats\": _head_dir + \"MODELS/FIGURES_STATS/\",\n",
        "  \"logs\": _head_dir + \"logs/\",\n",
        "  \"model_inputs\": _head_dir + \"MODELS/INPUTS/\",\n",
        "  \"model_outputs\": _head_dir + \"MODELS/OUTPUTS/\",\n",
        "  \"trained_models\": _head_dir + \"MODELS/TRAINED/\",\n",
        "  \"scripts\": _head_dir + \"MODELS/SCRIPTS/\"\n",
        "}\n",
        "\n",
        "# DNA sequence FASTA files\n",
        "_fastas = {\n",
        "  \"sc\": {\n",
        "    \"chrI\": \"SC_chr01.fsa\",\n",
        "    \"chrII\": \"SC_chr02.fsa\",\n",
        "    \"chrIII\": \"SC_chr03.fsa\",\n",
        "    \"chrIV\": \"SC_chr04.fsa\",\n",
        "    \"chrV\": \"SC_chr05.fsa\",\n",
        "    \"chrVI\": \"SC_chr06.fsa\",\n",
        "    \"chrVII\": \"SC_chr07.fsa\",\n",
        "    \"chrVIII\": \"SC_chr08.fsa\",\n",
        "    \"chrIX\": \"SC_chr09.fsa\",\n",
        "    \"chrX\": \"SC_chr10.fsa\",\n",
        "    \"chrXI\": \"SC_chr11.fsa\",\n",
        "    \"chrXII\": \"SC_chr12.fsa\",\n",
        "    \"chrXIII\": \"SC_chr13.fsa\",\n",
        "    \"chrXIV\": \"SC_chr14.fsa\",\n",
        "    \"chrXV\": \"SC_chr15.fsa\",\n",
        "    \"chrXVI\": \"SC_chr16.fsa\"\n",
        "  },\n",
        "  \"kl\": {\n",
        "    \"chrI\": \"KL_chr01.fsa\",\n",
        "    \"chrII\": \"KL_chr02.fsa\",\n",
        "    \"chrIII\": \"KL_chr03.fsa\",\n",
        "    \"chrIV\": \"KL_chr04.fsa\",\n",
        "    \"chrV\": \"KL_chr05.fsa\",\n",
        "    \"chrVI\": \"KL_chr06.fsa\"\n",
        "  }\n",
        "}\n",
        "\n",
        "# thresholds and region boundary considerations\n",
        "_heuristics = {\n",
        "  \"AUG_width\": 2048,\n",
        "  \"AUG_us_offset\": 64,\n",
        "  \"ol_reg_width\": 512,\n",
        "  \"TSS_width\": 1024,\n",
        "  \"TSS_seq_signal_percentile\": 95,\n",
        "  \"epsilon\": 1e-7,\n",
        "  \"training_regions\": 16\n",
        "}\n",
        "\n",
        "# S. cerevisiae cell states\n",
        "_sc_N_states = [\n",
        "  \"N_wt\",\n",
        "  \"N_ko_chd1\",\n",
        "  \"N_ko_ino80\",\n",
        "  \"N_ko_isw1\",\n",
        "  \"N_ko_isw2\",\n",
        "  \"N_ko_rsc\",\n",
        "  \"N_ko_spt6-1004\",\n",
        "  \"N_ko_spt6-YW\"\n",
        "]\n",
        "\n",
        "# strand conversions\n",
        "_strandedness = {\n",
        "  \"+\": \"pos\", \"pos\": \"+\",\n",
        "  \"-\": \"neg\", \"neg\": \"-\"\n",
        "}\n",
        "\n",
        "print(\"backend variables are loaded\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "cellView": "form",
        "id": "liGliVgyrgkd",
        "colab": {}
      },
      "source": [
        "#@title <-- custom functions\n",
        "\n",
        "def _np_5_to_3_kl(gene, data, width=1024):\n",
        "  \"\"\"orient extracted data in 5' --> 3' direction\"\"\"\n",
        "  # get genomic and sequencing data at gene location\n",
        "  seqid, strand = gene[\"seqid\"], gene[\"strand\"]\n",
        "  start, end = gene[\"TSS_start\"], gene[\"TSS_end\"]\n",
        "  locus = _np_get_locus_kl_sp(seqid, start, end, strand, data)\n",
        "  # ancillary information at gene location\n",
        "  ancillary_df = pd.DataFrame({\n",
        "    \"strand\": np.repeat(gene[\"strand\"], width),\n",
        "    \"ID\": np.repeat(gene[\"ID\"], width),\n",
        "    \"seqid\": np.repeat(gene[\"seqid\"], width)      \n",
        "  })\n",
        "  # return 5' --> 3' orientation genomic and sequencing data\n",
        "  if gene[\"strand\"] is \"+\":\n",
        "    df = pd.concat(\n",
        "      objs=(ancillary_df, locus.rename(columns={\n",
        "        \"NS_pos\": \"NS_sense\",\n",
        "        \"NS_neg\": \"NS_antisense\"\n",
        "      })), axis=1)\n",
        "  else: # gene[\"strand\"] is \"-\"\n",
        "    # reverse negative strand data\n",
        "    for col in [\n",
        "      \"DS_A\", \"DS_T\", \"DS_G\", \"DS_C\",\n",
        "      \"MS_1\", \"MS_2\",\n",
        "      \"NS_pos\", \"NS_neg\"\n",
        "    ]: locus[col] = locus[col].values[::-1]\n",
        "    # flip negative strand data\n",
        "    df = pd.concat(\n",
        "      objs=(ancillary_df, locus.rename(columns={\n",
        "        \"NS_pos\": \"NS_antisense\",\n",
        "        \"NS_neg\": \"NS_sense\",\n",
        "        \"DS_A\": \"DS_T\",\n",
        "        \"DS_T\": \"DS_A\",\n",
        "        \"DS_G\": \"DS_C\",\n",
        "        \"DS_C\": \"DS_G\"\n",
        "      })), axis=1)\n",
        "  return df[[\n",
        "    \"ID\", \"strand\", \"seqid\", \"position\", \n",
        "    \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "    \"NS_sense\", \"NS_antisense\",\n",
        "    \"MS_1\", \"MS_2\"]]\n",
        "\n",
        "\n",
        "def _np_5_to_3_sc(gene, approach, data, width=1024):\n",
        "  \"\"\"orient extracted data in 5' --> 3' direction\"\"\"\n",
        "  # get genomic and sequencing data at gene location\n",
        "  seqid, strand = gene[\"seqid\"], gene[\"strand\"]\n",
        "  start, end = gene[approach + \"_start\"], gene[approach + \"_end\"]\n",
        "  # random boundaries for AUG - phasing\n",
        "  if approach == \"AUG\":\n",
        "    start, end = _np_random_start_end(start, end, strand, width)\n",
        "  locus = _np_get_locus_sc(approach, seqid, start, end, strand, data)\n",
        "  # ancillary information at gene location\n",
        "  ancillary_df = pd.DataFrame({\n",
        "    \"strand\": np.repeat(gene[\"strand\"], width),\n",
        "    \"ID\": np.repeat(gene[\"ID\"], width),\n",
        "    \"seqid\": np.repeat(gene[\"seqid\"], width)      \n",
        "  })\n",
        "  # return 5' --> 3' orientation genomic and sequencing data\n",
        "  if approach == \"TSS\":\n",
        "    if gene[\"strand\"] is \"+\":\n",
        "      df = pd.concat(\n",
        "        objs=(ancillary_df, locus.rename(columns={\n",
        "          \"NS_wt_pos\": \"NS_wt_sense\",\n",
        "          \"NS_wt_neg\": \"NS_wt_antisense\",\n",
        "          \"NS_ko_chd1_pos\": \"NS_ko_chd1_sense\",\n",
        "          \"NS_ko_chd1_neg\": \"NS_ko_chd1_antisense\",\n",
        "          \"NS_ko_ino80_pos\": \"NS_ko_ino80_sense\",\n",
        "          \"NS_ko_ino80_neg\": \"NS_ko_ino80_antisense\",\n",
        "          \"NS_ko_isw1_pos\": \"NS_ko_isw1_sense\",\n",
        "          \"NS_ko_isw1_neg\": \"NS_ko_isw1_antisense\",\n",
        "          \"NS_ko_isw2_pos\": \"NS_ko_isw2_sense\",\n",
        "          \"NS_ko_isw2_neg\": \"NS_ko_isw2_antisense\",\n",
        "          \"NS_ko_rsc_pos\": \"NS_ko_rsc_sense\",\n",
        "          \"NS_ko_rsc_neg\": \"NS_ko_rsc_antisense\",\n",
        "          \"NS_ko_spt6-1004_pos\": \"NS_ko_spt6-1004_sense\",\n",
        "          \"NS_ko_spt6-1004_neg\": \"NS_ko_spt6-1004_antisense\",\n",
        "          \"NS_ko_spt6-YW_pos\": \"NS_ko_spt6-YW_sense\",\n",
        "          \"NS_ko_spt6-YW_neg\": \"NS_ko_spt6-YW_antisense\"\n",
        "        })), axis=1)\n",
        "    else: # gene[\"strand\"] is \"-\"\n",
        "      # reverse negative strand data\n",
        "      for col in [\n",
        "        \"DS_A\", \"DS_T\", \"DS_G\", \"DS_C\",\n",
        "        \"NS_wt_pos\", \"NS_wt_neg\",\n",
        "        \"NS_ko_chd1_pos\", \"NS_ko_chd1_neg\",\n",
        "        \"NS_ko_ino80_pos\", \"NS_ko_ino80_neg\",\n",
        "        \"NS_ko_isw1_pos\", \"NS_ko_isw1_neg\",\n",
        "        \"NS_ko_isw2_pos\", \"NS_ko_isw2_neg\",\n",
        "        \"NS_ko_rsc_pos\", \"NS_ko_rsc_neg\",\n",
        "        \"NS_ko_spt6-1004_pos\", \"NS_ko_spt6-1004_neg\",\n",
        "        \"NS_ko_spt6-YW_pos\", \"NS_ko_spt6-YW_neg\",\n",
        "        \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\",\n",
        "        \"MS_ko_spt6-1004\", \"MS_ko_spt6-YW\"\n",
        "      ]: locus[col] = locus[col].values[::-1]\n",
        "      # flip negative strand data\n",
        "      df = pd.concat(\n",
        "        objs=(ancillary_df, locus.rename(columns={\n",
        "          \"DS_A\": \"DS_T\",\n",
        "          \"DS_T\": \"DS_A\",\n",
        "          \"DS_G\": \"DS_C\",\n",
        "          \"DS_C\": \"DS_G\",\n",
        "          \"NS_wt_pos\": \"NS_wt_antisense\",\n",
        "          \"NS_wt_neg\": \"NS_wt_sense\",\n",
        "          \"NS_ko_chd1_pos\": \"NS_ko_chd1_antisense\",\n",
        "          \"NS_ko_chd1_neg\": \"NS_ko_chd1_sense\",\n",
        "          \"NS_ko_ino80_pos\": \"NS_ko_ino80_antisense\",\n",
        "          \"NS_ko_ino80_neg\": \"NS_ko_ino80_sense\",\n",
        "          \"NS_ko_isw1_pos\": \"NS_ko_isw1_antisense\",\n",
        "          \"NS_ko_isw1_neg\": \"NS_ko_isw1_sense\",\n",
        "          \"NS_ko_isw2_pos\": \"NS_ko_isw2_antisense\",\n",
        "          \"NS_ko_isw2_neg\": \"NS_ko_isw2_sense\",\n",
        "          \"NS_ko_rsc_pos\": \"NS_ko_rsc_antisense\",\n",
        "          \"NS_ko_rsc_neg\": \"NS_ko_rsc_sense\",\n",
        "          \"NS_ko_spt6-1004_pos\": \"NS_ko_spt6-1004_sense\",\n",
        "          \"NS_ko_spt6-1004_neg\": \"NS_ko_spt6-1004_antisense\",\n",
        "          \"NS_ko_spt6-YW_pos\": \"NS_ko_spt6-YW_sense\",\n",
        "          \"NS_ko_spt6-YW_neg\": \"NS_ko_spt6-YW_antisense\"\n",
        "        })), axis=1)\n",
        "    return df[[\n",
        "      \"ID\", \"strand\", \"seqid\", \"position\", \n",
        "      \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "      \"NS_wt_sense\", \"NS_wt_antisense\",\n",
        "      \"NS_ko_chd1_sense\", \"NS_ko_chd1_antisense\",\n",
        "      \"NS_ko_ino80_sense\", \"NS_ko_ino80_antisense\",\n",
        "      \"NS_ko_isw1_sense\", \"NS_ko_isw1_antisense\",\n",
        "      \"NS_ko_isw2_sense\", \"NS_ko_isw2_antisense\",\n",
        "      \"NS_ko_rsc_sense\", \"NS_ko_rsc_antisense\",\n",
        "      \"NS_ko_spt6-1004_sense\", \"NS_ko_spt6-1004_antisense\",\n",
        "      \"NS_ko_spt6-YW_sense\", \"NS_ko_spt6-YW_antisense\",\n",
        "      \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\",\n",
        "      \"MS_ko_spt6-1004\", \"MS_ko_spt6-YW\"\n",
        "    ]]\n",
        "  else: # approach == \"AUG\"\n",
        "    if gene[\"strand\"] is \"+\":\n",
        "      df = pd.concat(\n",
        "        objs=(ancillary_df, locus.rename(columns={\n",
        "          \"NS_wt_pos\": \"NS_wt_sense\",\n",
        "          \"NS_wt_neg\": \"NS_wt_antisense\"\n",
        "        })), axis=1)\n",
        "    else: # gene[\"strand\"] is \"-\"\n",
        "      # reverse negative strand data\n",
        "      for col in [\n",
        "        \"DS_A\", \"DS_T\", \"DS_G\", \"DS_C\",\n",
        "        \"NS_wt_pos\", \"NS_wt_neg\",\n",
        "        \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\"\n",
        "      ]: locus[col] = locus[col].values[::-1]\n",
        "      # flip negative strand data\n",
        "      df = pd.concat(\n",
        "        objs=(ancillary_df, locus.rename(columns={\n",
        "          \"DS_A\": \"DS_T\",\n",
        "          \"DS_T\": \"DS_A\",\n",
        "          \"DS_G\": \"DS_C\",\n",
        "          \"DS_C\": \"DS_G\",\n",
        "          \"NS_wt_pos\": \"NS_wt_antisense\",\n",
        "          \"NS_wt_neg\": \"NS_wt_sense\"\n",
        "        })), axis=1)\n",
        "    return df[[\n",
        "      \"ID\", \"strand\", \"seqid\", \"position\", \n",
        "      \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "      \"NS_wt_sense\", \"NS_wt_antisense\",\n",
        "      \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\"]]\n",
        "\n",
        "\n",
        "def _np_avg_var_df(df, reqs=[\"\"]):\n",
        "  \"\"\"average and variance of species / cell-state MNase-seq predictions\n",
        "  df := collated predictions dataframe\n",
        "  reqs := list of necessary strings\"\"\"\n",
        "  return {\n",
        "    \"avg\": np.average(_np_get_seq_data([\n",
        "      col for col in df.columns.values\n",
        "        if all([y in col for y in reqs + [\"mean\"]])\n",
        "    ], df), 1),\n",
        "    \"var\": np.var(_np_get_seq_data([\n",
        "      col for col in df.columns.values\n",
        "        if all([y in col for y in reqs + [\"variance\"]])\n",
        "    ], df), 1)\n",
        "  }\n",
        "\n",
        "\n",
        "def _np_bedgraph_to_df(seq, species, rep=None, strand=None, state=None):\n",
        "  \"\"\"single vector representation of bedgraph pandas DataFrame\"\"\"\n",
        "  if seq is \"MS\":\n",
        "    DFs = _np_get_ms_bg(species, rep, state).groupby(\"seqid\")\n",
        "    if species != \"sc\": column = \"_\".join([\"MS\", str(rep)])\n",
        "    elif state == None: column = \"_\".join([\"MS\", \"wt\", str(rep)])\n",
        "    else: column = \"_\".join([\"MS\", \"ko\", state])\n",
        "  else: # seq is \"NS\"\n",
        "    DFs = _np_get_ns_bg(species, strand, state).groupby(\"seqid\")\n",
        "    if species != \"sc\": column = \"_\".join([\"NS\", strand])\n",
        "    elif state == None: column = \"_\".join([\"NS\", \"wt\", strand])\n",
        "    else: column = \"_\".join([\"NS\", \"ko\", state, strand])\n",
        "  return pd.concat((\n",
        "    pd.DataFrame({\n",
        "      column: np.concatenate(\n",
        "        [np.repeat(c, diff)\n",
        "          for (c, diff) in zip(\n",
        "            DFs.get_group(seqid)[\"count\"].astype(\"int\").values,\n",
        "            DFs.get_group(seqid)[\"b\"].values - DFs.get_group(seqid)[\"a\"].values\n",
        "        )])}).reset_index().drop(\"index\", axis=1) for seqid in DFs.groups\n",
        "    ), ignore_index=True).reset_index().drop(\"index\", axis=1)\n",
        "\n",
        "\n",
        "def _np_ds_meta_repr(data, columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"]):\n",
        "  \"\"\"construct A/T nucleotide meta representation\"\"\"\n",
        "  meta = np.average(_np_get_seq_data(columns, data).squeeze(), 0)\n",
        "  return 100 * np.sum([meta[0, :], meta[3, :]], 0)\n",
        "\n",
        "\n",
        "def _np_flatten(p, q, eps=(1 / (2 ** 18))):\n",
        "  \"\"\"convert (p, q) to computationally stable probability vectors\"\"\"\n",
        "  p = p.flatten(); q = q.flatten()\n",
        "  p /= np.sum(p) + eps; q /= np.sum(q) + eps\n",
        "  return p, q\n",
        "\n",
        "\n",
        "def _np_get_ds_fa(species):\n",
        "  \"\"\"one hot encode fasta --> DataFrame\"\"\"\n",
        "  return pd.concat((\n",
        "    _np_one_hot_encode_fa(\n",
        "      _dirs[\"raw\"] + _fastas[species][seqid]\n",
        "    ) for seqid in _fastas[species]\n",
        "  ), ignore_index=True).reset_index()\n",
        "\n",
        "\n",
        "def _np_get_gene_ID(attr_str):\n",
        "  \"\"\"get gene ID from attribute string\"\"\"\n",
        "  gene_ID = re.compile(r'ID=gene(?P<id>.+?);')\n",
        "  result = gene_ID.search(attr_str).group(\"id\")\n",
        "  return result if result[0] is not \":\" else result[1:]\n",
        "\n",
        "\n",
        "def _np_get_gene_biotype(attr_str):\n",
        "  \"\"\"get gene_biotype / biotype (product) for annotated genes\"\"\"\n",
        "  gene_biotype = re.compile(r\"biotype=(?P<biotype>.+?);\")\n",
        "  return gene_biotype.search(attr_str).group(\"biotype\")\n",
        "\n",
        "\n",
        "def _np_get_gene_Name(attr_str):\n",
        "  \"\"\"get gene Name from attribute string\"\"\"\n",
        "  try: return re.compile(\";Name=(?P<Name>.+?);\").search(attr_str).group(\"Name\")\n",
        "  except: return \"\"\n",
        "\n",
        "\n",
        "def _np_get_kl_data(columns, data):\n",
        "  \"\"\"get column-specific sequencing data\"\"\"\n",
        "  return np.rot90(\n",
        "    data[columns].values.reshape((-1, 1024, len(columns))), axes=(1, 2))\n",
        "\n",
        "\n",
        "def _np_get_locus_kl(seqid, start, end, strand, data, cols=[\n",
        "  \"position\", \"seqid\",\n",
        "  \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "  \"NS_pos\", \"NS_neg\",\n",
        "  \"MS_1\", \"MS_2\"]):\n",
        "  \"\"\"extract DNA sequence, NET-seq and MNase-seq data at given locus\"\"\"\n",
        "  return data[\n",
        "    (data[\"seqid\"] == seqid) &\n",
        "    (data[\"position\"] >= start) &\n",
        "    (data[\"position\"] < end)\n",
        "  ][cols].reset_index(drop=True) if strand is \"+\" else data[\n",
        "      (data[\"seqid\"] == seqid) &\n",
        "      (data[\"position\"] > start) &\n",
        "      (data[\"position\"] <= end)\n",
        "    ][cols].reset_index(drop=True)\n",
        "\n",
        "\n",
        "def _np_get_locus_sc(approach, seqid, start, end, strand, data,\n",
        "  cols={\n",
        "    \"AUG\": [\n",
        "      \"position\", \"seqid\",\n",
        "      \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "      \"NS_wt_pos\", \"NS_wt_neg\",\n",
        "      \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\"\n",
        "    ],\n",
        "    \"TSS\": [\n",
        "      \"position\", \"seqid\",\n",
        "      \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "      \"NS_wt_pos\", \"NS_wt_neg\",\n",
        "      \"NS_ko_chd1_pos\", \"NS_ko_chd1_neg\",\n",
        "      \"NS_ko_ino80_pos\", \"NS_ko_ino80_neg\",\n",
        "      \"NS_ko_isw1_pos\", \"NS_ko_isw1_neg\",\n",
        "      \"NS_ko_isw2_pos\", \"NS_ko_isw2_neg\",\n",
        "      \"NS_ko_rsc_pos\", \"NS_ko_rsc_neg\",\n",
        "      \"NS_ko_spt6-1004_pos\", \"NS_ko_spt6-1004_neg\",\n",
        "      \"NS_ko_spt6-YW_pos\", \"NS_ko_spt6-YW_neg\",\n",
        "      \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\",\n",
        "      \"MS_ko_spt6-1004\", \"MS_ko_spt6-YW\"\n",
        "    ]\n",
        "  }):\n",
        "  \"\"\"extract DNA sequence, NET-seq and MNase-seq data at given locus\"\"\"\n",
        "  return data[\n",
        "    (data[\"seqid\"] == seqid) &\n",
        "    (data[\"position\"] >= start) &\n",
        "    (data[\"position\"] < end)\n",
        "  ][cols[approach]].reset_index(drop=True) if strand is \"+\" else data[\n",
        "      (data[\"seqid\"] == seqid) &\n",
        "      (data[\"position\"] > start) &\n",
        "      (data[\"position\"] <= end)\n",
        "    ][cols[approach]].reset_index(drop=True)\n",
        "\n",
        "\n",
        "def _np_get_ms_bg(species, rep=None, state=None):\n",
        "  \"\"\"(species, (rep)licate, state)-specific MNase-seq bedgraph to DataFrame\"\"\"\n",
        "  if rep != None:\n",
        "    df = pd.read_csv(\n",
        "      _dirs[\"raw\"] + _bedgraphs[\"ms\"][species][\"wt\"][rep],\n",
        "      sep=\"\\t\",\n",
        "      header=None,\n",
        "      names=[\"seqid\", \"a\", \"b\", \"count\"],\n",
        "      comment=\"t\"\n",
        "    )\n",
        "  if state != None:\n",
        "    df = pd.read_csv(\n",
        "      _dirs[\"raw\"] + _bedgraphs[\"ms\"][species][\"ko\"][state],\n",
        "      sep=\"\\t\",\n",
        "      header=None,\n",
        "      names=[\"seqid\", \"a\", \"b\", \"count\"],\n",
        "      comment=\"t\"\n",
        "    )\n",
        "  df = df[df.seqid != \"chrM\"].reset_index().drop(\"index\", axis=1)\n",
        "  df[\"seqid\"] = df[\"seqid\"].map(_a2r)\n",
        "  return df.sort_values(by=[\"seqid\", \"a\", \"b\"])\n",
        "\n",
        "\n",
        "def _np_get_ns_bg(species, strand, state=None):\n",
        "  \"\"\"(species, strand)-specific NET-seq bedgraph to DataFrame\"\"\"\n",
        "  if species != \"sc\":\n",
        "    df = pd.read_csv(\n",
        "      _dirs[\"raw\"] + _bedgraphs[\"ns\"][species][strand],\n",
        "      sep=\"\\t\",\n",
        "      header=None,\n",
        "      names=[\"seqid\", \"a\", \"b\", \"count\"],\n",
        "      comment=\"t\"\n",
        "    )\n",
        "  elif state == None: # species == \"sc\", wildtype condition\n",
        "    df = pd.read_csv(\n",
        "      _dirs[\"raw\"] + _bedgraphs[\"ns\"][species][\"wt\"][strand],\n",
        "      sep=\"\\t\",\n",
        "      header=None,\n",
        "      names=[\"seqid\", \"a\", \"b\", \"count\"],\n",
        "      comment=\"t\"\n",
        "    )\n",
        "  else: # species == \"sc\", knockout condition\n",
        "    df = pd.read_csv(\n",
        "      _dirs[\"raw\"] + _bedgraphs[\"ns\"][species][\"ko\"][state][strand],\n",
        "      sep=\"\\t\",\n",
        "      header=None,\n",
        "      names=[\"seqid\", \"a\", \"b\", \"count\"],\n",
        "      comment=\"t\"\n",
        "    )\n",
        "  df = df[df.seqid != \"chrM\"].reset_index().drop(\"index\", axis=1)\n",
        "  df[\"seqid\"] = df[\"seqid\"].map(_a2r)\n",
        "  df_seqids = df.sort_values(by=[\"seqid\", \"a\", \"b\"]).groupby(\"seqid\")\n",
        "  return pd.concat(\n",
        "    objs=(pd.concat(\n",
        "      objs=(df_seqids.get_group(seqid),\n",
        "        pd.DataFrame({\n",
        "          \"seqid\": seqid,\n",
        "          \"a\": df_seqids.get_group(seqid).iloc[-1][\"b\"],\n",
        "          \"b\": _chrom_lengths[species][_a2r[seqid]], \n",
        "          \"count\": 0.0\n",
        "        }, index=np.arange(1))\n",
        "      ), ignore_index=True).reset_index().drop(\"index\", axis=1)\n",
        "        for seqid in df_seqids.groups),\n",
        "    ignore_index=True).reset_index().drop(\"index\", axis=1)\n",
        "\n",
        "\n",
        "def _np_get_seq_data(columns, data):\n",
        "  \"\"\"get column-specific sequencing data according to approach\"\"\"\n",
        "  return np.rot90(\n",
        "    data[columns].values.reshape((-1, 1024, len(columns))), axes=(1, 2))\n",
        "\n",
        "\n",
        "def _np_get_sc_data(columns, seqid, data):\n",
        "  \"\"\"get column-specific sequencing data according to seqid\"\"\"\n",
        "  return np.rot90(\n",
        "    data[data[\"seqid\"] == seqid][\n",
        "      columns\n",
        "    ].values.reshape((-1, 1024, len(columns))), axes=(1, 2))\n",
        "\n",
        "\n",
        "def _np_jsd(p, q, eps=1e-7):\n",
        "  \"\"\"< jensen-shannon divergence > - information radius\"\"\"\n",
        "  p = p / (np.sum(p) + eps)\n",
        "  q = q / (np.sum(q) + eps)\n",
        "  M = (p + q) / 2\n",
        "  return np.sum(special.kl_div(p, M) / 2 + special.kl_div(q, M) / 2)\n",
        "\n",
        "\n",
        "def _np_mean_ci(v, z=1.96):\n",
        "  \"\"\"(mean - conf. interval, mean, mean + conf. interval) of vector (v)\"\"\"\n",
        "  mean, std = np.average(v, 0), np.std(v, 0)\n",
        "  return {\n",
        "    \"lwr\": mean - z * std / np.sqrt(v.shape[0]),\n",
        "    \"cen\": mean,\n",
        "    \"upr\": mean + z * std / np.sqrt(v.shape[0])\n",
        "  }\n",
        "\n",
        "\n",
        "def _np_ms_bedgraph_to_df_kl(rep):\n",
        "  \"\"\"K. lactis MNase-seq bedgraph to collated pandas dataframe\"\"\"\n",
        "  return pd.read_csv(\n",
        "    _dirs[\"raw\"] + _bedgraphs[\"ms\"][\"kl\"],\n",
        "    names=[\"MS_\" + str(rep)]\n",
        "  ).astype(\"int\")\n",
        "\n",
        "\n",
        "def _np_ms_meta_extrema(v, width=4):\n",
        "  \"\"\"calculate relative extrema of vector\n",
        "  v := vector\n",
        "  width := \"peak\" width lwr-bound\"\"\"\n",
        "  return sig.find_peaks(x=v, width=width)[0]\n",
        "\n",
        "\n",
        "def _np_ms_meta_repr(vectors, eps=_heuristics[\"epsilon\"]):\n",
        "  \"\"\"return smoothed, normalized, z-score transformed avg MNase-seq profile\n",
        "    vectors := matrix ( # regions, 1024 )\"\"\"\n",
        "  # sum gaussian smoothed\n",
        "  smooth = np.sum([\n",
        "    _np_smooth_gaussian(i / (np.sum(i) + eps))\n",
        "      for i in vectors\n",
        "  ], 0)\n",
        "  # z-score normalize \n",
        "  z_norm = (smooth - np.average(smooth)) / np.std(smooth)\n",
        "  # normalize by magnitude\n",
        "  z_norm_s_norm = z_norm / (np.sum(np.abs(z_norm)) + eps)\n",
        "  # translate to positive domain\n",
        "  return z_norm_s_norm + np.abs(np.min(z_norm_s_norm))\n",
        "\n",
        "\n",
        "def _np_ms_transform(vector, eps=_heuristics[\"epsilon\"]):\n",
        "  \"\"\"return smoothed, normalized, z-score transformed avg MNase-seq vector\n",
        "    vector := ( 1024 )\"\"\"\n",
        "  # no signal\n",
        "  if np.sum(vector) == 0: return np.zeros(1024)\n",
        "  # sum gaussian smoothed, normalized vectors\n",
        "  smooth = _np_smooth_gaussian(vector / (np.sum(vector) + eps))\n",
        "  # z-score normalize \n",
        "  z_norm = (smooth - np.average(smooth)) / np.std(smooth)\n",
        "  # translate to positive domain\n",
        "  z_norm_pos = z_norm + np.abs(np.min(z_norm))\n",
        "  # normalize by magnitude\n",
        "  return z_norm_pos / (np.sum(np.abs(z_norm_pos)) + eps)\n",
        "\n",
        "\n",
        "def _np_ns_bedgraph_to_df_kl(_):\n",
        "  \"\"\"K. lactis NET-seq bedgraph to collated pandas dataframe\"\"\"\n",
        "  # read in NET-seq data to dataframe\n",
        "  kl_ns = pd.read_csv(\n",
        "    _dirs[\"raw\"] + _bedgraphs[\"ns\"][\"kl\"],\n",
        "    header=None, sep=\"\\t\", names=[\"seqid\", \"start\", \"end\", \"count\"]\n",
        "  )\n",
        "  # homogenize chromosome nomenclature\n",
        "  kl_ns[\"seqid\"] = kl_ns[\"seqid\"].map(_chrom_convert).map(_a2r)\n",
        "  # homogenize data types\n",
        "  kl_ns = kl_ns.astype(\"int\")\n",
        "  # instantiate NET-seq arrays\n",
        "  kl_ns_arr = {\n",
        "    seqid: {\n",
        "      \"pos\": np.zeros(_chrom_lengths[\"kl\"][_a2r[seqid]]),\n",
        "      \"neg\": np.zeros(_chrom_lengths[\"kl\"][_a2r[seqid]])\n",
        "    } for seqid in range(1, len(_chrom_lengths[\"kl\"]) + 1)\n",
        "  }\n",
        "  # assign read values to NET-seq arrays\n",
        "  for row in kl_ns.iterrows():\n",
        "    seqid = row[1][\"seqid\"]\n",
        "    start = row[1][\"start\"]\n",
        "    end = row[1][\"end\"]\n",
        "    count = row[1][\"count\"]\n",
        "    # NOTE: bedgraph from GEO has strands flipped...\n",
        "      # strand specific\n",
        "    if count > 0: kl_ns_arr[seqid][\"neg\"][start:end] += count\n",
        "    else: kl_ns_arr[seqid][\"pos\"][start:end] += abs(count)\n",
        "  # collate NET-seq arrays to dataframe format\n",
        "  return pd.DataFrame({\n",
        "    \"NS_pos\": np.concatenate([\n",
        "      kl_ns_arr[seqid][\"pos\"]\n",
        "        for seqid in kl_ns_arr\n",
        "    ]).astype(\"int\"),\n",
        "    \"NS_neg\": np.concatenate([\n",
        "      kl_ns_arr[seqid][\"neg\"]\n",
        "        for seqid in kl_ns_arr\n",
        "    ]).astype(\"int\")\n",
        "  })\n",
        "\n",
        "\n",
        "def _np_ns_transform(vector, w=4, eps=_heuristics[\"epsilon\"]):\n",
        "  \"\"\"return smoothed NET-seq vector (same format as FIDDLE's conv. input)\n",
        "    w := 4 base pairs, FIDDLE conv. filter size\n",
        "    vector := ( 1024 )\"\"\"\n",
        "  if np.sum(vector) == 0: return np.zeros(1024) # no signal\n",
        "  return _np_smooth_gaussian(vector, w=w)\n",
        "\n",
        "\n",
        "def _np_one_hot_encode_fa(fa):\n",
        "  \"\"\"fasta --> string --> numpy array, 4 x sequence length\"\"\"\n",
        "  with open(fa, \"r\") as f:\n",
        "    # assumes single header line exists\n",
        "    _ = f.readline()\n",
        "    str_DNA = \"\".join(f.read().split(\"\\n\")).lower()\n",
        "    mapping = {\n",
        "      \"a\": [1, 0, 0, 0],\n",
        "      \"c\": [0, 1, 0, 0],\n",
        "      \"g\": [0, 0, 1, 0],\n",
        "      \"t\": [0, 0, 0, 1],\n",
        "      \"n\": [0.25] * 4\n",
        "    }\n",
        "    return pd.DataFrame(\n",
        "      np.array([mapping[nt] for nt in str_DNA]),\n",
        "      columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"])\n",
        "\n",
        "\n",
        "def _np_random_start_end(start, end, strand, width):\n",
        "  if strand == \"+\":\n",
        "    start = np.random.randint(start, start + width)\n",
        "    return start, start + width\n",
        "  else: # strand == \"-\"\n",
        "    end = np.random.randint(end - width, end)\n",
        "    return end - width, end\n",
        "\n",
        "\n",
        "def _np_smooth_gaussian(x, mode=\"same\", w=128):\n",
        "  \"\"\"smooth signal with gaussian window\"\"\"\n",
        "  vector = np.r_[x[(w - 1):0:-1], x, x[-2:-(w - 1):-1]]\n",
        "  conv = sig.windows.gaussian(w, w / 6)\n",
        "  return np.convolve(conv / conv.sum(), vector, mode=mode)[(w - 2):-(w - 2)]\n",
        "\n",
        "\n",
        "def _np_shuffle(M, shuff):\n",
        "  \"\"\"shuffle data if indicated\"\"\"\n",
        "  if shuff: return M[:, :, np.random.permutation(np.arange(M.shape[2])), :]\n",
        "  else: return M\n",
        "\n",
        "\n",
        "def _np_split(cols, seqids, split, data, w=1024):\n",
        "  \"\"\"split \"cols\"-specified data given seqids and split\"\"\"\n",
        "  if split == \"train\":\n",
        "    return np.expand_dims(np.rot90(\n",
        "      data[cols][\n",
        "        (data[\"approach\"] == \"AUG\") |\n",
        "        (data[\"seqid\"] == seqids[0]) |\n",
        "        (data[\"seqid\"] == seqids[1]) |\n",
        "        (data[\"seqid\"] == seqids[2]) |\n",
        "        (data[\"seqid\"] == seqids[3]) |\n",
        "        (data[\"seqid\"] == seqids[4]) |\n",
        "        (data[\"seqid\"] == seqids[5]) |\n",
        "        (data[\"seqid\"] == seqids[6]) |\n",
        "        (data[\"seqid\"] == seqids[7]) |\n",
        "        (data[\"seqid\"] == seqids[8]) |\n",
        "        (data[\"seqid\"] == seqids[9]) |\n",
        "        (data[\"seqid\"] == seqids[10]) |\n",
        "        (data[\"seqid\"] == seqids[11]) |\n",
        "        (data[\"seqid\"] == seqids[12]) |\n",
        "        (data[\"seqid\"] == seqids[13])\n",
        "      ].values.reshape((-1, w, len(cols))), axes=(1, 2)), axis=3)\n",
        "  else: # split == \"valid\" | split == \"test\"\n",
        "    return np.expand_dims(np.rot90(\n",
        "      data[cols][\n",
        "        (data[\"approach\"] == \"TSS\") &\n",
        "        (data[\"seqid\"] == seqids)\n",
        "      ].values.reshape((-1, w, len(cols))), axes=(1, 2)), axis=3)\n",
        "\n",
        "\n",
        "print(\"custom functions are loaded\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b7wyKDcJP31y",
        "colab_type": "text"
      },
      "source": [
        "# COLLATE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eYT-RNR_gGGL",
        "colab_type": "text"
      },
      "source": [
        "## S. _cerevisiae_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tGklhEqLdkl2",
        "colab_type": "text"
      },
      "source": [
        "### load / munge meta & TSS annotations\n",
        "\n",
        "1. grab protein-coding genes from GTF [here](http://fungi.ensembl.org/Saccharomyces_cerevisiae/Info/Index).\n",
        "2. combine with annotated transcription start sites (TSSs) [here](https://github.com/winston-lab/genomefiles-cerevisiae/tree/master/annotations)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-YAz-2NvP49p",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# load genome annotation file\n",
        "meta_sc = pd.read_csv(\n",
        "  os.path.join(_dirs[\"raw\"], _annotations[\"sc\"][\"gff\"]),\n",
        "  sep=\"\\t\",\n",
        "  comment=\"#\",\n",
        "  header=None,\n",
        "  names=[\n",
        "    \"seqid\",\n",
        "    \"source\",\n",
        "    \"type\",\n",
        "    \"start\",\n",
        "    \"end\",\n",
        "    \"score\",\n",
        "    \"strand\",                     \n",
        "    \"phase\",\n",
        "    \"attributes\"\n",
        "  ],\n",
        "  engine=\"python\"\n",
        ")\n",
        "# grab genes\n",
        "meta_sc = meta_sc[meta_sc.type == \"gene\"]\n",
        "# homogenize chromosome nomenclature\n",
        "meta_sc[\"seqid\"] = meta_sc[\"seqid\"].map(_chrom_convert)\n",
        "# remove unwanted chromosomes\n",
        "meta_sc = meta_sc[\n",
        "  meta_sc.seqid != \"remove\"\n",
        "].reset_index(drop=True)\n",
        "# gene ID from attributes\n",
        "meta_sc[\"ID\"] = meta_sc[\"attributes\"].apply(_np_get_gene_ID)\n",
        "# gene biotype from attributes\n",
        "meta_sc[\"biotype\"] = meta_sc[\"attributes\"].apply(_np_get_gene_biotype)\n",
        "# only protein coding genes\n",
        "meta_sc = meta_sc[meta_sc.biotype == \"protein_coding\"]\n",
        "# gene Name from attributes\n",
        "meta_sc[\"Name\"] = meta_sc[\"attributes\"].apply(_np_get_gene_Name)\n",
        "# zero-index all positions (python format)\n",
        "meta_sc[\"start\"] = meta_sc[\"start\"].apply(lambda x: x - 1)\n",
        "meta_sc[\"end\"] = meta_sc[\"end\"].apply(lambda x: x - 1)\n",
        "# read in Transcription Start Sites\n",
        "TSSs_sc = pd.read_csv(\n",
        "  os.path.join(_dirs[\"raw\"], _annotations[\"sc\"][\"tss\"]),\n",
        "  header=None,\n",
        "  sep=\"\\t\",\n",
        "  names=[\n",
        "    \"seqid\",\n",
        "    \"start\",\n",
        "    \"end\",\n",
        "    \"name\",\n",
        "    \"count\",\n",
        "    \"strand\"\n",
        "  ]\n",
        ")\n",
        "# remove unwanted chromosomes\n",
        "TSSs_sc = TSSs_sc[TSSs_sc[\"seqid\"] != \"chrM\"]\n",
        "# homogenize chromosome nomenclature\n",
        "_ph = TSSs_sc[\"seqid\"].map(_a2r)\n",
        "TSSs_sc[\"seqid\"] = _ph; del _ph\n",
        "TSSs_sc.sort_values([\"seqid\", \"start\", \"end\"], inplace=True)\n",
        "# assign TSS column to meta DataFrame\n",
        "tss_locs_sc = []\n",
        "for row in meta_sc.iterrows():\n",
        "  # from meta_sc\n",
        "  ID, Name = row[1][\"ID\"], row[1][\"Name\"]\n",
        "  strand = row[1][\"strand\"]\n",
        "  start, end = row[1][\"start\"], row[1][\"end\"]\n",
        "  # from TSS annotations\n",
        "  _id = TSSs_sc[TSSs_sc[\"name\"] == ID]\n",
        "  _name = TSSs_sc[TSSs_sc[\"name\"] == Name]\n",
        "  # gene in meta_sc but not in TSSs_sc\n",
        "  if len(_id) + len(_name) == 0:\n",
        "    tss_locs_sc.append(start) if strand is \"+\" \\\n",
        "      else tss_locs_sc.append(end)\n",
        "  # gene in meta_sc and TSSs_sc\n",
        "  else:\n",
        "    tss_source = _id if len(_id) == 1 else _name\n",
        "    tss_locs_sc.append(tss_source[\"start\"].values[0]) if strand is \"+\" \\\n",
        "      else tss_locs_sc.append(tss_source[\"end\"].values[0])\n",
        "# append TSS to meta_sc DataFrame\n",
        "meta_sc[\"TSS\"] = np.array(tss_locs_sc)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7_dw__HfP4kK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "meta_sc.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PyyE4IMtP4iM",
        "colab_type": "text"
      },
      "source": [
        "### define regions of consideration\n",
        "\n",
        "**_training regions semi-phased wrt AUG_**\n",
        "\n",
        "**_validation / testing regions phased wrt TSS_**\n",
        "\n",
        "1. training regions sampled from 1,024 * 2 bp (2,048 bp) regions centered on proxy TSS (64 bp upstream AUG)\n",
        "2. validation & testing regions are 1,024 bp-widthed windows centered on annotated TSS\n",
        "3. define out of bounds regions (any boundaries below zero or above chromosome length)\n",
        "4. define overlapping regions (any validation & testing regions with overlapping boundaries)\n",
        "5. clean-up dataframe (remove unnecessary columns / rows with out of bounds or overlapping regions)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q2jFwDVmP4gS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# training loci boundaries phased wrt AUG\n",
        "_meta_sc_strands = meta_sc.groupby(\"strand\")\n",
        "meta_sc = pd.concat((\n",
        "  _meta_sc_strands.get_group(\"+\").assign(\n",
        "    AUG_start = lambda x: \\\n",
        "      x.start - _heuristics[\"AUG_us_offset\"] - _heuristics[\"AUG_width\"] // 2\n",
        "        ).assign(AUG_end = lambda x: \\\n",
        "      x.start - _heuristics[\"AUG_us_offset\"] + _heuristics[\"AUG_width\"] // 2),\n",
        "  _meta_sc_strands.get_group(\"-\").assign(\n",
        "    AUG_start = lambda x: \\\n",
        "      x.end + _heuristics[\"AUG_us_offset\"] - _heuristics[\"AUG_width\"] // 2\n",
        "        ).assign(AUG_end = lambda x: \\\n",
        "      x.end + _heuristics[\"AUG_us_offset\"] + _heuristics[\"AUG_width\"] // 2)),\n",
        "  ignore_index=True).reset_index().drop(\"index\", axis=1)\n",
        "# validation / testing loci phased wrt TSS\n",
        "meta_sc = meta_sc.assign(\n",
        "  TSS_start = lambda x: x.TSS - _heuristics[\"TSS_width\"] // 2).assign(\n",
        "  TSS_end = lambda x: x.TSS + _heuristics[\"TSS_width\"] // 2)\n",
        "# define out of bounds (oob) loci wrt chromosome boundaries\n",
        "_seqids_sc = meta_sc.groupby(\"seqid\")\n",
        "oobs_sc = []\n",
        "for seqid in _chrom_lengths[\"sc\"]:\n",
        "  # oob := any (start, end), negative or greater than seqid length\n",
        "  oobs_sc.append(_seqids_sc.get_group(seqid).assign(\n",
        "    oob = \\\n",
        "      (_seqids_sc.get_group(seqid).TSS_start < 0) |\n",
        "      (_seqids_sc.get_group(seqid).TSS_end < 0) |\n",
        "      (_seqids_sc.get_group(seqid).AUG_start < 0) |\n",
        "      (_seqids_sc.get_group(seqid).AUG_end < 0) |\n",
        "      (_seqids_sc.get_group(seqid).TSS_start > _chrom_lengths[\"sc\"][seqid]) |\n",
        "      (_seqids_sc.get_group(seqid).TSS_end > _chrom_lengths[\"sc\"][seqid]) |\n",
        "      (_seqids_sc.get_group(seqid).AUG_start > _chrom_lengths[\"sc\"][seqid]) |\n",
        "      (_seqids_sc.get_group(seqid).AUG_end > _chrom_lengths[\"sc\"][seqid])\n",
        "  ))\n",
        "# concatenate to single DataFrame\n",
        "meta_sc = pd.concat(oobs_sc, ignore_index=True)\n",
        "# define overlapping (ol) loci wrt TSS_start / TSS_end\n",
        "  # overlapping base pairs (integer)\n",
        "_ol_int_sc = {\n",
        "  seqid: np.zeros(_chrom_lengths[\"sc\"][seqid])\n",
        "    for seqid in _chrom_lengths[\"sc\"]\n",
        "}\n",
        "  # overlapping loci (boolean)\n",
        "_ol_bool_sc = np.zeros(len(meta_sc)).astype(\"bool\")\n",
        "# number of loci at each base pair of each seqid\n",
        "for x in range(len(meta_sc)):\n",
        "  _ol_int_sc[meta_sc[\"seqid\"].iat[x]][\n",
        "    meta_sc[\"TSS_start\"].iat[x]:meta_sc[\"TSS_end\"].iat[x]\n",
        "  ] += 1\n",
        "# loci is overlapping w/ another loci\n",
        "for x in range(len(meta_sc)):\n",
        "  if any(_ol_int_sc[meta_sc[\"seqid\"].iat[x]][\n",
        "    meta_sc[\"TSS_start\"].iat[x]:meta_sc[\"TSS_end\"].iat[x]\n",
        "  ] > 1): _ol_bool_sc[x] = True\n",
        "# assign\n",
        "meta_sc[\"ol\"] = _ol_bool_sc\n",
        "# clean DataFrame - remove unnecessary columns, oob and ol genes\n",
        "meta_sc = meta_sc[\n",
        "  (meta_sc.oob == False) & \n",
        "  (meta_sc.ol == False)\n",
        "]\n",
        "meta_sc.drop(\n",
        "  labels=[\n",
        "    \"attributes\", \"biotype\", \"phase\", \"score\", \"source\", \"type\", \"ol\", \"oob\",\n",
        "    \"Name\"\n",
        "  ], axis=1, inplace=True\n",
        ")\n",
        "# homogenize chromosome nomenclature to integer format\n",
        "meta_sc[\"seqid\"] = meta_sc[\"seqid\"].map(_a2r)\n",
        "# reset index\n",
        "meta_sc = meta_sc.reset_index().drop(\"index\", axis=1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8L8T7s9gP4Xn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "meta_sc.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4kX4qKe_Pfon",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# save\n",
        "meta_sc_path = _dirs[\"collated\"] + \"meta_sc.csv\"\n",
        "if os.path.isfile(meta_sc_path):\n",
        "  if input(f\"overwrite {meta_sc_path}? ... yes, no\\n\") == \"yes\":\n",
        "    print(\"overwriting old meta_sc with new meta_sc\")\n",
        "    meta_sc.to_csv(meta_sc_path)\n",
        "else:\n",
        "  print(\"saving meta_sc\")\n",
        "  meta_sc.to_csv(meta_sc_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qARBUY0Sn8Xh",
        "colab_type": "text"
      },
      "source": [
        "### define genomic & sequencing data\n",
        "\n",
        "1. read in raw genomic data from fastas into dataframes\n",
        "2. read in sequencing data from bedgraphs into dataframes\n",
        "3. combine"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I5HZS3FErU5U",
        "colab_type": "code",
        "collapsed": true,
        "colab": {}
      },
      "source": [
        "# load data if previously checkpointed\n",
        "meta_sc = pd.read_csv(os.path.join(_dirs[\"collated\"], \"meta_sc.csv\"))\n",
        "nuisance = \"Unnamed: 0\"\n",
        "if nuisance in meta_sc.columns.values:\n",
        "  meta_sc.drop(nuisance, axis=1, inplace=True)\n",
        "meta_sc.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w13dnuwPppls",
        "colab_type": "code",
        "collapsed": true,
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# collate to single DataFrame\n",
        "raw_genomic_sc = pd.concat((\n",
        "  pd.DataFrame({\n",
        "    \"seqid\": np.concatenate([\n",
        "      np.repeat(_a2r[seqid], _chrom_lengths[\"sc\"][seqid])\n",
        "        for seqid in _chrom_lengths[\"sc\"]]),\n",
        "    \"position\": np.concatenate([\n",
        "      np.arange(_chrom_lengths[\"sc\"][seqid])\n",
        "        for seqid in _chrom_lengths[\"sc\"]])\n",
        "  }), _np_get_ds_fa(\"sc\")),\n",
        "  axis=1).reset_index().drop(\"index\", axis=1).drop(\"level_0\", axis=1)\n",
        "# combine genomic data with sequencing data\n",
        "gen_seq_sc = pd.concat((\n",
        "  raw_genomic_sc,\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"chd1\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"chd1\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"ino80\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"ino80\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"isw1\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"isw1\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"isw2\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"isw2\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"rsc\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"rsc\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"spt6-1004\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"spt6-1004\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"pos\", state=\"spt6-YW\"),\n",
        "  _np_bedgraph_to_df(\"NS\", \"sc\", strand=\"neg\", state=\"spt6-YW\"),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", rep=1),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", rep=2),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", rep=3),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", rep=4),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", state=\"spt6-1004\"),\n",
        "  _np_bedgraph_to_df(\"MS\", \"sc\", state=\"spt6-YW\")\n",
        "  ), axis=1).reset_index().drop(\"index\", axis=1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TO8H52HJ8MkK",
        "colab_type": "code",
        "collapsed": true,
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "gen_seq_sc.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vC3GsLB7Cyc5",
        "colab_type": "code",
        "collapsed": true,
        "colab": {}
      },
      "source": [
        "# save\n",
        "gen_seq_sc_path = _dirs[\"collated\"] + \"gen_seq_sc.csv\"\n",
        "if os.path.isfile(gen_seq_sc_path):\n",
        "  if input(f\"overwrite {gen_seq_sc_path}? ... yes, no\\n\") == \"yes\":\n",
        "    print(\"overwriting old gen_seq_sc with new gen_seq_sc\")\n",
        "    gen_seq_sc.to_csv(gen_seq_sc_path)\n",
        "else:\n",
        "  print(\"saving gen_seq_sc\")\n",
        "  gen_seq_sc.to_csv(gen_seq_sc_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ldOTIObHi_A9",
        "colab_type": "text"
      },
      "source": [
        "## K. _lactis_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "xzzTYPYJjP3e"
      },
      "source": [
        "### load / munge meta & TSS annotations\n",
        "\n",
        "1. grab protein-coding genes from GTF [here](https://fungi.ensembl.org/Kluyveromyces_lactis_gca_000002515/Info/Index)\n",
        "2. combine with annotated transcription start sites (TSSs) [here](http://www.yeastss.org/download/)\n",
        "3. define non-overlapping promoter regions wr.t. each gene's AUGf\n",
        "4. define TSS - locating region boundaries wr.t. each gene's AUG\n",
        "5. assign each gene's TSS from raw TSS-seq signal using approach similar to published methods"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "K-r34VXDjP3f",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# load genome annotation file\n",
        "meta_kl = pd.read_csv(\n",
        "  os.path.join(_dirs[\"raw\"], _annotations[\"kl\"][\"gff\"]),\n",
        "  sep=\"\\t\",\n",
        "  comment=\"#\",\n",
        "  header=None,\n",
        "  names=[\n",
        "    \"seqid\",\n",
        "    \"source\",\n",
        "    \"type\",\n",
        "    \"start\",\n",
        "    \"end\",\n",
        "    \"score\",\n",
        "    \"strand\",                     \n",
        "    \"phase\",\n",
        "    \"attributes\"\n",
        "  ],\n",
        "  engine=\"python\"\n",
        ")\n",
        "# grab genes\n",
        "meta_kl = meta_kl[meta_kl.type == \"gene\"]\n",
        "# homogenize chromosome nomenclature\n",
        "meta_kl[\"seqid\"] = meta_kl[\"seqid\"].map(_chrom_convert)\n",
        "# remove unwanted chromosomes\n",
        "meta_kl = meta_kl[\n",
        "  meta_kl.seqid != \"remove\"\n",
        "].reset_index(drop=True)\n",
        "# gene ID from attributes\n",
        "meta_kl[\"ID\"] = meta_kl[\"attributes\"].apply(_np_get_gene_ID)\n",
        "# gene biotype from attributes\n",
        "meta_kl[\"biotype\"] = meta_kl[\"attributes\"].apply(_np_get_gene_biotype)\n",
        "# only protein coding genes\n",
        "meta_kl = meta_kl[meta_kl.biotype == \"protein_coding\"]\n",
        "# gene Name from attributes\n",
        "meta_kl[\"Name\"] = meta_kl[\"attributes\"].apply(_np_get_gene_Name)\n",
        "# zero-index all positions (python format)\n",
        "meta_kl[\"start\"] = meta_kl[\"start\"].apply(lambda x: x - 1)\n",
        "meta_kl[\"end\"] = meta_kl[\"end\"].apply(lambda x: x - 1)\n",
        "# loci boundaries defined wrt AUG\n",
        "_meta_kl_strands = meta_kl.groupby(\"strand\")\n",
        "meta_kl = pd.concat((\n",
        "  _meta_kl_strands.get_group(\"+\").assign(\n",
        "    ol_reg_start = lambda x: x.start - _heuristics[\"ol_reg_width\"]).assign(\n",
        "    ol_reg_end = lambda x: x.start + _heuristics[\"ol_reg_width\"]).assign(\n",
        "    find_tss_start = lambda x: x.start - _heuristics[\"ol_reg_width\"] // 2).assign(\n",
        "    find_tss_end = lambda x: x.start + _heuristics[\"ol_reg_width\"] // 2),\n",
        "  _meta_kl_strands.get_group(\"-\").assign(\n",
        "    ol_reg_start = lambda x: x.end - _heuristics[\"ol_reg_width\"]).assign(\n",
        "    ol_reg_end = lambda x: x.end + _heuristics[\"ol_reg_width\"]).assign(\n",
        "    find_tss_start = lambda x: x.end - _heuristics[\"ol_reg_width\"] // 2).assign(\n",
        "    find_tss_end = lambda x: x.end + _heuristics[\"ol_reg_width\"] // 2),\n",
        "  ), ignore_index=True).reset_index().drop(\"index\", axis=1)\n",
        "# define out of bounds (oob) loci wrt chromosome boundaries\n",
        "_seqids_kl = meta_kl.groupby(\"seqid\")\n",
        "oobs_kl = []\n",
        "for seqid in _chrom_lengths[\"kl\"]:\n",
        "  # oob := any (start, end), negative or greater than seqid length\n",
        "  oobs_kl.append(_seqids_kl.get_group(seqid).assign(\n",
        "    oob = \\\n",
        "      (_seqids_kl.get_group(seqid).ol_reg_start < 0) |\n",
        "      (_seqids_kl.get_group(seqid).ol_reg_end < 0) |\n",
        "      (_seqids_kl.get_group(seqid).ol_reg_start > _chrom_lengths[\"kl\"][seqid]) |\n",
        "      (_seqids_kl.get_group(seqid).ol_reg_end > _chrom_lengths[\"kl\"][seqid])\n",
        "  ))\n",
        "# concatenate to single DataFrame\n",
        "meta_kl = pd.concat(oobs_kl, ignore_index=True)\n",
        "# define overlapping (ol) loci wrt AUG_start / AUG_end\n",
        "  # overlapping base pairs (integer)\n",
        "_ol_int_kl = {\n",
        "  seqid: np.zeros(_chrom_lengths[\"kl\"][seqid])\n",
        "    for seqid in _chrom_lengths[\"kl\"]\n",
        "}\n",
        "  # overlapping loci (boolean)\n",
        "_ol_bool_kl = np.zeros(len(meta_kl)).astype(\"bool\")\n",
        "# number of loci at each base pair of each seqid\n",
        "for x in range(len(meta_kl)):\n",
        "  _ol_int_kl[meta_kl[\"seqid\"].iat[x]][\n",
        "    meta_kl[\"ol_reg_start\"].iat[x]:meta_kl[\"ol_reg_end\"].iat[x]\n",
        "  ] += 1\n",
        "# loci is overlapping w/ another loci\n",
        "for x in range(len(meta_kl)):\n",
        "  if any(_ol_int_kl[meta_kl[\"seqid\"].iat[x]][\n",
        "    meta_kl[\"ol_reg_start\"].iat[x]:meta_kl[\"ol_reg_end\"].iat[x]\n",
        "  ] > 1): _ol_bool_kl[x] = True\n",
        "# assign\n",
        "meta_kl[\"ol\"] = _ol_bool_kl\n",
        "# clean DataFrame - remove unnecessary columns, oob and ol genes\n",
        "meta_kl = meta_kl[\n",
        "  (meta_kl.oob == False) & \n",
        "  (meta_kl.ol == False)\n",
        "]\n",
        "meta_kl.drop(\n",
        "  labels=[\n",
        "    \"attributes\", \"biotype\", \"phase\", \"score\", \"source\", \"type\", \"ol\", \"oob\",\n",
        "    \"Name\"\n",
        "  ], axis=1, inplace=True\n",
        ")\n",
        "# homogenize chromosome nomenclature to integer format\n",
        "meta_kl[\"seqid\"] = meta_kl[\"seqid\"].map(_a2r)\n",
        "# reset index\n",
        "meta_kl.reset_index(drop=True, inplace=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cgLUdswabIT9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# read in raw TSS-seq data to dataframe\n",
        "TSSs_kl = {\n",
        "  1: pd.read_csv(\n",
        "    os.path.join(_dirs[\"raw\"], _annotations[\"kl\"][\"tss\"][1]),\n",
        "    sep=\"\\t\", header=None, names=[\"seqid\", \"position\", \"strand\", \"count\"]\n",
        "  ),\n",
        "  2: pd.read_csv(\n",
        "    os.path.join(_dirs[\"raw\"], _annotations[\"kl\"][\"tss\"][2]),\n",
        "    sep=\"\\t\", header=None, names=[\"seqid\", \"position\", \"strand\", \"count\"]\n",
        "  )\n",
        "}\n",
        "# homogenize chromosome nomenclature, remove unwanted chromosomes\n",
        "for rep in [1, 2]:\n",
        "  TSSs_kl[rep][\"seqid\"] = TSSs_kl[rep][\"seqid\"].map(_chrom_convert)\n",
        "  TSSs_kl[rep] = TSSs_kl[rep][\n",
        "    TSSs_kl[rep].seqid != \"remove\"\n",
        "  ].reset_index(drop=True)\n",
        "# convert TSS-seq dataframe to TSS-seq array\n",
        "  # instantiate arrays\n",
        "_TSSs_kl_arr = {\n",
        "  rep: {\n",
        "    seqid: {\n",
        "      \"pos\": np.zeros(_chrom_lengths[\"kl\"][seqid]),\n",
        "      \"neg\": np.zeros(_chrom_lengths[\"kl\"][seqid])\n",
        "    } for seqid in _chrom_lengths[\"kl\"]\n",
        "  } for rep in [1, 2]\n",
        "}\n",
        "  # assign counts to arrays\n",
        "for rep in [1, 2]:\n",
        "  for row in TSSs_kl[rep].iterrows():\n",
        "    _TSSs_kl_arr[rep][\n",
        "      row[1][\"seqid\"]][\n",
        "        _strandedness[row[1][\"strand\"]]][\n",
        "          int(row[1][\"position\"])\n",
        "    ] += int(row[1][\"count\"])\n",
        "  # combined\n",
        "TSSs_kl_arr = {\n",
        "  seqid: {\n",
        "    \"pos\": \\\n",
        "      _TSSs_kl_arr[1][seqid][\"pos\"] + _TSSs_kl_arr[2][seqid][\"pos\"],\n",
        "    \"neg\": \\\n",
        "      _TSSs_kl_arr[1][seqid][\"neg\"] + _TSSs_kl_arr[2][seqid][\"neg\"]\n",
        "  } for seqid in _chrom_lengths[\"kl\"]\n",
        "}\n",
        "  # threshold arrays\n",
        "for seqid in _chrom_lengths[\"kl\"]:\n",
        "  for strand in [\"pos\", \"neg\"]:\n",
        "    TSSs_kl_arr[seqid][strand][\n",
        "      TSSs_kl_arr[seqid][strand] <\n",
        "        np.percentile(\n",
        "          a=TSSs_kl_arr[seqid][strand],\n",
        "          q=_heuristics[\"TSS_seq_signal_percentile\"]\n",
        "        )] = 0\n",
        "# obtain TSS for each gene\n",
        "TSS_annot = []\n",
        "for row in meta_kl.iterrows():  \n",
        "  # get data\n",
        "  seqid, strand = _a2r[row[1][\"seqid\"]], row[1][\"strand\"]\n",
        "  start, end = int(row[1][\"start\"]), int(row[1][\"end\"])\n",
        "  find_tss_start = int(row[1][\"find_tss_start\"])\n",
        "  find_tss_end = int(row[1][\"find_tss_end\"])\n",
        "  # find TSS-seq maxima\n",
        "  TSS = np.argmax(\n",
        "    TSSs_kl_arr[seqid][_strandedness[strand]][\n",
        "      find_tss_start:find_tss_end\n",
        "    ])\n",
        "  if (strand == \"+\") and (TSS == 0): TSS = start\n",
        "  elif (strand == \"-\") and (TSS == 0): TSS = end\n",
        "  else: TSS += find_tss_start\n",
        "  # append\n",
        "  TSS_annot.append(TSS)\n",
        "# append TSS annotation to meta dataframe\n",
        "meta_kl[\"TSS\"] = np.array(TSS_annot)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XS2NTeLYXNfz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "meta_kl.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "88e0wDN4XSTw",
        "colab_type": "text"
      },
      "source": [
        "### define regions of consideration\n",
        "\n",
        "**_1,024 bp-widthed windows centered on annotated TSS_**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PhzIGlk4XSlP",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# drop unnecessary columns\n",
        "meta_kl.drop(\n",
        "  [\"ol_reg_start\", \"ol_reg_end\", \"find_tss_start\", \"find_tss_end\"],\n",
        "  axis=1, inplace=True\n",
        ")\n",
        "# homogenize chromosome nomenclature\n",
        "meta_kl[\"seqid\"] = meta_kl[\"seqid\"].map(_a2r)\n",
        "# define testing region indices\n",
        "meta_kl = meta_kl.assign(\n",
        "  TSS_start = lambda x: x.TSS - _heuristics[\"TSS_width\"] // 2).assign(\n",
        "  TSS_end = lambda x: x.TSS + _heuristics[\"TSS_width\"] // 2)\n",
        "# define out of bounds (oob) loci wrt chromosome boundaries\n",
        "_seqids_kl = meta_kl.groupby(\"seqid\")\n",
        "oobs_kl = []\n",
        "for seqid in _chrom_lengths[\"kl\"]:\n",
        "  # oob := any (start, end), negative or greater than seqid length\n",
        "  oobs_kl.append(_seqids_kl.get_group(seqid).assign(\n",
        "    oob = \\\n",
        "      (_seqids_kl.get_group(seqid).TSS_start < 0) |\n",
        "      (_seqids_kl.get_group(seqid).TSS_end < 0) |\n",
        "      (_seqids_kl.get_group(seqid).TSS_start > _chrom_lengths[\"kl\"][seqid]) |\n",
        "      (_seqids_kl.get_group(seqid).TSS_end > _chrom_lengths[\"kl\"][seqid])\n",
        "  ))\n",
        "# concatenate to single DataFrame\n",
        "meta_kl = pd.concat(oobs_kl, ignore_index=True)\n",
        "# define overlapping (ol) loci wrt TSS_start / TSS_end\n",
        "  # overlapping base pairs (integer)\n",
        "_ol_int_kl = {\n",
        "  seqid: np.zeros(_chrom_lengths[\"kl\"][seqid])\n",
        "    for seqid in _chrom_lengths[\"kl\"]\n",
        "}\n",
        "  # overlapping loci (boolean)\n",
        "_ol_bool_kl = np.zeros(len(meta_kl)).astype(\"bool\")\n",
        "# number of loci at each base pair of each seqid\n",
        "for x in range(len(meta_kl)):\n",
        "  _ol_int_kl[meta_kl[\"seqid\"].iat[x]][\n",
        "    meta_kl[\"TSS_start\"].iat[x]:meta_kl[\"TSS_end\"].iat[x]\n",
        "  ] += 1\n",
        "# loci is overlapping w/ another loci\n",
        "for x in range(len(meta_kl)):\n",
        "  if any(_ol_int_kl[meta_kl[\"seqid\"].iat[x]][\n",
        "    meta_kl[\"TSS_start\"].iat[x]:meta_kl[\"TSS_end\"].iat[x]\n",
        "  ] > 1): _ol_bool_kl[x] = True\n",
        "# assign\n",
        "meta_kl[\"ol\"] = _ol_bool_kl\n",
        "# clean DataFrame - remove unnecessary columns, oob and ol genes\n",
        "meta_kl = meta_kl[\n",
        "  (meta_kl.oob == False) & \n",
        "  (meta_kl.ol == False)\n",
        "]\n",
        "meta_kl.drop(labels=[\"oob\", \"ol\"], axis=1, inplace=True)\n",
        "# homogenize chromosome nomenclature to integer format\n",
        "meta_kl[\"seqid\"] = meta_kl[\"seqid\"].map(_a2r)\n",
        "# reset index\n",
        "meta_kl = meta_kl.reset_index().drop(\"index\", axis=1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QwmMVDqfee6p",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# save\n",
        "meta_kl_path = _dirs[\"collated\"] + \"meta_kl.csv\"\n",
        "if os.path.isfile(meta_kl_path):\n",
        "  if input(f\"overwrite {meta_kl_path}? ... yes, no\\n\") == \"yes\":\n",
        "    print(\"overwriting old meta_kl with new meta_kl\")\n",
        "    meta_kl.to_csv(meta_kl_path)\n",
        "else:\n",
        "  print(\"saving meta_kl\")\n",
        "  meta_kl.to_csv(meta_kl_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "x3-4W64keKrx"
      },
      "source": [
        "### define genomic & sequencing data\n",
        "\n",
        "1. read in meta information of loci of interest \n",
        "2. read in raw genomic data from fastas into dataframes\n",
        "3. combine raw genomic dataframe with sequencing data from bedgraphs"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "7UsqSoCZeKry",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# collate to single DataFrame\n",
        "raw_genomic_kl = pd.concat((\n",
        "  pd.DataFrame({\n",
        "    \"seqid\": np.concatenate([\n",
        "      np.repeat(_a2r[seqid], _chrom_lengths[\"kl\"][seqid])\n",
        "        for seqid in _chrom_lengths[\"kl\"]]),\n",
        "    \"position\": np.concatenate([\n",
        "      np.arange(_chrom_lengths[\"kl\"][seqid])\n",
        "        for seqid in _chrom_lengths[\"kl\"]])\n",
        "  }), _np_get_ds_fa(\"kl\")),\n",
        "  axis=1\n",
        ").reset_index().drop(\"index\", axis=1).drop(\"level_0\", axis=1).astype(\"int\")\n",
        "# combine genomic data with sequencing data\n",
        "gen_seq_kl = pd.concat((\n",
        "  raw_genomic_kl,\n",
        "  _np_ns_bedgraph_to_df_kl(None),\n",
        "  _np_ms_bedgraph_to_df_kl(1),\n",
        "  _np_ms_bedgraph_to_df_kl(2)\n",
        "), axis=1).reset_index().drop(\"index\", axis=1).astype(\"int\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Gz4FkWhiuY6W",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# save\n",
        "gen_seq_kl_path = _dirs[\"collated\"] + \"gen_seq_kl.csv\"\n",
        "if os.path.isfile(gen_seq_kl_path):\n",
        "  if input(f\"overwrite {gen_seq_kl_path}? ... yes, no\\n\") == \"yes\":\n",
        "    print(\"overwriting old gen_seq_kl with new gen_seq_kl\")\n",
        "    gen_seq_kl.to_csv(gen_seq_kl_path)\n",
        "else:\n",
        "  print(\"saving gen_seq_kl\")\n",
        "  gen_seq_kl.to_csv(gen_seq_kl_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TtAM-2bCq_lS",
        "colab_type": "text"
      },
      "source": [
        "## prepare data for FIDDLE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-LSZGQqvf5jE",
        "colab_type": "text"
      },
      "source": [
        "### K. _lactis_\n",
        "\n",
        "1. collate to pandas dataframe, save\n",
        "2. collate to dictionary format,  amenable to model testing\n",
        "\n",
        "```\n",
        "(species)_test_df = pd.DataFrame({\n",
        "    ID, strand, seqid, position,\n",
        "    DS_A, DS_C, DS_G, DS_T,\n",
        "    NS_sense, NS_antisense,\n",
        "    MS_1, MS_2\n",
        "})\n",
        "\n",
        "(species)_test_dict = {\n",
        "    (species)_test_ds: ...,\n",
        "    (species)_test_ns: ...\n",
        "}\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_-wlvP2WgNgB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# load data\n",
        "data_kl = {\n",
        "  \"meta_kl\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"meta_kl.csv\")),\n",
        "  \"gen_seq_kl\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"gen_seq_kl.csv\"))\n",
        "}\n",
        "# remove unnecessary columns\n",
        "nuisance = \"Unnamed: 0\"\n",
        "for df in list(data_kl.keys()):\n",
        "  if nuisance in data_kl[df].columns.values:\n",
        "    data_kl[df].drop(nuisance, axis=1, inplace=True)\n",
        "# collate to dataframe\n",
        "test_dfs = {\n",
        "  \"kl\": pd.concat([\n",
        "    _np_5_to_3_kl_sp(gene[1], data_kl[\"gen_seq_kl\"])\n",
        "      for gene in data_kl[\"meta_kl\"].iterrows()\n",
        "  ])\n",
        "}\n",
        "# drop extra seqid columns\n",
        "for species in [\"kl\"]:\n",
        "  if test_dfs[species][\"seqid\"].values.shape[1] == 2:\n",
        "    test_dfs[species][\"SEQID\"] = test_dfs[species][\"seqid\"].values[:, 0]\n",
        "    test_dfs[species].drop(\"seqid\", axis=1, inplace=True)\n",
        "    test_dfs[species][\"seqid\"] = test_dfs[species][\"SEQID\"].values\n",
        "    test_dfs[species] = test_dfs[species][[\n",
        "      \"ID\", \"strand\", \"seqid\", \"position\",\n",
        "      \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "      \"NS_sense\", \"NS_antisense\",\n",
        "      \"MS_1\"]]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PKfL7fhCzSCG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# save test dataframes\n",
        "test_dfs[\"kl\"].to_csv(_dirs[\"collated\"] + \"kl_test_df.csv\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "smNPFNAqJgUJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# extract test arrays to dictionary\n",
        "test_dicts = {\n",
        "  species: {\n",
        "    species + \"_test_ds\": np.expand_dims(\n",
        "      a=_np_get_kl_sp_data(\n",
        "        columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"],\n",
        "        data=test_dfs[species]),\n",
        "      axis=3),\n",
        "    species + \"_test_ns\": np.expand_dims(\n",
        "      a=_np_get_kl_sp_data(\n",
        "        columns=[\"NS_antisense\", \"NS_sense\"],\n",
        "        data=test_dfs[species]),\n",
        "      axis=3),\n",
        "    species + \"_test_ms\": np.expand_dims(\n",
        "      a=_np_get_kl_sp_data(\n",
        "        columns=[\"MS_1\"],\n",
        "        data=test_dfs[species]),\n",
        "      axis=3)\n",
        "  } for species in [\"kl\"]\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8HWayqG7VWj_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# save test dictionaries\n",
        "np.save(_dirs[\"model_inputs\"] + \"kl_test_dict\", test_dicts[\"kl\"])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MMOoeTdH33oQ",
        "colab_type": "text"
      },
      "source": [
        "### S. _cerevisiae_\n",
        "\n",
        "1. collate to pandas dataframe\n",
        "2. define 16 randomly sampled AUG-phased regions for training, one TSS-phased region for validation and testing\n",
        "  - **note: training only done on wildtype**\n",
        "2. collate to dictionary format,  amenable to model testing\n",
        "\n",
        "```\n",
        "sc_train_df = pd.DataFrame({\n",
        "    ID, strand, seqid, position,\n",
        "    DS_A, DS_C, DS_G, DS_T,\n",
        "    NS_wt_pos, NS_wt_neg,\n",
        "    MS_wt_1\n",
        "})\n",
        "sc_train_dict = {\n",
        "    seqid (#): {\n",
        "      sc_train_ds: ...,\n",
        "      sc_train_ns: ...,\n",
        "      sc_train_ms: ...\n",
        "    } for seqid in S. cerevisiae\n",
        "}\n",
        "sc_valid_test_df = pd.DataFrame({\n",
        "    ID, strand, seqid, position,\n",
        "    DS_A, DS_C, DS_G, DS_T,\n",
        "    NS_wt_pos, NS_wt_neg,\n",
        "    NS_ko_chd1_pos, NS_ko_chd1_neg,\n",
        "    NS_ko_ino80_pos, NS_ko_ino80_neg,\n",
        "    NS_ko_isw1_pos, NS_ko_isw1_neg,\n",
        "    NS_ko_isw2_pos, NS_ko_isw2_neg,\n",
        "    NS_ko_rsc_pos, NS_ko_rsc_neg,\n",
        "    NS_ko_spt6-1004_pos, NS_ko_spt6-1004_neg,\n",
        "    NS_ko_spt6-YW_pos, NS_ko_spt6-YW_neg,\n",
        "    MS_wt_1, MS_wt_2, MS_wt_3, MS_wt_4,\n",
        "    MS_ko_spt6-1004, MS_ko_spt6-YW\n",
        "})\n",
        "sc_valid_test_dict = {\n",
        "    seqid (#): {\n",
        "      sc_valid_test_ds: ...,\n",
        "      sc_valid_test_ns_wt: ...,\n",
        "      sc_valid_test_ns_ko_chd1: ...,\n",
        "      sc_valid_test_ns_ko_ino80: ...,\n",
        "      sc_valid_test_ns_ko_isw1: ...,\n",
        "      sc_valid_test_ns_ko_isw2: ...,\n",
        "      sc_valid_test_ns_ko_rsc: ...,\n",
        "      sc_valid_test_ns_ko_spt6-1004: ...,\n",
        "      sc_valid_test_ns_ko_spt6-YW: ...,\n",
        "      sc_valid_test_ms: ...\n",
        "    } for seqid in S. cerevisiae\n",
        "}\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Agxyr5Cx_ZGg",
        "colab_type": "text"
      },
      "source": [
        "#### load meta and genomic / sequencing dataframes"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ceb2dDDbm3wj",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# load\n",
        "data_sc = {\n",
        "  \"meta_sc\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"meta_sc.csv\")),\n",
        "  \"gen_seq_sc\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"gen_seq_sc.csv\"))\n",
        "}\n",
        "# remove unnecessary columns\n",
        "nuisance = \"Unnamed: 0\"\n",
        "for df in list(data_sc.keys()):\n",
        "  if nuisance in data_sc[df].columns.values:\n",
        "    data_sc[df].drop(nuisance, axis=1, inplace=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JZn8NI0_nwuR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "data_sc[\"meta_sc\"].head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7Kxe7QgHFVbS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# sanity check\n",
        "data_sc[\"gen_seq_sc\"].head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uGwEmm16_a_H",
        "colab_type": "text"
      },
      "source": [
        "#### train dataframe via CoLab"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-2FtTx7wfME3",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# collate to dataframe\n",
        "sc_train_df = pd.concat([\n",
        "  pd.concat([\n",
        "    _np_5_to_3_sc(gene[1], \"AUG\", data_sc[\"gen_seq_sc\"])\n",
        "      for gene in data_sc[\"meta_sc\"].iterrows()\n",
        "    ], ignore_index=True)\n",
        "      for _ in range(_heuristics[\"training_regions\"])\n",
        "  ], ignore_index=True)\n",
        "# drop extra seqid columns\n",
        "sc_train_df[\"SEQID\"] = sc_train_df[\"seqid\"].values[:, 0]\n",
        "sc_train_df.drop(\"seqid\", axis=1, inplace=True)\n",
        "sc_train_df[\"seqid\"] = sc_train_df[\"SEQID\"].values\n",
        "sc_train_df = sc_train_df[[\n",
        "  \"ID\", \"strand\", \"seqid\", \"position\",\n",
        "  \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "  \"NS_wt_sense\", \"NS_wt_antisense\",\n",
        "  \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\"\n",
        "]]\n",
        "# save\n",
        "sc_train_df.to_csv(_dirs[\"collated\"] + \"sc_train_df.csv\")\n",
        "# extract train set arrays to dictionary\n",
        "sc_train_dict = {\n",
        "  seqid: {\n",
        "    \"sc_train_ds\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"],\n",
        "        seqid=seqid, data=sc_train_df),\n",
        "      axis=3),\n",
        "    \"sc_train_ns\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_wt_antisense\", \"NS_wt_sense\"],\n",
        "        seqid=seqid, data=sc_train_df),\n",
        "      axis=3),\n",
        "    \"sc_train_ms\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"MS_wt_1\"],\n",
        "        seqid=seqid, data=sc_train_df),\n",
        "      axis=3),\n",
        "  } for seqid in np.unique(sc_train_df[\"seqid\"].values)\n",
        "}\n",
        "# save\n",
        "np.save(_dirs[\"model_inputs\"] + \"sc_train_dict\", sc_train_dict)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X7aQLpfcE-sE",
        "colab_type": "text"
      },
      "source": [
        "#### train dataframe via script"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "LqJSjRM3FMDm",
        "colab": {}
      },
      "source": [
        "# define sbatch submissions wrt partitions\n",
        "submissions = [\n",
        "  f'module load gcc/6.2.0 python/3.6.0 cuda/9.0; source \\\n",
        "  ~/colab_venv/bin/activate; sbatch -p short -n 1 --mem 32G -t 0-02:00 \\\n",
        "  --mail-type=\"NONE\" --mail-user=\"' + _slurm[\"email\"] + '\" \\\n",
        "  -o /n/groups/churchman/dmm57/COLAB/logs/sc_train_split_{run}_out.txt \\\n",
        "  -e /n/groups/churchman/dmm57/COLAB/logs/sc_train_split_{run}_error.txt \\\n",
        "  -J {run} --wrap=\"python3.6 \\\n",
        "  /n/groups/churchman/dmm57/COLAB/MODELS/SCRIPTS/sc_train_split.py \\\n",
        "  {run}\"' for run in range(1, 17)\n",
        "]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "zhG8ztDTFMDo",
        "colab": {}
      },
      "source": [
        "# partition data\n",
        "for submission in submissions: subprocess.run(submission, shell=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xIEhmaYK_gDM",
        "colab_type": "text"
      },
      "source": [
        "#### valid / test dataframe"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8-vbLN8TxWrz",
        "colab_type": "code",
        "cellView": "both",
        "colab": {}
      },
      "source": [
        "%%time\n",
        "# collate to dataframe\n",
        "sc_valid_test_df = pd.concat([\n",
        "  _np_5_to_3_sc(gene[1], \"TSS\", data_sc[\"gen_seq_sc\"])\n",
        "      for gene in data_sc[\"meta_sc\"].iterrows()\n",
        "  ], ignore_index=True)\n",
        "# drop extra seqid columns\n",
        "sc_valid_test_df[\"SEQID\"] = sc_valid_test_df[\"seqid\"].values[:, 0]\n",
        "sc_valid_test_df.drop(\"seqid\", axis=1, inplace=True)\n",
        "sc_valid_test_df[\"seqid\"] = sc_valid_test_df[\"SEQID\"].values\n",
        "sc_valid_test_df = sc_valid_test_df[[\n",
        "  \"ID\", \"strand\", \"seqid\", \"position\",\n",
        "  \"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\",\n",
        "  \"NS_wt_sense\", \"NS_wt_antisense\",\n",
        "  \"NS_ko_chd1_sense\", \"NS_ko_chd1_antisense\",\n",
        "  \"NS_ko_ino80_sense\", \"NS_ko_ino80_antisense\",\n",
        "  \"NS_ko_isw1_sense\", \"NS_ko_isw1_antisense\",\n",
        "  \"NS_ko_isw2_sense\", \"NS_ko_isw2_antisense\",\n",
        "  \"NS_ko_rsc_sense\", \"NS_ko_rsc_antisense\",\n",
        "  \"NS_ko_spt6-1004_sense\", \"NS_ko_spt6-1004_antisense\",\n",
        "  \"NS_ko_spt6-YW_sense\", \"NS_ko_spt6-YW_antisense\",\n",
        "  \"MS_wt_1\", \"MS_wt_2\", \"MS_wt_3\", \"MS_wt_4\",\n",
        "  \"MS_ko_spt6-1004\", \"MS_ko_spt6-YW\"\n",
        "]]\n",
        "# save\n",
        "sc_valid_test_df.to_csv(_dirs[\"collated\"] + \"sc_valid_test_df.csv\")\n",
        "# extract valid / test set arrays to dictionary\n",
        "sc_valid_test_dict = {\n",
        "  seqid: {\n",
        "    \"sc_valid_test_ds\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_wt\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_wt_antisense\", \"NS_wt_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_chd1\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_chd1_antisense\", \"NS_ko_chd1_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_ino80\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_ino80_antisense\", \"NS_ko_ino80_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_isw1\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_isw1_antisense\", \"NS_ko_isw1_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_isw2\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_isw2_antisense\", \"NS_ko_isw2_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_rsc\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_rsc_antisense\", \"NS_ko_rsc_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_spt6-1004\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_spt6-1004_antisense\", \"NS_ko_spt6-1004_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ns_ko_spt6-YW\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"NS_ko_spt6-YW_antisense\", \"NS_ko_spt6-YW_sense\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "    \"sc_valid_test_ms\": np.expand_dims(\n",
        "      a=_np_get_sc_data(\n",
        "        columns=[\"MS_wt_1\"],\n",
        "        seqid=seqid, data=sc_valid_test_df),\n",
        "      axis=3),\n",
        "  } for seqid in np.unique(sc_valid_test_df[\"seqid\"].values)\n",
        "}\n",
        "# save\n",
        "np.save(_dirs[\"model_inputs\"] + \"sc_valid_test_dict\", sc_valid_test_dict)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vBYh84bRhQRs",
        "colab_type": "text"
      },
      "source": [
        "# MODEL"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "91pmO_OOpMou",
        "colab_type": "text"
      },
      "source": [
        "## high throughput model generation\n",
        "\n",
        "- mode (mode, parameters) defined as:\n",
        "  > < DS _ { native / shuffled / none } _ NS _ { native / shuffled / none } _ validSeqid _ { # } _ testSeqid { # } _ run _ { # } >\n",
        "\n",
        "- 8 DNA sequence / NET-seq modes, 15 validation chromosomes, 16 test chromosomes, 16 computational replicates\n",
        "  > < input modalities >\n",
        "    - DS_native_NS_native\n",
        "    - DS_native_NS_none\n",
        "    - DS_native_NS_shuffled\n",
        "    - DS_none_NS_native\n",
        "    - DS_none_NS_shuffled\n",
        "    - DS_shuffled_NS_native\n",
        "    - DS_shuffled_NS_none\n",
        "    - DS_shuffled_NS_shuffled\n",
        "    - note: DS_none_NS_none is **_not_** included for obvious reasons\n",
        "\n",
        "  > < amounts to: 30,720 models >"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ecX2FGkvFZy2",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title define mode parameters\n",
        "mode_params = [\n",
        "  mode_param for mode_param in [\n",
        "    \"_\".join([\n",
        "      \"DS\", ds_form, \"NS\", ns_form, \"validSeqid\", str(vs), \"testSeqid\", str(ts),\n",
        "      \"run\", str(run)\n",
        "    ]) for ds_form in [\"native\", \"shuffled\", \"none\"]\n",
        "      for ns_form in [\"native\", \"shuffled\", \"none\"]\n",
        "        for vs in np.arange(1, 17)\n",
        "          for ts in np.setdiff1d(np.arange(1, 17), vs)\n",
        "            for run in range(1, 17)\n",
        "  ] if \"DS_none_NS_none\" not in mode_param\n",
        "]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iaJQiCzxE_0w",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title define models from scratch\n",
        "submissions = [\n",
        "  f'module load gcc/6.2.0 python/3.6.0 cuda/9.0; source ' + \n",
        "  _slurm[\"virtual_env\"] + 'bin/activate; sbatch -p gpu -n 1 --mem 8G -t \\\n",
        "  0-00:16 --gres=gpu:1 --mail-type=\\\"NONE\\\" --mail-user=' + _slurm[\"email\"] +\n",
        "  ' -o ' + _dirs[\"logs\"] + mode_param + '_output_out.txt -e ' + \n",
        "  _dirs[\"logs\"] + mode_param + '_output_error.txt -J ' + mode_param + \n",
        "  ' --wrap=\"python3.6 ' + _dirs[\"scripts\"] + 'fiddle.py ' + mode_param + '\"'\n",
        "    for mode_param in unfinished[:6400]\n",
        "]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z4IAQHVQIzse",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# run model on all partitioned data\n",
        "for submission in submissions: subprocess.run(submission, shell=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dk4oJHoEEvTA",
        "colab_type": "text"
      },
      "source": [
        "## reload models\n",
        "###### _for runs that did not complete_"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "97ePZN3CYstn",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title define named models\n",
        "named_models = {\n",
        "  x: [h5 for h5 in os.listdir(_dirs[\"trained_models\"]) if y in h5]\n",
        "    for x, y in {\n",
        "      \"integrated\": \"DS_native_NS_native\",\n",
        "      \"DNA\": \"DS_native_NS_none\",\n",
        "      \"NET\": \"DS_none_NS_native\",\n",
        "      \"ctrl_ctrl\": \"DS_shuffled_NS_shuffled\"\n",
        "    }.items()\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XBIrtzS-R8OE",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title define reloaded model submissions\n",
        "named_model = \"integrated\" #@param [\"integrated\", \"DNA\", \"NET\", \"ctrl_ctrl\"]\n",
        "submissions = [\n",
        "  f'module load gcc/6.2.0 python/3.6.0 cuda/9.0; source ' + \n",
        "  _slurm[\"virtual_env\"] + 'bin/activate; sbatch -p gpu -n 1 --mem 4G -t \\\n",
        "  0-00:08 --gres=gpu:1 --mail-type=\\\"NONE\\\" --mail-user=' + _slurm[\"email\"] +\n",
        "  ' -o ' + _dirs[\"logs\"] + model + '_reload_output_out.txt -e ' + \n",
        "  _dirs[\"logs\"] + model + '_reload_output_error.txt -J ' + model + \n",
        "  ' --wrap=\"python3.6 ' + _dirs[\"scripts\"] + 'reload_fiddle.py ' + model + '\"'\n",
        "    for model in named_models[named_model]\n",
        "]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zACo-2wMFADl",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# run model on all partitioned data\n",
        "for submission in submissions: subprocess.run(submission, shell=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NPcWRSFcxJzW",
        "colab_type": "text"
      },
      "source": [
        "# ANALYZE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jUdA9dZ8eL4e",
        "colab_type": "text"
      },
      "source": [
        "## condense outputs"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wNZ0s3LD1nP1",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@markdown input modality\n",
        "mode_param = \"DS_native_NS_native\" #@param [\"DS_native_NS_native\", \"DS_native_NS_none\", \"DS_native_NS_shuffled\", \"DS_none_NS_native\", \"DS_none_NS_shuffled\", \"DS_shuffled_NS_native\", \"DS_shuffled_NS_none\", \"DS_shuffled_NS_shuffled\"]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xWz2ts5cBVU_",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title check input modality file existence\n",
        "print(\n",
        "  mode_param,\n",
        "  \"output dataframes:\\n\",\n",
        "  [x for x in \n",
        "    os.listdir(_dirs[\"collated\"])\n",
        "      if (mode_param in x) and (\"output_df.csv\" in x)\n",
        "  ]\n",
        ")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nWXJDxNzK16O",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title combine high throughput predictions\n",
        "# initialize output dataframes\n",
        "test_dfs = {\n",
        "  \"kl\": pd.read_csv(\n",
        "    _dirs[\"collated\"] + \"kl_test_df.csv\"\n",
        "  )[[\"ID\", \"strand\", \"seqid\", \"position\"]],\n",
        "  \"sc\": pd.read_csv(\n",
        "    _dirs[\"collated\"] + \"sc_valid_test_df.csv\"\n",
        "  )[[\"ID\", \"strand\", \"seqid\", \"position\"]]\n",
        "}\n",
        "  # define paths\n",
        "test_df_paths = [\n",
        "  _dirs[\"collated\"] + \"_\".join([species, mode_param, \"output_df.csv\"])\n",
        "    for species in [\"kl\", \"sc\"]\n",
        "]\n",
        "  # save\n",
        "for df in test_df_paths:\n",
        "  if os.path.isfile(df):\n",
        "    if input(f\"overwrite {df}? ... yes, no\\n\") == \"yes\":\n",
        "      print(f\"overwriting old {df}\")\n",
        "      test_dfs[df.split(\"/\")[-1].split(\"_\")[0]].to_csv(df)\n",
        "  else: test_dfs[df.split(\"/\")[-1].split(\"_\")[0]].to_csv(df)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Qr_OUGt-Ltuz",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title define sbatch high throughput submissions\n",
        "submissions = [\n",
        "  'module load gcc/6.2.0 python/3.6.0 cuda/9.0; source ' + \\\n",
        "  _slurm[\"virtual_env\"] + 'bin/activate; sbatch -p short -n 1 --mem 16G' + \\\n",
        "  ' -t 0-00:30 --mail-type=\\\"NONE\\\" --mail-user=' + _slurm[\"email\"] + ' -o ' + \\\n",
        "  _dirs[\"logs\"] + \"_\".join([mode_param, run_m]) + '.error.txt -e ' + \\\n",
        "  _dirs[\"logs\"] + \"_\".join([mode_param, run_m]) + '.error.txt -J ' + mode_param + \\\n",
        "  ' --dependency=singleton --wrap=\"python3.6 ' + _dirs[\"scripts\"] + \\\n",
        "  'collate_outputs.py ' + _dirs[\"collated\"] + ' ' + _dirs[\"model_outputs\"] + \\\n",
        "  ' ' + mode_param + ' ' + run_m + '\"'\n",
        "    for run_m in [\n",
        "      \"_\".join([\"run\", str(r), \"pred\"])\n",
        "        for r in range(1, 17)\n",
        "    ]\n",
        "]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y_zaQlf-238B",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# submit to command line\n",
        "for submission in submissions: subprocess.run(submission, shell=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "911AoZz24Ai1",
        "colab_type": "text"
      },
      "source": [
        "## prepare data for figures and related"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lU68UwgyAH8e",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title load measured and predicted data\n",
        "#@markdown < all models > or < named models (further considered) >\n",
        "data_context = \"named\" #@param [\"all\", \"named\"]\n",
        "  # mapping\n",
        "_data_context = {\n",
        "  \"all\": [\n",
        "    \"DS_native_NS_native\",\n",
        "    \"DS_native_NS_none\",\n",
        "    \"DS_native_NS_shuffled\",\n",
        "    \"DS_none_NS_native\",\n",
        "    \"DS_none_NS_shuffled\",\n",
        "    \"DS_shuffled_NS_native\",\n",
        "    \"DS_shuffled_NS_none\",\n",
        "    \"DS_shuffled_NS_shuffled\"\n",
        "  ],\n",
        "  \"named\": [\n",
        "    \"DS_native_NS_native\",\n",
        "    \"DS_native_NS_none\",\n",
        "    \"DS_none_NS_native\",\n",
        "    \"DS_shuffled_NS_shuffled\"\n",
        "  ]\n",
        "}\n",
        "  # load measured data\n",
        "meas_dfs = {\n",
        "  \"kl\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"kl_test_df.csv\")),\n",
        "  \"sc\": pd.read_csv(os.path.join(_dirs[\"collated\"], \"sc_valid_test_df.csv\"))\n",
        "}\n",
        "  # load collated predictions\n",
        "pred_dfs = {\n",
        "  \"_\".join([species, mode_param, \"output_df\"]): pd.read_csv(\n",
        "    _dirs[\"collated\"] + \"_\".join([species, mode_param, \"output_df.csv\"])\n",
        "  ) for species in [\"kl\", \"sc\"]\n",
        "    for mode_param in _data_context[data_context]\n",
        "}\n",
        "  # tidy up dataframes\n",
        "for df in meas_dfs.keys():\n",
        "  meas_dfs[df] = meas_dfs[df][[\n",
        "    col for col in meas_dfs[df].columns.values\n",
        "      if \"Unnamed\" not in col\n",
        "  ]]\n",
        "for df in pred_dfs.keys():\n",
        "  pred_dfs[df] = pred_dfs[df][[\n",
        "    col for col in pred_dfs[df].columns.values\n",
        "      if \"Unnamed\" not in col\n",
        "  ]]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r8VIvKR9WP-9",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title load meta / gen_seq dataframes\n",
        "  # load meta data\n",
        "meta_dfs = {\n",
        "  species: pd.read_csv(\n",
        "    os.path.join(_dirs[\"collated\"], \"meta_\" + species + \".csv\")\n",
        "  ) for species in [\"kl\", \"sc\"]\n",
        "}\n",
        "  # load genomic sequencing data\n",
        "gen_seq_dfs = {\n",
        "  species: pd.read_csv(\n",
        "    os.path.join(_dirs[\"collated\"], \"gen_seq_\" + species + \".csv\")\n",
        "  ) for species in [\"kl\", \"sc\"]\n",
        "}\n",
        "  # tidy up dataframes\n",
        "for df in meta_dfs.keys():\n",
        "  meta_dfs[df] = meta_dfs[df][[\n",
        "    col for col in meta_dfs[df].columns.values\n",
        "      if \"Unnamed\" not in col\n",
        "  ]]\n",
        "for df in gen_seq_dfs.keys():\n",
        "  gen_seq_dfs[df] = gen_seq_dfs[df][[\n",
        "    col for col in gen_seq_dfs[df].columns.values\n",
        "      if \"Unnamed\" not in col\n",
        "  ]]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5svljAaa4P-X",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title local (dict) and meta (df) MNase-seq\n",
        "  # local - ( # regions, 1024 )\n",
        "agg_ms_data = {\n",
        "  \"kl\": {\n",
        "    \"meas\": {\n",
        "      \"MS_1\": _np_get_seq_data([\"MS_1\"], meas_dfs[\"kl\"]).squeeze()\n",
        "    },\n",
        "    \"pred\": {\n",
        "      \"_\".join(x.split(\"_\")[1:5]): _np_avg_var_df(pred_dfs[x])\n",
        "        for x in pred_dfs.keys() if \"kl_\" in x\n",
        "    }\n",
        "  },\n",
        "  \"sc\": {\n",
        "    \"meas\": {\n",
        "      \"MS_wt_1\": _np_get_seq_data([\"MS_wt_1\"], meas_dfs[\"sc\"]).squeeze(),\n",
        "      \"MS_wt_2\": np.sum(np.concatenate([\n",
        "        _np_get_seq_data([\"MS_wt_2\"], meas_dfs[\"sc\"]),\n",
        "        _np_get_seq_data([\"MS_wt_3\"], meas_dfs[\"sc\"])\n",
        "      ], 1), 1),\n",
        "      \"MS_wt_4\": _np_get_seq_data([\"MS_wt_4\"], meas_dfs[\"sc\"]).squeeze()\n",
        "    },\n",
        "    \"pred\": {\n",
        "      state: {\n",
        "        \"_\".join(x.split(\"_\")[1:5]): _np_avg_var_df(pred_dfs[x], [state])\n",
        "          for x in pred_dfs.keys() if \"sc_\" in x\n",
        "      } for state in _sc_N_states\n",
        "    }\n",
        "  }\n",
        "}\n",
        "  # global / meta - ( 1024 )\n",
        "agg_ms_meta = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"kl_MS_1\": _np_ms_meta_repr(agg_ms_data[\"kl\"][\"meas\"][\"MS_1\"], \"meas\"),\n",
        "    \"sc_MS_wt_1\": _np_ms_meta_repr(agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_1\"], \"meas\"),\n",
        "    \"sc_MS_wt_2\": _np_ms_meta_repr(agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_2\"], \"meas\"),\n",
        "    \"sc_MS_wt_4\": _np_ms_meta_repr(agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_4\"], \"meas\")\n",
        "  }),\n",
        "  pd.DataFrame({\n",
        "    \"_\".join([\"kl\", x]): _np_ms_meta_repr(\n",
        "      agg_ms_data[\"kl\"][\"pred\"][x][\"avg\"], \"meas\"\n",
        "    ) for x in _data_context[data_context]\n",
        "  }),\n",
        "  pd.DataFrame({\n",
        "    \"_\".join([\"sc\", state, x]): _np_ms_meta_repr(\n",
        "      agg_ms_data[\"sc\"][\"pred\"][state][x][\"avg\"], \"meas\"\n",
        "    ) for x in _data_context[data_context]\n",
        "      for state in _sc_N_states\n",
        "  })\n",
        "], 1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XdWDh1KJpPeO",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title nucleosome locations\n",
        "  # local maxima\n",
        "nuc_data = {\n",
        "  meta: _np_ms_meta_extrema(agg_ms_meta[meta].values)\n",
        "    for meta in [\n",
        "      m for m in agg_ms_meta.columns.values\n",
        "        if any([\n",
        "          (\"DS_shuffled_NS_shuffled\" in m) or \n",
        "          (\"DS_none_NS_native\" in m) or \n",
        "          (\"DS_native_NS_none\" in m) or \n",
        "          (\"DS_native_NS_native\" in m) or \n",
        "          (\"MS\" in m)\n",
        "        ])\n",
        "    ]\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XSkjmWu13niS",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title local (dict) and meta (dict) NET-seq\n",
        "  # local - ( # regions, 2, 1024 )\n",
        "agg_ns_data = {\n",
        "  \"kl\": _np_get_seq_data(\n",
        "    columns=[\"NS_antisense\", \"NS_sense\"],\n",
        "    data=meas_dfs[\"kl\"]).squeeze(),\n",
        "  \"sc_wt\": _np_get_seq_data(\n",
        "    columns=[\"NS_wt_antisense\", \"NS_wt_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_chd1\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_chd1_antisense\", \"NS_ko_chd1_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze()#,\n",
        "  \"sc_ko_ino80\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_ino80_antisense\", \"NS_ko_ino80_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_isw1\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_isw1_antisense\", \"NS_ko_isw1_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_isw2\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_isw2_antisense\", \"NS_ko_isw2_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_rsc\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_rsc_antisense\", \"NS_ko_rsc_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_spt6-1004\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_spt6-1004_antisense\", \"NS_ko_spt6-1004_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze(),\n",
        "  \"sc_ko_spt6-YW\": _np_get_seq_data(\n",
        "    columns=[\"NS_ko_spt6-YW_antisense\", \"NS_ko_spt6-YW_sense\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze()\n",
        "}\n",
        "  # meta - ( 1024 )\n",
        "# total NET-seq signal for each species / cell state\n",
        "ns_sums = {\n",
        "  species_state: np.sum(agg_ns_data[species_state], (1, 2))\n",
        "    for species_state in agg_ns_data.keys()\n",
        "}\n",
        "# NET-seq signal indices for each species / cell state g.t. 95th %\n",
        "ns_viz_idx = {\n",
        "  species_state: np.where(\n",
        "    ns_sums[species_state] < np.percentile(ns_sums[species_state], 95)\n",
        "  )[0] for species_state in ns_sums.keys()\n",
        "}\n",
        "# NET-seq meta profiles\n",
        "ns_meta = {\n",
        "  species_state: {\n",
        "    \"sense\": np.array([\n",
        "      _np_ns_transform(x) for x in agg_ns_data[species_state][\n",
        "        ns_viz_idx[species_state], 0, :\n",
        "    ]]),\n",
        "    \"antisense\": np.array([\n",
        "      _np_ns_transform(x) for x in agg_ns_data[species_state][\n",
        "        ns_viz_idx[species_state], 1, :\n",
        "    ]])\n",
        "  } for species_state in ns_viz_idx.keys()\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GZ-scnr3PyNq",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title local (dict) and meta (df) DNA sequence\n",
        "  # local - ( # regions, 4, 1024 )\n",
        "agg_ds_data = {\n",
        "  \"kl\": _np_get_seq_data(\n",
        "    columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"],\n",
        "    data=meas_dfs[\"kl\"]).squeeze(),\n",
        "  \"sc\": _np_get_seq_data(\n",
        "    columns=[\"DS_A\", \"DS_C\", \"DS_G\", \"DS_T\"],\n",
        "    data=meas_dfs[\"sc\"]).squeeze()\n",
        "}\n",
        "  # global / meta - ( 1024 )\n",
        "agg_ds_meta = pd.DataFrame({\n",
        "  species: _np_ds_meta_repr(meas_dfs[species])\n",
        "    for species in meas_dfs.keys()\n",
        "})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V1hUloAywgvN",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title sequence identity comparison\n",
        "  # K. lactis vs K. lactis\n",
        "kl_kl_seq_identity = np.array([\n",
        "  100 * (np.sum([np.dot(i[nt, :], j[nt, :]) for nt in range(4)]) / 1024)\n",
        "      for i in agg_ds_data[\"kl\"] for j in agg_ds_data[\"kl\"]\n",
        "]).reshape((agg_ds_data[\"kl\"].shape[0], agg_ds_data[\"kl\"].shape[0]))\n",
        "  # S. cerevisiae vs S. cerevisiae\n",
        "sc_sc_seq_identity = np.array([\n",
        "  100 * (np.sum([np.dot(i[nt, :], j[nt, :]) for nt in range(4)]) / 1024)\n",
        "      for i in agg_ds_data[\"sc\"] for j in agg_ds_data[\"sc\"]\n",
        "]).reshape((agg_ds_data[\"sc\"].shape[0], agg_ds_data[\"sc\"].shape[0]))\n",
        "  # K. lactis vs S. cerevisiae\n",
        "kl_sc_seq_identity = np.array([\n",
        "  100 * (np.sum([np.dot(i[nt, :], j[nt, :]) for nt in range(4)]) / 1024)\n",
        "      for i in agg_ds_data[\"kl\"] for j in agg_ds_data[\"sc\"]\n",
        "]).reshape((agg_ds_data[\"kl\"].shape[0], agg_ds_data[\"sc\"].shape[0]))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AdqveKCx3XWV",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title loss function, stats (df)\n",
        "\n",
        "# S. cerevisiae Wildtype\n",
        "  # calculate losses to wide form\n",
        "sc_wt_loss_df_wide = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"rep1_rep2\": np.array([_np_jsd(\n",
        "      _np_ms_transform(p),\n",
        "      _np_ms_transform(q)\n",
        "    ) for p, q in zip(\n",
        "      agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_1\"],\n",
        "      agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_2\"]\n",
        "    )])\n",
        "  }),\n",
        "  pd.DataFrame({\n",
        "    \"rep1_CM\": np.array([_np_jsd(\n",
        "      _np_ms_transform(p),\n",
        "      _np_ms_transform(q)\n",
        "    ) for p, q in zip(\n",
        "      agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_1\"],\n",
        "      agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_4\"]\n",
        "    )])\n",
        "  }),\n",
        "  pd.DataFrame({\n",
        "    x: np.array([_np_jsd(\n",
        "      _np_ms_transform(p),\n",
        "      _np_ms_transform(q)\n",
        "    ) for p, q in zip(\n",
        "      agg_ms_data[\"sc\"][\"meas\"][\"MS_wt_1\"],\n",
        "      agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"][x][\"avg\"]\n",
        "    )]) for x in agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"].keys()\n",
        "  })\n",
        "], 1)\n",
        "  # reorder columns by median loss\n",
        "sc_wt_loss_df_wide = sc_wt_loss_df_wide[\n",
        "  list(sc_wt_loss_df_wide.median().sort_values().index.values[::-1])\n",
        "]\n",
        "  # collate losses to long form\n",
        "sc_wt_loss_df_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"p_vs_q\": np.repeat(col, len(sc_wt_loss_df_wide)),\n",
        "    \"jsd\": sc_wt_loss_df_wide[col].values\n",
        "  }) for col in sc_wt_loss_df_wide.columns.values\n",
        "], 0)\n",
        "\n",
        "# K. lactis\n",
        "  # calculate losses to wide form\n",
        "kl_loss_df_wide = pd.DataFrame({\n",
        "  x: np.array([_np_jsd(\n",
        "    _np_ms_transform(p),\n",
        "    _np_ms_transform(q)\n",
        "  ) for p, q in zip(\n",
        "    agg_ms_data[\"kl\"][\"meas\"][\"MS_1\"],\n",
        "    agg_ms_data[\"kl\"][\"pred\"][x][\"avg\"]\n",
        "  )]) for x in agg_ms_data[\"kl\"][\"pred\"].keys()\n",
        "})\n",
        "  # reorder columns by median loss\n",
        "kl_loss_df_wide = kl_loss_df_wide[\n",
        "  list(kl_loss_df_wide.median().sort_values().index.values[::-1])\n",
        "]\n",
        "  # collate losses to long form\n",
        "kl_loss_df_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"p_vs_q\": np.repeat(col, len(kl_loss_df_wide)),\n",
        "    \"jsd\": kl_loss_df_wide[col].values\n",
        "  }) for col in kl_loss_df_wide.columns.values\n",
        "], 0)\n",
        "# statistical significance via KS-test\n",
        "  # S. cerevisiae Wildtype\n",
        "sc_wt_stat_sig = pd.DataFrame({\n",
        "  col_i: np.array([\n",
        "    stats.ks_2samp(\n",
        "      sc_wt_loss_df_wide[col_i].values,\n",
        "      sc_wt_loss_df_wide[col_j].values\n",
        "    )[1] for col_j in sc_wt_loss_df_wide.columns\n",
        "  ]) for col_i in sc_wt_loss_df_wide.columns\n",
        "})\n",
        "  # K. lactis\n",
        "kl_stat_sig = pd.DataFrame({\n",
        "  col_i: np.array([\n",
        "    stats.ks_2samp(\n",
        "      kl_loss_df_wide[col_i].values,\n",
        "      kl_loss_df_wide[col_j].values\n",
        "    )[1] for col_j in kl_loss_df_wide.columns\n",
        "  ]) for col_i in kl_loss_df_wide.columns\n",
        "})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XM95zaCzuO6v",
        "colab_type": "text"
      },
      "source": [
        "## plot aggregate considerations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dOUMAc0hsQsD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# plotting style\n",
        "sns.set_style(\"ticks\")\n",
        "# sns.set(style=\"white\", rc={\"axes.facecolor\": (0, 0, 0, 0)})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e-9uvmZzCHkn",
        "colab_type": "text"
      },
      "source": [
        "### loss results"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wg0_NijLyUT1",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Wildtype loss function distribution\n",
        "plot_sc_wt_loss_distr = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_loss_distr = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# plot\n",
        "if plot_sc_wt_loss_distr:\n",
        "  # loss plot details\n",
        "  _sc_wt_loss = {\n",
        "    \"x\": \"p_vs_q\",\n",
        "    \"y\": \"jsd\",\n",
        "    \"data\": sc_wt_loss_df_long,\n",
        "    \"palette\": {\n",
        "      \"DS_none_NS_shuffled\": \"#FFFFFF\",\n",
        "      \"DS_shuffled_NS_none\": \"#FFFFFF\",\n",
        "      \"DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "      \"DS_shuffled_NS_native\": \"#FFFFFF\",\n",
        "      \"DS_none_NS_native\": \"#284292\",\n",
        "      \"DS_native_NS_shuffled\": \"#FFFFFF\",\n",
        "      \"DS_native_NS_none\": \"#CA2128\",\n",
        "      \"DS_native_NS_native\": \"#BB478D\",\n",
        "      \"rep1_rep2\": \"#3BAC49\"\n",
        "    }\n",
        "  }\n",
        "  _sc_wt_loss_box = {\n",
        "    \"fliersize\": 0,\n",
        "    \"linewidth\": 2.5,\n",
        "    \"width\": 0.55\n",
        "  }\n",
        "  # plot\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  sns.boxplot(**{**_sc_wt_loss, **_sc_wt_loss_box}, ax=A)\n",
        "  A.set_xticks([])\n",
        "  A.set_xlabel(None)\n",
        "  A.set_ylabel(None)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"right\", \"bottom\", \"top\"]]\n",
        "  A.ticklabel_format(axis=\"y\", style=\"sci\", useOffset=True)\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if sc_wt_save_loss_distr:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_loss_distr.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_loss_distr.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DyFxpAnVtAkf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Summary Statistics of Loss Against Rep. 2\n",
        "sc_wt_loss_df_wide.describe()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vi5QC9PU5u7Y",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Wildtype loss function distribution\n",
        "plot_sc_wt_loss_distr = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_loss_distr = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# plot\n",
        "if plot_sc_wt_loss_distr:\n",
        "  # loss plot details\n",
        "  _sc_wt_loss = {\n",
        "    \"x\": \"p_vs_q\",\n",
        "    \"y\": \"jsd\",\n",
        "    \"data\": sc_wt_loss_df_long,\n",
        "    \"palette\": {\n",
        "      \"DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "      \"DS_none_NS_native\": \"#284292\",\n",
        "      \"DS_native_NS_none\": \"#CA2128\",\n",
        "      \"DS_native_NS_native\": \"#BB478D\",\n",
        "      \"rep1_rep2\": \"#3BAC49\",\n",
        "      \"rep1_CM\": \"#e6b800\"\n",
        "    }\n",
        "  }\n",
        "  _sc_wt_loss_box = {\n",
        "    \"fliersize\": 0,\n",
        "    \"linewidth\": 2.5,\n",
        "    \"width\": 0.55\n",
        "  }\n",
        "  # plot\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  sns.boxplot(**{**_sc_wt_loss, **_sc_wt_loss_box}, ax=A)\n",
        "  A.set_xticks([])\n",
        "  A.set_ylim([0, 0.25])\n",
        "  A.set_xlabel(None)\n",
        "  A.set_ylabel(None)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"right\", \"bottom\", \"top\"]]\n",
        "  A.ticklabel_format(axis=\"y\", style=\"sci\", useOffset=True)\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_sc_wt_loss_distr:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_loss_distr.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_loss_distr.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O_2P6HkC5u_U",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Summary Statistics of Loss Against Rep. 2\n",
        "sc_wt_loss_df_wide.describe()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kT6t0ML9WbwN",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Wildtype binary statistical significance\n",
        "plot_sc_wt_stat_sig = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# plot\n",
        "if plot_sc_wt_stat_sig:\n",
        "  F, A = plt.subplots(1, 1, figsize=(9, 7))\n",
        "  # heatmap details\n",
        "  _sc_wt_stat_sig = {\n",
        "    \"lbls\": [\n",
        "      \"\\n\".join([x.split(\"_\")[1], x.split(\"_\")[3]])\n",
        "        for x in sc_wt_loss_df_wide.columns.values[:-1]\n",
        "    ] + [\"rep1_rep2\"]\n",
        "  }\n",
        "  # heatmaps\n",
        "  sns.heatmap(\n",
        "    data=sc_wt_stat_sig > 0.05,\n",
        "    cbar=False,\n",
        "    linewidths=2,\n",
        "    xticklabels=_sc_wt_stat_sig[\"lbls\"],\n",
        "    yticklabels=_sc_wt_stat_sig[\"lbls\"],\n",
        "    ax=A\n",
        "  )\n",
        "  # show\n",
        "  plt.tight_layout()\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oXUQ-RA62oqf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# all vs all loss fcn statistical significance\n",
        "sc_wt_stat_sig"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fu015uP0z9qS",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis loss function distribution\n",
        "plot_kl_loss_distr = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_loss_distr = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# plot\n",
        "if plot_kl_loss_distr:\n",
        "  # loss plot details\n",
        "  _kl_loss = {\n",
        "    \"x\": \"p_vs_q\",\n",
        "    \"y\": \"jsd\",\n",
        "    \"data\": kl_loss_df_long,\n",
        "    \"palette\": {\n",
        "      \"DS_none_NS_shuffled\": \"#FFFFFF\",\n",
        "      \"DS_shuffled_NS_none\": \"#FFFFFF\",\n",
        "      \"DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "      \"DS_shuffled_NS_native\": \"#FFFFFF\",\n",
        "      \"DS_none_NS_native\": \"#284292\",\n",
        "      \"DS_native_NS_shuffled\": \"#FFFFFF\",\n",
        "      \"DS_native_NS_none\": \"#CA2128\",\n",
        "      \"DS_native_NS_native\": \"#BB478D\"\n",
        "    },\n",
        "    \"fliersize\": 0,\n",
        "    \"linewidth\": 2.5,\n",
        "    \"width\": 0.55   \n",
        "  }\n",
        "  # plot\n",
        "  F, A = plt.subplots(1, 1, figsize=(8, 7))\n",
        "  sns.boxplot(**_kl_loss, ax=A)\n",
        "  A.set_xticks([])\n",
        "  A.set_ylim([0, 0.15])\n",
        "  A.set_yticks([0, 0.05, 0.10, 0.15])\n",
        "  A.set_xlabel(None)\n",
        "  A.set_ylabel(None)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"right\", \"bottom\", \"top\"]]\n",
        "  A.ticklabel_format(axis=\"y\", style=\"sci\", useOffset=True)\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_kl_loss_distr:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"kl_loss_distr.pdf\", format=\"pdf\")\n",
        "    print(\"kl_loss_distr.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3xNnweuKzQCv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "kl_loss_df_wide.describe()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "cellView": "form",
        "id": "7cxNyf-O5Ps7",
        "colab": {}
      },
      "source": [
        "#@title K. lactis binary statistical significance\n",
        "plot_kl_stat_sig = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# plot\n",
        "if plot_kl_stat_sig:\n",
        "  F, A = plt.subplots(1, 1, figsize=(9, 7))\n",
        "  # heatmap details\n",
        "  _kl_stat_sig = {\n",
        "    \"lbls\": [\n",
        "      \"\\n\".join([x.split(\"_\")[1], x.split(\"_\")[3]])\n",
        "        for x in kl_loss_df_wide.columns.values[:-1]\n",
        "    ]\n",
        "  }\n",
        "  # heatmaps\n",
        "  sns.heatmap(\n",
        "    data=kl_stat_sig > 0.05,\n",
        "    cbar=False,\n",
        "    linewidths=2,\n",
        "    xticklabels=_kl_stat_sig[\"lbls\"],\n",
        "    yticklabels=_kl_stat_sig[\"lbls\"],\n",
        "    ax=A\n",
        "  )\n",
        "  # show\n",
        "  plt.tight_layout()\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "YsunjVV05PtB",
        "colab": {}
      },
      "source": [
        "# all vs all loss fcn statistical significance\n",
        "kl_stat_sig"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6foHpV5LzY1E",
        "colab_type": "text"
      },
      "source": [
        "### MNase-seq"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WENdY0K8xwUy",
        "colab_type": "text"
      },
      "source": [
        "#### canonical meta"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P6EcST8lf7ik",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis\n",
        "plot_kl_ms_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "show_kl_neg_ctrls = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_ms_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "_kl_ms_meta = {\n",
        "  \"kl_MS_1\": {\n",
        "    \"color\": \"#47A677\",\n",
        "    \"label\": \"rep. 1\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.8,\n",
        "    \"alpha\": 0.65\n",
        "  },\n",
        "  \"kl_DS_shuffled_NS_native\": {\n",
        "    \"color\": \"#B49FFA\",\n",
        "    \"label\": \"shuffled native\",\n",
        "    \"linestyle\": (0, (5, 1.5)),\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"kl_DS_none_NS_native\": {\n",
        "    \"color\": \"#6060FD\",\n",
        "    \"label\": \"absent native\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"kl_DS_native_NS_none\": {\n",
        "    \"color\": \"#F5271C\",\n",
        "    \"label\": \"native absent\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"kl_DS_native_NS_shuffled\": {\n",
        "    \"color\": \"#ec1358\",\n",
        "    \"label\": \"native shuffled\",\n",
        "    \"linestyle\": (0, (5, 1.5)),\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"kl_DS_native_NS_native\": {\n",
        "    \"color\": \"#8D05B7\",\n",
        "    \"label\": \"native native\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 3.0,\n",
        "    \"alpha\": 0.65\n",
        "  }\n",
        "}\n",
        "# plot negative controls\n",
        "if show_kl_neg_ctrls:\n",
        "  for neg_ctrl in [\n",
        "    \"kl_DS_shuffled_NS_shuffled\",\n",
        "    \"kl_DS_shuffled_NS_none\",\n",
        "    \"kl_DS_none_NS_shuffled\"\n",
        "  ]: _kl_ms_meta.update({\n",
        "    neg_ctrl: {\n",
        "      \"color\": \"k\",\n",
        "      \"label\": \" \".join([neg_ctrl.split(\"_\")[2], neg_ctrl.split(\"_\")[4]]),\n",
        "      \"linestyle\": \"-\",\n",
        "      \"linewidth\": 0.8,\n",
        "      \"alpha\": 0.5\n",
        "    }\n",
        "  })\n",
        "# plot\n",
        "if plot_kl_ms_meta:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  # measured and predicted\n",
        "  for i, details in _kl_ms_meta.items():\n",
        "    A.plot(agg_ms_meta[i].values, **details)\n",
        "  # ancillary\n",
        "  A.set_yticks([])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticks(np.linspace(0, 1024, 9))\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\"]]\n",
        "  plt.grid(axis=\"x\")\n",
        "  plt.legend(loc=2)\n",
        "  plt.tight_layout()\n",
        "  if save_kl_ms_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"kl_ms_meta.pdf\", format=\"pdf\")\n",
        "    print(\"kl_ms_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "cellView": "form",
        "id": "eDsVytXTOGaa",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae wt\n",
        "plot_sc_wt_ms_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "show_sc_wt_neg_ctrls = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_ms_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "_sc_wt_ms_meta = {\n",
        "  \"sc_MS_wt_1\": {\n",
        "    \"color\": \"#47A677\",\n",
        "    \"label\": \"rep. 1\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.8,\n",
        "    \"alpha\": 0.65\n",
        "  },\n",
        "  \"sc_MS_wt_2\": {\n",
        "    \"color\": \"#10A105\",\n",
        "    \"label\": \"rep. 2\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.8,\n",
        "    \"alpha\": 0.65\n",
        "  },\n",
        "  \"sc_MS_wt_4\": {\n",
        "    \"color\": \"#e6b800\",\n",
        "    \"label\": \"chemical mapping\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.8,\n",
        "    \"alpha\": 0.65\n",
        "  },\n",
        "  \"sc_N_wt_DS_shuffled_NS_native\": {\n",
        "    \"color\": \"#B49FFA\",\n",
        "    \"label\": \"shuffled native\",\n",
        "    \"linestyle\": (0, (5, 1.5)),\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"sc_N_wt_DS_none_NS_native\": {\n",
        "    \"color\": \"#6060FD\",\n",
        "    \"label\": \"absent native\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"sc_N_wt_DS_native_NS_none\": {\n",
        "    \"color\": \"#F5271C\",\n",
        "    \"label\": \"native absent\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"sc_N_wt_DS_native_NS_shuffled\": {\n",
        "    \"color\": \"#ec1358\",\n",
        "    \"label\": \"native shuffled\",\n",
        "    \"linestyle\": (0, (5, 1.5)),\n",
        "    \"linewidth\": 1.6,\n",
        "    \"alpha\": 0.8\n",
        "  },\n",
        "  \"sc_N_wt_DS_native_NS_native\": {\n",
        "    \"color\": \"#8D05B7\",\n",
        "    \"label\": \"native native\",\n",
        "    \"linestyle\": \"-\",\n",
        "    \"linewidth\": 3.0,\n",
        "    \"alpha\": 0.65\n",
        "  }\n",
        "}\n",
        "# plot negative controls\n",
        "if show_sc_wt_neg_ctrls:\n",
        "  for neg_ctrl in [\n",
        "    \"sc_N_wt_DS_shuffled_NS_shuffled\",\n",
        "    \"sc_N_wt_DS_shuffled_NS_none\",\n",
        "    \"sc_N_wt_DS_none_NS_shuffled\"\n",
        "  ]: _sc_wt_ms_meta.update({\n",
        "    neg_ctrl: {\n",
        "      \"color\": \"k\",\n",
        "      \"label\": \" \".join([neg_ctrl.split(\"_\")[4], neg_ctrl.split(\"_\")[6]]),\n",
        "      \"linestyle\": \"-\",\n",
        "      \"linewidth\": 0.8,\n",
        "      \"alpha\": 0.5\n",
        "    }\n",
        "  })\n",
        "# plot\n",
        "if plot_sc_wt_ms_meta:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  # measured and predicted\n",
        "  for i, details in _sc_wt_ms_meta.items():\n",
        "    A.plot(agg_ms_meta[i].values, **details)\n",
        "  # ancillary\n",
        "  A.set_yticks([])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticks(np.linspace(0, 1024, 9))\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\"]]\n",
        "  plt.grid(axis=\"x\")\n",
        "  plt.legend(loc=2)\n",
        "  plt.tight_layout()\n",
        "  if save_sc_wt_ms_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_ms_meta.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_ms_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t_wnvJc4Yovp",
        "colab_type": "text"
      },
      "source": [
        "#### ridge-line meta"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1LLh0O2zOhmp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sns.set(style=\"white\", rc={\"axes.facecolor\": (0, 0, 0, 0)})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5eDgKfiOmKk0",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae wt\n",
        "plot_sc_wt_ms_ridge_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_wt_ms_meta_ridge_column_palettes = {\n",
        "  \"sc_N_wt_DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "  \"sc_MS_wt_1\": \"#63a463\",\n",
        "  \"sc_MS_wt_2\": \"#39ac39\",\n",
        "  \"sc_MS_wt_4\": \"#e6b800\",\n",
        "  \"sc_N_wt_DS_none_NS_native\": \"#284292\",\n",
        "  \"sc_N_wt_DS_native_NS_none\": \"#CA2128\",\n",
        "  \"sc_N_wt_DS_native_NS_native\": \"#BB478D\"\n",
        "}\n",
        "sc_wt_ms_meta_ridge_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"col\": np.repeat(col, 1024),\n",
        "    \"x\": np.arange(1024),\n",
        "    \"val\": agg_ms_meta[col].values,\n",
        "  }) for col in sc_wt_ms_meta_ridge_column_palettes.keys()\n",
        "], 0)\n",
        "# plot\n",
        "if plot_sc_wt_ms_ridge_meta:\n",
        "  # initialize\n",
        "  g = sns.FacetGrid(\n",
        "    data=sc_wt_ms_meta_ridge_long,\n",
        "    row=\"col\",\n",
        "    hue=\"col\",\n",
        "    aspect=14,\n",
        "    sharex=True,\n",
        "    sharey=True,\n",
        "    xlim=(0, 1024),\n",
        "    ylim=(0, np.max(sc_wt_ms_meta_ridge_long[\"val\"].values)),\n",
        "    height=1.1,\n",
        "    palette=sc_wt_ms_meta_ridge_column_palettes.values()\n",
        "  )\n",
        "  # map to facet grid\n",
        "  g.map(plt.fill_between, \"x\", \"val\", alpha=0.66, clip_on=False)\n",
        "  g.map(plt.plot, \"val\", linewidth=1.9, color=\"k\", clip_on=False)\n",
        "  g.map(plt.axhline, y=0, linewidth=0.5, alpha=0.66, clip_on=False)\n",
        "  # ancillary\n",
        "  g.fig.subplots_adjust(hspace=-.7)\n",
        "  g.set_titles(\"\")\n",
        "  g.set(xticks=[])\n",
        "  g.set(yticks=[])\n",
        "  g.despine(bottom=True, left=True);\n",
        "  if save_sc_wt_ms_ridge_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_ms_ridge_meta.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_ms_ridge_meta.pdf - saved\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rwnr_toAGbWt",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae wt nuc maxima\n",
        "plot_sc_wt_nuc_maxima = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_wt_nuc_maxima = {\n",
        "  \"DS_shuffled_NS_shuffled\": {\n",
        "    \"val\": nuc_data[\"sc_N_wt_DS_shuffled_NS_shuffled\"],\n",
        "    \"clr\": \"#808080\"\n",
        "  },\n",
        "  \"MS_1\": {\n",
        "    \"val\": nuc_data[\"sc_MS_1\"],\n",
        "    \"clr\": \"#63a463\"\n",
        "  },\n",
        "  \"MS_2\": {\n",
        "    \"val\": nuc_data[\"sc_MS_2\"],\n",
        "    \"clr\": \"#39ac39\"\n",
        "  },\n",
        "  \"DS_none_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_wt_DS_none_NS_native\"],\n",
        "    \"clr\": \"#284292\"\n",
        "  },\n",
        "  \"DS_native_NS_none\": {\n",
        "    \"val\": nuc_data[\"sc_N_wt_DS_native_NS_none\"],\n",
        "    \"clr\": \"#CA2128\"\n",
        "  },\n",
        "  \"DS_native_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_wt_DS_native_NS_native\"],\n",
        "    \"clr\": \"#BB478D\"\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_sc_wt_nuc_maxima:\n",
        "  F, A = plt.subplots(1, 1, figsize=(12 * 1.2, 7))\n",
        "  for i, x in enumerate(sc_wt_nuc_maxima):\n",
        "    plt.scatter(\n",
        "      x=sc_wt_nuc_maxima[x][\"val\"],\n",
        "      y=(1 - 0.1 * i) * np.ones(5),\n",
        "      c=sc_wt_nuc_maxima[x][\"clr\"],\n",
        "      s=300,\n",
        "      marker=\"|\"\n",
        "    ) \n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  A.set_yticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\", \"bottom\"]]\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_sc_wt_nuc_maxima:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_nuc_maxima.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_nuc_maxima.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m3cEklNj-jxy",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ ISW1\n",
        "plot_sc_ko_isw1_ms_ridge_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_isw1_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "sys.exit(1)\n",
        "# meta plot details\n",
        "sc_ko_isw1_ms_meta_ridge_column_palettes = {\n",
        "  \"sc_N_ko_isw1_DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "  \"sc_N_ko_isw1_DS_none_NS_native\": \"#284292\",\n",
        "  \"sc_N_ko_isw1_DS_native_NS_none\": \"#CA2128\",\n",
        "  \"sc_N_ko_isw1_DS_native_NS_native\": \"#BB478D\"\n",
        "}\n",
        "sc_ko_isw1_ms_meta_ridge_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"col\": np.repeat(col, 1024),\n",
        "    \"x\": np.arange(1024),\n",
        "    \"val\": agg_ms_meta[col].values\n",
        "  }) for col in sc_ko_isw1_ms_meta_ridge_column_palettes.keys()\n",
        "], 0)\n",
        "# plot\n",
        "if plot_sc_ko_isw1_ms_ridge_meta:\n",
        "  # initialize\n",
        "  g = sns.FacetGrid(\n",
        "    data=sc_ko_isw1_ms_meta_ridge_long,\n",
        "    row=\"col\",\n",
        "    hue=\"col\",\n",
        "    aspect=12,\n",
        "    sharex=True,\n",
        "    sharey=True,\n",
        "    xlim=(0, 1024),\n",
        "    ylim=(0, np.max(sc_ko_isw1_ms_meta_ridge_long[\"val\"].values)),\n",
        "    height=1.2,\n",
        "    palette=sc_ko_isw1_ms_meta_ridge_column_palettes.values()\n",
        "  )\n",
        "  # map to facet grid\n",
        "  g.map(plt.fill_between, \"x\", \"val\", alpha=0.66, clip_on=False)\n",
        "  g.map(plt.plot, \"val\", linewidth=2.1, color=\"k\", clip_on=False)\n",
        "  g.map(plt.axhline, y=0, linewidth=0.5, alpha=0.66, clip_on=False)\n",
        "  # ancillary\n",
        "  g.fig.subplots_adjust(hspace=-.8)\n",
        "  g.set_titles(\"\")\n",
        "  g.set(xticks=[])\n",
        "  g.set(yticks=[])\n",
        "  g.despine(bottom=True, left=True);\n",
        "  if save_sc_ko_isw1_ms_ridge_meta:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"sc_ko_isw1_ms_ridge_meta.pdf\",\n",
        "      format=\"pdf\"\n",
        "    )\n",
        "    print(\"sc_ko_isw1_ms_ridge_meta.pdf - saved\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TYPRsdgaHpBJ",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ ISW1 nuc maxima\n",
        "plot_sc_ko_isw1_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_isw1_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_isw1_nuc_maxima = {\n",
        "  \"DS_shuffled_NS_shuffled\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_isw1_DS_shuffled_NS_shuffled\"],\n",
        "    \"clr\": \"#808080\"\n",
        "  },\n",
        "  \"DS_none_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_isw1_DS_none_NS_native\"],\n",
        "    \"clr\": \"#284292\"\n",
        "  },\n",
        "  \"DS_native_NS_none\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_isw1_DS_native_NS_none\"],\n",
        "    \"clr\": \"#CA2128\"\n",
        "  },\n",
        "  \"DS_native_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_isw1_DS_native_NS_native\"],\n",
        "    \"clr\": \"#BB478D\"\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_sc_ko_isw1_nuc_maxima:\n",
        "  F, A = plt.subplots(1, 1, figsize=(12 * 1.2, 7))\n",
        "  for i, x in enumerate(sc_ko_isw1_nuc_maxima):\n",
        "    plt.scatter(\n",
        "      x=sc_ko_isw1_nuc_maxima[x][\"val\"],\n",
        "      y=(1 - 0.1 * i) * np.ones(5),\n",
        "      c=sc_ko_isw1_nuc_maxima[x][\"clr\"],\n",
        "      s=300,\n",
        "      marker=\"|\"\n",
        "    ) \n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  A.set_yticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\", \"bottom\"]]\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_sc_ko_isw1_nuc_maxima:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_ko_isw1_nuc_maxima.pdf\", format=\"pdf\")\n",
        "    print(\"sc_ko_isw1_nuc_maxima.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GCjvYAwDEGpK",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ CHD1\n",
        "plot_sc_ko_chd1_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_chd1_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_chd1_ms_meta_ridge_column_palettes = {\n",
        "  \"sc_N_ko_chd1_DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "  \"sc_N_ko_chd1_DS_none_NS_native\": \"#284292\",\n",
        "  \"sc_N_ko_chd1_DS_native_NS_none\": \"#CA2128\",\n",
        "  \"sc_N_ko_chd1_DS_native_NS_native\": \"#BB478D\"\n",
        "}\n",
        "sc_ko_chd1_ms_meta_ridge_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"col\": np.repeat(col, 1024),\n",
        "    \"x\": np.arange(1024),\n",
        "    \"val\": agg_ms_meta[col].values\n",
        "  }) for col in sc_ko_chd1_ms_meta_ridge_column_palettes.keys()\n",
        "], 0)\n",
        "# plot\n",
        "if plot_sc_ko_chd1_ms_ridge_meta:\n",
        "  # initialize\n",
        "  g = sns.FacetGrid(\n",
        "    data=sc_ko_chd1_ms_meta_ridge_long,\n",
        "    row=\"col\",\n",
        "    hue=\"col\",\n",
        "    aspect=12,\n",
        "    sharex=True,\n",
        "    sharey=True,\n",
        "    xlim=(0, 1024),\n",
        "    ylim=(0, np.max(sc_ko_chd1_ms_meta_ridge_long[\"val\"].values)),\n",
        "    height=1.2,\n",
        "    palette=sc_ko_chd1_ms_meta_ridge_column_palettes.values()\n",
        "  )\n",
        "  # map to facet grid\n",
        "  g.map(plt.fill_between, \"x\", \"val\", alpha=0.66, clip_on=False)\n",
        "  g.map(plt.plot, \"val\", linewidth=2.1, color=\"k\", clip_on=False)\n",
        "  g.map(plt.axhline, y=0, linewidth=0.5, alpha=0.66, clip_on=False)\n",
        "  # ancillary\n",
        "  g.fig.subplots_adjust(hspace=-.8)\n",
        "  g.set_titles(\"\")\n",
        "  g.set(xticks=[])\n",
        "  g.set(yticks=[])\n",
        "  g.despine(bottom=True, left=True);\n",
        "  if save_sc_ko_chd1_ms_ridge_meta:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"sc_ko_chd1_ms_ridge_meta.pdf\",\n",
        "      format=\"pdf\"\n",
        "    )\n",
        "    print(\"sc_ko_chd1_ms_ridge_meta.pdf - saved\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "toBJ8rZkERfq",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ CHD1 nuc maxima\n",
        "plot_sc_ko_chd1_nuc_maxima = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_chd1_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_chd1_nuc_maxima = {\n",
        "  \"DS_shuffled_NS_shuffled\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_chd1_DS_shuffled_NS_shuffled\"],\n",
        "    \"clr\": \"#808080\"\n",
        "  },\n",
        "  \"DS_none_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_chd1_DS_none_NS_native\"],\n",
        "    \"clr\": \"#284292\"\n",
        "  },\n",
        "  \"DS_native_NS_none\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_chd1_DS_native_NS_none\"],\n",
        "    \"clr\": \"#CA2128\"\n",
        "  },\n",
        "  \"DS_native_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_chd1_DS_native_NS_native\"],\n",
        "    \"clr\": \"#BB478D\"\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_sc_ko_chd1_nuc_maxima:\n",
        "  F, A = plt.subplots(1, 1, figsize=(12 * 1.2, 7))\n",
        "  for i, x in enumerate(sc_ko_chd1_nuc_maxima):\n",
        "    plt.scatter(\n",
        "      x=sc_ko_chd1_nuc_maxima[x][\"val\"],\n",
        "      y=(1 - 0.1 * i) * np.ones(5),\n",
        "      c=sc_ko_chd1_nuc_maxima[x][\"clr\"],\n",
        "      s=300,\n",
        "      marker=\"|\"\n",
        "    ) \n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  A.set_yticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\", \"bottom\"]]\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_sc_ko_chd1_nuc_maxima:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_ko_chd1_nuc_maxima.pdf\", format=\"pdf\")\n",
        "    print(\"sc_ko_chd1_nuc_maxima.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9U75a0KvGoct",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ INO80\n",
        "plot_sc_ko_ino80_ms_ridge_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_ino80_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_ino80_ms_meta_ridge_column_palettes = {\n",
        "  \"sc_N_ko_ino80_DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "  \"sc_N_ko_ino80_DS_none_NS_native\": \"#284292\",\n",
        "  \"sc_N_ko_ino80_DS_native_NS_none\": \"#CA2128\",\n",
        "  \"sc_N_ko_ino80_DS_native_NS_native\": \"#BB478D\"\n",
        "}\n",
        "sc_ko_ino80_ms_meta_ridge_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"col\": np.repeat(col, 1024),\n",
        "    \"x\": np.arange(1024),\n",
        "    \"val\": agg_ms_meta[col].values\n",
        "  }) for col in sc_ko_ino80_ms_meta_ridge_column_palettes.keys()\n",
        "], 0)\n",
        "# plot\n",
        "if plot_sc_ko_ino80_ms_ridge_meta:\n",
        "  # initialize\n",
        "  g = sns.FacetGrid(\n",
        "    data=sc_ko_ino80_ms_meta_ridge_long,\n",
        "    row=\"col\",\n",
        "    hue=\"col\",\n",
        "    aspect=12,\n",
        "    sharex=True,\n",
        "    sharey=True,\n",
        "    xlim=(0, 1024),\n",
        "    ylim=(0, np.max(sc_ko_ino80_ms_meta_ridge_long[\"val\"].values)),\n",
        "    height=1.2,\n",
        "    palette=sc_ko_ino80_ms_meta_ridge_column_palettes.values()\n",
        "  )\n",
        "  # map to facet grid\n",
        "  g.map(plt.fill_between, \"x\", \"val\", alpha=0.66, clip_on=False)\n",
        "  g.map(plt.plot, \"val\", linewidth=2.1, color=\"k\", clip_on=False)\n",
        "  g.map(plt.axhline, y=0, linewidth=0.5, alpha=0.66, clip_on=False)\n",
        "  # ancillary\n",
        "  g.fig.subplots_adjust(hspace=-.8)\n",
        "  g.set_titles(\"\")\n",
        "  g.set(xticks=[])\n",
        "  g.set(yticks=[])\n",
        "  g.despine(bottom=True, left=True);\n",
        "  if save_sc_ko_ino80_ms_ridge_meta:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"sc_ko_ino80_ms_ridge_meta.pdf\",\n",
        "      format=\"pdf\"\n",
        "    )\n",
        "    print(\"sc_ko_ino80_ms_ridge_meta.pdf - saved\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l-yScUXJGojs",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae Δ INO80 nuc maxima\n",
        "plot_sc_ko_ino80_nuc_maxima = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_ino80_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_ino80_nuc_maxima = {\n",
        "  \"DS_shuffled_NS_shuffled\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_ino80_DS_shuffled_NS_shuffled\"],\n",
        "    \"clr\": \"#808080\"\n",
        "  },\n",
        "  \"DS_none_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_ino80_DS_none_NS_native\"],\n",
        "    \"clr\": \"#284292\"\n",
        "  },\n",
        "  \"DS_native_NS_none\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_ino80_DS_native_NS_none\"],\n",
        "    \"clr\": \"#CA2128\"\n",
        "  },\n",
        "  \"DS_native_NS_native\": {\n",
        "    \"val\": nuc_data[\"sc_N_ko_ino80_DS_native_NS_native\"],\n",
        "    \"clr\": \"#BB478D\"\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_sc_ko_ino80_nuc_maxima:\n",
        "  F, A = plt.subplots(1, 1, figsize=(12 * 1.2, 7))\n",
        "  for i, x in enumerate(sc_ko_ino80_nuc_maxima):\n",
        "    plt.scatter(\n",
        "      x=sc_ko_ino80_nuc_maxima[x][\"val\"],\n",
        "      y=(1 - 0.1 * i) * np.ones(5),\n",
        "      c=sc_ko_ino80_nuc_maxima[x][\"clr\"],\n",
        "      s=300,\n",
        "      marker=\"|\"\n",
        "    ) \n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  A.set_yticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\", \"bottom\"]]\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_sc_ko_ino80_nuc_maxima:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_ko_ino80_nuc_maxima.pdf\", format=\"pdf\")\n",
        "    print(\"sc_ko_ino80_nuc_maxima.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NnGFTF4pFhRb",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis\n",
        "plot_kl_ms_ridge_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_ms_ridge_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "kl_ms_meta_ridge_column_palettes = {\n",
        "  \"kl_DS_shuffled_NS_shuffled\": \"#808080\",\n",
        "  \"kl_MS_1\": \"#63a463\",\n",
        "  \"kl_DS_none_NS_native\": \"#284292\",\n",
        "  \"kl_DS_native_NS_none\": \"#CA2128\",\n",
        "  \"kl_DS_native_NS_native\": \"#BB478D\"\n",
        "}\n",
        "kl_ms_meta_ridge_long = pd.concat([\n",
        "  pd.DataFrame({\n",
        "    \"col\": np.repeat(col, 1024),\n",
        "    \"x\": np.arange(1024),\n",
        "    \"val\": agg_ms_meta[col].values,\n",
        "  }) for col in kl_ms_meta_ridge_column_palettes.keys()\n",
        "], 0)\n",
        "# plot\n",
        "if plot_kl_ms_ridge_meta:\n",
        "  # initialize\n",
        "  g = sns.FacetGrid(\n",
        "    data=kl_ms_meta_ridge_long,\n",
        "    row=\"col\",\n",
        "    hue=\"col\",\n",
        "    aspect=14,\n",
        "    sharex=True,\n",
        "    sharey=True,\n",
        "    xlim=(0, 1024),\n",
        "    ylim=(0, np.max(kl_ms_meta_ridge_long[\"val\"].values)),\n",
        "    height=1.1,\n",
        "    palette=kl_ms_meta_ridge_column_palettes.values()\n",
        "  )\n",
        "  # map to facet grid\n",
        "  g.map(plt.fill_between, \"x\", \"val\", alpha=0.66, clip_on=False)\n",
        "  g.map(plt.plot, \"val\", linewidth=1.9, color=\"k\", clip_on=False)\n",
        "  g.map(plt.axhline, y=0, linewidth=0.5, alpha=0.66, clip_on=False)\n",
        "  # ancillary\n",
        "  g.fig.subplots_adjust(hspace=-.7)\n",
        "  g.set_titles(\"\")\n",
        "  g.set(xticks=[])\n",
        "  g.set(yticks=[])\n",
        "  g.despine(bottom=True, left=True);\n",
        "  if save_kl_ms_ridge_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"kl_ms_ridge_meta.pdf\", format=\"pdf\")\n",
        "    print(\"kl_ms_ridge_meta.pdf - saved\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b1kf5YY1GZVE",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis nuc maxima\n",
        "plot_kl_nuc_maxima = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_nuc_maxima = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "kl_nuc_maxima = {\n",
        "  \"DS_shuffled_NS_shuffled\": {\n",
        "    \"val\": nuc_data[\"kl_DS_shuffled_NS_shuffled\"],\n",
        "    \"clr\": \"#808080\"\n",
        "  },\n",
        "  \"MS_1\": {\n",
        "    \"val\": nuc_data[\"kl_MS_1\"],\n",
        "    \"clr\": \"#63a463\"\n",
        "  },\n",
        "  \"DS_none_NS_native\": {\n",
        "    \"val\": nuc_data[\"kl_DS_none_NS_native\"],\n",
        "    \"clr\": \"#284292\"\n",
        "  },\n",
        "  \"DS_native_NS_none\": {\n",
        "    \"val\": nuc_data[\"kl_DS_native_NS_none\"],\n",
        "    \"clr\": \"#CA2128\"\n",
        "  },\n",
        "  \"DS_native_NS_native\": {\n",
        "    \"val\": nuc_data[\"kl_DS_native_NS_native\"],\n",
        "    \"clr\": \"#BB478D\"\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_kl_nuc_maxima:\n",
        "  F, A = plt.subplots(1, 1, figsize=(12 * 1.2, 7))\n",
        "  for i, x in enumerate(kl_nuc_maxima):\n",
        "    plt.scatter(\n",
        "      x=kl_nuc_maxima[x][\"val\"],\n",
        "      y=(1 - 0.1 * i) * np.ones(5),\n",
        "      c=kl_nuc_maxima[x][\"clr\"],\n",
        "      s=300,\n",
        "      marker=\"|\"\n",
        "    ) \n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_ylabel(None)\n",
        "  A.set_xticklabels([\"\"] * 9)\n",
        "  A.set_yticklabels([\"\"] * 9)\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\", \"bottom\"]]\n",
        "  plt.grid(axis=\"y\")\n",
        "  plt.tight_layout()\n",
        "  if save_kl_nuc_maxima:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"kl_nuc_maxima.pdf\", format=\"pdf\")\n",
        "    print(\"kl_nuc_maxima.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CkOUInXY-A2Z",
        "colab_type": "text"
      },
      "source": [
        "#### heatmap meta"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hWUJdJ7-L_Iu",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title Sc wt\n",
        "plot_sc_wt_ms_heatmap = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_ms_heatmap = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_wt_ms_heatmap_sort_idx = np.argsort(\n",
        "  np.sum(agg_ns_data[\"sc_wt\"][:, 0, :], 1)\n",
        ")[::-1]\n",
        "sc_wt_ms_heatmap_data = {\n",
        "  \"wt_DS_shuffled_NS_shuffled\": np.array([x\n",
        "    for x in agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"][\"DS_shuffled_NS_shuffled\"][\"avg\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :],\n",
        "  \"wt_DS_none_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"][\"DS_none_NS_native\"][\"avg\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :],\n",
        "  \"wt_DS_native_NS_none\": np.array([x\n",
        "    for x in agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"][\"DS_native_NS_none\"][\"avg\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :],\n",
        "  \"wt_DS_native_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\"sc\"][\"pred\"][\"N_wt\"][\"DS_native_NS_native\"][\"avg\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :],\n",
        "  \"wt_1\": np.array([_np_ms_transform(x)\n",
        "    for x in agg_ms_data[\"sc\"][\"meas\"][\"MS_1\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :],\n",
        "  \"wt_2\": np.array([_np_ms_transform(x)\n",
        "    for x in agg_ms_data[\"sc\"][\"meas\"][\"MS_2\"]\n",
        "  ])[sc_wt_ms_heatmap_sort_idx, :]\n",
        "}\n",
        "sc_wt_ms_heatmap_extrema = {\n",
        "  \"min\": np.min(np.concatenate([\n",
        "    x for x in sc_wt_ms_heatmap_data.values()\n",
        "  ])),\n",
        "  \"max\": np.max(np.concatenate([\n",
        "    x for x in sc_wt_ms_heatmap_data.values()\n",
        "  ]))\n",
        "}\n",
        "if plot_sc_wt_ms_heatmap:\n",
        "  F, A = plt.subplots(1, len(sc_wt_ms_heatmap_data), figsize=(12, 10))\n",
        "  for i, (m, m_data) in enumerate(sc_wt_ms_heatmap_data.items()):\n",
        "    print(i, m)\n",
        "    sns.heatmap(\n",
        "      data=m_data,\n",
        "      cbar=False,\n",
        "      vmin=sc_wt_ms_heatmap_extrema[\"min\"],\n",
        "      vmax=sc_wt_ms_heatmap_extrema[\"max\"],\n",
        "      cmap=\"ocean_r\",\n",
        "      ax=A[i]\n",
        "    )\n",
        "    A[i].set_xticks([]); A[i].set_yticks([])\n",
        "    A[i].set_title(m, fontdict={\"fontsize\": 8})\n",
        "  plt.tight_layout(w_pad=-0.3)\n",
        "  if save_sc_wt_ms_heatmap:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"sc_wt_ms_meta_heatmap.png\",\n",
        "      format=\"png\",\n",
        "      dpi=1024\n",
        "    )\n",
        "    print(\"sc_wt_ms_meta_heatmap.png - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B1drdydfkJnd",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title Sc ko\n",
        "ko = \"chd1\" #@param [\"isw1\", \"ino80\", \"chd1\"]\n",
        "plot_sc_ko_ms_heatmap = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_ms_heatmap = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "sc_ko_ms_heatmap_sort_idx = np.argsort(\n",
        "  np.sum(agg_ns_data[\"sc_ko_\" + ko][:, 0, :], 1)\n",
        ")[::-1]\n",
        "sc_ko_ms_heatmap_data = {\n",
        "  \"ko_DS_shuffled_NS_shuffled\": np.array([x\n",
        "    for x in agg_ms_data[\n",
        "      \"sc\"][\"pred\"][\"N_ko_\" + ko][\"DS_shuffled_NS_shuffled\"][\"avg\"]\n",
        "  ])[sc_ko_ms_heatmap_sort_idx, :],\n",
        "  \"ko_DS_none_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\n",
        "      \"sc\"][\"pred\"][\"N_ko_\" + ko][\"DS_none_NS_native\"][\"avg\"]\n",
        "  ])[sc_ko_ms_heatmap_sort_idx, :],\n",
        "  \"ko_DS_native_NS_none\": np.array([x\n",
        "    for x in agg_ms_data[\n",
        "      \"sc\"][\"pred\"][\"N_ko_\" + ko][\"DS_native_NS_none\"][\"avg\"]\n",
        "  ])[sc_ko_ms_heatmap_sort_idx, :],\n",
        "  \"ko_DS_native_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\n",
        "      \"sc\"][\"pred\"][\"N_ko_\" + ko][\"DS_native_NS_native\"][\"avg\"]\n",
        "  ])[sc_ko_ms_heatmap_sort_idx, :]\n",
        "}\n",
        "sc_ko_ms_heatmap_extrema = {\n",
        "  \"min\": np.min(np.concatenate([\n",
        "    x for x in sc_ko_ms_heatmap_data.values()\n",
        "  ])),\n",
        "  \"max\": np.max(np.concatenate([\n",
        "    x for x in sc_ko_ms_heatmap_data.values()\n",
        "  ]))\n",
        "}\n",
        "if plot_sc_ko_ms_heatmap:\n",
        "  F, A = plt.subplots(1, len(sc_ko_ms_heatmap_data), figsize=(8, 10))\n",
        "  for i, (m, m_data) in enumerate(sc_ko_ms_heatmap_data.items()):\n",
        "    print(i, m)\n",
        "    sns.heatmap(\n",
        "      data=m_data,\n",
        "      cbar=False,\n",
        "      vmin=sc_ko_ms_heatmap_extrema[\"min\"],\n",
        "      vmax=sc_ko_ms_heatmap_extrema[\"max\"],\n",
        "      cmap=\"ocean_r\",\n",
        "      ax=A[i]\n",
        "    )\n",
        "    A[i].set_xticks([]); A[i].set_yticks([])\n",
        "    A[i].set_title(m, fontdict={\"fontsize\": 8})\n",
        "  plt.tight_layout(w_pad=-0.3)\n",
        "  if save_sc_ko_ms_heatmap:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"sc_ko_\" + ko + \"_ms_meta_heatmap.png\",\n",
        "      format=\"png\",\n",
        "      dpi=1024\n",
        "    )\n",
        "    print(\"sc_ko_\" + ko + \"_ms_meta_heatmap.png - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eZkIrSfGSFLs",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis\n",
        "plot_kl_ms_heatmap = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_ms_heatmap = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "kl_ms_heatmap_sort_idx = np.argsort(\n",
        "  np.sum(agg_ns_data[\"kl\"][:, 0, :], 1)\n",
        ")[::-1]\n",
        "kl_ms_heatmap_data = {\n",
        "  \"kl_DS_shuffled_NS_shuffled\": np.array([x\n",
        "    for x in agg_ms_data[\"kl\"][\"pred\"][\"DS_shuffled_NS_shuffled\"][\"avg\"]\n",
        "  ])[kl_ms_heatmap_sort_idx, :],\n",
        "  \"kl_DS_none_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\"kl\"][\"pred\"][\"DS_none_NS_native\"][\"avg\"]\n",
        "  ])[kl_ms_heatmap_sort_idx, :],\n",
        "  \"kl_DS_native_NS_none\": np.array([x\n",
        "    for x in agg_ms_data[\"kl\"][\"pred\"][\"DS_native_NS_none\"][\"avg\"]\n",
        "  ])[kl_ms_heatmap_sort_idx, :],\n",
        "  \"kl_DS_native_NS_native\": np.array([x\n",
        "    for x in agg_ms_data[\"kl\"][\"pred\"][\"DS_native_NS_native\"][\"avg\"]\n",
        "  ])[kl_ms_heatmap_sort_idx, :],\n",
        "  \"kl_wt_1\": np.array([_np_ms_transform(x)\n",
        "    for x in agg_ms_data[\"kl\"][\"meas\"][\"MS_1\"]\n",
        "  ])[kl_ms_heatmap_sort_idx, :]\n",
        "}   \n",
        "kl_ms_heatmap_vmax = np.max(np.concatenate([\n",
        "  x for x in kl_ms_heatmap_data.values()\n",
        "]))\n",
        "# plot\n",
        "if plot_kl_ms_heatmap:\n",
        "  F, A = plt.subplots(1, len(kl_ms_heatmap_data), figsize=(10, 10))\n",
        "  for i, (m, m_data) in enumerate(kl_ms_heatmap_data.items()):\n",
        "    print(i, m)\n",
        "    sns.heatmap(\n",
        "      data=m_data,\n",
        "      cbar=False,\n",
        "      vmin=0,\n",
        "      vmax=kl_ms_heatmap_vmax,\n",
        "      cmap=\"ocean_r\",\n",
        "      ax=A[i]\n",
        "    )\n",
        "    A[i].set_xticks([]); A[i].set_yticks([])\n",
        "    A[i].set_title(m, fontdict={\"fontsize\": 4})\n",
        "  plt.tight_layout(w_pad=-0.3)\n",
        "  if save_kl_ms_heatmap:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"kl_ms_meta_heatmap.png\",\n",
        "      format=\"png\", dpi=1024\n",
        "    )\n",
        "    print(\"kl_ms_meta_heatmap.png - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tar6re3PzUwf",
        "colab_type": "text"
      },
      "source": [
        "#### nucleosome positioning"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SavSKjMUXlcc",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title nucleosome location and distances\n",
        "nuc_dist_all = {\n",
        "  x: {\n",
        "    \"N\": nuc_data[x],\n",
        "    \"D\": np.array([\n",
        "      nuc_data[x][i] - nuc_data[x][i - 1]\n",
        "        for i in range(1, len(nuc_data[x]))\n",
        "    ])\n",
        "  } for x in nuc_data.keys()\n",
        "}\n",
        "# absolute differences\n",
        "  # Sc wt wrt Sc meas 2\n",
        "sc_wt_wrt_MS_2 = {\n",
        "  x: {\n",
        "    \"N\": nuc_dist_all[x][\"N\"] - nuc_dist_all[\"sc_MS_2\"][\"N\"],\n",
        "    \"D\": nuc_dist_all[x][\"D\"] - nuc_dist_all[\"sc_MS_2\"][\"D\"]\n",
        "  } for x in [\n",
        "    \"sc_MS_1\",\n",
        "    \"sc_N_wt_DS_shuffled_NS_shuffled\",\n",
        "    \"sc_N_wt_DS_none_NS_native\",\n",
        "    \"sc_N_wt_DS_native_NS_none\",\n",
        "    \"sc_N_wt_DS_native_NS_native\"\n",
        "  ]\n",
        "}\n",
        "  # Kl wrt Kl meas 1\n",
        "kl_wrt_kl_MS_1 = {\n",
        "  x: {\n",
        "    \"N\": nuc_dist_all[x][\"N\"] - nuc_dist_all[\"kl_MS_1\"][\"N\"],\n",
        "    \"D\": nuc_dist_all[x][\"D\"] - nuc_dist_all[\"kl_MS_1\"][\"D\"]\n",
        "  } for x in [\n",
        "    \"kl_DS_shuffled_NS_shuffled\",\n",
        "    \"kl_DS_none_NS_native\",\n",
        "    \"kl_DS_native_NS_none\",\n",
        "    \"kl_DS_native_NS_native\"\n",
        "  ]\n",
        "}\n",
        "# relative differences\n",
        "  # Kl wrt Sc\n",
        "kl_wrt_sc = {\n",
        "  (x, y): {\n",
        "    \"N\": nuc_dist_all[x][\"N\"] - nuc_dist_all[y][\"N\"],\n",
        "    \"D\": nuc_dist_all[x][\"D\"] - nuc_dist_all[y][\"D\"]\n",
        "  } for x, y in {\n",
        "    \"kl_MS_1\": \"sc_MS_1\",\n",
        "    \"kl_DS_shuffled_NS_shuffled\": \"sc_N_wt_DS_shuffled_NS_shuffled\",\n",
        "    \"kl_DS_none_NS_native\": \"sc_N_wt_DS_none_NS_native\",\n",
        "    \"kl_DS_native_NS_none\": \"sc_N_wt_DS_native_NS_none\",\n",
        "    \"kl_DS_native_NS_native\": \"sc_N_wt_DS_native_NS_native\"\n",
        "  }.items()\n",
        "}\n",
        "  # Sc ko wrt Sc wt\n",
        "sc_ko_wrt_sc_wt = {\n",
        "  ko: {\n",
        "    (x, y): {\n",
        "      \"N\": nuc_dist_all[x][\"N\"] - nuc_dist_all[y][\"N\"],\n",
        "      \"D\": nuc_dist_all[x][\"D\"] - nuc_dist_all[y][\"D\"]\n",
        "    } for x, y in ko_val.items()\n",
        "  } for ko, ko_val in {\n",
        "    ko: dict(zip(\n",
        "      [\"sc_N_ko_\" + ko + \"_DS_shuffled_NS_shuffled\",\n",
        "       \"sc_N_ko_\" + ko + \"_DS_none_NS_native\",\n",
        "       \"sc_N_ko_\" + ko + \"_DS_native_NS_none\",\n",
        "       \"sc_N_ko_\" + ko + \"_DS_native_NS_native\"],\n",
        "      [\"sc_N_wt_DS_shuffled_NS_shuffled\",\n",
        "       \"sc_N_wt_DS_none_NS_native\",\n",
        "       \"sc_N_wt_DS_native_NS_none\",\n",
        "       \"sc_N_wt_DS_native_NS_native\"]\n",
        "    )) for ko in [\"isw1\", \"chd1\", \"ino80\"]\n",
        "  }.items()\n",
        "}"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GdjASXm32ajX",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title absolute / relative differences\n",
        "  # absolute\n",
        "    # Sc\n",
        "print(\"K. lactis - predictions minus measured\")\n",
        "for x, y in sc_wt_wrt_MS_2.items():\n",
        "  sum_N = np.sum(np.abs(y[\"N\"]))\n",
        "  sum_D = np.sum(np.abs(y[\"D\"]))\n",
        "  print(\"  \", x, \"...\", sum_N + sum_D)\n",
        "  print(\"\\tnucs:\", y[\"N\"], \"...\", sum_N)\n",
        "  print(\"\\tdist:\", y[\"D\"], \"...\", sum_D)\n",
        "    # Kl\n",
        "print(\"\\nK. lactis - predictions minus measured\")\n",
        "for x, y in kl_wrt_kl_MS_1.items():\n",
        "  sum_N = np.sum(np.abs(y[\"N\"]))\n",
        "  sum_D = np.sum(np.abs(y[\"D\"]))\n",
        "  print(\"  \", x, \"...\", sum_N + sum_D)\n",
        "  print(\"\\tnucs:\", y[\"N\"], \"...\", sum_N)\n",
        "  print(\"\\tdist:\", y[\"D\"], \"...\", sum_D)\n",
        "  # relative\n",
        "    # Kl vs Sc\n",
        "print(\"\\nK. lactis (pred / meas) minus S. cerevisiae (pred / meas)\")\n",
        "for (x, y), ND in kl_wrt_sc.items():\n",
        "  print(\"  \", x, \"...\", y)\n",
        "  print(\"\\tdist:\", ND[\"D\"], \"...\", np.around(np.sum(ND[\"D\"],), 2))\n",
        "    # Sc - ko vs wt\n",
        "print(\"\\nS. cerevisiae - knockout minus wildtype\")\n",
        "for ko, ko_val in sc_ko_wrt_sc_wt.items():\n",
        "  print(\"   ko:\", ko)\n",
        "  for (x, y), ND in ko_val.items():\n",
        "    sum_N = np.sum(np.abs(ND[\"N\"]))\n",
        "    sum_D = np.sum(np.abs(ND[\"D\"]))\n",
        "    print(\"     \", (x, y), \"...\", sum_N + sum_D)\n",
        "    print(\"\\tnucs:\", ND[\"N\"], \"...\", sum_N)\n",
        "    print(\"\\tdist:\", ND[\"D\"], \"...\", sum_D)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TGqN3e7wx80Z",
        "colab_type": "text"
      },
      "source": [
        "### DNA sequence"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "neqGY4NQyokJ",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title A/T content\n",
        "plot_ds_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_ds_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "_ds_meta = {\n",
        "  \"cmap\": \"RdGy\", #\"plasma_r\",\n",
        "  \"vmin\": np.min(np.concatenate([\n",
        "    agg_ds_meta[species].values.flatten()\n",
        "      for species in agg_ds_meta.keys()\n",
        "  ])),\n",
        "  \"vmax\": np.max(np.concatenate([\n",
        "    agg_ds_meta[species].values.flatten()\n",
        "      for species in agg_ds_meta.keys()\n",
        "  ])),\n",
        "  \"cbar_kws\": {\n",
        "    \"pad\": 0.01,\n",
        "    \"shrink\": 0.9,\n",
        "    \"aspect\": 8,\n",
        "    \"ticks\": [44, 82]\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_ds_meta:\n",
        "  F, A = plt.subplots(agg_ds_meta.shape[1], 1, figsize=(11, 3))\n",
        "  for i, species in enumerate(agg_ds_meta.keys()):\n",
        "    sns.heatmap(\n",
        "      np.expand_dims(agg_ds_meta[species], 0),\n",
        "      **_ds_meta, ax=A[i]\n",
        "    )\n",
        "    A[i].set_xticks([])\n",
        "    A[i].set_xlim([0, 1024])\n",
        "    A[i].set_yticks([])\n",
        "  # ancillary\n",
        "  plt.tight_layout(pad=0.001)\n",
        "  if save_ds_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"ds_meta.pdf\", format=\"pdf\")\n",
        "    print(\"ds_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9ASlVFZqymVx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Summary Statistics of A/T content\n",
        "agg_ds_meta.describe()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UFYYnl1MlskB",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title Kl vs Sc seq. id comparison\n",
        "plot_kl_sc_ds_id = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_sc_ds_id = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# meta plot details\n",
        "_kl_sc_ds_id = {\n",
        "  \"data\": kl_sc_seq_identity.flatten(),\n",
        "  \"bins\": 64,\n",
        "  \"hist\": False,\n",
        "  \"kde\": True,\n",
        "  \"kde_kws\": {\n",
        "    \"shade\": True\n",
        "  }\n",
        "}\n",
        "# plot\n",
        "if plot_kl_sc_ds_id:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  sns.distplot(\n",
        "    a=_kl_sc_ds_id[\"data\"],\n",
        "    bins=_kl_sc_ds_id[\"bins\"],\n",
        "    hist=_kl_sc_ds_id[\"hist\"],\n",
        "    kde=_kl_sc_ds_id[\"kde\"],\n",
        "    kde_kws=_kl_sc_ds_id[\"kde_kws\"],\n",
        "    ax=A\n",
        "  )\n",
        "  A.set_xticks(list(np.arange(20, 60, 5)))\n",
        "  A.set_xlim([\n",
        "    np.floor(np.min(_kl_sc_ds_id[\"data\"])),\n",
        "    np.ceil(np.max(_kl_sc_ds_id[\"data\"]))\n",
        "  ])\n",
        "  A.tick_params(width=1.5, length=8)\n",
        "  # ancillary\n",
        "  plt.tight_layout(pad=0.01)\n",
        "  if save_kl_sc_ds_id:\n",
        "    plt.savefig(\n",
        "      fname=_dirs[\"figs_stats\"] + \"kl_sc_seq_identity.pdf\",\n",
        "      transparent=True,\n",
        "      format=\"pdf\"\n",
        "    )\n",
        "    print(\"kl_sc_seq_identity.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jTanQKTpPbj7",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title describe seq. identity\n",
        "for compare, mtx in {\n",
        "  \"kl_kl\": kl_kl_seq_identity,\n",
        "  \"sc_sc\": sc_sc_seq_identity,\n",
        "  \"kl_sc\": kl_sc_seq_identity,\n",
        "}.items():\n",
        "  print(compare, mtx.shape); mtx = mtx.flatten()\n",
        "    # min, max\n",
        "  print(\"min, max:\", np.floor(np.min(mtx)), np.ceil(np.max(mtx)))\n",
        "    # median, IQR\n",
        "  print(\n",
        "    \"25th, 50th, 75th %:\",\n",
        "    np.around(np.percentile(mtx, [25, 50, 75]), 1)\n",
        "  )\n",
        "    # mean, std dev.\n",
        "  print(\n",
        "    \"avg:\", np.around(np.average(mtx), 1),\n",
        "    \"std:\", np.around(np.std(mtx), 1), \"\\n\"\n",
        "  )"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cWyoUiZHyDzQ",
        "colab_type": "text"
      },
      "source": [
        "### NET-seq"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZvNrKCzjDWm8",
        "colab_type": "text"
      },
      "source": [
        "#### canonical meta"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8TfhjYne_fCf",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title K. lactis\n",
        "plot_kl_ns_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_kl_ns_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "_kl_ns_meta = {\n",
        "  \"color\": \"#2c3b96\",\n",
        "  \"linewidth\": 2.5,\n",
        "  \"ylim\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[\"kl\"][\"antisense\"], 0))),\n",
        "    np.ceil(np.max(np.sum(ns_meta[\"kl\"][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"yticks\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[\"kl\"][\"antisense\"], 0))),\n",
        "    0,\n",
        "    np.ceil(np.max(np.sum(ns_meta[\"kl\"][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"sense\": np.sum(ns_meta[\"kl\"][\"sense\"], 0),\n",
        "  \"antisense\": np.sum(ns_meta[\"kl\"][\"antisense\"], 0)\n",
        "}\n",
        "# plot\n",
        "if plot_kl_ns_meta:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  # sense signal\n",
        "  A.plot(\n",
        "    y=_kl_ns_meta[\"sense\"],\n",
        "    color=_kl_ns_meta[\"color\"],\n",
        "    linewidth=_kl_ns_meta[\"linewidth\"]\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=_kl_ns_meta[\"sense\"],\n",
        "    color=_kl_ns_meta[\"color\"],\n",
        "    linewidth=_kl_ns_meta[\"linewidth\"],\n",
        "    alpha=0.6\n",
        "  )\n",
        "  # antisense signal\n",
        "  A.plot(\n",
        "    y=-_kl_ns_meta[\"antisense\"],\n",
        "    color=_kl_ns_meta[\"color\"],\n",
        "    linewidth=_kl_ns_meta[\"linewidth\"],\n",
        "    alpha=0.8\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=-_kl_ns_meta[\"antisense\"],\n",
        "    y2=np.zeros(1024),\n",
        "    color=_kl_ns_meta[\"color\"],\n",
        "    linewidth=_kl_ns_meta[\"linewidth\"],\n",
        "    alpha=0.3\n",
        "  )\n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_xticks([])\n",
        "  A.set_ylim(_kl_ns_meta[\"ylim\"])\n",
        "  A.set_yticks(_kl_ns_meta[\"yticks\"])\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\"]]\n",
        "  plt.tight_layout()\n",
        "  # save\n",
        "  if save_kl_ns_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"kl_ns_meta.pdf\", format=\"pdf\")\n",
        "    print(\"kl_ns_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3ff6elvGxPVe",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae wt\n",
        "plot_sc_wt_ns_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_wt_ns_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "_sc_wt_ns_meta = {\n",
        "  \"color\": \"#2c3b96\",\n",
        "  \"linewidth\": 2.5,\n",
        "  \"ylim\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[\"sc_wt\"][\"antisense\"], 0))),\n",
        "    np.ceil(np.max(np.sum(ns_meta[\"sc_wt\"][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"yticks\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[\"sc_wt\"][\"antisense\"], 0))),\n",
        "    0,\n",
        "    np.ceil(np.max(np.sum(ns_meta[\"sc_wt\"][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"sense\": np.sum(ns_meta[\"sc_wt\"][\"sense\"], 0),\n",
        "  \"antisense\": np.sum(ns_meta[\"sc_wt\"][\"antisense\"], 0)\n",
        "}\n",
        "# plot\n",
        "if plot_sc_wt_ns_meta:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  # sense signal\n",
        "  A.plot(\n",
        "    y=_sc_wt_ns_meta[\"sense\"],\n",
        "    color=_sc_wt_ns_meta[\"color\"],\n",
        "    linewidth=_sc_wt_ns_meta[\"linewidth\"]\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=_sc_wt_ns_meta[\"sense\"],\n",
        "    color=_sc_wt_ns_meta[\"color\"],\n",
        "    linewidth=_sc_wt_ns_meta[\"linewidth\"],\n",
        "    alpha=0.6\n",
        "  )\n",
        "  # antisense signal\n",
        "  A.plot(\n",
        "    y=-_sc_wt_ns_meta[\"antisense\"],\n",
        "    color=_sc_wt_ns_meta[\"color\"],\n",
        "    linewidth=_sc_wt_ns_meta[\"linewidth\"],\n",
        "    alpha=0.8\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=-_sc_wt_ns_meta[\"antisense\"],\n",
        "    y2=np.zeros(1024),\n",
        "    color=_sc_wt_ns_meta[\"color\"],\n",
        "    linewidth=_sc_wt_ns_meta[\"linewidth\"],\n",
        "    alpha=0.3\n",
        "  )\n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_xticks([])\n",
        "  A.set_ylim(_sc_wt_ns_meta[\"ylim\"])\n",
        "  A.set_yticks(_sc_wt_ns_meta[\"yticks\"])\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\"]]\n",
        "  plt.tight_layout()\n",
        "  # save\n",
        "  if save_sc_wt_ns_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"sc_wt_ns_meta.pdf\", format=\"pdf\")\n",
        "    print(\"sc_wt_ns_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k_hcat8uCoZd",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title S. cerevisiae KO\n",
        "ko = \"spt6-YW\" #@param [\"chd1\", \"ino80\", \"isw1\", \"isw2\", \"rsc\", \"spt6-1004\", \"spt6-YW\"]\n",
        "plot_sc_ko_ns_meta = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_sc_ko_ns_meta = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "_ko = \"sc_ko_\" + ko\n",
        "_sc_ko_ns_meta = {\n",
        "  \"color\": \"#2c3b96\",\n",
        "  \"linewidth\": 2.5,\n",
        "  \"ylim\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[_ko][\"antisense\"], 0))),\n",
        "    np.ceil(np.max(np.sum(ns_meta[_ko][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"yticks\": [\n",
        "    -np.ceil(np.max(np.sum(ns_meta[_ko][\"antisense\"], 0))),\n",
        "    0,\n",
        "    np.ceil(np.max(np.sum(ns_meta[_ko][\"sense\"], 0)))\n",
        "  ],\n",
        "  \"sense\": np.sum(ns_meta[_ko][\"sense\"], 0),\n",
        "  \"antisense\": np.sum(ns_meta[_ko][\"antisense\"], 0)\n",
        "}\n",
        "# plot\n",
        "if plot_sc_ko_ns_meta:\n",
        "  F, A = plt.subplots(1, 1, figsize=(11, 7))\n",
        "  # sense signal\n",
        "  A.plot(\n",
        "    y=_sc_ko_ns_meta[\"sense\"],\n",
        "    color=_sc_ko_ns_meta[\"color\"],\n",
        "    linewidth=_sc_ko_ns_meta[\"linewidth\"]\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=_sc_ko_ns_meta[\"sense\"],\n",
        "    color=_sc_ko_ns_meta[\"color\"],\n",
        "    linewidth=_sc_ko_ns_meta[\"linewidth\"],\n",
        "    alpha=0.6\n",
        "  )\n",
        "  # antisense signal\n",
        "  A.plot(\n",
        "    y=-_sc_ko_ns_meta[\"antisense\"],\n",
        "    color=_sc_ko_ns_meta[\"color\"],\n",
        "    linewidth=_sc_ko_ns_meta[\"linewidth\"],\n",
        "    alpha=0.8\n",
        "  )\n",
        "  A.fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=-_sc_ko_ns_meta[\"antisense\"],\n",
        "    y2=np.zeros(1024),\n",
        "    color=_sc_ko_ns_meta[\"color\"],\n",
        "    linewidth=_sc_ko_ns_meta[\"linewidth\"],\n",
        "    alpha=0.3\n",
        "  )\n",
        "  # ancillary\n",
        "  A.set_xlim([0, 1024])\n",
        "  A.set_xticks([])\n",
        "  A.set_ylim(_sc_ko_ns_meta[\"ylim\"])\n",
        "  A.set_yticks(_sc_ko_ns_meta[\"yticks\"])\n",
        "  [A.spines[s].set_color(None) for s in [\"left\", \"top\", \"right\"]]\n",
        "  plt.tight_layout()\n",
        "  # save\n",
        "  if save_sc_ko_ns_meta:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + _ko + \"_ns_meta.pdf\", format=\"pdf\")\n",
        "    print(_ko + \"_ns_meta.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e_v4KYJIBox6",
        "colab_type": "text"
      },
      "source": [
        "#### heatmap meta"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "80eVJcrpOMKV",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title all species / cell states\n",
        "plot_ns_heatmap = True #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "save_ns_heatmap = False #@param [\"False\", \"True\"] {type:\"raw\"}\n",
        "# details\n",
        "_ns_heatmap = {\n",
        "  species_state: np.arcsinh(ns_meta[species_state][\"sense\"])[\n",
        "    np.argsort(np.sum(ns_meta[species_state][\"sense\"], 1))[::-1]]\n",
        "      for species_state in ns_meta.keys()\n",
        "}\n",
        "if plot_ns_heatmap:\n",
        "  F, A = plt.subplots(1, len(_ns_heatmap), figsize=(len(_ns_heatmap) * 2, 8))\n",
        "  for i, (n, n_data) in enumerate(_ns_heatmap.items()):\n",
        "    print(i, n,\n",
        "      \"\\n   arcsinh max:\", np.around(np.max(n_data), 1),\n",
        "      \"\\n   normal max:\", np.around(np.max(ns_meta[n][\"sense\"]), 1)\n",
        "    )\n",
        "    sns.heatmap(n_data, 0, np.max(n_data), \"twilight\", cbar=False, ax=A[i])\n",
        "    A[i].set_xticks([]); A[i].set_yticks([])\n",
        "    A[i].set_title(n, fontdict={\"fontsize\": 4})\n",
        "  plt.tight_layout(w_pad=-0.15)\n",
        "  if save_ns_heatmap:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"ns_heatmap.png\", format=\"png\", dpi=1024)\n",
        "    print(\"ns_heatmap.png - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2fcBZkAHlm0-",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title all species / cell states - cbar\n",
        "if plot_ns_heatmap:\n",
        "  F, A = plt.subplots(1, 1)\n",
        "  sns.heatmap(np.expand_dims(n_data[0], 1), 0, np.max(n_data), \"twilight\", ax=A)\n",
        "  if save_ns_heatmap:\n",
        "    plt.savefig(_dirs[\"figs_stats\"] + \"ns_heatmap_cbar.pdf\", format=\"pdf\")\n",
        "    print(\"ns_heatmap_cbar.pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ki9lBEBjbA3C",
        "colab_type": "text"
      },
      "source": [
        "## plot individual considerations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zM5Vq0EzNO4I",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# plotting style\n",
        "sns.set_style(\"ticks\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cXmPQna5YJcf",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title extract index to gene ID\n",
        "# index : gene ID\n",
        "  # S. cerevisiae Wildtype\n",
        "sc_wt_ex_idx_ID = dict(zip(\n",
        "  np.arange(len(np.unique(meas_dfs[\"sc\"][\"ID\"]))),\n",
        "  np.unique(meas_dfs[\"sc\"][\"ID\"])\n",
        "))\n",
        "  # S. cerevisiae Δ CHD1\n",
        "sc_ko_chd1_ex_idx_ID = dict(zip(\n",
        "  np.arange(len(np.unique(meas_dfs[\"sc\"][\"ID\"]))),\n",
        "  np.unique(meas_dfs[\"sc\"][\"ID\"])\n",
        "))\n",
        "  # K. lactis\n",
        "kl_ex_idx_ID = dict(zip(\n",
        "  np.arange(len(np.unique(meas_dfs[\"kl\"][\"ID\"]))),\n",
        "  np.unique(meas_dfs[\"kl\"][\"ID\"])\n",
        "))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V2vvWamoCrIf",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title Sc WT\n",
        "idx += 1\n",
        "sc_wt_ex_ID = sc_wt_ex_idx_ID[idx]\n",
        "# get data\n",
        "  # raw dataframes\n",
        "sc_wt_ex_dfs = {\n",
        "  \"meas\": meas_dfs[\"sc\"][\n",
        "    meas_dfs[\"sc\"][\"ID\"] == sc_wt_ex_ID\n",
        "  ],\n",
        "  \"pred\": {\n",
        "    x: pred_dfs[y][pred_dfs[y][\"ID\"] == sc_wt_ex_ID]\n",
        "      for x, y in {\n",
        "        \"ctrl_ctrl\": \"sc_DS_shuffled_NS_shuffled_output_df\",\n",
        "        \"NET\": \"sc_DS_none_NS_native_output_df\",\n",
        "        \"DNA\": \"sc_DS_native_NS_none_output_df\",\n",
        "        \"integrated\": \"sc_DS_native_NS_native_output_df\"\n",
        "      }.items()\n",
        "  }\n",
        "}\n",
        "  # raw MNase-seq data\n",
        "sc_wt_ex_raw = {\n",
        "  \"meas\": {\n",
        "    \"MS_wt_1\": sc_wt_ex_dfs[\"meas\"][\"MS_wt_1\"].values,\n",
        "    \"MS_wt_2\": sc_wt_ex_dfs[\"meas\"][\"MS_wt_2\"].values\n",
        "  },\n",
        "  \"pred\": {\n",
        "    x: {\n",
        "      \"avg\": np.sum(sc_wt_ex_dfs[\"pred\"][x][[\n",
        "        col for col in sc_wt_ex_dfs[\"pred\"][x].columns\n",
        "          if (\"N_wt_DS_\" in col) and (\"pred_mean\" in col)\n",
        "      ]].values, 1),\n",
        "      \"var\": np.sum(sc_wt_ex_dfs[\"pred\"][x][[\n",
        "        col for col in sc_wt_ex_dfs[\"pred\"][x].columns\n",
        "          if (\"N_wt_DS_\" in col) and (\"pred_variance\" in col)\n",
        "      ]].values, 1)\n",
        "    } for x in sc_wt_ex_dfs[\"pred\"].keys()\n",
        "  }\n",
        "}\n",
        "  # smoothed MNase-seq data\n",
        "sc_wt_ex_smooth = {\n",
        "  \"meas\": {\n",
        "    \"MS_wt_1\": _np_ms_transform(sc_wt_ex_raw[\"meas\"][\"MS_wt_1\"]),\n",
        "    \"MS_wt_2\": _np_ms_transform(sc_wt_ex_raw[\"meas\"][\"MS_wt_2\"])\n",
        "  },\n",
        "  \"pred\": {\n",
        "    x: {\n",
        "      \"avg\": _np_ms_transform(sc_wt_ex_raw[\"pred\"][x][\"avg\"]),\n",
        "      \"var\": _np_ms_transform(sc_wt_ex_raw[\"pred\"][x][\"var\"])\n",
        "    } for x in sc_wt_ex_raw[\"pred\"].keys()\n",
        "  }\n",
        "}\n",
        "# plot details\n",
        "_ex_details = {\n",
        "  \"clr_pred\": {\n",
        "    \"ctrl_ctrl\": \"#e5e6e6\",\n",
        "    \"NET\": \"#c3cfef\",\n",
        "    \"DNA\": \"#f4bdbf\",\n",
        "    \"integrated\": \"#bb478d\"\n",
        "  },\n",
        "  \"clr_meas\": {\n",
        "    \"MS_wt_1\": \"#3cac49\",\n",
        "    \"MS_wt_2\": \"#b3ffb3\"\n",
        "  },\n",
        "  \"lw_min\": 1.0,\n",
        "  \"lw_max\": 2.3,\n",
        "  \"a_min\": 0.5,\n",
        "  \"a_max\": 0.9,\n",
        "  \"_fname_suff\": \"sc_wt_\" + sc_wt_ex_ID\n",
        "}\n",
        "# plot\n",
        "print(\"plotting \" + sc_wt_ex_ID)\n",
        "  # set-up\n",
        "F, A = plt.subplots(4, 1, sharex=True, figsize=(11, 10))\n",
        "# measured data\n",
        "for ii, (x, clr) in enumerate(_ex_details[\"clr_meas\"].items()):\n",
        "    # raw\n",
        "  A[ii].plot(\n",
        "    sc_wt_ex_raw[\"meas\"][x],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_min\"],\n",
        "    color=_ex_details[\"clr_meas\"][x]\n",
        "  )\n",
        "  A[ii].fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=sc_wt_ex_raw[\"meas\"][x],\n",
        "    alpha=_ex_details[\"a_min\"],\n",
        "    lw=_ex_details[\"lw_min\"],\n",
        "    color=_ex_details[\"clr_meas\"][x]\n",
        "  )\n",
        "    # smoothed\n",
        "  A[3].plot(\n",
        "    sc_wt_ex_smooth[\"meas\"][x],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_max\"],\n",
        "    color=_ex_details[\"clr_meas\"][x]\n",
        "  )\n",
        "# predicted data\n",
        "for x, clr in _ex_details[\"clr_pred\"].items():\n",
        "    # raw\n",
        "  A[2].plot(\n",
        "    sc_wt_ex_raw[\"pred\"][x][\"avg\"],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_min\"],\n",
        "    color=clr\n",
        "  )\n",
        "    # smoothed\n",
        "  A[3].plot(\n",
        "    sc_wt_ex_smooth[\"pred\"][x][\"avg\"],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_max\"],\n",
        "    color=clr\n",
        "  )\n",
        "# plot ancillaries\n",
        "A[0].set_ylim([0, np.max(sc_wt_ex_raw[\"meas\"][\"MS_wt_1\"])])\n",
        "A[0].set_yticks([0, np.max(sc_wt_ex_raw[\"meas\"][\"MS_wt_1\"])])\n",
        "A[1].set_ylim([0, np.max(sc_wt_ex_raw[\"meas\"][\"MS_wt_2\"])])\n",
        "A[1].set_yticks([0, np.max(sc_wt_ex_raw[\"meas\"][\"MS_wt_2\"])])\n",
        "A[2].set_yticks([])\n",
        "A[3].set_yticks([])\n",
        "for i in range(4):\n",
        "  A[i].set_xlim([-2, 1026])\n",
        "  A[i].set_xticks([])\n",
        "  for s in [\"bottom\", \"left\", \"top\", \"right\"]:\n",
        "    A[i].spines[s].set_linewidth(1)\n",
        "    A[i].spines[s].set_color(\"k\")\n",
        "plt.tight_layout(pad=0.2, h_pad=-0.2)\n",
        "plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ts7C33ytYWBc",
        "colab_type": "code",
        "cellView": "form",
        "colab": {}
      },
      "source": [
        "#@title Kl\n",
        "%%capture\n",
        "for idx, kl_ex_ID in kl_ex_idx_ID.items():\n",
        "  # get data\n",
        "    # raw dataframes\n",
        "  kl_ex_dfs = {\n",
        "    \"meas\": meas_dfs[\"kl\"][\n",
        "      meas_dfs[\"kl\"][\"ID\"] == kl_ex_ID\n",
        "    ],\n",
        "    \"pred\": {\n",
        "      x: pred_dfs[y][pred_dfs[y][\"ID\"] == kl_ex_ID]\n",
        "        for x, y in {\n",
        "          \"ctrl_ctrl\": \"kl_DS_shuffled_NS_shuffled_output_df\",\n",
        "          \"NET\": \"kl_DS_none_NS_native_output_df\",\n",
        "          \"DNA\": \"kl_DS_native_NS_none_output_df\",\n",
        "          \"integrated\": \"kl_DS_native_NS_native_output_df\"\n",
        "        }.items()\n",
        "    }\n",
        "  }\n",
        "    # raw MNase-seq data\n",
        "  kl_ex_raw = {\n",
        "    \"meas\": kl_ex_dfs[\"meas\"][\"MS_1\"].values,\n",
        "    \"pred\": {\n",
        "      x: {\n",
        "        \"avg\": np.sum(kl_ex_dfs[\"pred\"][x][[\n",
        "          col for col in kl_ex_dfs[\"pred\"][x].columns\n",
        "            if \"pred_mean\" in col\n",
        "        ]].values, 1),\n",
        "        \"var\": np.sum(kl_ex_dfs[\"pred\"][x][[\n",
        "          col for col in kl_ex_dfs[\"pred\"][x].columns\n",
        "            if \"pred_variance\" in col\n",
        "        ]].values, 1)\n",
        "      } for x in kl_ex_dfs[\"pred\"].keys()\n",
        "    }\n",
        "  }\n",
        "    # smoothed MNase-seq data\n",
        "  kl_ex_smooth = {\n",
        "    \"meas\": _np_ms_transform(kl_ex_raw[\"meas\"]),\n",
        "    \"pred\": {\n",
        "      x: {\n",
        "        \"avg\": _np_ms_transform(kl_ex_raw[\"pred\"][x][\"avg\"]),\n",
        "        \"var\": _np_ms_transform(kl_ex_raw[\"pred\"][x][\"var\"])\n",
        "      } for x in kl_ex_raw[\"pred\"].keys()\n",
        "    }\n",
        "  }\n",
        "  # plot details\n",
        "  _ex_details = {\n",
        "    \"clr_pred\": {\n",
        "      \"ctrl_ctrl\": \"#d1d3d4\",\n",
        "      \"NET\": \"#284492\",\n",
        "      \"DNA\": \"#cb2128\",\n",
        "      \"integrated\": \"#bb478d\"\n",
        "    },\n",
        "    \"clr_meas\": \"#3cac49\",\n",
        "    \"lw_min\": 1.4,\n",
        "    \"lw_max\": 2.6,\n",
        "    \"a_min\": 0.4,\n",
        "    \"a_max\": 0.7,\n",
        "    \"_fname_suff\": \"kl_\" + str(kl_ex_ID)\n",
        "  }\n",
        "  # plot\n",
        "  print(\"plotting \" + str(kl_ex_ID))\n",
        "  # set-up\n",
        "  F, A = plt.subplots(3, 1, sharex=True, figsize=(11, 11))\n",
        "  # measured data\n",
        "    # raw\n",
        "  A[0].plot(\n",
        "    kl_ex_raw[\"meas\"],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_min\"],\n",
        "    color=_ex_details[\"clr_meas\"]\n",
        "  )\n",
        "  A[0].fill_between(\n",
        "    x=np.arange(1024),\n",
        "    y1=kl_ex_raw[\"meas\"],\n",
        "    alpha=_ex_details[\"a_min\"],\n",
        "    lw=_ex_details[\"lw_min\"],\n",
        "    color=_ex_details[\"clr_meas\"]\n",
        "  )\n",
        "    # smoothed\n",
        "  A[2].plot(\n",
        "    kl_ex_smooth[\"meas\"],\n",
        "    alpha=_ex_details[\"a_max\"],\n",
        "    lw=_ex_details[\"lw_max\"],\n",
        "    color=_ex_details[\"clr_meas\"]\n",
        "  )\n",
        "  # predicted data\n",
        "  for x, clr in _ex_details[\"clr_pred\"].items():\n",
        "      # raw\n",
        "    A[1].plot(\n",
        "      kl_ex_raw[\"pred\"][x][\"avg\"],\n",
        "      alpha=_ex_details[\"a_max\"],\n",
        "      lw=_ex_details[\"lw_max\"],\n",
        "      color=clr\n",
        "    )\n",
        "      # smoothed\n",
        "    A[2].plot(\n",
        "      kl_ex_smooth[\"pred\"][x][\"avg\"],\n",
        "      alpha=_ex_details[\"a_max\"],\n",
        "      lw=_ex_details[\"lw_max\"],\n",
        "      color=clr\n",
        "    )\n",
        "  # plot ancillaries\n",
        "  A[0].set_ylim([0, np.max(kl_ex_raw[\"meas\"])])\n",
        "  A[1].set_yticks([])\n",
        "  A[2].set_yticks([])\n",
        "  for i in range(3):\n",
        "    A[i].set_xlim([-2, 1026])\n",
        "    A[i].set_xticks([])\n",
        "    for s in [\"bottom\", \"left\", \"top\", \"right\"]:\n",
        "      A[i].spines[s].set_linewidth(1)\n",
        "      A[i].spines[s].set_color(\"k\")\n",
        "  plt.tight_layout(pad=0.2, h_pad=-0.2)\n",
        "  # save\n",
        "  plt.savefig(\n",
        "    fname=\"/\".join([\n",
        "      _dirs[\"figs_stats\"], \"example_loci\", _ex_details[\"_fname_suff\"] + \".pdf\"\n",
        "    ]),\n",
        "    bbox_inches=\"tight\",\n",
        "    format=\"pdf\"\n",
        "  )\n",
        "  print(_ex_details[\"_fname_suff\"] + \".pdf - saved\")\n",
        "  plt.show();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7y9xKGTflGui",
        "colab_type": "text"
      },
      "source": [
        "# scratch"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gJiElRm14o3z",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4vSGLPNmleBT",
        "colab_type": "text"
      },
      "source": [
        "# backend"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1uaf2Jq5nbpE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%reset -f"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8bc7ZLxPv4Wi",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%whos"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tgJxdJBclx3K",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!cat /proc/meminfo"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0UCMzqutpMii",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip freeze"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}