hidden_markov_package.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "## Install Libraries"
      ],
      "metadata": {
        "id": "xay-hjRKnLsO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 0: Student Information\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "- Student Name = Taraneh Khosrojerdi\n",
        "- Student Number = 400104929\n"
      ],
      "metadata": {
        "id": "K3LRemwqjrk1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install hidden_markov"
      ],
      "metadata": {
        "id": "G7TptWBCcmti",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4f627d0f-3c6c-4da9-b798-5373fc03b8c4"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting hidden_markov\n",
            "  Downloading hidden_markov-0.3.2.tar.gz (6.6 kB)\n",
            "  Preparing metadata (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from hidden_markov) (1.23.5)\n",
            "Building wheels for collected packages: hidden_markov\n",
            "  Building wheel for hidden_markov (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for hidden_markov: filename=hidden_markov-0.3.2-py3-none-any.whl size=6887 sha256=54e8d53bc6b7050b68e98c2133ea03fa83225871329259a941e8a632b172a1dc\n",
            "  Stored in directory: /root/.cache/pip/wheels/3c/56/8e/207019e308cefe1e24152109a944d0b8b7ef9f9863ddf5f47d\n",
            "Successfully built hidden_markov\n",
            "Installing collected packages: hidden_markov\n",
            "Successfully installed hidden_markov-0.3.2\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from hidden_markov import hmm"
      ],
      "metadata": {
        "id": "0E50GDimcmrt"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "In the human genome cytosine (C) is typically methylated. There is a relatively high change of mutation\n",
        "of this methyl-C into a thymine(T). As a result, in general CpG dinucleotides are rare in the genome than\n",
        "would be expected from the independent probabilities of C and G. [1]"
      ],
      "metadata": {
        "id": "tGA4_CNoeINa"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "<img src=\"https://i.ibb.co/9w4vG7t/Screen-Shot-2023-04-15-at-18-12-01.png\" alt=\"Screen-Shot-2023-04-15-at-18-12-01\" width=\"800\" border=\"0\">"
      ],
      "metadata": {
        "id": "0uPn2d8XefV7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 1: Creating HMM Model (5 Mark)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "- Define states and observations\n",
        "- Define the probabilities\n",
        "- Create the model"
      ],
      "metadata": {
        "id": "_AaFshr3ey0k"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "states = ('non_island', 'cpg_island')\n",
        "observations = ('A', 'C', 'T', 'G')\n",
        "\n",
        "start_probs = np.matrix('0.5 0.5')\n",
        "trans_probs = np.matrix('0.95 0.05 ; 0.1 0.9')\n",
        "emission_probs = np.matrix('0.27 0.24 0.26 0.23 ; 0.15 0.33 0.16 0.36')\n",
        "\n",
        "model = hmm(states, observations, start_probs, trans_probs, emission_probs)"
      ],
      "metadata": {
        "id": "K1wzHT1Qcmp_"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 2: Forward Algorithm (2 Mark)\n",
        "\n",
        "By using Forward algorithm, calculate the probability of all possible paths which give rise to the\n",
        "sequence **CGCG**."
      ],
      "metadata": {
        "id": "U6evvJ9XeCLW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Use forward algorithm from hidden_markov library\n",
        "obs_sequence = ('C', 'G', 'C', 'G')\n",
        "model.forward_algo(obs_sequence)"
      ],
      "metadata": {
        "id": "aff362LicmoL",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "d4e1be82-1883-41de-826c-a6e78270360e"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.007838160412500001"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 3: Viterbi Algorithm (2 Mark)\n",
        "\n",
        "Using the HMM model shown in Figure 2, predict the most probable path by Viterbi algorithm which\n",
        "give rise to the seqeunce CGCG."
      ],
      "metadata": {
        "id": "zReVSUQEfccb"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Use viterbi algorithm from hidden_markov library\n",
        "model.viterbi(obs_sequence)"
      ],
      "metadata": {
        "id": "eXvkAqdPihKZ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "d4e8da0e-f707-4253-cb39-acea2c87d2ca"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']"
            ]
          },
          "metadata": {},
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 4: Find K Best Possible Sequences (5 Mark)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "uB__n4Jgfhmo"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "- Complete find_k_obs function with following information below.\n",
        "- **Using find_k_obs** function with k = 5, print:\n",
        "  - Sequences\n",
        "  - Likelihoods\n",
        "  - Most probable hidden states by Viterbi algorithm for specific sequence\n",
        "\n",
        "**NOTE:** Take the length of the possible observation sequences as 4."
      ],
      "metadata": {
        "id": "S45N6CBghTJU"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import itertools\n",
        "\n",
        "def generate_strings(alphabets, k):\n",
        "    return [''.join(combination) for combination in itertools.product(alphabets, repeat=k)]"
      ],
      "metadata": {
        "id": "L22KHqgdEqkN"
      },
      "execution_count": 23,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def find_k_obs(hmm_model, possible_observation, k_seq = 5):\n",
        "  prob_dict = {}\n",
        "  best_seqs, best_likelihoods = [], []\n",
        "  all_seqs = generate_strings(possible_observation, 4)\n",
        "  for seq in all_seqs:\n",
        "    prob = hmm_model.forward_algo(seq)\n",
        "    prob_dict[seq] = prob\n",
        "  sorted_prob_dict = sorted(prob_dict.items(), key=lambda x: x[1], reverse=True)\n",
        "  for key, value in sorted_prob_dict:\n",
        "    best_seqs.append(list(key))\n",
        "    best_likelihoods.append(value)\n",
        "  return best_seqs, best_likelihoods\n",
        "\n",
        "  \"\"\"\n",
        "  The function returns the k most likely sequences and their likelihoods\n",
        "\n",
        "  Arguments:\n",
        "  hmm_model: hmm model object to perform forward and viterbi algorithm\n",
        "  possible observations: Unique set of possible observations\n",
        "  k_seq: The number of most likely observation sequences to return\n",
        "\n",
        "  Return:\n",
        "    best_seq (list): Includes k most likely sequences\n",
        "    best_likelihood (list): Includes likelihoods of k most likely sequences\n",
        "\n",
        "  Example:\n",
        "\n",
        "  possible_observations = ['X','Y']\n",
        "  best_obs_seq, best_likelihood = find_k_obs(hmm_model, possible_observation, k_seq = 2)\n",
        "  best_obs_seq: [['X','X'], ['Y','X']]\n",
        "  best_likelihood: [0.09, 0.08]\n",
        "\n",
        "\n",
        "  \"\"\""
      ],
      "metadata": {
        "id": "YHzVDoJTcmme"
      },
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Task 5: Print Sequences, Likelihoods, Most Probable Paths (1 Mark)\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "PH1x1Q03i7Jq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "k = 5\n",
        "seqs, probs = find_k_obs(model, observations, k)\n",
        "for i in range(k):\n",
        "  print(seqs[i], end=\"\\t\")\n",
        "  print(probs[i], end='\\t')\n",
        "  print(model.viterbi(seqs[i]))"
      ],
      "metadata": {
        "id": "7XlQypApcmkw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "5768a403-a2d2-4c7f-fa35-9c1e781f2b5b"
      },
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "['G', 'G', 'G', 'G']\t0.008813228509375\t['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']\n",
            "['G', 'G', 'G', 'C']\t0.008353077309375\t['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']\n",
            "['G', 'G', 'C', 'G']\t0.00831223985625\t['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']\n",
            "['C', 'G', 'G', 'G']\t0.008294296668750001\t['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']\n",
            "['G', 'C', 'G', 'G']\t0.00829376510625\t['cpg_island', 'cpg_island', 'cpg_island', 'cpg_island']\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# References\n",
        "\n",
        "\n",
        "\n",
        "[1] Richard Durbin, Sean R Eddy, Anders Krogh, and Graeme Mitchison. Biological sequence analysis:\n",
        "probabilistic models of proteins and nucleic acids. Cambridge university press, 1998."
      ],
      "metadata": {
        "id": "9rNzCxR3ixQq"
      }
    }
  ]
}