get_hobohm_PSSM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering sequence weighting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from time import time\n",
    "import math\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data_utils():\n",
    "    # Load the alphabet file\n",
    "    alphabet_file = data_dir + \"Matrices/alphabet\"\n",
    "\n",
    "    alphabet = np.loadtxt(alphabet_file, dtype=str)\n",
    "\n",
    "    # Load the matrices file \n",
    "    blosum62_file = data_dir + \"Matrices/blosum62.freq_rownorm\"\n",
    "    _blosum62 = np.loadtxt(blosum62_file, dtype=float).T\n",
    "\n",
    "    blosum62 = {}\n",
    "\n",
    "    for i, letter_1 in enumerate(alphabet):\n",
    "\n",
    "        blosum62[letter_1] = {}\n",
    "\n",
    "        for j, letter_2 in enumerate(alphabet):\n",
    "\n",
    "            blosum62[letter_1][letter_2] = _blosum62[i, j]\n",
    "            \n",
    "    # Load the bg information\n",
    "    bg_file = data_dir + \"Matrices/bg.freq.fmt\"\n",
    "    _bg = np.loadtxt(bg_file, dtype=float)\n",
    "\n",
    "    bg = {}\n",
    "    for i in range(0, len(alphabet)):\n",
    "        bg[alphabet[i]] = _bg[i]\n",
    "\n",
    "    return alphabet, blosum62, bg\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_matrix(peptide_len, alphabet):\n",
    "    init_matrix = [0]*peptide_len\n",
    "\n",
    "    for i in range(0, peptide_len):\n",
    "        row = {}\n",
    "        for letter in alphabet: \n",
    "            row[letter] = 0.0\n",
    "        init_matrix[i] = row\n",
    "        \n",
    "    return init_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def c_matrix(peptides, peptide_len, alphabet):\n",
    "    \n",
    "    cmatrix = initialize_matrix(peptide_len, alphabet)\n",
    "\n",
    "    for position in range(0, peptide_len):\n",
    "        for peptide in peptides:\n",
    "            cmatrix[position][peptide[position]] += 1\n",
    "    return cmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clustering_weighting(peptides, cmatrix):\n",
    "\n",
    "    weights = {}\n",
    "    for peptide in peptides:\n",
    "        # apply sequence weighting\n",
    "        weights[peptide] = 1\n",
    "        \n",
    "    return weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_matrix(peptides, peptide_len, alphabet):\n",
    "    \n",
    "    fmatrix = initialize_matrix(peptide_len, alphabet)\n",
    "\n",
    "    for position in range(0, peptide_len):\n",
    "        n = 0;\n",
    "\n",
    "        for peptide in peptides:\n",
    "            fmatrix[position][peptide[position]] += weights[peptide]\n",
    "            n += weights[peptide]\n",
    "\n",
    "        for letter in alphabet: \n",
    "            fmatrix[position][letter] = fmatrix[position][letter]/n\n",
    "\n",
    "    return fmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g_matrix(f_matrix, peptide_len, alphabet, blosum62):\n",
    "    \n",
    "    gmatrix = initialize_matrix(peptide_len, alphabet)\n",
    "\n",
    "    for position in range(0, peptide_len):\n",
    "        for letter_1 in alphabet:\n",
    "            for letter_2 in alphabet:\n",
    "              gmatrix[position][letter_1] += fmatrix[position][letter_2] * blosum62[letter_1][letter_2]\n",
    "\n",
    "    return gmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_matrix(f_matrix, g_matrix, peptide_len, alphabet, beta, neff):\n",
    "    \n",
    "    pmatrix = initialize_matrix(peptide_len, alphabet)\n",
    "    alpha = neff - 1\n",
    "\n",
    "    for position in range(0, peptide_len):\n",
    "        for a in alphabet:\n",
    "            pmatrix[position][a] = (alpha*fmatrix[position][a] + beta*gmatrix[position][a]) / (alpha + beta)\n",
    "\n",
    "    return pmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def w_matrix(p_matrix, bg, peptide_len, alphabet):\n",
    "    wmatrix = initialize_matrix(peptide_len, alphabet)\n",
    "\n",
    "    for position in range(0, peptide_len):\n",
    "\n",
    "        for letter in alphabet:\n",
    "            if pmatrix[position][letter] > 0:\n",
    "                wmatrix[position][letter] = 2 * math.log(pmatrix[position][letter]/bg[letter])/math.log(2)\n",
    "            else:\n",
    "                wmatrix[position][letter] = -999.9\n",
    "\n",
    "    return wmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_psi_blast(matrix):\n",
    "\n",
    "    header = [\"\", \"A\", \"R\", \"N\", \"D\", \"C\", \"Q\", \"E\", \"G\", \"H\", \"I\", \"L\", \"K\", \"M\", \"F\", \"P\", \"S\", \"T\", \"W\", \"Y\", \"V\"]\n",
    "    print ('{:>4} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8}'.format(*header)) \n",
    "    letter_order = [\"A\", \"R\", \"N\", \"D\", \"C\", \"Q\", \"E\", \"G\", \"H\", \"I\", \"L\", \"K\", \"M\", \"F\", \"P\", \"S\", \"T\", \"W\", \"Y\", \"V\"]\n",
    "\n",
    "    for i, row in enumerate(matrix):\n",
    "        scores = []\n",
    "        scores.append(str(i+1) + \" A\")\n",
    "\n",
    "        for letter in letter_order:\n",
    "            score = row[letter]\n",
    "            scores.append(round(score, 4))\n",
    "\n",
    "        print('{:>4} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8}'.format(*scores)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_psi_blast_file(matrix, file_name):\n",
    "    \n",
    "    with open(file_name, 'w') as file:\n",
    "        header = [\"\", \"A\", \"R\", \"N\", \"D\", \"C\", \"Q\", \"E\", \"G\", \"H\", \"I\", \"L\", \"K\", \"M\", \"F\", \"P\", \"S\", \"T\", \"W\", \"Y\", \"V\"]\n",
    "        file.write ('{:>4} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8}\\n'.format(*header)) \n",
    "        letter_order = [\"A\", \"R\", \"N\", \"D\", \"C\", \"Q\", \"E\", \"G\", \"H\", \"I\", \"L\", \"K\", \"M\", \"F\", \"P\", \"S\", \"T\", \"W\", \"Y\", \"V\"]\n",
    "\n",
    "        for i, row in enumerate(matrix):\n",
    "            scores = []\n",
    "            scores.append(str(i+1) + \" A\")\n",
    "\n",
    "            for letter in letter_order:\n",
    "                score = row[letter]\n",
    "                scores.append(round(score, 4))\n",
    "\n",
    "            file.write('{:>4} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8} {:>8}\\n'.format(*scores)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"/Users/laurasansc/github/algorithms_bioinformatics/data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphabet, blosum62, bg = load_data_utils()\n",
    "\n",
    "sequence_weighting = True\n",
    "beta = 50\n",
    "peptide_len = 9\n",
    "train_dir = data_dir + 'train/'\n",
    "\n",
    "# Iterate through data directory to find train files \n",
    "for peptides_file in os.listdir(train_dir):\n",
    "    if peptides_file.endswith(\"_train_cluster_centroids.csv\"): \n",
    "        allele = peptides_file.replace(\"_train_cluster_centroids.csv\",\"\")\n",
    "        # Read file \n",
    "        peptides_df = pd.read_csv(train_dir+peptides_file)\n",
    "        raw_peptides = peptides_df[\"sequence\"].to_list()\n",
    "\n",
    "        peptides = []\n",
    "        for i in range(0, len(raw_peptides)):\n",
    "            if len(raw_peptides[i]) == peptide_len:          \n",
    "                peptides.append(raw_peptides[i])\n",
    "            #else:\n",
    "             #   print(\"Peptide length too short discard\", raw_peptides[i])\n",
    "        neff = len(peptides)\n",
    "        ### get count matrices\n",
    "        cmatrix = c_matrix(peptides, peptide_len, alphabet)\n",
    "        ### weights \n",
    "        weights = clustering_weighting(peptides, cmatrix)\n",
    "        ### f, g, p and w matrices\n",
    "        fmatrix = f_matrix(peptides, peptide_len, alphabet)\n",
    "        gmatrix = g_matrix(fmatrix, peptide_len, alphabet, blosum62)\n",
    "        pmatrix = p_matrix(fmatrix, g_matrix, peptide_len, alphabet, beta, neff)\n",
    "        wmatrix = w_matrix(pmatrix, bg, peptide_len, alphabet)\n",
    "        \n",
    "        # Write out PSSM in Psi-Blast format to file\n",
    "        file_name = data_dir + \"hobohm_PSSM/w_matrix_\" + allele + '.csv'\n",
    "        to_psi_blast_file(wmatrix, file_name)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}