spaced_kmer_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate binding scores based on a best-configuration spaced k-mer model\n",
    "### (Figure S6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import random\n",
    "import pandas as pd\n",
    "from Bio import SeqIO\n",
    "import itertools\n",
    "from sklearn.metrics import roc_auc_score, average_precision_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kmer_probs(sequences, k):\n",
    "    kmer_inx = generate_kmer_inx(l=k)\n",
    "    kmer_inx6 = generate_kmer_inx(l=6)\n",
    "    \n",
    "    kmer_counts = np.zeros(4**6)\n",
    "    \n",
    "    for seq in sequences:\n",
    "        kmers = get_kmers_sequence(seq, k=k)\n",
    "        for kmer in kmers:\n",
    "            kmer_counts[kmer_inx6[kmer[:3]+kmer[-3:]]] += 1\n",
    "            \n",
    "    #add pseudo count\n",
    "    kmer_counts = kmer_counts+1\n",
    "    \n",
    "    #normalize\n",
    "    kmer_probs = kmer_counts/np.sum(kmer_counts)\n",
    "    return kmer_probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kmer_ll(ps_seqs, bg_seqs, k):\n",
    "    \n",
    "    ps_probs = get_kmer_probs(ps_seqs, k=k)\n",
    "    bg_probs = get_kmer_probs(bg_seqs, k=k)\n",
    "    \n",
    "    return np.log(ps_probs)-np.log(bg_probs)    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_best_ll(sequence, ll, kmer_inx6, k):\n",
    "    kmers = get_kmers_sequence(sequence, k)\n",
    "    ll_values = [ll[kmer_inx6[kmer[:3]+kmer[-3:]]] for kmer in kmers]\n",
    "    return max(ll_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rank_sequences_by_ll(sequences, ll, k):\n",
    "    kmer_inx6 = generate_kmer_inx(l=6)\n",
    "    ll_by_sequence = [get_best_ll(sequence, ll, kmer_inx6, k) for sequence in sequences]\n",
    "    return ll_by_sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(6,12):\n",
    "    file_table = pd.read_csv('../../../cbscratch/selex_taipale/ps_and_bg_file_table_first_and_last.txt', sep='\\t', header=None)\n",
    "    fasta_path = '../../../cbscratch/selex_taipale/fasta_split_first_and_last/'\n",
    "    pred_dir = f'param/selex_nodist0_1000/selex_spaced_{k}mers'\n",
    "    !mkdir -p $pred_dir\n",
    "\n",
    "    for i in range(len(file_table)):\n",
    "        factor = file_table.iloc[i,0]\n",
    "\n",
    "        ps_file = os.path.join(fasta_path, f'{factor}_train_ps.fasta')\n",
    "        ps_seqs = parse_fasta(ps_file)\n",
    "\n",
    "        bg_file = os.path.join(fasta_path, f'{factor}_train_bg.fasta')\n",
    "        bg_seqs = parse_fasta(bg_file)\n",
    "\n",
    "        ps_file_test = os.path.join(fasta_path, f'{factor}_test_ps.fasta')\n",
    "        ps_seqs_test = parse_fasta(ps_file_test)\n",
    "\n",
    "        bg_file_test = os.path.join(fasta_path, f'{factor}_test_bg.fasta')\n",
    "        bg_seqs_test = parse_fasta(bg_file_test)\n",
    "\n",
    "        ll = get_kmer_ll(ps_seqs,bg_seqs,k=k)\n",
    "\n",
    "        ps_ll = rank_sequences_by_ll(ps_seqs_test, ll, k=k)\n",
    "        bg_ll = rank_sequences_by_ll(bg_seqs_test, ll, k=k)\n",
    "\n",
    "        y_scores = np.hstack([ps_ll, bg_ll])\n",
    "\n",
    "        output_prefix = os.path.join(pred_dir, factor)\n",
    "        np.savetxt(fname=f'{output_prefix}.predictions', X=y_scores)"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}