sigE_discovery_peak_summits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SigE\n",
    "Analysis of chip-seq for SigE in EtOH stres vs nostres conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import re\n",
    "from Bio import SeqIO\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from IPython.display import IFrame, display, Image\n",
    "import shutil\n",
    "import subprocess\n",
    "import functools\n",
    "import operator\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import time\n",
    "import seaborn as sns\n",
    "import sys\n",
    "from Bio import SeqUtils\n",
    "from itertools import islice      \n",
    "\n",
    "sns.set()\n",
    "\n",
    "def chunker(seq, l):\n",
    "    for i in range(0, len(seq), l):\n",
    "        yield seq[i:i+l]\n",
    "        \n",
    "def window(seq, n=2):\n",
    "    \"Returns a sliding window (of width n) over data from the iterable\"\n",
    "    \"   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   \"\n",
    "    it = iter(seq)\n",
    "    result = tuple(islice(it, n))\n",
    "    if len(result) == n:\n",
    "        yield result\n",
    "    for elem in it:\n",
    "        result = result[1:] + (elem,)\n",
    "        yield result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = 'peak_summit_motif_discovery'\n",
    "if os.path.exists(temp):\n",
    "    shutil.rmtree(temp)\n",
    "os.makedirs(temp, exist_ok=True)\n",
    "os.environ['JUPYTER_TEMP'] = temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Streptomyces genome accession AL645882.2\n",
    "SC_ACCESSION = \"AL645882.2\"\n",
    "\n",
    "sc_genome = [\n",
    "    s for s in SeqIO.parse('data/GCA_000203835.1_ASM20383v1_genomic.fna', format='fasta') if s.id == SC_ACCESSION\n",
    "][0]\n",
    "\n",
    "# write fasta file\n",
    "SeqIO.write([sc_genome], os.path.join(temp, 'AL645882.2.fasta'), format='fasta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[Kcessed: 100.0%"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "fasta-get-markov -m 1 $JUPYTER_TEMP/AL645882.2.fasta > $JUPYTER_TEMP/AL645882.2_m1.markov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare sequences\n",
    "---\n",
    "Extract sigE peak sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the SigE peaks\n",
    "df_e = pd.read_excel('data/sige_macs_peaks_curated.xlsx', sheet_name='sigE_peaks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>length</th>\n",
       "      <th>abs_summit</th>\n",
       "      <th>pileup</th>\n",
       "      <th>neglog10pvalue</th>\n",
       "      <th>fold_enrichment</th>\n",
       "      <th>neglog10qvalue</th>\n",
       "      <th>norm/EtOH</th>\n",
       "      <th>sequence_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8662666</td>\n",
       "      <td>8663025</td>\n",
       "      <td>360</td>\n",
       "      <td>8662837</td>\n",
       "      <td>4159.26</td>\n",
       "      <td>399.36166</td>\n",
       "      <td>2.10648</td>\n",
       "      <td>397.12149</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8483428</td>\n",
       "      <td>8484106</td>\n",
       "      <td>679</td>\n",
       "      <td>8483707</td>\n",
       "      <td>9107.42</td>\n",
       "      <td>2952.46606</td>\n",
       "      <td>4.61191</td>\n",
       "      <td>2948.32202</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8038581</td>\n",
       "      <td>8038938</td>\n",
       "      <td>358</td>\n",
       "      <td>8038755</td>\n",
       "      <td>6019.20</td>\n",
       "      <td>1160.16821</td>\n",
       "      <td>3.04824</td>\n",
       "      <td>1157.51782</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>8005927</td>\n",
       "      <td>8006952</td>\n",
       "      <td>1026</td>\n",
       "      <td>8006657</td>\n",
       "      <td>9276.34</td>\n",
       "      <td>1427.76892</td>\n",
       "      <td>2.65750</td>\n",
       "      <td>1425.00598</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>7885747</td>\n",
       "      <td>7886423</td>\n",
       "      <td>677</td>\n",
       "      <td>7886038</td>\n",
       "      <td>8795.63</td>\n",
       "      <td>2617.73657</td>\n",
       "      <td>4.26193</td>\n",
       "      <td>2614.10645</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>188</th>\n",
       "      <td>6107752</td>\n",
       "      <td>6107920</td>\n",
       "      <td>169</td>\n",
       "      <td>6107856</td>\n",
       "      <td>7593.41</td>\n",
       "      <td>1068.68542</td>\n",
       "      <td>2.52894</td>\n",
       "      <td>1066.07446</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_188</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>189</th>\n",
       "      <td>5999380</td>\n",
       "      <td>5999532</td>\n",
       "      <td>153</td>\n",
       "      <td>5999460</td>\n",
       "      <td>6709.26</td>\n",
       "      <td>833.43805</td>\n",
       "      <td>2.37112</td>\n",
       "      <td>830.94513</td>\n",
       "      <td>x</td>\n",
       "      <td>sigE_189</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>190</th>\n",
       "      <td>4336288</td>\n",
       "      <td>4336472</td>\n",
       "      <td>185</td>\n",
       "      <td>4336360</td>\n",
       "      <td>7307.68</td>\n",
       "      <td>841.24341</td>\n",
       "      <td>2.28611</td>\n",
       "      <td>838.74652</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>191</th>\n",
       "      <td>1869170</td>\n",
       "      <td>1869322</td>\n",
       "      <td>153</td>\n",
       "      <td>1869226</td>\n",
       "      <td>9659.11</td>\n",
       "      <td>3326.14575</td>\n",
       "      <td>4.89125</td>\n",
       "      <td>3321.38403</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>192</th>\n",
       "      <td>3958080</td>\n",
       "      <td>3958352</td>\n",
       "      <td>273</td>\n",
       "      <td>3958160</td>\n",
       "      <td>10132.63</td>\n",
       "      <td>2805.60889</td>\n",
       "      <td>4.00254</td>\n",
       "      <td>2801.70435</td>\n",
       "      <td>1</td>\n",
       "      <td>sigE_192</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>193 rows × 10 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       start      end  length  abs_summit    pileup  neglog10pvalue  \\\n",
       "0    8662666  8663025     360     8662837   4159.26       399.36166   \n",
       "1    8483428  8484106     679     8483707   9107.42      2952.46606   \n",
       "2    8038581  8038938     358     8038755   6019.20      1160.16821   \n",
       "3    8005927  8006952    1026     8006657   9276.34      1427.76892   \n",
       "4    7885747  7886423     677     7886038   8795.63      2617.73657   \n",
       "..       ...      ...     ...         ...       ...             ...   \n",
       "188  6107752  6107920     169     6107856   7593.41      1068.68542   \n",
       "189  5999380  5999532     153     5999460   6709.26       833.43805   \n",
       "190  4336288  4336472     185     4336360   7307.68       841.24341   \n",
       "191  1869170  1869322     153     1869226   9659.11      3326.14575   \n",
       "192  3958080  3958352     273     3958160  10132.63      2805.60889   \n",
       "\n",
       "     fold_enrichment  neglog10qvalue norm/EtOH sequence_name  \n",
       "0            2.10648       397.12149         1        sigE_0  \n",
       "1            4.61191      2948.32202         1        sigE_1  \n",
       "2            3.04824      1157.51782         1        sigE_2  \n",
       "3            2.65750      1425.00598         1        sigE_3  \n",
       "4            4.26193      2614.10645         1        sigE_4  \n",
       "..               ...             ...       ...           ...  \n",
       "188          2.52894      1066.07446         1      sigE_188  \n",
       "189          2.37112       830.94513         x      sigE_189  \n",
       "190          2.28611       838.74652         1      sigE_190  \n",
       "191          4.89125      3321.38403         1      sigE_191  \n",
       "192          4.00254      2801.70435         1      sigE_192  \n",
       "\n",
       "[193 rows x 10 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = pd.read_table(\n",
    "    'data/GCA_000203835.1_ASM20383v1_genomic.gff', comment=\"#\", names=['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes']\n",
    ")\n",
    "\n",
    "features = features.loc[features.seqid == SC_ACCESSION]\n",
    "\n",
    "def gff_attributes_mapper(gff_attr):\n",
    "    d = dict()\n",
    "    for attr in gff_attr.split(';'):\n",
    "        k,v = attr.split('=')\n",
    "        d[k] = v\n",
    "    return d\n",
    "\n",
    "features = pd.concat([features, pd.DataFrame(map(lambda x:x, features.attributes.apply(gff_attributes_mapper)))], axis=1)\n",
    "features = features.drop(columns=['attributes'])\n",
    "\n",
    "name_tag = 'gene'\n",
    "\n",
    "f_fw = {k: v for k, v in features.loc[(features.strand == '+') & (~features[name_tag].isna()) & (features.type.isin({'CDS', 'gene'})), [name_tag, 'start']].values}\n",
    "f_rc = {k: v for k, v in features.loc[(features.strand == '-') & (~features[name_tag].isna()) & (features.type.isin({'CDS', 'gene'})), [name_tag, 'end']].values}\n",
    "\n",
    "# copy data set - retrieve data from source, set on new\n",
    "df_sigE = df_e.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run motif search\n",
    "---\n",
    "\n",
    "- using FIMO from MEME suite\n",
    "- we use backgroud - first-order markov model from whole genome (AL645882.2)\n",
    "- the individual hits will be processed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write peak sequences\n",
    "\n",
    "# use peak summit\n",
    "summ_delta = 50\n",
    "peak_seqs = []\n",
    "for id, summ in df_sigE.loc[:, ['sequence_name', 'abs_summit']].values:\n",
    "    s = summ-summ_delta\n",
    "    e = summ+summ_delta\n",
    "    peak_seqs.append(\n",
    "        '>{}-{}-{}\\n{}\\n'.format(id, s, e, str(sc_genome.seq[s:e]))\n",
    "    )\n",
    "    \n",
    "with open(os.path.join(temp, 'sigE_peaks_50.fasta'), 'w') as o:\n",
    "    o.write('\\n'.join(peak_seqs))\n",
    "    \n",
    "    \n",
    "peak_seqs2 = []\n",
    "for id, s, e in df_sigE.loc[:, ['sequence_name', 'start', 'end']].values:\n",
    "    peak_seqs2.append(\n",
    "        '>{}-{}-{}\\n{}\\n'.format(id, s, e-s+1, str(sc_genome.seq[s:e]))\n",
    "   )\n",
    "with open(os.path.join(temp, 'sigE_peaks_whole.fasta'), 'w') as o:\n",
    "    o.write('\\n'.join(peak_seqs2))\n",
    "    \n",
    "    \n",
    "    \n",
    "# use peak summit\n",
    "summ_delta = 100\n",
    "peak_seqs_100 = []\n",
    "for id, summ in df_sigE.loc[:, ['sequence_name', 'abs_summit']].values:\n",
    "    s = summ-summ_delta\n",
    "    e = summ+summ_delta\n",
    "    peak_seqs_100.append(\n",
    "        '>{}-{}-{}\\n{}\\n'.format(id, s, e, str(sc_genome.seq[s:e]))\n",
    "    )\n",
    "with open(os.path.join(temp, 'sigE_peaks_100.fasta'), 'w') as o:\n",
    "    o.write('\\n'.join(peak_seqs_100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MEME\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MEME zoops with background, delta 50 nts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Writing results to output directory 'peak_summit_motif_discovery/meme_zoops_test'.\n",
      "BACKGROUND: using background model of order 1\n",
      "PRIMARY (classic): n 193 p0 193 p1 0 p2 0\n",
      "SEQUENCE GROUP USAGE-- Starts/EM: p0; Trim: p0; pvalue: p0; nsites: p0,p1,p2\n",
      "SEEDS: maxwords 19300 highwater mark: seq 193 pos 96\n",
      "BALANCE: samples 193 chars 19300 nodes 4 chars/node 4825\n",
      "Initializing the motif probability tables for 2 to 193 sites...\n",
      "nsites = 192\n",
      "Done initializing.\n",
      "all widths from min to max\n",
      "\n",
      "seqs=   193, min_w= 100, max_w=  100, total_size=    19300\n",
      "\n",
      "motif=1\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=2\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=3\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=4\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=5\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=6\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=7\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=8\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=9\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n",
      "motif=10\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 4 5 6 7\n",
      "em: w=   7, psites=  32, iter=  40 \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "meme $JUPYTER_TEMP/sigE_peaks_50.fasta -revcomp -bfile $JUPYTER_TEMP/AL645882.2_m1.markov -dna -mod zoops -minw 4 -maxw 7 -allw -oc $JUPYTER_TEMP/meme_zoops_test -nmotifs 10 -p 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MEME zoops, background, delta 100 nts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Writing results to output directory 'peak_summit_motif_discovery/meme_zoops_test_100'.\n",
      "BACKGROUND: using background model of order 1\n",
      "PRIMARY (classic): n 193 p0 193 p1 0 p2 0\n",
      "SEQUENCE GROUP USAGE-- Starts/EM: p0; Trim: p0; pvalue: p0; nsites: p0,p1,p2\n",
      "SEEDS: maxwords 38600 highwater mark: seq 193 pos 194\n",
      "BALANCE: samples 193 chars 38600 nodes 4 chars/node 9650\n",
      "Initializing the motif probability tables for 2 to 193 sites...\n",
      "nsites = 192\n",
      "Done initializing.\n",
      "all widths from min to max\n",
      "\n",
      "seqs=   193, min_w= 200, max_w=  200, total_size=    38600\n",
      "\n",
      "motif=1\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=2\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=3\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=4\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=5\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=6\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=7\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=8\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=9\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=10\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "meme $JUPYTER_TEMP/sigE_peaks_100.fasta -revcomp -bfile $JUPYTER_TEMP/AL645882.2_m1.markov -dna -mod zoops -minw 6 -maxw 10 -allw -oc $JUPYTER_TEMP/meme_zoops_test_100 -nmotifs 10 -p 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MEME zoops, background, delta 100 nts, central distance as objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Writing results to output directory 'peak_summit_motif_discovery/meme_zoops_test_100_cd'.\n",
      "CD: cefrac 0.25 length 200 central region 50\n",
      "BACKGROUND: using background model of order 1\n",
      "PRIMARY (cd): n 193 p0 96 p1 48 p2 49\n",
      "SEQUENCE GROUP USAGE-- Starts/EM: p0; Trim: p1; pvalue: p2; nsites: p0,p1,p2\n",
      "SEEDS: maxwords 38600 highwater mark: seq 193 pos 194\n",
      "BALANCE: samples 193 chars 38600 nodes 4 chars/node 9650\n",
      "all widths from min to max\n",
      "\n",
      "seqs=   193, min_w= 200, max_w=  200, total_size=    38600\n",
      "\n",
      "motif=1\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  30 \n",
      "motif=2\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=3\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=4\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=5\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=6\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=7\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=8\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=9\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=10\n",
      "SEED DEPTHS: 2 4 8 16 32 64 96\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "meme $JUPYTER_TEMP/sigE_peaks_100.fasta -objfun cd -revcomp -bfile $JUPYTER_TEMP/AL645882.2_m1.markov -dna -mod zoops -minw 6 -maxw 10 -allw -oc $JUPYTER_TEMP/meme_zoops_test_100_cd -nmotifs 10 -p 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MEME zoops, background, delta 50, motifsw 6-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Writing results to output directory 'peak_summit_motif_discovery/meme_zoops_test_6-10'.\n",
      "BACKGROUND: using background model of order 1\n",
      "PRIMARY (classic): n 193 p0 193 p1 0 p2 0\n",
      "SEQUENCE GROUP USAGE-- Starts/EM: p0; Trim: p0; pvalue: p0; nsites: p0,p1,p2\n",
      "SEEDS: maxwords 19300 highwater mark: seq 193 pos 94\n",
      "BALANCE: samples 193 chars 19300 nodes 4 chars/node 4825\n",
      "Initializing the motif probability tables for 2 to 193 sites...\n",
      "nsites = 192\n",
      "Done initializing.\n",
      "all widths from min to max\n",
      "\n",
      "seqs=   193, min_w= 100, max_w=  100, total_size=    19300\n",
      "\n",
      "motif=1\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=2\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=3\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=4\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=5\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=6\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=7\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=8\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=9\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n",
      "motif=10\n",
      "SEED DEPTHS: 2 4 8 16 32 64 128 193\n",
      "SEED WIDTHS: 6 7 8 9 10\n",
      "em: w=  10, psites=  32, iter=  40 \n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "meme $JUPYTER_TEMP/sigE_peaks_50.fasta -revcomp -bfile $JUPYTER_TEMP/AL645882.2_m1.markov -dna -mod zoops -minw 6 -maxw 10 -allw -oc $JUPYTER_TEMP/meme_zoops_test_6-10 -nmotifs 10 -p 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import meme_txt_parser\n",
    "selected_motif = meme_txt_parser.read_file(os.path.join(temp, 'meme_zoops_test', 'meme.txt'))\n",
    "\n",
    "# select motif MEME-1\n",
    "selected_motif = selected_motif.loc[selected_motif.motif_id == \"MEME-1\", :]\n",
    "seq_name = selected_motif.sequence_name.str.split('-')\n",
    "selected_motif['sequence_name'] = [i[0] for i in seq_name]\n",
    "selected_motif['sequence_offset'] = [int(i[1]) for i in seq_name]\n",
    "selected_motif['stop'] = selected_motif.start + selected_motif.site.str.len()\n",
    "selected_motif.drop(columns=['_site_left', '_site_right'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>motif_id</th>\n",
       "      <th>sequence_name</th>\n",
       "      <th>strand</th>\n",
       "      <th>start</th>\n",
       "      <th>P-value</th>\n",
       "      <th>site</th>\n",
       "      <th>sequence_offset</th>\n",
       "      <th>stop</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_181</td>\n",
       "      <td>+</td>\n",
       "      <td>16</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>2088569</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_171</td>\n",
       "      <td>+</td>\n",
       "      <td>53</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>3659392</td>\n",
       "      <td>60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_168</td>\n",
       "      <td>+</td>\n",
       "      <td>81</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>4074167</td>\n",
       "      <td>88</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_156</td>\n",
       "      <td>-</td>\n",
       "      <td>43</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>6883700</td>\n",
       "      <td>50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_154</td>\n",
       "      <td>+</td>\n",
       "      <td>51</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>7631695</td>\n",
       "      <td>58</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  motif_id sequence_name strand  start   P-value     site  sequence_offset  \\\n",
       "0   MEME-1      sigE_181      +     16  0.000046  GAAGACG          2088569   \n",
       "1   MEME-1      sigE_171      +     53  0.000046  GAAGACG          3659392   \n",
       "2   MEME-1      sigE_168      +     81  0.000046  GAAGACG          4074167   \n",
       "3   MEME-1      sigE_156      -     43  0.000046  GAAGACG          6883700   \n",
       "4   MEME-1      sigE_154      +     51  0.000046  GAAGACG          7631695   \n",
       "\n",
       "   stop  \n",
       "0    23  \n",
       "1    60  \n",
       "2    88  \n",
       "3    50  \n",
       "4    58  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected_motif.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find direct downstream genes\n",
    "\n",
    "features = pd.read_table(\n",
    "    'data/GCA_000203835.1_ASM20383v1_genomic.gff', comment=\"#\", names=['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes']\n",
    ")\n",
    "\n",
    "features = features.loc[features.seqid == SC_ACCESSION]\n",
    "\n",
    "def gff_attributes_mapper(gff_attr):\n",
    "    d = dict()\n",
    "    for attr in gff_attr.split(';'):\n",
    "        k,v = attr.split('=')\n",
    "        d[k] = v\n",
    "    return d\n",
    "\n",
    "features = pd.concat([features, pd.DataFrame(map(lambda x:x, features.attributes.apply(gff_attributes_mapper)))], axis=1)\n",
    "features = features.drop(columns=['attributes'])\n",
    "\n",
    "features_fw = features.loc[features.type.isin(['gene', 'CDS']) & (features.strand == '+'), :]\n",
    "features_rc = features.loc[features.type.isin(['gene', 'CDS']) & (features.strand == '-'), :]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "downstream_features = []\n",
    "downstream_feature_distances = []\n",
    "for r in selected_motif.itertuples():\n",
    "    if r.strand == '+':\n",
    "        feature_index = features_fw.start >= r.stop + r.sequence_offset\n",
    "        if feature_index.any():\n",
    "            downstream_feature = features_fw.loc[features_fw.loc[feature_index].start.idxmin()].to_dict()\n",
    "            distance = downstream_feature['start'] - (r.stop + r.sequence_offset)\n",
    "        else:\n",
    "            downstream_feature = dict()\n",
    "            distance = None\n",
    "    elif r.strand == '-':\n",
    "        feature_index = features_rc.end <= r.start + r.sequence_offset\n",
    "        if feature_index.any():\n",
    "            downstream_feature = features_rc.loc[features_rc.loc[feature_index].end.idxmax()].to_dict()\n",
    "            distance = (r.start + r.sequence_offset) - downstream_feature['end']\n",
    "        else:\n",
    "            downstream_feature = dict()\n",
    "            distance = None\n",
    "    else:\n",
    "        raise ValueError('Invalid strand')\n",
    "    downstream_features.append(downstream_feature)\n",
    "    downstream_feature_distances.append(distance)\n",
    "\n",
    "downstream_features = pd.DataFrame(downstream_features)\n",
    "downstream_features['distance'] = downstream_feature_distances\n",
    "downstream_features.dropna(axis=1, how='all', inplace=True)\n",
    "\n",
    "selected_motif_with_downstream = pd.concat([selected_motif, downstream_features.loc[:, ['gene', 'distance']]], axis=1)\n",
    "selected_motif_with_downstream.rename(columns={'gene': 'downstream_gene', 'distance': 'distance_to_downstream'}, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>motif_id</th>\n",
       "      <th>sequence_name</th>\n",
       "      <th>strand</th>\n",
       "      <th>start</th>\n",
       "      <th>P-value</th>\n",
       "      <th>site</th>\n",
       "      <th>sequence_offset</th>\n",
       "      <th>stop</th>\n",
       "      <th>downstream_gene</th>\n",
       "      <th>distance_to_downstream</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_181</td>\n",
       "      <td>+</td>\n",
       "      <td>16</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>2088569</td>\n",
       "      <td>23</td>\n",
       "      <td>SCO1954</td>\n",
       "      <td>1617.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_171</td>\n",
       "      <td>+</td>\n",
       "      <td>53</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>3659392</td>\n",
       "      <td>60</td>\n",
       "      <td>SCO3310</td>\n",
       "      <td>2285.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_168</td>\n",
       "      <td>+</td>\n",
       "      <td>81</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>4074167</td>\n",
       "      <td>88</td>\n",
       "      <td>SCO3696</td>\n",
       "      <td>1683.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_156</td>\n",
       "      <td>-</td>\n",
       "      <td>43</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>6883700</td>\n",
       "      <td>50</td>\n",
       "      <td>SCO6255</td>\n",
       "      <td>5568.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_154</td>\n",
       "      <td>+</td>\n",
       "      <td>51</td>\n",
       "      <td>0.000046</td>\n",
       "      <td>GAAGACG</td>\n",
       "      <td>7631695</td>\n",
       "      <td>58</td>\n",
       "      <td>SCO6885</td>\n",
       "      <td>22393.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_116</td>\n",
       "      <td>-</td>\n",
       "      <td>42</td>\n",
       "      <td>0.000339</td>\n",
       "      <td>ATAGACG</td>\n",
       "      <td>3064610</td>\n",
       "      <td>49</td>\n",
       "      <td>SCO2806</td>\n",
       "      <td>88.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>83</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_178</td>\n",
       "      <td>-</td>\n",
       "      <td>53</td>\n",
       "      <td>0.000394</td>\n",
       "      <td>GTAGATG</td>\n",
       "      <td>2645931</td>\n",
       "      <td>60</td>\n",
       "      <td>SCO2458</td>\n",
       "      <td>563.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>84</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_158</td>\n",
       "      <td>+</td>\n",
       "      <td>49</td>\n",
       "      <td>0.000394</td>\n",
       "      <td>AAATAAG</td>\n",
       "      <td>6509941</td>\n",
       "      <td>56</td>\n",
       "      <td>SCO5943</td>\n",
       "      <td>28.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>85</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_187</td>\n",
       "      <td>-</td>\n",
       "      <td>36</td>\n",
       "      <td>0.000404</td>\n",
       "      <td>ATAGAAG</td>\n",
       "      <td>490356</td>\n",
       "      <td>43</td>\n",
       "      <td>SCO0468</td>\n",
       "      <td>42.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>86</th>\n",
       "      <td>MEME-1</td>\n",
       "      <td>sigE_115</td>\n",
       "      <td>+</td>\n",
       "      <td>75</td>\n",
       "      <td>0.000696</td>\n",
       "      <td>ATATAAG</td>\n",
       "      <td>3095306</td>\n",
       "      <td>82</td>\n",
       "      <td>SCO2836</td>\n",
       "      <td>118.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>87 rows × 10 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   motif_id sequence_name strand  start   P-value     site  sequence_offset  \\\n",
       "0    MEME-1      sigE_181      +     16  0.000046  GAAGACG          2088569   \n",
       "1    MEME-1      sigE_171      +     53  0.000046  GAAGACG          3659392   \n",
       "2    MEME-1      sigE_168      +     81  0.000046  GAAGACG          4074167   \n",
       "3    MEME-1      sigE_156      -     43  0.000046  GAAGACG          6883700   \n",
       "4    MEME-1      sigE_154      +     51  0.000046  GAAGACG          7631695   \n",
       "..      ...           ...    ...    ...       ...      ...              ...   \n",
       "82   MEME-1      sigE_116      -     42  0.000339  ATAGACG          3064610   \n",
       "83   MEME-1      sigE_178      -     53  0.000394  GTAGATG          2645931   \n",
       "84   MEME-1      sigE_158      +     49  0.000394  AAATAAG          6509941   \n",
       "85   MEME-1      sigE_187      -     36  0.000404  ATAGAAG           490356   \n",
       "86   MEME-1      sigE_115      +     75  0.000696  ATATAAG          3095306   \n",
       "\n",
       "    stop downstream_gene  distance_to_downstream  \n",
       "0     23         SCO1954                  1617.0  \n",
       "1     60         SCO3310                  2285.0  \n",
       "2     88         SCO3696                  1683.0  \n",
       "3     50         SCO6255                  5568.0  \n",
       "4     58         SCO6885                 22393.0  \n",
       "..   ...             ...                     ...  \n",
       "82    49         SCO2806                    88.0  \n",
       "83    60         SCO2458                   563.0  \n",
       "84    56         SCO5943                    28.0  \n",
       "85    43         SCO0468                    42.0  \n",
       "86    82         SCO2836                   118.0  \n",
       "\n",
       "[87 rows x 10 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selected_motif_with_downstream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pd.ExcelWriter(os.path.join(temp, 'sigE_motif_discovery.xlsx'), engine='xlsxwriter') as writer:\n",
    "    desc = [\n",
    "        \"Motif sites for MEME-1 motif discovered when searching for conserved 4-7nt mers in 100bp region of sigE peaks centered at abs_summit\",\n",
    "        \"first order markov model of S.c. genome was used as background.\",\n",
    "        \"Columns:\",\n",
    "        \"  motif_id: meme assigned motif id\",\n",
    "        \"  sequence_name: peak sequence id\",\n",
    "        \"  strand: strand on which motif is found\",\n",
    "        \"  start: motif site start position\",\n",
    "        \"  P-value: motif site P-value (reported by meme)\",\n",
    "        \"  site: discovered motif site sequence 5'->3'\",\n",
    "        \"  sequence_offset: sequence start position in genome\",\n",
    "        \"  stop: motif site stop position\",\n",
    "        \"  downstream_gene: first gene downstream from motif stop\",\n",
    "        \"  distance_to_downstream: distance from motif site to downstream gene [nt]\"\n",
    "    ]\n",
    "    \n",
    "    sheetname = 'sigE-motif'\n",
    "    selected_motif_with_downstream.to_excel(writer, sheet_name=sheetname, startrow=len(desc) + 1, index=False)\n",
    "    for i, row in enumerate(desc):\n",
    "        writer.sheets[sheetname].write(i, 0, row)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.11"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}