BA10E -profile HMM_week2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 10: Hidden Markov Models 2: Profile HMM construction\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying to align some protein of interest to an alignment of a family of proteins...\n",
    "\n",
    "$Alignment$ contains many space symbols and likely do not represent meaningful characteristics of the family. Ignore columns for which the fraction of space symbols is ≥ a **column removal threshold θ** in the profile alignment. Column removal results in a 5 × 8 **seed alignment**:\n",
    "\n",
    "Given seed alignment $Alignment∗$, build an HMM that models the propensities of symbols in $Alignment∗$ in the **profile matrix** $Profile(Alignment∗)$ , as shown below. \n",
    "\n",
    "Compute the **probability that the HMM emits Text**. If the HMM is designed well, then the more similar Text is to the strings in Alignment∗, the more likely it will be emitted by the HMM.\n",
    "\n",
    "<br>\n",
    "<br> <img src = http://bioinformaticsalgorithms.com/images/HMM/multiple_alignment_2.png width = 500px> \n",
    "<br>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a simple HMM that treats the columns of $Alignment$ as $k$ sequentially linked states called **match states** (adding a panel to the figure from the previous step) denoted $Match(1), . . . , Match(k)$. \n",
    "\n",
    "When the HMM enters state $Match(i)$, it emits symbol $x_i$ with probability equal to the frequency of this symbol in the $i$-th column of $Profile(Alignment∗)$. The HMM then moves into state $Match(i + 1)$ with transition probability equal to 1.\n",
    "\n",
    "The **similarity score** between $Alignment_*$ and $Text$ is the probability $Pr(Text)$ that the HMM for $Alignment∗$ emits $Text$. This score is equal to the product of frequencies in $Profile(Alignment∗)$ corresponding to each symbol of $Text$. For example, the probability that the HMM in the figure below emits ADDAFFDF is\n",
    "\n",
    "$ 1 · (1/4) · (3/4) · (1/5) · 1 · (1/5) · (3/4) · (3/5) = 0.003375 $\n",
    "\n",
    "STOP and Think: What are the limitations of the HMM below?\n",
    "\n",
    "\n",
    "    \n",
    "<br> <img src = http://bioinformaticsalgorithms.com/images/HMM/multiple_alignment_3.png width = 500px> <br>\n",
    "\n",
    "Figure: Adding the diagram of a simple HMM that models the profile matrix introduced on the previous step. The match states Match(i) are abbreviated as Mi . The HMM only has one possible path; it is initially in state Match(1), the transition probability from state Match(i) to state Match(i+1) is equal to 1 for all i, and all other transitions are forbidden. Emission probabilities are equal to frequencies in the profile, e.g., emission probabilities for M2 are 0 for A, 2/4 for C, 1/4 for D, 0 for E, and 1/4 for F."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Profile HMM Problem: Construct a profile HMM from a multiple alignment.\n",
    "\n",
    "**Input**: A multiple alignment Alignment and a threshold θ.\n",
    "**Output**: HMM(Alignment, θ), in the form of transition and emission matrices.\n",
    "\n",
    "\n",
    "Note: Your matrices can be either space-separated or tab-separated.\n",
    "\n",
    "\n",
    "        Sample Input:\n",
    "\n",
    "        0.289\n",
    "        --------\n",
    "        A B C D E\n",
    "        --------\n",
    "        EBA\n",
    "        E-D\n",
    "        EB-\n",
    "        EED\n",
    "        EBD\n",
    "        EBE\n",
    "        E-D\n",
    "        E-D\n",
    "        Sample Output:\n",
    "\n",
    "            S\tI0\tM1\tD1\tI1\tM2\tD2\tI2\tE\t\n",
    "        S\t0\t0\t1.0\t0\t0\t0\t0\t0\t0\n",
    "        I0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n",
    "        M1\t0\t0\t0\t0\t0.625\t0.375\t0\t0\t0\n",
    "        D1\t0\t0\t0\t0\t0\t0\t0\t0\t0\n",
    "        I1\t0\t0\t0\t0\t0\t0.8\t0.2\t0\t0\n",
    "        M2\t0\t0\t0\t0\t0\t0\t0\t0\t1.0\n",
    "        D2\t0\t0\t0\t0\t0\t0\t0\t0\t1.0\n",
    "        I2\t0\t0\t0\t0\t0\t0\t0\t0\t0\n",
    "        E\t0\t0\t0\t0\t0\t0\t0\t0\t0\n",
    "        --------\n",
    "            A\tB\tC\tD\tE\n",
    "        S\t0\t0\t0\t0\t0\n",
    "        I0\t0\t0\t0\t0\t0\n",
    "        M1\t0\t0\t0\t0\t1.0\n",
    "        D1\t0\t0\t0\t0\t0\n",
    "        I1\t0\t0.8\t0\t0\t0.2\n",
    "        M2\t0.143\t0\t0\t0.714\t0.143\n",
    "        D2\t0\t0\t0\t0\t0\n",
    "        I2\t0\t0\t0\t0\t0\n",
    "        E\t0\t0\t0\t0\t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "def parse_HMM(file):\n",
    "    lines = f.readlines()\n",
    "    theta = float(lines[0].strip())\n",
    "    alpha = lines[2].strip().split()\n",
    "    alpha = dict(zip(alpha, range(len(alpha))))\n",
    "    alignment = [line.strip() for line in lines[4:]]\n",
    "    print(\"Alignment:\")\n",
    "    for a in alignment:\n",
    "        print('\\t', a)\n",
    "    print(\"Alpha:\", alpha)\n",
    "    print(\"Theta:\", theta)\n",
    "    return alignment, alpha, theta\n",
    "\n",
    "### algorithm for building profile HMM from an alignment|theta\n",
    "def profileHMM(alignment, alpha, theta):\n",
    "\n",
    "    # # 'seed' indicates which columns pass threshold theta for align*\n",
    "    def get_seed_alignment(alignment, theta, alphabet):\n",
    "        \"\"\"'seed': True if gaps < theta; included in align*\n",
    "            returns seed alignment columns for HMM\"\"\"\n",
    "        k = len(alignment[0])\n",
    "        a = len(alphabet.keys())\n",
    "        print(\"Length of alignment(k) = \",k)\n",
    "        print(\"size of alphabet = \",k,'\\n')\n",
    "        freq = np.zeros(shape=(a+1, k)) \n",
    "\n",
    "        for seq in alignment:\n",
    "            for i in range(k):\n",
    "                if seq[i] == '-':\n",
    "#                     print('\\tgap', seq[i])\n",
    "                    freq[a][i] += 1\n",
    "                else:\n",
    "                    freq[alphabet[seq[i]]][i] +=1\n",
    "#                     print('\\tbase', seq[i])\n",
    "        n = len(alignment)\n",
    "        print('freq matrix:\\n', freq)\n",
    "        seed = [x/n < theta for x in freq[a]]\n",
    "        return seed\n",
    "    \n",
    "    def normalize_matrices(T,E):\n",
    "        \"\"\"Normalize Transmission and Emission frequencies to sum of row for state\"\"\"\n",
    "        for state in range(len(S)):\n",
    "            if sum(T[state]) > 0:\n",
    "                T[state] =  T[state]/sum(T[state])\n",
    "            if sum(E[state]) >0:\n",
    "                E[state] =  E[state]/sum(E[state])\n",
    "        return T, E\n",
    "    \n",
    "    def state_transition(T, prev, kind, S):\n",
    "        \"\"\" find next Match, Del, Ins, or End state in [states] \"\"\"\n",
    "        x=0\n",
    "        if S[prev][0] == 'M':\n",
    "            x=1\n",
    "        for nxt in range(prev+1+x, len(T)):\n",
    "            if S[nxt][0] == kind[0]:\n",
    "                T[prev][nxt] += 1\n",
    "#                 print('  >Transition added:',((S[prev],prev),(S[nxt], nxt)))\n",
    "                return T, nxt\n",
    "    \n",
    "    ### Walk each sequence to count transitions, emissions\n",
    "    seed = get_seed_alignment(alignment, theta, alpha)\n",
    "    print('seed:',seed)\n",
    "    n = len(alignment)\n",
    "    k = len(alignment[0])\n",
    "    S = ['S','I0']+list(c+str(n) for n in range(1,sum(seed)+1) for c in \"MDI\")+['E']\n",
    "    E = np.zeros(shape=(len(S), len(alpha.keys())))\n",
    "    T = np.zeros(shape=(len(S), len(S)))\n",
    "\n",
    "    for seq in alignment:\n",
    "        # 'state' is the hidden state col/row; i is pos in alignment\n",
    "        state = 0; i = 0\n",
    "        while i < k:\n",
    "            # seed: either match or del as hidden state\n",
    "            if seed[i]:\n",
    "                if seq[i] in alpha:\n",
    "                    T, state = state_transition(T, state, 'Match', S)\n",
    "                    E[state][alpha[seq[i]]] +=1\n",
    "                else:\n",
    "                    T, state = state_transition(T, state, 'Deletion', S)\n",
    "            # not seed: either insert or nothing\n",
    "            else:\n",
    "                # count any emissions before next seed column\n",
    "                emits = []\n",
    "                while not seed[i]:\n",
    "                    if seq[i] in alpha:\n",
    "                        emits.append(seq[i])\n",
    "                    i += 1\n",
    "                    if i == k: \n",
    "                        break  # hit end of sequence\n",
    "                i-=1                \n",
    "                if len(emits)>0:\n",
    "                    # count the transition(state, insertion)\n",
    "                    T, state = state_transition(T, state, 'Insert', S)\n",
    "                    # get all emissions(state) in insertion\n",
    "                    for symbol in emits:\n",
    "                        E[state][alpha[symbol]] += 1\n",
    "                    # count all symbols as cyclic transition t(ins_x, ins_x)\n",
    "                    if len(emits) >1:\n",
    "                        T[state][state] += len(emits)-1\n",
    "                # else do nothing, just gaps in align\n",
    "            i += 1\n",
    "        # from last state to 'End'    \n",
    "        T, state = state_transition(T, state, 'End', S)\n",
    "\n",
    "    T,E = normalize_matrices(T,E)\n",
    "    return S, T, E\n",
    "\n",
    "def outputHMM(T,E,S,alpha):\n",
    "    lists = []\n",
    "    lists.append(['']+S)\n",
    "    for i in range(len(S)):\n",
    "        # use {0:g} format to remove trailing zeros\n",
    "        lists.append([S[i]]+['{0:g}'.format(x) for x in T[i]])\n",
    "    lists.append('--------')\n",
    "    lists.append([''] + alpha)\n",
    "    for i in range(len(S)):\n",
    "        lists.append([S[i]]+['{0:g}'.format(x) for x in E[i]])\n",
    "    return lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alignment:\n",
      "\t DCDABACED\n",
      "\t DCCA--CA-\n",
      "\t DCDAB-CA-\n",
      "\t BCDA---A-\n",
      "\t BC-ABE-AE\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4}\n",
      "Theta: 0.252\n",
      "Length of alignment(k) =  9\n",
      "size of alphabet =  9 \n",
      "\n",
      "freq matrix:\n",
      " [[0. 0. 0. 5. 0. 1. 0. 4. 0.]\n",
      " [2. 0. 0. 0. 3. 0. 0. 0. 0.]\n",
      " [0. 5. 1. 0. 0. 0. 3. 0. 0.]\n",
      " [3. 0. 3. 0. 0. 0. 0. 0. 1.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 1. 1.]\n",
      " [0. 0. 1. 0. 2. 3. 2. 0. 3.]]\n",
      "seed: [True, True, True, True, False, False, False, True, False]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(True,\n",
       " True,\n",
       "        0     1     2     3     4     5     6     7     8     9    10    11  \\\n",
       " 0   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 1   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 2   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 3   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 4   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 5   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 6   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 7   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 8   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 9   True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 10  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 11  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 12  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 13  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 14  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 15  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 16  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " 17  True  True  True  True  True  True  True  True  True  True  True  True   \n",
       " \n",
       "       12    13    14    15    16    17  \n",
       " 0   True  True  True  True  True  True  \n",
       " 1   True  True  True  True  True  True  \n",
       " 2   True  True  True  True  True  True  \n",
       " 3   True  True  True  True  True  True  \n",
       " 4   True  True  True  True  True  True  \n",
       " 5   True  True  True  True  True  True  \n",
       " 6   True  True  True  True  True  True  \n",
       " 7   True  True  True  True  True  True  \n",
       " 8   True  True  True  True  True  True  \n",
       " 9   True  True  True  True  True  True  \n",
       " 10  True  True  True  True  True  True  \n",
       " 11  True  True  True  True  True  True  \n",
       " 12  True  True  True  True  True  True  \n",
       " 13  True  True  True  True  True  True  \n",
       " 14  True  True  True  True  True  True  \n",
       " 15  True  True  True  True  True  True  \n",
       " 16  True  True  True  True  True  True  \n",
       " 17  True  True  True  True  True  True  )"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### Extra test data 1:\n",
    "\n",
    "\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/10e.extra.txt') as f:\n",
    "    alignment, alpha, theta = parse_HMM(f)\n",
    "\n",
    "S, T, E = profileHMM(alignment, alpha, theta)\n",
    "T = np.around(T,3)\n",
    "E = np.around(E,3)\n",
    "alpha = list(alpha.keys())\n",
    "lists =outputHMM(T,E,S, alpha)\n",
    "with open('10e.out.txt','w') as file:\n",
    "    for line in lists[:-1]:\n",
    "        if line[0] != '-':\n",
    "            file.write('\\t'.join(line))\n",
    "        else:\n",
    "            file.write(line)\n",
    "        file.write('\\n')\n",
    "    file.write('\\t'.join(lists[-1]))\n",
    "\n",
    "    \n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/10e.extra.E.csv', 'r') as f:\n",
    "    e = pd.read_csv(f, header=None)\n",
    "    solE = np.array(e)\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/10e.extra.T.csv', 'r') as f:\n",
    "    t = pd.read_csv(f, header=None)\n",
    "    solT = np.array(t)\n",
    "np.all(T== solT),np.all(E== solE), pd.DataFrame(T == solT)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Challenge dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "could not convert string to float: 'S\\tI0\\tM1\\tD1\\tI1\\tM2\\tD2\\tI2\\tM3\\tD3\\tI3\\tM4\\tD4\\tI4\\tE'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-87-ea27f3e625f5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/rosalind_ba10e.txt'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0malignment\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtheta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparse_HMM\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprofileHMM\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malignment\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtheta\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-85-4184144fc89b>\u001b[0m in \u001b[0;36mparse_HMM\u001b[0;34m(file)\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mparse_HMM\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mlines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreadlines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m     \u001b[0mtheta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlines\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m     \u001b[0malpha\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlines\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0malpha\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: could not convert string to float: 'S\\tI0\\tM1\\tD1\\tI1\\tM2\\tD2\\tI2\\tM3\\tD3\\tI3\\tM4\\tD4\\tI4\\tE'"
     ]
    }
   ],
   "source": [
    "## rosalind data\n",
    "\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/rosalind_ba10e.txt') as f:\n",
    "    alignment, alpha, theta = parse_HMM(f)\n",
    "\n",
    "S, T, E = profileHMM(alignment, alpha, theta)\n",
    "T = np.around(T,3)\n",
    "E = np.around(E,3)\n",
    "alpha = list(alpha.keys())\n",
    "lists =outputHMM(T,E,S, alpha)\n",
    "with open(\"/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/rosalind_ba10e.output.txt\",'w') as file:\n",
    "    for line in lists[:-1]:\n",
    "        if line[0] != '-':\n",
    "            file.write('\\t'.join(line))\n",
    "        else:\n",
    "            file.write(line)\n",
    "        file.write('\\n')\n",
    "    file.write('\\t'.join(lists[-1]))\n",
    "\n",
    "print(\"\\nT:\", pd.DataFrame(T), \"\\nE:\", pd.DataFrame(E))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alignment:\n",
      "\t BEA-E-AD-\n",
      "\t DC-DDCAEC\n",
      "\t BE-CE-ABC\n",
      "\t -EAC-CADE\n",
      "\t B-ACECEDE\n",
      "\t EEA--CADC\n",
      "\t BED--CA-C\n",
      "\t BEA-E-ADC\n",
      "Alpha: {'A': 0, 'B': 1, 'C': 2, 'D': 3, 'E': 4}\n",
      "Theta: 0.382\n",
      "Length of alignment(k) =  9\n",
      "size of alphabet =  9 \n",
      "\n",
      "freq matrix:\n",
      " [[0. 0. 5. 0. 0. 0. 7. 0. 0.]\n",
      " [5. 0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 1. 0. 3. 0. 5. 0. 0. 5.]\n",
      " [1. 0. 1. 1. 1. 0. 0. 5. 0.]\n",
      " [1. 6. 0. 0. 4. 0. 1. 1. 2.]\n",
      " [1. 1. 2. 4. 3. 3. 0. 1. 1.]]\n",
      "seed: [True, True, True, False, True, True, True, True, True]\n",
      "\n",
      "\n",
      "seq: BEA-E-AD- \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "***no insert in this bitch ***, state is still M3\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 4 adding to ? trans from state M3 8\n",
      "  >findtrans S[prev]: M3 8 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('M4', 11))\n",
      "new state is M4 11\n",
      "\n",
      "pos.i: 5 adding to ? trans from state M4 11\n",
      "  >findtrans S[prev]: M4 11 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M4', 11), ('D5', 15))\n",
      "new state is D5 15\n",
      "\n",
      "pos.i: 6 adding to ? trans from state D5 15\n",
      "  >findtrans S[prev]: D5 15 to Match\n",
      "  >Transition added: (('D5', 15), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('D8', 24))\n",
      "new state is D8 24\n",
      "  >findtrans S[prev]: D8 24 to End\n",
      "  >Transition added: (('D8', 24), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: BEA-E-AD- \n",
      " S - M1 - M2 - M3 - M4 - D5 - M6 - M7 - D8 - E\n",
      "\n",
      "\n",
      "seq: DC-DDCAEC \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('D3', 9))\n",
      "new state is D3 9\n",
      "\n",
      "pos.i: 3 adding to ? trans from state D3 9\n",
      "not seed, check insert:\n",
      "*I* > emits ['D'] insert @ state: (9, 'D3'), in col i:3 \n",
      "  >findtrans S[prev]: D3 9 to Insert\n",
      "  >Transition added: (('D3', 9), ('I3', 10))\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2 - D3\n",
      "new state is I3 10\n",
      "\n",
      "pos.i: 4 adding to ? trans from state I3 10\n",
      "  >findtrans S[prev]: I3 10 to Match\n",
      "  >Transition added: (('I3', 10), ('M4', 11))\n",
      "new state is M4 11\n",
      "\n",
      "pos.i: 5 adding to ? trans from state M4 11\n",
      "  >findtrans S[prev]: M4 11 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M4', 11), ('M5', 14))\n",
      "new state is M5 14\n",
      "\n",
      "pos.i: 6 adding to ? trans from state M5 14\n",
      "  >findtrans S[prev]: M5 14 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M5', 14), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: DC-DDCAEC \n",
      " S - M1 - M2 - D3 - I3 - M4 - M5 - M6 - M7 - M8 - E\n",
      "\n",
      "\n",
      "seq: BE-CE-ABC \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('D3', 9))\n",
      "new state is D3 9\n",
      "\n",
      "pos.i: 3 adding to ? trans from state D3 9\n",
      "not seed, check insert:\n",
      "*I* > emits ['C'] insert @ state: (9, 'D3'), in col i:3 \n",
      "  >findtrans S[prev]: D3 9 to Insert\n",
      "  >Transition added: (('D3', 9), ('I3', 10))\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2 - D3\n",
      "new state is I3 10\n",
      "\n",
      "pos.i: 4 adding to ? trans from state I3 10\n",
      "  >findtrans S[prev]: I3 10 to Match\n",
      "  >Transition added: (('I3', 10), ('M4', 11))\n",
      "new state is M4 11\n",
      "\n",
      "pos.i: 5 adding to ? trans from state M4 11\n",
      "  >findtrans S[prev]: M4 11 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M4', 11), ('D5', 15))\n",
      "new state is D5 15\n",
      "\n",
      "pos.i: 6 adding to ? trans from state D5 15\n",
      "  >findtrans S[prev]: D5 15 to Match\n",
      "  >Transition added: (('D5', 15), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: BE-CE-ABC \n",
      " S - M1 - M2 - D3 - I3 - M4 - D5 - M6 - M7 - M8 - E\n",
      "\n",
      "\n",
      "seq: -EAC-CADE \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Deletion\n",
      "  >Transition added: (('S', 0), ('D1', 3))\n",
      "new state is D1 3\n",
      "\n",
      "pos.i: 1 adding to ? trans from state D1 3\n",
      "  >findtrans S[prev]: D1 3 to Match\n",
      "  >Transition added: (('D1', 3), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "*I* > emits ['C'] insert @ state: (8, 'M3'), in col i:3 \n",
      "  >findtrans S[prev]: M3 8 to Insert\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('I3', 10))\n",
      "  > hidden path after ins:\n",
      "\t S - D1 - M2 - M3\n",
      "new state is I3 10\n",
      "\n",
      "pos.i: 4 adding to ? trans from state I3 10\n",
      "  >findtrans S[prev]: I3 10 to Deletion\n",
      "  >Transition added: (('I3', 10), ('D4', 12))\n",
      "new state is D4 12\n",
      "\n",
      "pos.i: 5 adding to ? trans from state D4 12\n",
      "  >findtrans S[prev]: D4 12 to Match\n",
      "  >Transition added: (('D4', 12), ('M5', 14))\n",
      "new state is M5 14\n",
      "\n",
      "pos.i: 6 adding to ? trans from state M5 14\n",
      "  >findtrans S[prev]: M5 14 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M5', 14), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: -EAC-CADE \n",
      " S - D1 - M2 - M3 - I3 - D4 - M5 - M6 - M7 - M8 - E\n",
      "\n",
      "\n",
      "seq: B-ACECEDE \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('D2', 6))\n",
      "new state is D2 6\n",
      "\n",
      "pos.i: 2 adding to ? trans from state D2 6\n",
      "  >findtrans S[prev]: D2 6 to Match\n",
      "  >Transition added: (('D2', 6), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "*I* > emits ['C'] insert @ state: (8, 'M3'), in col i:3 \n",
      "  >findtrans S[prev]: M3 8 to Insert\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('I3', 10))\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - D2 - M3\n",
      "new state is I3 10\n",
      "\n",
      "pos.i: 4 adding to ? trans from state I3 10\n",
      "  >findtrans S[prev]: I3 10 to Match\n",
      "  >Transition added: (('I3', 10), ('M4', 11))\n",
      "new state is M4 11\n",
      "\n",
      "pos.i: 5 adding to ? trans from state M4 11\n",
      "  >findtrans S[prev]: M4 11 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M4', 11), ('M5', 14))\n",
      "new state is M5 14\n",
      "\n",
      "pos.i: 6 adding to ? trans from state M5 14\n",
      "  >findtrans S[prev]: M5 14 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M5', 14), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: B-ACECEDE \n",
      " S - M1 - D2 - M3 - I3 - M4 - M5 - M6 - M7 - M8 - E\n",
      "\n",
      "\n",
      "seq: EEA--CADC \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "***no insert in this bitch ***, state is still M3\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 4 adding to ? trans from state M3 8\n",
      "  >findtrans S[prev]: M3 8 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('D4', 12))\n",
      "new state is D4 12\n",
      "\n",
      "pos.i: 5 adding to ? trans from state D4 12\n",
      "  >findtrans S[prev]: D4 12 to Match\n",
      "  >Transition added: (('D4', 12), ('M5', 14))\n",
      "new state is M5 14\n",
      "\n",
      "pos.i: 6 adding to ? trans from state M5 14\n",
      "  >findtrans S[prev]: M5 14 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M5', 14), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: EEA--CADC \n",
      " S - M1 - M2 - M3 - D4 - M5 - M6 - M7 - M8 - E\n",
      "\n",
      "\n",
      "seq: BED--CA-C \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "***no insert in this bitch ***, state is still M3\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 4 adding to ? trans from state M3 8\n",
      "  >findtrans S[prev]: M3 8 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('D4', 12))\n",
      "new state is D4 12\n",
      "\n",
      "pos.i: 5 adding to ? trans from state D4 12\n",
      "  >findtrans S[prev]: D4 12 to Match\n",
      "  >Transition added: (('D4', 12), ('M5', 14))\n",
      "new state is M5 14\n",
      "\n",
      "pos.i: 6 adding to ? trans from state M5 14\n",
      "  >findtrans S[prev]: M5 14 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M5', 14), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('D7', 21))\n",
      "new state is D7 21\n",
      "\n",
      "pos.i: 8 adding to ? trans from state D7 21\n",
      "  >findtrans S[prev]: D7 21 to Match\n",
      "  >Transition added: (('D7', 21), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: BED--CA-C \n",
      " S - M1 - M2 - M3 - D4 - M5 - M6 - D7 - M8 - E\n",
      "\n",
      "\n",
      "seq: BEA-E-ADC \n",
      "\n",
      "-------------------------------\n",
      "\n",
      "\n",
      "pos.i: 0 adding to ? trans from state S 0\n",
      "  >findtrans S[prev]: S 0 to Match\n",
      "  >Transition added: (('S', 0), ('M1', 2))\n",
      "new state is M1 2\n",
      "\n",
      "pos.i: 1 adding to ? trans from state M1 2\n",
      "  >findtrans S[prev]: M1 2 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M1', 2), ('M2', 5))\n",
      "new state is M2 5\n",
      "\n",
      "pos.i: 2 adding to ? trans from state M2 5\n",
      "  >findtrans S[prev]: M2 5 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M2', 5), ('M3', 8))\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 3 adding to ? trans from state M3 8\n",
      "not seed, check insert:\n",
      "***no insert in this bitch ***, state is still M3\n",
      "  > hidden path after ins:\n",
      "\t S - M1 - M2\n",
      "new state is M3 8\n",
      "\n",
      "pos.i: 4 adding to ? trans from state M3 8\n",
      "  >findtrans S[prev]: M3 8 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M3', 8), ('M4', 11))\n",
      "new state is M4 11\n",
      "\n",
      "pos.i: 5 adding to ? trans from state M4 11\n",
      "  >findtrans S[prev]: M4 11 to Deletion\n",
      "   >offset =1\n",
      "  >Transition added: (('M4', 11), ('D5', 15))\n",
      "new state is D5 15\n",
      "\n",
      "pos.i: 6 adding to ? trans from state D5 15\n",
      "  >findtrans S[prev]: D5 15 to Match\n",
      "  >Transition added: (('D5', 15), ('M6', 17))\n",
      "new state is M6 17\n",
      "\n",
      "pos.i: 7 adding to ? trans from state M6 17\n",
      "  >findtrans S[prev]: M6 17 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M6', 17), ('M7', 20))\n",
      "new state is M7 20\n",
      "\n",
      "pos.i: 8 adding to ? trans from state M7 20\n",
      "  >findtrans S[prev]: M7 20 to Match\n",
      "   >offset =1\n",
      "  >Transition added: (('M7', 20), ('M8', 23))\n",
      "new state is M8 23\n",
      "  >findtrans S[prev]: M8 23 to End\n",
      "   >offset =1\n",
      "  >Transition added: (('M8', 23), ('E', 26))\n",
      "HIDDEN PATH FOR SEQ: BEA-E-ADC \n",
      " S - M1 - M2 - M3 - M4 - D5 - M6 - M7 - M8 - E\n"
     ]
    }
   ],
   "source": [
    "## dataset_26258_15\n",
    "\n",
    "with open('/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data/dataset_26258_15.txt') as f:\n",
    "    alignment, alpha, theta = parse_HMM(f)\n",
    "\n",
    "S, T, E = profileHMM(alignment, alpha, theta)\n",
    "T = np.around(T,3)\n",
    "E = np.around(E,3)\n",
    "alpha = list(alpha.keys())\n",
    "lists =outputHMM(T,E,S, alpha)\n",
    "with open('10e.finalout.txt','w') as file:\n",
    "    for line in lists[:-1]:\n",
    "        if line[0] != '-':\n",
    "            file.write('\\t'.join(line))\n",
    "        else:\n",
    "            file.write(line)\n",
    "        file.write('\\n')\n",
    "    file.write('\\t'.join(lists[-1]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "## FORMAT OUTPUT ~ the major issue that cost me hours of wasted efforts\n",
    "\n",
    "# Column of row titles for output (S, I0, M1, D1, I1,....In, E)\n",
    "\n",
    "rows = np.array(S, dtype='|S20')[:, np.newaxis]\n",
    "# use hstack to put header row(s) on top of array\n",
    "with open('test.csv', 'w') as f:\n",
    "    np.savetxt(f, np.hstack((rows, data)), delimiter=', ', fmt='%s')\n",
    "\n",
    "# how to drop all trailing 0's from float\n",
    "# '{0:g}'.format(somefloat)\n",
    "\n",
    "# how to round floats in np.array( np.around(matrix, digits)\n",
    "T = np.around(T,3)\n",
    "E = np.around(E,3)\n"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}